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 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(edgeTangents)::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(faceTanU)::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(sideNormals)::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(faceNormals)::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::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:518
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::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:372
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:715
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::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:1546
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:688
Intrepid2::FunctorCellTools::F_refEdgeTangent
Definition Intrepid2_CellToolsDefNodeInfo.hpp:428
Intrepid2::FunctorCellTools::F_refFaceTangents
Definition Intrepid2_CellToolsDefNodeInfo.hpp:584
Generated by doxygen 1.9.8