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
24namespace 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 using ViewType = Kokkos::DynRankView< typename decltype(refFaceNormal)::value_type, DeviceType >;
359 ViewType refFaceTanU = Impl::createMatchingView<ViewType>(refFaceNormal,"CellTools::getReferenceFaceNormal::refFaceTanU", dim);
360 ViewType refFaceTanV = Impl::createMatchingView<ViewType>(refFaceNormal,"CellTools::getReferenceFaceNormal::refFaceTanV", dim);
361 getReferenceFaceTangents(refFaceTanU, refFaceTanV, faceOrd, parentCell);
362
363 RealSpaceTools<DeviceType>::vecprod(refFaceNormal, refFaceTanU, refFaceTanV);
364 }
365
366 template<typename DeviceType>
367 template<typename edgeTangentValueType, class ...edgeTangentProperties,
368 typename worksetJacobianValueType, class ...worksetJacobianProperties>
369 void
370 CellTools<DeviceType>::
371 getPhysicalEdgeTangents( Kokkos::DynRankView<edgeTangentValueType,edgeTangentProperties...> edgeTangents,
372 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
373 const ordinal_type worksetEdgeOrd,
374 const shards::CellTopology parentCell ) {
375#ifdef HAVE_INTREPID2_DEBUG
376 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3 &&
377 parentCell.getDimension() != 2, std::invalid_argument,
378 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required." );
379
380 // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
381 INTREPID2_TEST_FOR_EXCEPTION( rank(edgeTangents) != 3, std::invalid_argument,
382 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents requires rank 3." );
383 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(2) != 2 &&
384 edgeTangents.extent(2) != 3, std::invalid_argument,
385 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension(2) must be 2 or 3." );
386
387 // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
388 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
389 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians requires rank 4." );
390 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 2 &&
391 worksetJacobians.extent(2) != 3, std::invalid_argument,
392 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) must be 2 or 3." );
393 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
394 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) and (3) must match each other." );
395
396 // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
397 for (auto i=0;i<3;++i) {
398 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
399 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension (i) does not match to worksetJacobians dimension(i)." );
400 }
401#endif
402 constexpr bool are_accessible =
403 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
404 typename decltype(edgeTangents)::memory_space>::accessible &&
405 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
406 typename decltype(worksetJacobians)::memory_space>::accessible;
407 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalEdgeTangents(..): input/output views' memory spaces are not compatible with DeviceType");
408
409
410 // Storage for constant reference edge tangent: rank-1 (D) arrays
411 const auto dim = parentCell.getDimension();
412 using ViewType = Kokkos::DynRankView< typename decltype(edgeTangents)::value_type, DeviceType >;
413 ViewType refEdgeTan = Impl::createMatchingView<ViewType>(edgeTangents,"CellTools::getPhysicalEdgeTangents::refEdgeTan", dim);
414 getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
415
416 RealSpaceTools<DeviceType>::matvec(edgeTangents, worksetJacobians, refEdgeTan);
417 }
418
419
420 namespace FunctorCellTools {
421
422 template<typename tangentViewType,
423 typename faceOrdinalViewType,
424 typename parametrizationViewType
425 >
426 struct F_refEdgeTangent {
427 tangentViewType refEdgeTan_;
428 const faceOrdinalViewType edgeOrdView_;
429 const parametrizationViewType edgeParametrization_;
430
431 KOKKOS_INLINE_FUNCTION
432 F_refEdgeTangent( tangentViewType refEdgeTan,
433 const faceOrdinalViewType edgeOrdView,
434 const parametrizationViewType edgeParametrization)
435 : refEdgeTan_(refEdgeTan), edgeOrdView_(edgeOrdView), edgeParametrization_(edgeParametrization){};
436
437 KOKKOS_INLINE_FUNCTION
438 void operator()(const size_type ic) const {
439 for (size_type d=0; d<refEdgeTan_.extent(1); d++) {
440 refEdgeTan_(ic,d) = edgeParametrization_(edgeOrdView_(ic), d, 1);
441 }
442 }
443 };
444 }
445
446 template<typename DeviceType>
447 template<typename edgeTangentValueType, class ...edgeTangentProperties,
448 typename worksetJacobianValueType, class ...worksetJacobianProperties,
449 typename edgeOrdValueType, class ...edgeOrdProperties>
450 void
451 CellTools<DeviceType>::
452 getPhysicalEdgeTangents( Kokkos::DynRankView<edgeTangentValueType,edgeTangentProperties...> edgeTangents,
453 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
454 const Kokkos::DynRankView<edgeOrdValueType,edgeOrdProperties...> worksetEdgeOrds,
455 const shards::CellTopology parentCell ) {
456#ifdef HAVE_INTREPID2_DEBUG
457 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3 &&
458 parentCell.getDimension() != 2, std::invalid_argument,
459 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required." );
460
461 // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
462 INTREPID2_TEST_FOR_EXCEPTION( rank(edgeTangents) != 3, std::invalid_argument,
463 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents requires rank 3." );
464 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(2) != 2 &&
465 edgeTangents.extent(2) != 3, std::invalid_argument,
466 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension(2) must be 2 or 3." );
467
468 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(0) != worksetEdgeOrds.extent(0), std::invalid_argument,
469 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetEdgeOrds extent 0 should match that of edgeTangents." );
470
471 // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
472 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
473 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians requires rank 4." );
474 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 2 &&
475 worksetJacobians.extent(2) != 3, std::invalid_argument,
476 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) must be 2 or 3." );
477 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
478 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) and (3) must match each other." );
479
480 // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
481 for (auto i=0;i<3;++i) {
482 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
483 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension (i) does not match to worksetJacobians dimension(i)." );
484 }
485#endif
486 constexpr bool are_accessible =
487 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
488 typename decltype(edgeTangents)::memory_space>::accessible &&
489 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
490 typename decltype(worksetJacobians)::memory_space>::accessible &&
491 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
492 typename decltype(worksetEdgeOrds)::memory_space>::accessible;
493 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalEdgeTangents(..): input/output views' memory spaces are not compatible with DeviceType");
494
495
496 // Storage for constant reference edge tangent: rank-1 (D) arrays
497 const ordinal_type dim = parentCell.getDimension();
498
499 using ViewType = Kokkos::DynRankView< typename decltype(edgeTangents)::value_type, DeviceType >;
500 ViewType refEdgeTan = Impl::createMatchingView<ViewType>(edgeTangents,"CellTools::getPhysicalEdgeTangents::refEdgeTan", edgeTangents.extent(0), dim);
501
502 const auto edgeMap = RefSubcellParametrization<DeviceType>::get(1, parentCell.getKey());
503
504 using FunctorType = FunctorCellTools::F_refEdgeTangent<decltype(refEdgeTan),decltype(worksetEdgeOrds),decltype(edgeMap)>;
505 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,refEdgeTan.extent(0)), FunctorType(refEdgeTan, worksetEdgeOrds, edgeMap) );
506
507 typename DeviceType::execution_space().fence();
508 RealSpaceTools<DeviceType>::matvec(edgeTangents, worksetJacobians, refEdgeTan);
509 }
510
511 template<typename DeviceType>
512 template<typename faceTanValueType, class ...faceTanProperties,
513 typename worksetJacobianValueType, class ...worksetJacobianProperties>
514 void
515 CellTools<DeviceType>::
516 getPhysicalFaceTangents( Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanU,
517 Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanV,
518 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
519 const ordinal_type worksetFaceOrd,
520 const shards::CellTopology parentCell ) {
521#ifdef HAVE_INTREPID2_DEBUG
522 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
523 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
524
525 // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
526 INTREPID2_TEST_FOR_EXCEPTION( rank(faceTanU) != 3 ||
527 rank(faceTanV) != 3, std::invalid_argument,
528 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V must have rank 3." );
529
530 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(2) != 3 ||
531 faceTanV.extent(2) != 3, std::invalid_argument,
532 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (2) must be 3." );
533
534 for (auto i=0;i<3;++i) {
535 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != faceTanV.extent(i), std::invalid_argument,
536 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (i) must match each other." );
537 }
538
539 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
540 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
541 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians must have rank 4." );
542
543 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
544 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) must be 3." );
545
546 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
547 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) and dimension(3) must match." );
548
549 // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
550 for (auto i=0;i<3;++i) {
551 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
552 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(i) and faceTan dimension (i) must match." );
553 }
554#endif
555 constexpr bool are_accessible =
556 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
557 typename decltype(faceTanU)::memory_space>::accessible &&
558 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
559 typename decltype(worksetJacobians)::memory_space>::accessible;
560 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceTangents(..): input/output views' memory spaces are not compatible with DeviceType");
561
562 // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
563 const auto dim = parentCell.getDimension();
564
565 using ViewType = Kokkos::DynRankView< typename decltype(faceTanU)::value_type, DeviceType >;
566 ViewType refFaceTanU = Impl::createMatchingView<ViewType>(faceTanU,"CellTools::getPhysicalFaceTangents::refFaceTanU", dim);
567 ViewType refFaceTanV = Impl::createMatchingView<ViewType>(faceTanV,"CellTools::getPhysicalFaceTangents::refFaceTanV", dim);
568
569 getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
570
571 RealSpaceTools<DeviceType>::matvec(faceTanU, worksetJacobians, refFaceTanU);
572 RealSpaceTools<DeviceType>::matvec(faceTanV, worksetJacobians, refFaceTanV);
573 }
574
575 namespace FunctorCellTools {
576
577 template<typename tangentsViewType,
578 typename faceOrdinalViewType,
579 typename parametrizationViewType
580 >
581 struct F_refFaceTangents {
582 tangentsViewType refFaceTanU_;
583 tangentsViewType refFaceTanV_;
584 const faceOrdinalViewType faceOrdView_;
585 const parametrizationViewType faceParametrization_;
586
587 KOKKOS_INLINE_FUNCTION
588 F_refFaceTangents( tangentsViewType refFaceTanU,
589 tangentsViewType refFaceTanV,
590 const faceOrdinalViewType faceOrdView,
591 const parametrizationViewType faceParametrization)
592 : refFaceTanU_(refFaceTanU), refFaceTanV_(refFaceTanV), faceOrdView_(faceOrdView), faceParametrization_(faceParametrization){};
593
594 KOKKOS_INLINE_FUNCTION
595 void operator()(const size_type ic) const {
596 for (size_type d=0; d<refFaceTanU_.extent(1); d++) {
597 refFaceTanU_(ic,d) = faceParametrization_(faceOrdView_(ic), d, 1);
598 refFaceTanV_(ic,d) = faceParametrization_(faceOrdView_(ic), d, 2);
599 }
600 }
601 };
602 }
603
604
605 template<typename DeviceType>
606 template<typename faceTanValueType, class ...faceTanProperties,
607 typename worksetJacobianValueType, class ...worksetJacobianProperties,
608 typename faceOrdValueType, class ...faceOrdProperties>
609 void
610 CellTools<DeviceType>::
611 getPhysicalFaceTangents( Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanU,
612 Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanV,
613 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
614 const Kokkos::DynRankView<faceOrdValueType,faceOrdProperties...> worksetFaceOrds,
615 const shards::CellTopology parentCell ) {
616#ifdef HAVE_INTREPID2_DEBUG
617 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
618 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
619
620 // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
621 INTREPID2_TEST_FOR_EXCEPTION( rank(faceTanU) != 3 ||
622 rank(faceTanV) != 3, std::invalid_argument,
623 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V must have rank 3." );
624
625 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(2) != 3 ||
626 faceTanV.extent(2) != 3, std::invalid_argument,
627 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (2) must be 3." );
628
629 for (auto i=0;i<3;++i) {
630 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != faceTanV.extent(i), std::invalid_argument,
631 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (i) must match each other." );
632 }
633
634 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(0) != worksetFaceOrds.extent(0), std::invalid_argument,
635 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetFaceOrds extent 0 should match that of faceTanU." );
636
637
638 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
639 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
640 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians must have rank 4." );
641
642 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
643 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) must be 3." );
644
645 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
646 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) and dimension(3) must match." );
647
648 // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
649 for (auto i=0;i<3;++i) {
650 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
651 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(i) and faceTan dimension (i) must match." );
652 }
653#endif
654 constexpr bool are_accessible =
655 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
656 typename decltype(faceTanU)::memory_space>::accessible &&
657 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
658 typename decltype(worksetJacobians)::memory_space>::accessible &&
659 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
660 typename decltype(worksetFaceOrds)::memory_space>::accessible;
661 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceTangents(..): input/output views' memory spaces are not compatible with DeviceType");
662
663 // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
664 const ordinal_type dim = parentCell.getDimension();
665
666 using ViewType = Kokkos::DynRankView< typename decltype(faceTanU)::value_type, DeviceType >;
667 ViewType refFaceTanU = Impl::createMatchingView<ViewType>(faceTanU,"CellTools::getPhysicalFaceTangents::refFaceTanU", faceTanU.extent(0), dim);
668 ViewType refFaceTanV = Impl::createMatchingView<ViewType>(faceTanV,"CellTools::getPhysicalFaceTangents::refFaceTanV", faceTanV.extent(0), dim);
669
670
671 const auto faceMap = RefSubcellParametrization<DeviceType>::get(2, parentCell.getKey());
672
673 using FunctorType = FunctorCellTools::F_refFaceTangents<decltype(refFaceTanU),decltype(worksetFaceOrds),decltype(faceMap)>;
674 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,refFaceTanU.extent(0)), FunctorType(refFaceTanU, refFaceTanV, worksetFaceOrds, faceMap) );
675
676 typename DeviceType::execution_space().fence();
677 RealSpaceTools<DeviceType>::matvec(faceTanU, worksetJacobians, refFaceTanU);
678 RealSpaceTools<DeviceType>::matvec(faceTanV, worksetJacobians, refFaceTanV);
679 }
680
681 namespace FunctorCellTools {
682
683 template<typename normalsViewType,
684 typename tangentsViewType>
685 struct F_edgeNormalsFromTangents {
686 normalsViewType edgeNormals_;
687 const tangentsViewType edgeTangents_;
688
689 KOKKOS_INLINE_FUNCTION
690 F_edgeNormalsFromTangents( normalsViewType edgeNormals,
691 const tangentsViewType refEdgeTangents)
692 : edgeNormals_(edgeNormals), edgeTangents_(refEdgeTangents){};
693
694 KOKKOS_INLINE_FUNCTION
695 void operator()(const size_type iter) const {
696 size_type cell, pt;
697 unrollIndex( cell, pt, edgeTangents_.extent(0),
698 edgeTangents_.extent(1), iter );
699 edgeNormals_(cell,pt,0) = edgeTangents_(cell,pt,1);
700 edgeNormals_(cell,pt,1) = -edgeTangents_(cell,pt,0);
701 }
702 };
703 }
704
705
706
707 template<typename DeviceType>
708 template<typename sideNormalValueType, class ...sideNormalProperties,
709 typename worksetJacobianValueType, class ...worksetJacobianProperties>
710 void
711 CellTools<DeviceType>::
712 getPhysicalSideNormals( Kokkos::DynRankView<sideNormalValueType,sideNormalProperties...> sideNormals,
713 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
714 const ordinal_type worksetSideOrd,
715 const shards::CellTopology parentCell ) {
716#ifdef HAVE_INTREPID2_DEBUG
717 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
718 parentCell.getDimension() != 3, std::invalid_argument,
719 ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
720
721 // Check side ordinal: by definition side is subcell whose dimension = parentCell.getDimension()-1
722 INTREPID2_TEST_FOR_EXCEPTION( worksetSideOrd < 0 ||
723 worksetSideOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(parentCell.getDimension() - 1)), std::invalid_argument,
724 ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): side ordinal out of bounds");
725#endif
726 constexpr bool are_accessible =
727 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
728 typename decltype(sideNormals)::memory_space>::accessible &&
729 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
730 typename decltype(worksetJacobians)::memory_space>::accessible;
731 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalSideNormals(..): input/output views' memory spaces are not compatible with DeviceType");
732
733 const auto dim = parentCell.getDimension();
734
735 if (dim == 2) {
736 // compute edge tangents and rotate it
737 using ViewType = Kokkos::DynRankView< typename decltype(sideNormals)::value_type, DeviceType >;
738 ViewType edgeTangents = Impl::createMatchingView<ViewType>(sideNormals,"CellTools::getPhysicalSideNormals::edgeTan", sideNormals.extent(0), sideNormals.extent(1), sideNormals.extent(2));
739 getPhysicalEdgeTangents(edgeTangents, worksetJacobians, worksetSideOrd, parentCell);
740
741 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
742 using FunctorType = FunctorCellTools::F_edgeNormalsFromTangents<decltype(sideNormals), decltype(edgeTangents)>;
743 const auto loopSize = edgeTangents.extent(0)*edgeTangents.extent(1);
744 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
745 Kokkos::parallel_for( policy, FunctorType(sideNormals, edgeTangents) );
746
747 } else {
748 getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
749 }
750 }
751
752
753 template<typename DeviceType>
754 template<typename sideNormalValueType, class ...sideNormalProperties,
755 typename worksetJacobianValueType, class ...worksetJacobianProperties,
756 typename edgeOrdValueType, class ...edgeOrdProperties>
757 void
758 CellTools<DeviceType>::
759 getPhysicalSideNormals( Kokkos::DynRankView<sideNormalValueType,sideNormalProperties...> sideNormals,
760 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
761 const Kokkos::DynRankView<edgeOrdValueType,edgeOrdProperties...> worksetSideOrds,
762 const shards::CellTopology parentCell ) {
763#ifdef HAVE_INTREPID2_DEBUG
764 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
765 parentCell.getDimension() != 3, std::invalid_argument,
766 ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
767#endif
768 constexpr bool are_accessible =
769 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
770 typename decltype(sideNormals)::memory_space>::accessible &&
771 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
772 typename decltype(worksetJacobians)::memory_space>::accessible &&
773 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
774 typename decltype(worksetSideOrds)::memory_space>::accessible;
775 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalSideNormals(..): input/output views' memory spaces are not compatible with DeviceType");
776
777 const auto dim = parentCell.getDimension();
778
779 if (dim == 2) {
780 using ViewType = Kokkos::DynRankView< typename decltype(sideNormals)::value_type, DeviceType >;
781 ViewType edgeTangents = Impl::createMatchingView<ViewType>(sideNormals,"CellTools::getPhysicalSideNormals::edgeTan", sideNormals.extent(0), sideNormals.extent(1), sideNormals.extent(2));
782 getPhysicalEdgeTangents(edgeTangents, worksetJacobians, worksetSideOrds, parentCell);
783
784 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
785 using FunctorType = FunctorCellTools::F_edgeNormalsFromTangents<decltype(sideNormals), decltype(edgeTangents)>;
786 const auto loopSize = edgeTangents.extent(0)*edgeTangents.extent(1);
787 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
788 Kokkos::parallel_for( policy, FunctorType(sideNormals, edgeTangents) );
789
790 } else {
791 getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrds, parentCell);
792 }
793 }
794
795
796 template<typename DeviceType>
797 template<typename faceNormalValueType, class ...faceNormalProperties,
798 typename worksetJacobianValueType, class ...worksetJacobianProperties>
799 void
800 CellTools<DeviceType>::
801 getPhysicalFaceNormals( Kokkos::DynRankView<faceNormalValueType,faceNormalProperties...> faceNormals,
802 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
803 const ordinal_type worksetFaceOrd,
804 const shards::CellTopology parentCell ) {
805#ifdef HAVE_INTREPID2_DEBUG
806 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
807 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required." );
808
809 // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
810 INTREPID2_TEST_FOR_EXCEPTION( rank(faceNormals) != 3, std::invalid_argument,
811 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals must have a rank 3." );
812 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(2) != 3, std::invalid_argument,
813 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (2) must be 3." );
814
815 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
816 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
817 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians must have a rank 4." );
818 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
819 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must be 3." );
820 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
821 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must match to dimension (3)." );
822
823 // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
824 for (auto i=0;i<3;++i) {
825 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
826 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (i) must match to worksetJacobians dimension (i)." );
827 }
828#endif
829 constexpr bool are_accessible =
830 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
831 typename decltype(faceNormals)::memory_space>::accessible &&
832 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
833 typename decltype(worksetJacobians)::memory_space>::accessible;
834 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceNormals(..): input/output views' memory spaces are not compatible with DeviceType");
835
836
837 // this should be provided from users
838 // Storage for physical face tangents: rank-3 (C,P,D) arrays
839 const auto worksetSize = worksetJacobians.extent(0);
840 const auto facePtCount = worksetJacobians.extent(1);
841 const auto dim = parentCell.getDimension();
842
843 using ViewType = Kokkos::DynRankView< typename decltype(faceNormals)::value_type, DeviceType >;
844 ViewType faceTanU = Impl::createMatchingView<ViewType>(faceNormals, "CellTools::getPhysicalFaceNormals::faceTanU", worksetSize, facePtCount, dim);
845 ViewType faceTanV = Impl::createMatchingView<ViewType>(faceNormals, "CellTools::getPhysicalFaceNormals::faceTanV", worksetSize, facePtCount, dim);
846
847 getPhysicalFaceTangents(faceTanU, faceTanV,
848 worksetJacobians,
849 worksetFaceOrd,
850 parentCell);
851
852 typename DeviceType::execution_space().fence();
853 RealSpaceTools<DeviceType>::vecprod(faceNormals, faceTanU, faceTanV);
854 }
855
856 template<typename DeviceType>
857 template<typename faceNormalValueType, class ...faceNormalProperties,
858 typename worksetJacobianValueType, class ...worksetJacobianProperties,
859 typename faceOrdValueType, class ...faceOrdProperties>
860 void
861 CellTools<DeviceType>::
862 getPhysicalFaceNormals( Kokkos::DynRankView<faceNormalValueType,faceNormalProperties...> faceNormals,
863 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
864 const Kokkos::DynRankView<faceOrdValueType,faceOrdProperties...> worksetFaceOrds,
865 const shards::CellTopology parentCell ) {
866#ifdef HAVE_INTREPID2_DEBUG
867 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
868 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required." );
869
870 // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
871 INTREPID2_TEST_FOR_EXCEPTION( rank(faceNormals) != 3, std::invalid_argument,
872 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals must have a rank 3." );
873 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(2) != 3, std::invalid_argument,
874 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (2) must be 3." );
875 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(0) != worksetFaceOrds.extent(0), std::invalid_argument,
876 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetFaceOrds extent 0 should match that of faceNormals." );
877
878 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
879 INTREPID2_TEST_FOR_EXCEPTION( rank(worksetJacobians) != 4, std::invalid_argument,
880 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians must have a rank 4." );
881 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
882 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must be 3." );
883 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
884 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must match to dimension (3)." );
885
886 // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
887 for (auto i=0;i<3;++i) {
888 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
889 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (i) must match to worksetJacobians dimension (i)." );
890 }
891#endif
892 constexpr bool are_accessible =
893 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
894 typename decltype(faceNormals)::memory_space>::accessible &&
895 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
896 typename decltype(worksetJacobians)::memory_space>::accessible &&
897 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
898 typename decltype(worksetFaceOrds)::memory_space>::accessible;
899 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceNormals(..): input/output views' memory spaces are not compatible with DeviceType");
900
901
902 // this should be provided from users
903 // Storage for physical face tangents: rank-3 (C,P,D) arrays
904 const auto worksetSize = worksetJacobians.extent(0);
905 const auto facePtCount = worksetJacobians.extent(1);
906 const auto dim = parentCell.getDimension();
907
908 using ViewType = Kokkos::DynRankView< typename decltype(faceNormals)::value_type, DeviceType >;
909 ViewType faceTanU = Impl::createMatchingView<ViewType>(faceNormals,"CellTools::getPhysicalFaceNormals::faceTanU", worksetSize, facePtCount, dim);
910 ViewType faceTanV = Impl::createMatchingView<ViewType>(faceNormals,"CellTools::getPhysicalFaceNormals::faceTanV", worksetSize, facePtCount, dim);
911
912 getPhysicalFaceTangents(faceTanU, faceTanV,
913 worksetJacobians,
914 worksetFaceOrds,
915 parentCell);
916
917 typename DeviceType::execution_space().fence();
918 RealSpaceTools<DeviceType>::vecprod(faceNormals, faceTanU, faceTanV);
919 }
920}
921
922#endif
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::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:516
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:801
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::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:371
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::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::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::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::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::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::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
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:712
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:1406
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:1596
Intrepid2::RefCellCenter::get
static ConstViewType get(const unsigned cellTopoKey)
Retrieves the Cartesian coordinates of a reference cell barycenter.
Definition Intrepid2_CellDataDef.hpp:734
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::FunctorCellTools::F_edgeNormalsFromTangents
Definition Intrepid2_CellToolsDefNodeInfo.hpp:685
Intrepid2::FunctorCellTools::F_refEdgeTangent
Definition Intrepid2_CellToolsDefNodeInfo.hpp:426
Intrepid2::FunctorCellTools::F_refFaceTangents
Definition Intrepid2_CellToolsDefNodeInfo.hpp:581
Generated by doxygen 1.9.8