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