Intrepid2
Intrepid2_CellToolsDefNodeInfo.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid2 Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38 // Mauro Perego (mperego@sandia.gov)
39 //
40 // ************************************************************************
41 // @HEADER
42 
43 
49 #ifndef __INTREPID2_CELLTOOLS_DEF_NODE_INFO_HPP__
50 #define __INTREPID2_CELLTOOLS_DEF_NODE_INFO_HPP__
51 
52 // disable clang warnings
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
55 #endif
56 
57 namespace Intrepid2 {
58 
59  //============================================================================================//
60  // //
61  // Reference nodes //
62  // //
63  //============================================================================================//
64 
65  template<typename DeviceType>
66  template<typename cellCenterValueType, class ...cellCenterProperties>
67  void
68  CellTools<DeviceType>::
69  getReferenceCellCenter( Kokkos::DynRankView<cellCenterValueType,cellCenterProperties...> cellCenter,
70  const shards::CellTopology cell ) {
71 #ifdef HAVE_INTREPID2_DEBUG
72  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(cell), std::invalid_argument,
73  ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): the specified cell topology does not have a reference cell." );
74 
75  INTREPID2_TEST_FOR_EXCEPTION( cellCenter.rank() != 1, std::invalid_argument,
76  ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): cellCenter must have rank 1." );
77 
78  INTREPID2_TEST_FOR_EXCEPTION( cellCenter.extent(0) < cell.getDimension(), std::invalid_argument,
79  ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): cellCenter must have dimension at least as large as cell.getDimension()." );
80 #endif
81 
82  constexpr bool is_accessible =
83  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
84  typename decltype(cellCenter)::memory_space>::accessible;
85  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceCellCenter(..): output view's memory space is not compatible with DeviceType");
86 
87  const ordinal_type dim = cell.getDimension();
88 
89  const auto refCellCenter = RefCellCenter<DeviceType>::get(cell.getKey());
90 
91  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,dim),
92  KOKKOS_LAMBDA (const int &i) {cellCenter(i) = refCellCenter(i);}
93  );
94  }
95 
96 
97  template<typename DeviceType>
98  template<typename cellVertexValueType, class ...cellVertexProperties>
99  void
100  CellTools<DeviceType>::
101  getReferenceVertex( Kokkos::DynRankView<cellVertexValueType,cellVertexProperties...> cellVertex,
102  const shards::CellTopology cell,
103  const ordinal_type vertexOrd ) {
104 #ifdef HAVE_INTREPID2_DEBUG
105  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(cell), std::invalid_argument,
106  ">>> ERROR (Intrepid2::CellTools::getReferenceVertex): the specified cell topology does not have a reference cell." );
107 
108  INTREPID2_TEST_FOR_EXCEPTION( (vertexOrd < 0) || vertexOrd > static_cast<ordinal_type>(cell.getVertexCount()), std::invalid_argument,
109  ">>> ERROR (Intrepid2::CellTools::getReferenceVertex): invalid node ordinal for the specified cell topology." );
110 
111  INTREPID2_TEST_FOR_EXCEPTION( cellVertex.rank() != 1, std::invalid_argument,
112  ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNodes must have rank 1." );
113 
114  INTREPID2_TEST_FOR_EXCEPTION( cellVertex.extent(0) < cell.getDimension(), std::invalid_argument,
115  ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNodes must have dimension at least as large as cell.getDimension()." );
116 #endif
117 
118  constexpr bool is_accessible =
119  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
120  typename decltype(cellVertex)::memory_space>::accessible;
121  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceVertex(..): output view's memory space is not compatible with DeviceType");
122 
123  const auto refNodes = RefCellNodes<DeviceType>::get(cell.getKey());
124 
125  ordinal_type dim = cell.getDimension();
126  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,dim),
127  KOKKOS_LAMBDA (const int &i) {cellVertex(i) = refNodes(vertexOrd,i);}
128  );
129  }
130 
131 
132  template<typename DeviceType>
133  template<typename subcellVertexValueType, class ...subcellVertexProperties>
134  void
135  CellTools<DeviceType>::
136  getReferenceSubcellVertices( Kokkos::DynRankView<subcellVertexValueType,subcellVertexProperties...> subcellVertices,
137  const ordinal_type subcellDim,
138  const ordinal_type subcellOrd,
139  const shards::CellTopology parentCell ) {
140 #ifdef HAVE_INTREPID2_DEBUG
141  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
142  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): the specified cell topology does not have a reference cell." );
143 
144  INTREPID2_TEST_FOR_EXCEPTION( subcellDim > static_cast<ordinal_type>(parentCell.getDimension()), std::invalid_argument,
145  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcell dimension cannot exceed cell dimension." );
146 
147  INTREPID2_TEST_FOR_EXCEPTION( subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
148  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcell ordinal cannot exceed subcell count." );
149 
150  // Verify subcellVertices rank and dimensions
151  INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.rank() != 2, std::invalid_argument,
152  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces must have rank 2." );
153 
154  // need to match to node count as it invokes getReferenceSubcellNodes
155  INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.extent(0) != parentCell.getNodeCount(subcellDim, subcellOrd), std::invalid_argument,
156  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces dimension(0) must match to parent cell vertex count." );
157 
158  INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.extent(1) != parentCell.getDimension(), std::invalid_argument,
159  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces dimension(1) must match to parent cell dimension." );
160 #endif
161  getReferenceSubcellNodes(subcellVertices,
162  subcellDim,
163  subcellOrd,
164  parentCell);
165  }
166 
167 
168  template<typename DeviceType>
169  template<typename cellNodeValueType, class ...cellNodeProperties>
170  void
171  CellTools<DeviceType>::
172  getReferenceNode( Kokkos::DynRankView<cellNodeValueType,cellNodeProperties...> cellNode,
173  const shards::CellTopology cell,
174  const ordinal_type nodeOrd ) {
175 #ifdef HAVE_INTREPID2_DEBUG
176  INTREPID2_TEST_FOR_EXCEPTION( nodeOrd >= static_cast<ordinal_type>(cell.getNodeCount()), std::invalid_argument,
177  ">>> ERROR (Intrepid2::CellTools::getReferenceNode): invalid node ordinal for the specified cell topology." );
178 
179  INTREPID2_TEST_FOR_EXCEPTION( cellNode.rank() != 1, std::invalid_argument,
180  ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNode must have rank 1." );
181 
182  INTREPID2_TEST_FOR_EXCEPTION( cellNode.extent(0) < cell.getDimension(), std::invalid_argument,
183  ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNode must have dimension at least as large as cell.getDimension()." );
184 #endif
185 
186  constexpr bool is_accessible =
187  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
188  typename decltype(cellNode)::memory_space>::accessible;
189  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceNode(..): output view's memory space is not compatible with DeviceType");
190 
191  const auto refNodes = RefCellNodes<DeviceType>::get(cell.getKey());
192 
193  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,cell.getDimension()),
194  KOKKOS_LAMBDA (const int &i) {cellNode(i) = refNodes(nodeOrd,i);}
195  );
196  }
197 
198  template<typename DeviceType>
199  template<typename subcellNodeValueType, class ...subcellNodeProperties>
200  void
201  CellTools<DeviceType>::
202  getReferenceSubcellNodes( Kokkos::DynRankView<subcellNodeValueType,subcellNodeProperties...> subcellNodes,
203  const ordinal_type subcellDim,
204  const ordinal_type subcellOrd,
205  const shards::CellTopology parentCell ) {
206 #ifdef HAVE_INTREPID2_DEBUG
207  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
208  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): the specified cell topology does not have a reference cell.");
209 
210  INTREPID2_TEST_FOR_EXCEPTION( subcellDim > static_cast<ordinal_type>(parentCell.getDimension()), std::invalid_argument,
211  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcell dimension out of range.");
212 
213  INTREPID2_TEST_FOR_EXCEPTION( subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
214  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcell ordinal out of range.");
215 
216  // Verify subcellNodes rank and dimensions
217  INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.rank() != 2, std::invalid_argument,
218  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes must have rank 2.");
219 
220  INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.extent(0) != parentCell.getNodeCount(subcellDim, subcellOrd), std::invalid_argument,
221  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes dimension (0) must match to node count in cell.");
222 
223  INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.extent(1) != parentCell.getDimension(), std::invalid_argument,
224  ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes dimension (1) must match to cell dimension.");
225 #endif
226 
227  // Find how many cellWorkset does the specified subcell have.
228  const auto subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
229 
230  // Loop over subcell cellWorkset
231  for (size_type subcNodeOrd=0;subcNodeOrd<subcNodeCount;++subcNodeOrd) {
232  // Get the node number relative to the parent reference cell
233  const auto cellNodeOrd = parentCell.getNodeMap(subcellDim, subcellOrd, subcNodeOrd);
234 
235  auto dst = Kokkos::subdynrankview(subcellNodes, subcNodeOrd, Kokkos::ALL());
236  getReferenceNode(dst, parentCell, cellNodeOrd);
237  }
238  }
239 
240  template<typename DeviceType>
241  template<typename refEdgeTangentValueType, class ...refEdgeTangentProperties>
242  void
243  CellTools<DeviceType>::
244  getReferenceEdgeTangent( Kokkos::DynRankView<refEdgeTangentValueType,refEdgeTangentProperties...> refEdgeTangent,
245  const ordinal_type edgeOrd,
246  const shards::CellTopology parentCell ) {
247 #ifdef HAVE_INTREPID2_DEBUG
248  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
249  parentCell.getDimension() != 3, std::invalid_argument,
250  ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): two or three-dimensional parent cell required");
251 
252  INTREPID2_TEST_FOR_EXCEPTION( refEdgeTangent.rank() != 1, std::invalid_argument,
253  ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): rank = 1 required for output arrays");
254 
255  INTREPID2_TEST_FOR_EXCEPTION( refEdgeTangent.extent(0) != parentCell.getDimension(), std::invalid_argument,
256  ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): output array size is required to match space dimension");
257 
258  INTREPID2_TEST_FOR_EXCEPTION( edgeOrd < 0 ||
259  edgeOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(1)), std::invalid_argument,
260  ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): edge ordinal out of bounds");
261 
262 #endif
263  constexpr bool is_accessible =
264  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
265  typename decltype(refEdgeTangent)::memory_space>::accessible;
266  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceEdgeTangent(..): output view's memory space is not compatible with DeviceType");
267 
268  const auto edgeMap = RefSubcellParametrization<DeviceType>::get(1, parentCell.getKey());
269 
270  // All ref. edge maps have affine coordinate functions: f_dim(u) = C_0(dim) + C_1(dim)*u,
271  // => edge Tangent: -> C_1(*)
272  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,parentCell.getDimension()),
273  KOKKOS_LAMBDA (const int &i) {refEdgeTangent(i) = edgeMap(edgeOrd, i, 1);}
274  );
275  }
276 
277 
278  template<typename DeviceType>
279  template<typename refFaceTanValueType, class ...refFaceTanProperties>
280  void
281  CellTools<DeviceType>::
282  getReferenceFaceTangents( Kokkos::DynRankView<refFaceTanValueType,refFaceTanProperties...> refFaceTanU,
283  Kokkos::DynRankView<refFaceTanValueType,refFaceTanProperties...> refFaceTanV,
284  const ordinal_type faceOrd,
285  const shards::CellTopology parentCell ) {
286 #ifdef HAVE_INTREPID2_DEBUG
287  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
288  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): three-dimensional parent cell required");
289 
290  INTREPID2_TEST_FOR_EXCEPTION( faceOrd < 0 || faceOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(2)), std::invalid_argument,
291  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): face ordinal out of bounds");
292 
293  INTREPID2_TEST_FOR_EXCEPTION( refFaceTanU.rank() != 1 || refFaceTanV.rank() != 1, std::invalid_argument,
294  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): rank = 1 required for output arrays");
295 
296  INTREPID2_TEST_FOR_EXCEPTION( refFaceTanU.extent(0) != parentCell.getDimension(), std::invalid_argument,
297  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
298 
299  INTREPID2_TEST_FOR_EXCEPTION( refFaceTanV.extent(0) != parentCell.getDimension(), std::invalid_argument,
300  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
301 #endif
302  constexpr bool is_accessible =
303  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
304  typename decltype(refFaceTanU)::memory_space>::accessible;
305  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceFaceTangents(..): output views' memory spaces are not compatible with DeviceType");
306 
307 
308  const auto faceMap = RefSubcellParametrization<DeviceType>::get(2, parentCell.getKey());
309 
310  /* All ref. face maps have affine coordinate functions: f_dim(u,v) = C_0(dim) + C_1(dim)*u + C_2(dim)*v
311  * ` => Tangent vectors are: refFaceTanU -> C_1(*); refFaceTanV -> C_2(*)
312  */
313 
314  // set refFaceTanU -> C_1(*)
315  // set refFaceTanV -> C_2(*)
316  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,parentCell.getDimension()),
317  KOKKOS_LAMBDA (const int &i) {
318  refFaceTanU(i) = faceMap(faceOrd, i, 1);
319  refFaceTanV(i) = faceMap(faceOrd, i, 2);
320  });
321  }
322 
323  template<typename DeviceType>
324  template<typename refSideNormalValueType, class ...refSideNormalProperties>
325  void
326  CellTools<DeviceType>::
327  getReferenceSideNormal( Kokkos::DynRankView<refSideNormalValueType,refSideNormalProperties...> refSideNormal,
328  const ordinal_type sideOrd,
329  const shards::CellTopology parentCell ) {
330 #ifdef HAVE_INTREPID2_DEBUG
331  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
332  parentCell.getDimension() != 3, std::invalid_argument,
333  ">>> ERROR (Intrepid2::CellTools::getReferenceSideNormal): two or three-dimensional parent cell required");
334 
335  // Check side ordinal: by definition side is subcell whose dimension = parentCell.getDimension()-1
336  INTREPID2_TEST_FOR_EXCEPTION( sideOrd < 0 || sideOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(parentCell.getDimension() - 1)), std::invalid_argument,
337  ">>> ERROR (Intrepid2::CellTools::getReferenceSideNormal): side ordinal out of bounds");
338 #endif
339  constexpr bool is_accessible = Kokkos::Impl::MemorySpaceAccess<
340  MemSpaceType,
341  typename decltype(refSideNormal)::memory_space>::accessible;
342  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceSideNormal(..): output view's memory space is not compatible with DeviceType");
343 
344  const auto dim = parentCell.getDimension();
345  if (dim == 2) {
346  // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
347  getReferenceEdgeTangent(refSideNormal, sideOrd, parentCell);
348 
349  // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
350  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,1),
351  KOKKOS_LAMBDA (const int &) {
352  refSideNormalValueType tmp = refSideNormal(0);
353  refSideNormal(0) = refSideNormal(1);
354  refSideNormal(1) = -tmp;
355  });
356  } else {
357  // 3D parent cell: side = 2D subcell (face), call the face normal method.
358  getReferenceFaceNormal(refSideNormal, sideOrd, parentCell);
359  }
360  }
361 
362 
363  template<typename DeviceType>
364  template<typename refFaceNormalValueType, class ...refFaceNormalProperties>
365  void
366  CellTools<DeviceType>::
367  getReferenceFaceNormal( Kokkos::DynRankView<refFaceNormalValueType,refFaceNormalProperties...> refFaceNormal,
368  const ordinal_type faceOrd,
369  const shards::CellTopology parentCell ) {
370 #ifdef HAVE_INTREPID2_DEBUG
371  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
372  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): three-dimensional parent cell required");
373 
374  INTREPID2_TEST_FOR_EXCEPTION( faceOrd < 0 || faceOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(2)), std::invalid_argument,
375  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): face ordinal out of bounds");
376 
377  INTREPID2_TEST_FOR_EXCEPTION( refFaceNormal.rank() != 1, std::invalid_argument,
378  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): rank = 1 required for output array");
379 
380  INTREPID2_TEST_FOR_EXCEPTION( refFaceNormal.extent(0) != parentCell.getDimension(), std::invalid_argument,
381  ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): dim0 (spatial dim) must match that of parent cell");
382 #endif
383  constexpr bool is_accessible = Kokkos::Impl::MemorySpaceAccess<
384  MemSpaceType,
385  typename decltype(refFaceNormal)::memory_space>::accessible;
386  static_assert(is_accessible, "CellTools<DeviceType>::getReferenceFaceNormal(..): output view's memory space is not compatible with DeviceType");
387 
388  // Reference face normal = vector product of reference face tangents. Allocate temp FC storage:
389  const auto dim = parentCell.getDimension();
390  auto vcprop = Kokkos::common_view_alloc_prop(refFaceNormal);
391  using common_value_type = typename decltype(vcprop)::value_type;
392  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanU ( Kokkos::view_alloc("CellTools::getReferenceFaceNormal::refFaceTanU", vcprop), dim );
393  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanV ( Kokkos::view_alloc("CellTools::getReferenceFaceNormal::refFaceTanV", vcprop), dim );
394  getReferenceFaceTangents(refFaceTanU, refFaceTanV, faceOrd, parentCell);
395 
396  RealSpaceTools<DeviceType>::vecprod(refFaceNormal, refFaceTanU, refFaceTanV);
397  }
398 
399  template<typename DeviceType>
400  template<typename edgeTangentValueType, class ...edgeTangentProperties,
401  typename worksetJacobianValueType, class ...worksetJacobianProperties>
402  void
403  CellTools<DeviceType>::
404  getPhysicalEdgeTangents( Kokkos::DynRankView<edgeTangentValueType,edgeTangentProperties...> edgeTangents,
405  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
406  const ordinal_type worksetEdgeOrd,
407  const shards::CellTopology parentCell ) {
408 #ifdef HAVE_INTREPID2_DEBUG
409  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3 &&
410  parentCell.getDimension() != 2, std::invalid_argument,
411  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required." );
412 
413  // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
414  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.rank() != 3, std::invalid_argument,
415  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents requires rank 3." );
416  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(2) != 2 &&
417  edgeTangents.extent(2) != 3, std::invalid_argument,
418  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension(2) must be 2 or 3." );
419 
420  // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
421  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
422  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians requires rank 4." );
423  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 2 &&
424  worksetJacobians.extent(2) != 3, std::invalid_argument,
425  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) must be 2 or 3." );
426  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
427  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) and (3) must match each other." );
428 
429  // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
430  for (auto i=0;i<3;++i) {
431  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
432  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension (i) does not match to worksetJacobians dimension(i)." );
433  }
434 #endif
435  constexpr bool are_accessible =
436  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
437  typename decltype(edgeTangents)::memory_space>::accessible &&
438  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
439  typename decltype(worksetJacobians)::memory_space>::accessible;
440  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalEdgeTangents(..): input/output views' memory spaces are not compatible with DeviceType");
441 
442 
443  // Storage for constant reference edge tangent: rank-1 (D) arrays
444  const auto dim = parentCell.getDimension();
445  auto vcprop = Kokkos::common_view_alloc_prop(edgeTangents);
446  using common_value_type = typename decltype(vcprop)::value_type;
447  Kokkos::DynRankView< common_value_type, DeviceType > refEdgeTan ( Kokkos::view_alloc("CellTools::getPhysicalEdgeTangents::refEdgeTan", vcprop), dim);
448  getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
449 
450  RealSpaceTools<DeviceType>::matvec(edgeTangents, worksetJacobians, refEdgeTan);
451  }
452 
453 
454  namespace FunctorCellTools {
455 
456  template<typename tangentViewType,
457  typename faceOrdinalViewType,
458  typename parametrizationViewType
459  >
460  struct F_refEdgeTangent {
461  tangentViewType refEdgeTan_;
462  const faceOrdinalViewType edgeOrdView_;
463  const parametrizationViewType edgeParametrization_;
464 
465  KOKKOS_INLINE_FUNCTION
466  F_refEdgeTangent( tangentViewType refEdgeTan,
467  const faceOrdinalViewType edgeOrdView,
468  const parametrizationViewType edgeParametrization)
469  : refEdgeTan_(refEdgeTan), edgeOrdView_(edgeOrdView), edgeParametrization_(edgeParametrization){};
470 
471  KOKKOS_INLINE_FUNCTION
472  void operator()(const size_type ic) const {
473  for (size_type d=0; d<refEdgeTan_.extent(1); d++) {
474  refEdgeTan_(ic,d) = edgeParametrization_(edgeOrdView_(ic), d, 1);
475  }
476  }
477  };
478  }
479 
480  template<typename DeviceType>
481  template<typename edgeTangentValueType, class ...edgeTangentProperties,
482  typename worksetJacobianValueType, class ...worksetJacobianProperties,
483  typename edgeOrdValueType, class ...edgeOrdProperties>
484  void
485  CellTools<DeviceType>::
486  getPhysicalEdgeTangents( Kokkos::DynRankView<edgeTangentValueType,edgeTangentProperties...> edgeTangents,
487  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
488  const Kokkos::DynRankView<edgeOrdValueType,edgeOrdProperties...> worksetEdgeOrds,
489  const shards::CellTopology parentCell ) {
490 #ifdef HAVE_INTREPID2_DEBUG
491  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3 &&
492  parentCell.getDimension() != 2, std::invalid_argument,
493  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required." );
494 
495  // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
496  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.rank() != 3, std::invalid_argument,
497  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents requires rank 3." );
498  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(2) != 2 &&
499  edgeTangents.extent(2) != 3, std::invalid_argument,
500  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension(2) must be 2 or 3." );
501 
502  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(0) != worksetEdgeOrds.extent(0), std::invalid_argument,
503  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetEdgeOrds extent 0 should match that of edgeTangents." );
504 
505  // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
506  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
507  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians requires rank 4." );
508  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 2 &&
509  worksetJacobians.extent(2) != 3, std::invalid_argument,
510  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) must be 2 or 3." );
511  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
512  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) and (3) must match each other." );
513 
514  // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
515  for (auto i=0;i<3;++i) {
516  INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
517  ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension (i) does not match to worksetJacobians dimension(i)." );
518  }
519 #endif
520  constexpr bool are_accessible =
521  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
522  typename decltype(edgeTangents)::memory_space>::accessible &&
523  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
524  typename decltype(worksetJacobians)::memory_space>::accessible &&
525  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
526  typename decltype(worksetEdgeOrds)::memory_space>::accessible;
527  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalEdgeTangents(..): input/output views' memory spaces are not compatible with DeviceType");
528 
529 
530  // Storage for constant reference edge tangent: rank-1 (D) arrays
531  const ordinal_type dim = parentCell.getDimension();
532  auto vcprop = Kokkos::common_view_alloc_prop(edgeTangents);
533  using common_value_type = typename decltype(vcprop)::value_type;
534  Kokkos::DynRankView< common_value_type, DeviceType > refEdgeTan ( Kokkos::view_alloc("CellTools::getPhysicalEdgeTangents::refEdgeTan", vcprop), edgeTangents.extent(0), dim);
535 
536  const auto edgeMap = RefSubcellParametrization<DeviceType>::get(1, parentCell.getKey());
537 
538  using FunctorType = FunctorCellTools::F_refEdgeTangent<decltype(refEdgeTan),decltype(worksetEdgeOrds),decltype(edgeMap)>;
539  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,refEdgeTan.extent(0)), FunctorType(refEdgeTan, worksetEdgeOrds, edgeMap) );
540 
541  typename DeviceType::execution_space().fence();
542  RealSpaceTools<DeviceType>::matvec(edgeTangents, worksetJacobians, refEdgeTan);
543  }
544 
545  template<typename DeviceType>
546  template<typename faceTanValueType, class ...faceTanProperties,
547  typename worksetJacobianValueType, class ...worksetJacobianProperties>
548  void
549  CellTools<DeviceType>::
550  getPhysicalFaceTangents( Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanU,
551  Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanV,
552  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
553  const ordinal_type worksetFaceOrd,
554  const shards::CellTopology parentCell ) {
555 #ifdef HAVE_INTREPID2_DEBUG
556  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
557  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
558 
559  // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
560  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.rank() != 3 ||
561  faceTanV.rank() != 3, std::invalid_argument,
562  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V must have rank 3." );
563 
564  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(2) != 3 ||
565  faceTanV.extent(2) != 3, std::invalid_argument,
566  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (2) must be 3." );
567 
568  for (auto i=0;i<3;++i) {
569  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != faceTanV.extent(i), std::invalid_argument,
570  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (i) must match each other." );
571  }
572 
573  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
574  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
575  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians must have rank 4." );
576 
577  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
578  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) must be 3." );
579 
580  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
581  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) and dimension(3) must match." );
582 
583  // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
584  for (auto i=0;i<3;++i) {
585  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
586  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(i) and faceTan dimension (i) must match." );
587  }
588 #endif
589  constexpr bool are_accessible =
590  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
591  typename decltype(faceTanU)::memory_space>::accessible &&
592  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
593  typename decltype(worksetJacobians)::memory_space>::accessible;
594  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceTangents(..): input/output views' memory spaces are not compatible with DeviceType");
595 
596  // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
597  const auto dim = parentCell.getDimension();
598 
599  auto vcprop = Kokkos::common_view_alloc_prop(faceTanU);
600  using common_value_type = typename decltype(vcprop)::value_type;
601  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanU", vcprop), dim);
602  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanV", vcprop), dim);
603 
604  getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
605 
606  RealSpaceTools<DeviceType>::matvec(faceTanU, worksetJacobians, refFaceTanU);
607  RealSpaceTools<DeviceType>::matvec(faceTanV, worksetJacobians, refFaceTanV);
608  }
609 
610  namespace FunctorCellTools {
611 
612  template<typename tangentsViewType,
613  typename faceOrdinalViewType,
614  typename parametrizationViewType
615  >
616  struct F_refFaceTangents {
617  tangentsViewType refFaceTanU_;
618  tangentsViewType refFaceTanV_;
619  const faceOrdinalViewType faceOrdView_;
620  const parametrizationViewType faceParametrization_;
621 
622  KOKKOS_INLINE_FUNCTION
623  F_refFaceTangents( tangentsViewType refFaceTanU,
624  tangentsViewType refFaceTanV,
625  const faceOrdinalViewType faceOrdView,
626  const parametrizationViewType faceParametrization)
627  : refFaceTanU_(refFaceTanU), refFaceTanV_(refFaceTanV), faceOrdView_(faceOrdView), faceParametrization_(faceParametrization){};
628 
629  KOKKOS_INLINE_FUNCTION
630  void operator()(const size_type ic) const {
631  for (size_type d=0; d<refFaceTanU_.extent(1); d++) {
632  refFaceTanU_(ic,d) = faceParametrization_(faceOrdView_(ic), d, 1);
633  refFaceTanV_(ic,d) = faceParametrization_(faceOrdView_(ic), d, 2);
634  }
635  }
636  };
637  }
638 
639 
640  template<typename DeviceType>
641  template<typename faceTanValueType, class ...faceTanProperties,
642  typename worksetJacobianValueType, class ...worksetJacobianProperties,
643  typename faceOrdValueType, class ...faceOrdProperties>
644  void
645  CellTools<DeviceType>::
646  getPhysicalFaceTangents( Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanU,
647  Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanV,
648  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
649  const Kokkos::DynRankView<faceOrdValueType,faceOrdProperties...> worksetFaceOrds,
650  const shards::CellTopology parentCell ) {
651 #ifdef HAVE_INTREPID2_DEBUG
652  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
653  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
654 
655  // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
656  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.rank() != 3 ||
657  faceTanV.rank() != 3, std::invalid_argument,
658  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V must have rank 3." );
659 
660  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(2) != 3 ||
661  faceTanV.extent(2) != 3, std::invalid_argument,
662  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (2) must be 3." );
663 
664  for (auto i=0;i<3;++i) {
665  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != faceTanV.extent(i), std::invalid_argument,
666  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (i) must match each other." );
667  }
668 
669  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(0) != worksetFaceOrds.extent(0), std::invalid_argument,
670  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetFaceOrds extent 0 should match that of faceTanU." );
671 
672 
673  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
674  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
675  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians must have rank 4." );
676 
677  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
678  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) must be 3." );
679 
680  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
681  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) and dimension(3) must match." );
682 
683  // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
684  for (auto i=0;i<3;++i) {
685  INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
686  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(i) and faceTan dimension (i) must match." );
687  }
688 #endif
689  constexpr bool are_accessible =
690  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
691  typename decltype(faceTanU)::memory_space>::accessible &&
692  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
693  typename decltype(worksetJacobians)::memory_space>::accessible &&
694  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
695  typename decltype(worksetFaceOrds)::memory_space>::accessible;
696  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceTangents(..): input/output views' memory spaces are not compatible with DeviceType");
697 
698  // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
699  const ordinal_type dim = parentCell.getDimension();
700 
701  auto vcprop = Kokkos::common_view_alloc_prop(faceTanU);
702  using common_value_type = typename decltype(vcprop)::value_type;
703  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanU", vcprop), faceTanU.extent(0), dim);
704  Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanV", vcprop), faceTanV.extent(0), dim);
705 
706  const auto faceMap = RefSubcellParametrization<DeviceType>::get(2, parentCell.getKey());
707 
708  using FunctorType = FunctorCellTools::F_refFaceTangents<decltype(refFaceTanU),decltype(worksetFaceOrds),decltype(faceMap)>;
709  Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,refFaceTanU.extent(0)), FunctorType(refFaceTanU, refFaceTanV, worksetFaceOrds, faceMap) );
710 
711  typename DeviceType::execution_space().fence();
712  RealSpaceTools<DeviceType>::matvec(faceTanU, worksetJacobians, refFaceTanU);
713  RealSpaceTools<DeviceType>::matvec(faceTanV, worksetJacobians, refFaceTanV);
714  }
715 
716  namespace FunctorCellTools {
717 
718  template<typename normalsViewType,
719  typename tangentsViewType>
720  struct F_edgeNormalsFromTangents {
721  normalsViewType edgeNormals_;
722  const tangentsViewType edgeTangents_;
723 
724  KOKKOS_INLINE_FUNCTION
725  F_edgeNormalsFromTangents( normalsViewType edgeNormals,
726  const tangentsViewType refEdgeTangents)
727  : edgeNormals_(edgeNormals), edgeTangents_(refEdgeTangents){};
728 
729  KOKKOS_INLINE_FUNCTION
730  void operator()(const size_type iter) const {
731  size_type cell, pt;
732  unrollIndex( cell, pt, edgeTangents_.extent(0),
733  edgeTangents_.extent(1), iter );
734  edgeNormals_(cell,pt,0) = edgeTangents_(cell,pt,1);
735  edgeNormals_(cell,pt,1) = -edgeTangents_(cell,pt,0);
736  }
737  };
738  }
739 
740 
741 
742  template<typename DeviceType>
743  template<typename sideNormalValueType, class ...sideNormalProperties,
744  typename worksetJacobianValueType, class ...worksetJacobianProperties>
745  void
746  CellTools<DeviceType>::
747  getPhysicalSideNormals( Kokkos::DynRankView<sideNormalValueType,sideNormalProperties...> sideNormals,
748  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
749  const ordinal_type worksetSideOrd,
750  const shards::CellTopology parentCell ) {
751 #ifdef HAVE_INTREPID2_DEBUG
752  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
753  parentCell.getDimension() != 3, std::invalid_argument,
754  ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
755 
756  // Check side ordinal: by definition side is subcell whose dimension = parentCell.getDimension()-1
757  INTREPID2_TEST_FOR_EXCEPTION( worksetSideOrd < 0 ||
758  worksetSideOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(parentCell.getDimension() - 1)), std::invalid_argument,
759  ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): side ordinal out of bounds");
760 #endif
761  constexpr bool are_accessible =
762  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
763  typename decltype(sideNormals)::memory_space>::accessible &&
764  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
765  typename decltype(worksetJacobians)::memory_space>::accessible;
766  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalSideNormals(..): input/output views' memory spaces are not compatible with DeviceType");
767 
768  const auto dim = parentCell.getDimension();
769 
770  if (dim == 2) {
771  // compute edge tangents and rotate it
772  auto vcprop = Kokkos::common_view_alloc_prop(sideNormals);
773  using common_value_type = typename decltype(vcprop)::value_type;
774  Kokkos::DynRankView< common_value_type, DeviceType > edgeTangents ( Kokkos::view_alloc("CellTools::getPhysicalSideNormals::edgeTan", vcprop),
775  sideNormals.extent(0),
776  sideNormals.extent(1),
777  sideNormals.extent(2));
778  getPhysicalEdgeTangents(edgeTangents, worksetJacobians, worksetSideOrd, parentCell);
779 
780  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
781  using FunctorType = FunctorCellTools::F_edgeNormalsFromTangents<decltype(sideNormals), decltype(edgeTangents)>;
782  const auto loopSize = edgeTangents.extent(0)*edgeTangents.extent(1);
783  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
784  Kokkos::parallel_for( policy, FunctorType(sideNormals, edgeTangents) );
785 
786  } else {
787  getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
788  }
789  }
790 
791 
792  template<typename DeviceType>
793  template<typename sideNormalValueType, class ...sideNormalProperties,
794  typename worksetJacobianValueType, class ...worksetJacobianProperties,
795  typename edgeOrdValueType, class ...edgeOrdProperties>
796  void
797  CellTools<DeviceType>::
798  getPhysicalSideNormals( Kokkos::DynRankView<sideNormalValueType,sideNormalProperties...> sideNormals,
799  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
800  const Kokkos::DynRankView<edgeOrdValueType,edgeOrdProperties...> worksetSideOrds,
801  const shards::CellTopology parentCell ) {
802 #ifdef HAVE_INTREPID2_DEBUG
803  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
804  parentCell.getDimension() != 3, std::invalid_argument,
805  ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
806 #endif
807  constexpr bool are_accessible =
808  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
809  typename decltype(sideNormals)::memory_space>::accessible &&
810  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
811  typename decltype(worksetJacobians)::memory_space>::accessible &&
812  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
813  typename decltype(worksetSideOrds)::memory_space>::accessible;
814  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalSideNormals(..): input/output views' memory spaces are not compatible with DeviceType");
815 
816  const auto dim = parentCell.getDimension();
817 
818  if (dim == 2) {
819  // compute edge tangents and rotate it
820  auto vcprop = Kokkos::common_view_alloc_prop(sideNormals);
821  using common_value_type = typename decltype(vcprop)::value_type;
822  Kokkos::DynRankView< common_value_type, DeviceType > edgeTangents ( Kokkos::view_alloc("CellTools::getPhysicalSideNormals::edgeTan", vcprop),
823  sideNormals.extent(0),
824  sideNormals.extent(1),
825  sideNormals.extent(2));
826  getPhysicalEdgeTangents(edgeTangents, worksetJacobians, worksetSideOrds, parentCell);
827 
828  //Note: this function has several template parameters and the compiler gets confused if using a lambda function
829  using FunctorType = FunctorCellTools::F_edgeNormalsFromTangents<decltype(sideNormals), decltype(edgeTangents)>;
830  const auto loopSize = edgeTangents.extent(0)*edgeTangents.extent(1);
831  Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
832  Kokkos::parallel_for( policy, FunctorType(sideNormals, edgeTangents) );
833 
834  } else {
835  getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrds, parentCell);
836  }
837  }
838 
839 
840  template<typename DeviceType>
841  template<typename faceNormalValueType, class ...faceNormalProperties,
842  typename worksetJacobianValueType, class ...worksetJacobianProperties>
843  void
844  CellTools<DeviceType>::
845  getPhysicalFaceNormals( Kokkos::DynRankView<faceNormalValueType,faceNormalProperties...> faceNormals,
846  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
847  const ordinal_type worksetFaceOrd,
848  const shards::CellTopology parentCell ) {
849 #ifdef HAVE_INTREPID2_DEBUG
850  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
851  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required." );
852 
853  // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
854  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.rank() != 3, std::invalid_argument,
855  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals must have a rank 3." );
856  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(2) != 3, std::invalid_argument,
857  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (2) must be 3." );
858 
859  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
860  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
861  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians must have a rank 4." );
862  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
863  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must be 3." );
864  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
865  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must match to dimension (3)." );
866 
867  // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
868  for (auto i=0;i<3;++i) {
869  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
870  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (i) must match to worksetJacobians dimension (i)." );
871  }
872 #endif
873  constexpr bool are_accessible =
874  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
875  typename decltype(faceNormals)::memory_space>::accessible &&
876  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
877  typename decltype(worksetJacobians)::memory_space>::accessible;
878  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceNormals(..): input/output views' memory spaces are not compatible with DeviceType");
879 
880 
881  // this should be provided from users
882  // Storage for physical face tangents: rank-3 (C,P,D) arrays
883  const auto worksetSize = worksetJacobians.extent(0);
884  const auto facePtCount = worksetJacobians.extent(1);
885  const auto dim = parentCell.getDimension();
886 
887  auto vcprop = Kokkos::common_view_alloc_prop(faceNormals);
888  using common_value_type = typename decltype(vcprop)::value_type;
889  Kokkos::DynRankView< common_value_type, DeviceType > faceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanU", vcprop), worksetSize, facePtCount, dim);
890  Kokkos::DynRankView< common_value_type, DeviceType > faceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanV", vcprop), worksetSize, facePtCount, dim);
891 
892  getPhysicalFaceTangents(faceTanU, faceTanV,
893  worksetJacobians,
894  worksetFaceOrd,
895  parentCell);
896 
897  typename DeviceType::execution_space().fence();
898  RealSpaceTools<DeviceType>::vecprod(faceNormals, faceTanU, faceTanV);
899  }
900 
901  template<typename DeviceType>
902  template<typename faceNormalValueType, class ...faceNormalProperties,
903  typename worksetJacobianValueType, class ...worksetJacobianProperties,
904  typename faceOrdValueType, class ...faceOrdProperties>
905  void
906  CellTools<DeviceType>::
907  getPhysicalFaceNormals( Kokkos::DynRankView<faceNormalValueType,faceNormalProperties...> faceNormals,
908  const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
909  const Kokkos::DynRankView<faceOrdValueType,faceOrdProperties...> worksetFaceOrds,
910  const shards::CellTopology parentCell ) {
911 #ifdef HAVE_INTREPID2_DEBUG
912  INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
913  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required." );
914 
915  // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
916  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.rank() != 3, std::invalid_argument,
917  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals must have a rank 3." );
918  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(2) != 3, std::invalid_argument,
919  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (2) must be 3." );
920  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(0) != worksetFaceOrds.extent(0), std::invalid_argument,
921  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetFaceOrds extent 0 should match that of faceNormals." );
922 
923  // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
924  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
925  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians must have a rank 4." );
926  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
927  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must be 3." );
928  INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
929  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must match to dimension (3)." );
930 
931  // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
932  for (auto i=0;i<3;++i) {
933  INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
934  ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (i) must match to worksetJacobians dimension (i)." );
935  }
936 #endif
937  constexpr bool are_accessible =
938  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
939  typename decltype(faceNormals)::memory_space>::accessible &&
940  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
941  typename decltype(worksetJacobians)::memory_space>::accessible &&
942  Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
943  typename decltype(worksetFaceOrds)::memory_space>::accessible;
944  static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceNormals(..): input/output views' memory spaces are not compatible with DeviceType");
945 
946 
947  // this should be provided from users
948  // Storage for physical face tangents: rank-3 (C,P,D) arrays
949  const auto worksetSize = worksetJacobians.extent(0);
950  const auto facePtCount = worksetJacobians.extent(1);
951  const auto dim = parentCell.getDimension();
952 
953  auto vcprop = Kokkos::common_view_alloc_prop(faceNormals);
954  using common_value_type = typename decltype(vcprop)::value_type;
955  Kokkos::DynRankView< common_value_type, DeviceType > faceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanU", vcprop), worksetSize, facePtCount, dim);
956  Kokkos::DynRankView< common_value_type, DeviceType > faceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanV", vcprop), worksetSize, facePtCount, dim);
957 
958  getPhysicalFaceTangents(faceTanU, faceTanV,
959  worksetJacobians,
960  worksetFaceOrds,
961  parentCell);
962 
963  typename DeviceType::execution_space().fence();
964  RealSpaceTools<DeviceType>::vecprod(faceNormals, faceTanU, faceTanV);
965  }
966 }
967 
968 #endif
Intrepid2::CellTools::getReferenceCellCenter
static void getReferenceCellCenter(Kokkos::DynRankView< cellCenterValueType, cellCenterProperties... > cellCenter, const shards::CellTopology cell)
Computes the Cartesian coordinates of reference cell barycenter.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:69
Intrepid2::CellTools::getPhysicalFaceTangents
static void getPhysicalFaceTangents(Kokkos::DynRankView< faceTanValueType, faceTanProperties... > faceTanU, Kokkos::DynRankView< faceTanValueType, faceTanProperties... > faceTanV, const Kokkos::DynRankView< worksetJacobianValueType, worksetJacobianProperties... > worksetJacobians, const ordinal_type worksetFaceOrd, const shards::CellTopology parentCell)
Computes non-normalized tangent vector pairs to physical faces in a face workset ; (see Subcell works...
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:550
Intrepid2::RefCellNodes::get
static ConstViewType get(const unsigned cellTopoKey)
Retrieves the Cartesian coordinates of reference cell nodes.
Definition: Intrepid2_CellDataDef.hpp:520
Intrepid2::RefSubcellParametrization::get
static ConstViewType get(const ordinal_type subcellDim, const unsigned parentCellKey)
Returns a Kokkos view with the coefficients of the parametrization maps for the edges or faces of a r...
Definition: Intrepid2_CellDataDef.hpp:103
Intrepid2::CellTools::getReferenceEdgeTangent
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties... > refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:244
Intrepid2::RealSpaceTools::matvec
static void matvec(Kokkos::DynRankView< matVecValueType, matVecProperties... > matVecs, const Kokkos::DynRankView< inMatValueType, inMatProperties... > inMats, const Kokkos::DynRankView< inVecValueType, inVecProperties... > inVecs)
Matrix-vector left multiply using multidimensional arrays: matVec = inMat * inVec.
Definition: Intrepid2_RealSpaceToolsDef.hpp:1389
Intrepid2::CellTools::getReferenceNode
static void getReferenceNode(Kokkos::DynRankView< cellNodeValueType, cellNodeProperties... > cellNode, const shards::CellTopology cell, const ordinal_type nodeOrd)
Retrieves the Cartesian coordinates of a reference cell node.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:172
Intrepid2::CellTools::getReferenceFaceTangents
static void getReferenceFaceTangents(Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanU, Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:282
Intrepid2::FunctorCellTools::F_edgeNormalsFromTangents
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:720
Intrepid2::FunctorCellTools::F_refFaceTangents
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:616
Intrepid2::CellTools::getPhysicalFaceNormals
static void getPhysicalFaceNormals(Kokkos::DynRankView< faceNormalValueType, faceNormalProperties... > faceNormals, const Kokkos::DynRankView< worksetJacobianValueType, worksetJacobianProperties... > worksetJacobians, const ordinal_type worksetFaceOrd, const shards::CellTopology parentCell)
Computes non-normalized normal vectors to physical faces in a face workset ; (see Subcell worksets fo...
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:845
Intrepid2::CellTools::getReferenceSubcellNodes
static void getReferenceSubcellNodes(Kokkos::DynRankView< subcellNodeValueType, subcellNodeProperties... > subcellNodes, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Retrieves the Cartesian coordinates of all nodes of a reference subcell.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:202
Intrepid2::CellTools::getReferenceSubcellVertices
static void getReferenceSubcellVertices(Kokkos::DynRankView< subcellVertexValueType, subcellVertexProperties... > subcellVertices, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Retrieves the Cartesian coordinates of all vertices of a reference subcell.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:136
Intrepid2::CellTools::getReferenceSideNormal
static void getReferenceSideNormal(Kokkos::DynRankView< refSideNormalValueType, refSideNormalProperties... > refSideNormal, const ordinal_type sideOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:327
Intrepid2::CellTools::getPhysicalEdgeTangents
static void getPhysicalEdgeTangents(Kokkos::DynRankView< edgeTangentValueType, edgeTangentProperties... > edgeTangents, const Kokkos::DynRankView< worksetJacobianValueType, worksetJacobianProperties... > worksetJacobians, const ordinal_type worksetEdgeOrd, const shards::CellTopology parentCell)
Computes non-normalized tangent vectors to physical edges in an edge workset ; (see Subcell worksets ...
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:404
Intrepid2::RefCellCenter::get
static ConstViewType get(const unsigned cellTopoKey)
Retrieves the Cartesian coordinates of a reference cell barycenter.
Definition: Intrepid2_CellDataDef.hpp:768
Intrepid2::CellTools::getReferenceVertex
static void getReferenceVertex(Kokkos::DynRankView< cellVertexValueType, cellVertexProperties... > cellVertex, const shards::CellTopology cell, const ordinal_type vertexOrd)
Retrieves the Cartesian coordinates of a reference cell vertex.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:101
Intrepid2::RealSpaceTools::vecprod
static void vecprod(Kokkos::DynRankView< vecProdValueType, vecProdProperties... > vecProd, const Kokkos::DynRankView< inLeftValueType, inLeftProperties... > inLeft, const Kokkos::DynRankView< inRightValueType, inRightProperties... > inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
Definition: Intrepid2_RealSpaceToolsDef.hpp:1573
Intrepid2::CellTools::getPhysicalSideNormals
static void getPhysicalSideNormals(Kokkos::DynRankView< sideNormalValueType, sideNormalProperties... > sideNormals, const Kokkos::DynRankView< worksetJacobianValueType, worksetJacobianProperties... > worksetJacobians, const ordinal_type worksetSideOrd, const shards::CellTopology parentCell)
Computes non-normalized normal vectors to physical sides in a side workset .
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:747
Intrepid2::CellTools::getReferenceFaceNormal
static void getReferenceFaceNormal(Kokkos::DynRankView< refFaceNormalValueType, refFaceNormalProperties... > refFaceNormal, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to faces of 3D reference cell.
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:367
Intrepid2::FunctorCellTools::F_refEdgeTangent
Definition: Intrepid2_CellToolsDefNodeInfo.hpp:460
Generated by   doxygen 1.8.14