Intrepid2
Intrepid2_CellToolsDefNodeInfo.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid2 Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Kyungjoo Kim (kyukim@sandia.gov), or
38// Mauro Perego (mperego@sandia.gov)
39//
40// ************************************************************************
41// @HEADER
42
43
49#ifndef __INTREPID2_CELLTOOLS_DEF_NODE_INFO_HPP__
50#define __INTREPID2_CELLTOOLS_DEF_NODE_INFO_HPP__
51
52// disable clang warnings
53#if defined (__clang__) && !defined (__INTEL_COMPILER)
54#pragma clang system_header
55#endif
56
57namespace Intrepid2 {
58
59 //============================================================================================//
60 // //
61 // Reference nodes //
62 // //
63 //============================================================================================//
64
65 template<typename DeviceType>
66 template<typename cellCenterValueType, class ...cellCenterProperties>
67 void
69 getReferenceCellCenter( Kokkos::DynRankView<cellCenterValueType,cellCenterProperties...> cellCenter,
70 const shards::CellTopology cell ) {
71#ifdef HAVE_INTREPID2_DEBUG
72 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(cell), std::invalid_argument,
73 ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): the specified cell topology does not have a reference cell." );
74
75 INTREPID2_TEST_FOR_EXCEPTION( cellCenter.rank() != 1, std::invalid_argument,
76 ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): cellCenter must have rank 1." );
77
78 INTREPID2_TEST_FOR_EXCEPTION( cellCenter.extent(0) < cell.getDimension(), std::invalid_argument,
79 ">>> ERROR (Intrepid2::CellTools::getReferenceCellCenter): cellCenter must have dimension at least as large as cell.getDimension()." );
80#endif
81
82 constexpr bool is_accessible =
83 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
84 typename decltype(cellCenter)::memory_space>::accessible;
85 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceCellCenter(..): output view's memory space is not compatible with DeviceType");
86
87 const ordinal_type dim = cell.getDimension();
88
89 const auto refCellCenter = RefCellCenter<DeviceType>::get(cell.getKey());
90
91 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,dim),
92 KOKKOS_LAMBDA (const int &i) {cellCenter(i) = refCellCenter(i);}
93 );
94 }
95
96
97 template<typename DeviceType>
98 template<typename cellVertexValueType, class ...cellVertexProperties>
99 void
101 getReferenceVertex( Kokkos::DynRankView<cellVertexValueType,cellVertexProperties...> cellVertex,
102 const shards::CellTopology cell,
103 const ordinal_type vertexOrd ) {
104#ifdef HAVE_INTREPID2_DEBUG
105 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(cell), std::invalid_argument,
106 ">>> ERROR (Intrepid2::CellTools::getReferenceVertex): the specified cell topology does not have a reference cell." );
107
108 INTREPID2_TEST_FOR_EXCEPTION( (vertexOrd < 0) || vertexOrd > static_cast<ordinal_type>(cell.getVertexCount()), std::invalid_argument,
109 ">>> ERROR (Intrepid2::CellTools::getReferenceVertex): invalid node ordinal for the specified cell topology." );
110
111 INTREPID2_TEST_FOR_EXCEPTION( cellVertex.rank() != 1, std::invalid_argument,
112 ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNodes must have rank 1." );
113
114 INTREPID2_TEST_FOR_EXCEPTION( cellVertex.extent(0) < cell.getDimension(), std::invalid_argument,
115 ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNodes must have dimension at least as large as cell.getDimension()." );
116#endif
117
118 constexpr bool is_accessible =
119 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
120 typename decltype(cellVertex)::memory_space>::accessible;
121 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceVertex(..): output view's memory space is not compatible with DeviceType");
122
123 const auto refNodes = RefCellNodes<DeviceType>::get(cell.getKey());
124
125 ordinal_type dim = cell.getDimension();
126 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,dim),
127 KOKKOS_LAMBDA (const int &i) {cellVertex(i) = refNodes(vertexOrd,i);}
128 );
129 }
130
131
132 template<typename DeviceType>
133 template<typename subcellVertexValueType, class ...subcellVertexProperties>
134 void
136 getReferenceSubcellVertices( Kokkos::DynRankView<subcellVertexValueType,subcellVertexProperties...> subcellVertices,
137 const ordinal_type subcellDim,
138 const ordinal_type subcellOrd,
139 const shards::CellTopology parentCell ) {
140#ifdef HAVE_INTREPID2_DEBUG
141 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
142 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): the specified cell topology does not have a reference cell." );
143
144 INTREPID2_TEST_FOR_EXCEPTION( subcellDim > static_cast<ordinal_type>(parentCell.getDimension()), std::invalid_argument,
145 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcell dimension cannot exceed cell dimension." );
146
147 INTREPID2_TEST_FOR_EXCEPTION( subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
148 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcell ordinal cannot exceed subcell count." );
149
150 // Verify subcellVertices rank and dimensions
151 INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.rank() != 2, std::invalid_argument,
152 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces must have rank 2." );
153
154 // need to match to node count as it invokes getReferenceSubcellNodes
155 INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.extent(0) != parentCell.getNodeCount(subcellDim, subcellOrd), std::invalid_argument,
156 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces dimension(0) must match to parent cell vertex count." );
157
158 INTREPID2_TEST_FOR_EXCEPTION( subcellVertices.extent(1) != parentCell.getDimension(), std::invalid_argument,
159 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellVertices): subcellVertieces dimension(1) must match to parent cell dimension." );
160#endif
161 getReferenceSubcellNodes(subcellVertices,
162 subcellDim,
163 subcellOrd,
164 parentCell);
165 }
166
167
168 template<typename DeviceType>
169 template<typename cellNodeValueType, class ...cellNodeProperties>
170 void
172 getReferenceNode( Kokkos::DynRankView<cellNodeValueType,cellNodeProperties...> cellNode,
173 const shards::CellTopology cell,
174 const ordinal_type nodeOrd ) {
175#ifdef HAVE_INTREPID2_DEBUG
176 INTREPID2_TEST_FOR_EXCEPTION( nodeOrd >= static_cast<ordinal_type>(cell.getNodeCount()), std::invalid_argument,
177 ">>> ERROR (Intrepid2::CellTools::getReferenceNode): invalid node ordinal for the specified cell topology." );
178
179 INTREPID2_TEST_FOR_EXCEPTION( cellNode.rank() != 1, std::invalid_argument,
180 ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNode must have rank 1." );
181
182 INTREPID2_TEST_FOR_EXCEPTION( cellNode.extent(0) < cell.getDimension(), std::invalid_argument,
183 ">>> ERROR (Intrepid2::CellTools::getReferenceNode): cellNode must have dimension at least as large as cell.getDimension()." );
184#endif
185
186 constexpr bool is_accessible =
187 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
188 typename decltype(cellNode)::memory_space>::accessible;
189 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceNode(..): output view's memory space is not compatible with DeviceType");
190
191 const auto refNodes = RefCellNodes<DeviceType>::get(cell.getKey());
192
193 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,cell.getDimension()),
194 KOKKOS_LAMBDA (const int &i) {cellNode(i) = refNodes(nodeOrd,i);}
195 );
196 }
197
198 template<typename DeviceType>
199 template<typename subcellNodeValueType, class ...subcellNodeProperties>
200 void
202 getReferenceSubcellNodes( Kokkos::DynRankView<subcellNodeValueType,subcellNodeProperties...> subcellNodes,
203 const ordinal_type subcellDim,
204 const ordinal_type subcellOrd,
205 const shards::CellTopology parentCell ) {
206#ifdef HAVE_INTREPID2_DEBUG
207 INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
208 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): the specified cell topology does not have a reference cell.");
209
210 INTREPID2_TEST_FOR_EXCEPTION( subcellDim > static_cast<ordinal_type>(parentCell.getDimension()), std::invalid_argument,
211 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcell dimension out of range.");
212
213 INTREPID2_TEST_FOR_EXCEPTION( subcellOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(subcellDim)), std::invalid_argument,
214 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcell ordinal out of range.");
215
216 // Verify subcellNodes rank and dimensions
217 INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.rank() != 2, std::invalid_argument,
218 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes must have rank 2.");
219
220 INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.extent(0) != parentCell.getNodeCount(subcellDim, subcellOrd), std::invalid_argument,
221 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes dimension (0) must match to node count in cell.");
222
223 INTREPID2_TEST_FOR_EXCEPTION( subcellNodes.extent(1) != parentCell.getDimension(), std::invalid_argument,
224 ">>> ERROR (Intrepid2::CellTools::getReferenceSubcellNodes): subcellNodes dimension (1) must match to cell dimension.");
225#endif
226
227 // Find how many cellWorkset does the specified subcell have.
228 const auto subcNodeCount = parentCell.getNodeCount(subcellDim, subcellOrd);
229
230 // Loop over subcell cellWorkset
231 for (size_type subcNodeOrd=0;subcNodeOrd<subcNodeCount;++subcNodeOrd) {
232 // Get the node number relative to the parent reference cell
233 const auto cellNodeOrd = parentCell.getNodeMap(subcellDim, subcellOrd, subcNodeOrd);
234
235 auto dst = Kokkos::subdynrankview(subcellNodes, subcNodeOrd, Kokkos::ALL());
236 getReferenceNode(dst, parentCell, cellNodeOrd);
237 }
238 }
239
240 template<typename DeviceType>
241 template<typename refEdgeTangentValueType, class ...refEdgeTangentProperties>
242 void
244 getReferenceEdgeTangent( Kokkos::DynRankView<refEdgeTangentValueType,refEdgeTangentProperties...> refEdgeTangent,
245 const ordinal_type edgeOrd,
246 const shards::CellTopology parentCell ) {
247#ifdef HAVE_INTREPID2_DEBUG
248 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
249 parentCell.getDimension() != 3, std::invalid_argument,
250 ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): two or three-dimensional parent cell required");
251
252 INTREPID2_TEST_FOR_EXCEPTION( refEdgeTangent.rank() != 1, std::invalid_argument,
253 ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): rank = 1 required for output arrays");
254
255 INTREPID2_TEST_FOR_EXCEPTION( refEdgeTangent.extent(0) != parentCell.getDimension(), std::invalid_argument,
256 ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): output array size is required to match space dimension");
257
258 INTREPID2_TEST_FOR_EXCEPTION( edgeOrd < 0 ||
259 edgeOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(1)), std::invalid_argument,
260 ">>> ERROR (Intrepid2::CellTools::getReferenceEdgeTangent): edge ordinal out of bounds");
261
262#endif
263 constexpr bool is_accessible =
264 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
265 typename decltype(refEdgeTangent)::memory_space>::accessible;
266 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceEdgeTangent(..): output view's memory space is not compatible with DeviceType");
267
268 const auto edgeMap = RefSubcellParametrization<DeviceType>::get(1, parentCell.getKey());
269
270 // All ref. edge maps have affine coordinate functions: f_dim(u) = C_0(dim) + C_1(dim)*u,
271 // => edge Tangent: -> C_1(*)
272 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,parentCell.getDimension()),
273 KOKKOS_LAMBDA (const int &i) {refEdgeTangent(i) = edgeMap(edgeOrd, i, 1);}
274 );
275 }
276
277
278 template<typename DeviceType>
279 template<typename refFaceTanValueType, class ...refFaceTanProperties>
280 void
282 getReferenceFaceTangents( Kokkos::DynRankView<refFaceTanValueType,refFaceTanProperties...> refFaceTanU,
283 Kokkos::DynRankView<refFaceTanValueType,refFaceTanProperties...> refFaceTanV,
284 const ordinal_type faceOrd,
285 const shards::CellTopology parentCell ) {
286#ifdef HAVE_INTREPID2_DEBUG
287 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
288 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): three-dimensional parent cell required");
289
290 INTREPID2_TEST_FOR_EXCEPTION( faceOrd < 0 || faceOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(2)), std::invalid_argument,
291 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): face ordinal out of bounds");
292
293 INTREPID2_TEST_FOR_EXCEPTION( refFaceTanU.rank() != 1 || refFaceTanV.rank() != 1, std::invalid_argument,
294 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): rank = 1 required for output arrays");
295
296 INTREPID2_TEST_FOR_EXCEPTION( refFaceTanU.extent(0) != parentCell.getDimension(), std::invalid_argument,
297 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
298
299 INTREPID2_TEST_FOR_EXCEPTION( refFaceTanV.extent(0) != parentCell.getDimension(), std::invalid_argument,
300 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceTangents): dim0 (spatial dim) must match that of parent cell");
301#endif
302 constexpr bool is_accessible =
303 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
304 typename decltype(refFaceTanU)::memory_space>::accessible;
305 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceFaceTangents(..): output views' memory spaces are not compatible with DeviceType");
306
307
308 const auto faceMap = RefSubcellParametrization<DeviceType>::get(2, parentCell.getKey());
309
310 /* All ref. face maps have affine coordinate functions: f_dim(u,v) = C_0(dim) + C_1(dim)*u + C_2(dim)*v
311 * ` => Tangent vectors are: refFaceTanU -> C_1(*); refFaceTanV -> C_2(*)
312 */
313
314 // set refFaceTanU -> C_1(*)
315 // set refFaceTanV -> C_2(*)
316 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,parentCell.getDimension()),
317 KOKKOS_LAMBDA (const int &i) {
318 refFaceTanU(i) = faceMap(faceOrd, i, 1);
319 refFaceTanV(i) = faceMap(faceOrd, i, 2);
320 });
321 }
322
323 template<typename DeviceType>
324 template<typename refSideNormalValueType, class ...refSideNormalProperties>
325 void
327 getReferenceSideNormal( Kokkos::DynRankView<refSideNormalValueType,refSideNormalProperties...> refSideNormal,
328 const ordinal_type sideOrd,
329 const shards::CellTopology parentCell ) {
330#ifdef HAVE_INTREPID2_DEBUG
331 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
332 parentCell.getDimension() != 3, std::invalid_argument,
333 ">>> ERROR (Intrepid2::CellTools::getReferenceSideNormal): two or three-dimensional parent cell required");
334
335 // Check side ordinal: by definition side is subcell whose dimension = parentCell.getDimension()-1
336 INTREPID2_TEST_FOR_EXCEPTION( sideOrd < 0 || sideOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(parentCell.getDimension() - 1)), std::invalid_argument,
337 ">>> ERROR (Intrepid2::CellTools::getReferenceSideNormal): side ordinal out of bounds");
338#endif
339 constexpr bool is_accessible = Kokkos::Impl::MemorySpaceAccess<
340 MemSpaceType,
341 typename decltype(refSideNormal)::memory_space>::accessible;
342 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceSideNormal(..): output view's memory space is not compatible with DeviceType");
343
344 const auto dim = parentCell.getDimension();
345 if (dim == 2) {
346 // 2D parent cells: side = 1D subcell (edge), call the edge tangent method and rotate tangents
347 getReferenceEdgeTangent(refSideNormal, sideOrd, parentCell);
348
349 // rotate t(t1, t2) to get n(t2, -t1) so that (n,t) is positively oriented: det(n1,n2/t1,t2)>0
350 Kokkos::parallel_for(Kokkos::RangePolicy<ExecSpaceType>(0,1),
351 KOKKOS_LAMBDA (const int &) {
352 refSideNormalValueType tmp = refSideNormal(0);
353 refSideNormal(0) = refSideNormal(1);
354 refSideNormal(1) = -tmp;
355 });
356 } else {
357 // 3D parent cell: side = 2D subcell (face), call the face normal method.
358 getReferenceFaceNormal(refSideNormal, sideOrd, parentCell);
359 }
360 }
361
362
363 template<typename DeviceType>
364 template<typename refFaceNormalValueType, class ...refFaceNormalProperties>
365 void
367 getReferenceFaceNormal( Kokkos::DynRankView<refFaceNormalValueType,refFaceNormalProperties...> refFaceNormal,
368 const ordinal_type faceOrd,
369 const shards::CellTopology parentCell ) {
370#ifdef HAVE_INTREPID2_DEBUG
371 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
372 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): three-dimensional parent cell required");
373
374 INTREPID2_TEST_FOR_EXCEPTION( faceOrd < 0 || faceOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(2)), std::invalid_argument,
375 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): face ordinal out of bounds");
376
377 INTREPID2_TEST_FOR_EXCEPTION( refFaceNormal.rank() != 1, std::invalid_argument,
378 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): rank = 1 required for output array");
379
380 INTREPID2_TEST_FOR_EXCEPTION( refFaceNormal.extent(0) != parentCell.getDimension(), std::invalid_argument,
381 ">>> ERROR (Intrepid2::CellTools::getReferenceFaceNormal): dim0 (spatial dim) must match that of parent cell");
382#endif
383 constexpr bool is_accessible = Kokkos::Impl::MemorySpaceAccess<
384 MemSpaceType,
385 typename decltype(refFaceNormal)::memory_space>::accessible;
386 static_assert(is_accessible, "CellTools<DeviceType>::getReferenceFaceNormal(..): output view's memory space is not compatible with DeviceType");
387
388 // Reference face normal = vector product of reference face tangents. Allocate temp FC storage:
389 const auto dim = parentCell.getDimension();
390 auto vcprop = Kokkos::common_view_alloc_prop(refFaceNormal);
391 using common_value_type = typename decltype(vcprop)::value_type;
392 Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanU ( Kokkos::view_alloc("CellTools::getReferenceFaceNormal::refFaceTanU", vcprop), dim );
393 Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanV ( Kokkos::view_alloc("CellTools::getReferenceFaceNormal::refFaceTanV", vcprop), dim );
394 getReferenceFaceTangents(refFaceTanU, refFaceTanV, faceOrd, parentCell);
395
396 RealSpaceTools<DeviceType>::vecprod(refFaceNormal, refFaceTanU, refFaceTanV);
397 }
398
399 template<typename DeviceType>
400 template<typename edgeTangentValueType, class ...edgeTangentProperties,
401 typename worksetJacobianValueType, class ...worksetJacobianProperties>
402 void
404 getPhysicalEdgeTangents( Kokkos::DynRankView<edgeTangentValueType,edgeTangentProperties...> edgeTangents,
405 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
406 const ordinal_type worksetEdgeOrd,
407 const shards::CellTopology parentCell ) {
408#ifdef HAVE_INTREPID2_DEBUG
409 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3 &&
410 parentCell.getDimension() != 2, std::invalid_argument,
411 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): 2D or 3D parent cell required." );
412
413 // (1) edgeTangents is rank-3 (C,P,D) and D=2, or 3 is required
414 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.rank() != 3, std::invalid_argument,
415 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents requires rank 3." );
416 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(2) != 2 &&
417 edgeTangents.extent(2) != 3, std::invalid_argument,
418 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension(2) must be 2 or 3." );
419
420 // (2) worksetJacobians in rank-4 (C,P,D,D) and D=2, or 3 is required
421 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
422 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians requires rank 4." );
423 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 2 &&
424 worksetJacobians.extent(2) != 3, std::invalid_argument,
425 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) must be 2 or 3." );
426 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
427 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): worksetJacobians dimension(2) and (3) must match each other." );
428
429 // (4) cross-check array dimensions: edgeTangents (C,P,D) vs. worksetJacobians (C,P,D,D)
430 for (auto i=0;i<3;++i) {
431 INTREPID2_TEST_FOR_EXCEPTION( edgeTangents.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
432 ">>> ERROR (Intrepid2::CellTools::getPhysicalEdgeTangents): edgeTangents dimension (i) does not match to worksetJacobians dimension(i)." );
433 }
434#endif
435 constexpr bool are_accessible =
436 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
437 typename decltype(edgeTangents)::memory_space>::accessible &&
438 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
439 typename decltype(worksetJacobians)::memory_space>::accessible;
440 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalEdgeTangents(..): input/output views' memory spaces are not compatible with DeviceType");
441
442
443 // Storage for constant reference edge tangent: rank-1 (D) arrays
444 const auto dim = parentCell.getDimension();
445 auto vcprop = Kokkos::common_view_alloc_prop(edgeTangents);
446 using common_value_type = typename decltype(vcprop)::value_type;
447 Kokkos::DynRankView< common_value_type, DeviceType > refEdgeTan ( Kokkos::view_alloc("CellTools::getPhysicalEdgeTangents::refEdgeTan", vcprop), dim);
448 getReferenceEdgeTangent(refEdgeTan, worksetEdgeOrd, parentCell);
449
450 RealSpaceTools<DeviceType>::matvec(edgeTangents, worksetJacobians, refEdgeTan);
451 }
452
453
454 template<typename DeviceType>
455 template<typename faceTanValueType, class ...faceTanProperties,
456 typename worksetJacobianValueType, class ...worksetJacobianProperties>
457 void
459 getPhysicalFaceTangents( Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanU,
460 Kokkos::DynRankView<faceTanValueType,faceTanProperties...> faceTanV,
461 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
462 const ordinal_type worksetFaceOrd,
463 const shards::CellTopology parentCell ) {
464#ifdef HAVE_INTREPID2_DEBUG
465 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
466 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): three-dimensional parent cell required");
467
468 // (1) faceTanU and faceTanV are rank-3 (C,P,D) and D=3 is required
469 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.rank() != 3 ||
470 faceTanV.rank() != 3, std::invalid_argument,
471 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V must have rank 3." );
472
473 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(2) != 3 ||
474 faceTanV.extent(2) != 3, std::invalid_argument,
475 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (2) must be 3." );
476
477 for (auto i=0;i<3;++i) {
478 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != faceTanV.extent(i), std::invalid_argument,
479 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): faceTan U,V dimension (i) must match each other." );
480 }
481
482 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
483 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
484 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians must have rank 4." );
485
486 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
487 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) must be 3." );
488
489 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
490 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(2) and dimension(3) must match." );
491
492 // (4) cross-check array dimensions: faceTanU (C,P,D) vs. worksetJacobians (C,P,D,D)
493 for (auto i=0;i<3;++i) {
494 INTREPID2_TEST_FOR_EXCEPTION( faceTanU.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
495 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceTangents): worksetJacobians dimension(i) and faceTan dimension (i) must match." );
496 }
497#endif
498 constexpr bool are_accessible =
499 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
500 typename decltype(faceTanU)::memory_space>::accessible &&
501 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
502 typename decltype(worksetJacobians)::memory_space>::accessible;
503 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceTangents(..): input/output views' memory spaces are not compatible with DeviceType");
504
505 // Temp storage for the pair of constant ref. face tangents: rank-1 (D) arrays
506 const auto dim = parentCell.getDimension();
507
508 auto vcprop = Kokkos::common_view_alloc_prop(faceTanU);
509 using common_value_type = typename decltype(vcprop)::value_type;
510 Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanU", vcprop), dim);
511 Kokkos::DynRankView< common_value_type, DeviceType > refFaceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceTangents::refFaceTanV", vcprop), dim);
512
513 getReferenceFaceTangents(refFaceTanU, refFaceTanV, worksetFaceOrd, parentCell);
514
515 RealSpaceTools<DeviceType>::matvec(faceTanU, worksetJacobians, refFaceTanU);
516 RealSpaceTools<DeviceType>::matvec(faceTanV, worksetJacobians, refFaceTanV);
517 }
518
519 namespace FunctorCellTools {
520 template<typename normalsViewType,
521 typename tangentsViewType>
523 normalsViewType edgeNormals_;
524 const tangentsViewType edgeTangents_;
525
526 KOKKOS_INLINE_FUNCTION
527 F_edgeNormalsFromTangents( normalsViewType edgeNormals,
528 const tangentsViewType refEdgeTangents)
529 : edgeNormals_(edgeNormals), edgeTangents_(refEdgeTangents){};
530
531 KOKKOS_INLINE_FUNCTION
532 void operator()(const size_type iter) const {
533 size_type cell, pt;
534 unrollIndex( cell, pt, edgeTangents_.extent(0),
535 edgeTangents_.extent(1), iter );
536 edgeNormals_(cell,pt,0) = edgeTangents_(cell,pt,1);
537 edgeNormals_(cell,pt,1) = -edgeTangents_(cell,pt,0);
538 }
539 };
540}
541
542 template<typename DeviceType>
543 template<typename sideNormalValueType, class ...sideNormalProperties,
544 typename worksetJacobianValueType, class ...worksetJacobianProperties>
545 void
547 getPhysicalSideNormals( Kokkos::DynRankView<sideNormalValueType,sideNormalProperties...> sideNormals,
548 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
549 const ordinal_type worksetSideOrd,
550 const shards::CellTopology parentCell ) {
551#ifdef HAVE_INTREPID2_DEBUG
552 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 2 &&
553 parentCell.getDimension() != 3, std::invalid_argument,
554 ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): two or three-dimensional parent cell required");
555
556 // Check side ordinal: by definition side is subcell whose dimension = parentCell.getDimension()-1
557 INTREPID2_TEST_FOR_EXCEPTION( worksetSideOrd < 0 ||
558 worksetSideOrd >= static_cast<ordinal_type>(parentCell.getSubcellCount(parentCell.getDimension() - 1)), std::invalid_argument,
559 ">>> ERROR (Intrepid2::CellTools::getPhysicalSideNormals): side ordinal out of bounds");
560#endif
561 constexpr bool are_accessible =
562 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
563 typename decltype(sideNormals)::memory_space>::accessible &&
564 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
565 typename decltype(worksetJacobians)::memory_space>::accessible;
566 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalSideNormals(..): input/output views' memory spaces are not compatible with DeviceType");
567
568 const auto dim = parentCell.getDimension();
569
570 if (dim == 2) {
571 // compute edge tangents and rotate it
572 auto vcprop = Kokkos::common_view_alloc_prop(sideNormals);
573 using common_value_type = typename decltype(vcprop)::value_type;
574 Kokkos::DynRankView< common_value_type, DeviceType > edgeTangents ( Kokkos::view_alloc("CellTools::getPhysicalSideNormals::edgeTan", vcprop),
575 sideNormals.extent(0),
576 sideNormals.extent(1),
577 sideNormals.extent(2));
578 getPhysicalEdgeTangents(edgeTangents, worksetJacobians, worksetSideOrd, parentCell);
579
580 //Note: this function has several template parameters and the compiler gets confused if using a lambda function
581 using FunctorType = FunctorCellTools::F_edgeNormalsFromTangents<decltype(sideNormals), decltype(edgeTangents)>;
582 const auto loopSize = edgeTangents.extent(0)*edgeTangents.extent(1);
583 Kokkos::RangePolicy<ExecSpaceType,Kokkos::Schedule<Kokkos::Static> > policy(0, loopSize);
584 Kokkos::parallel_for( policy, FunctorType(sideNormals, edgeTangents) );
585
586 } else {
587 getPhysicalFaceNormals(sideNormals, worksetJacobians, worksetSideOrd, parentCell);
588 }
589 }
590
591
592 template<typename DeviceType>
593 template<typename faceNormalValueType, class ...faceNormalProperties,
594 typename worksetJacobianValueType, class ...worksetJacobianProperties>
595 void
597 getPhysicalFaceNormals( Kokkos::DynRankView<faceNormalValueType,faceNormalProperties...> faceNormals,
598 const Kokkos::DynRankView<worksetJacobianValueType,worksetJacobianProperties...> worksetJacobians,
599 const ordinal_type worksetFaceOrd,
600 const shards::CellTopology parentCell ) {
601#ifdef HAVE_INTREPID2_DEBUG
602 INTREPID2_TEST_FOR_EXCEPTION( parentCell.getDimension() != 3, std::invalid_argument,
603 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): three-dimensional parent cell required." );
604
605 // (1) faceNormals is rank-3 (C,P,D) and D=3 is required
606 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.rank() != 3, std::invalid_argument,
607 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals must have a rank 3." );
608 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(2) != 3, std::invalid_argument,
609 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (2) must be 3." );
610
611 // (3) worksetJacobians in rank-4 (C,P,D,D) and D=3 is required
612 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.rank() != 4, std::invalid_argument,
613 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians must have a rank 4." );
614 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != 3, std::invalid_argument,
615 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must be 3." );
616 INTREPID2_TEST_FOR_EXCEPTION( worksetJacobians.extent(2) != worksetJacobians.extent(3), std::invalid_argument,
617 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): worksetJacobians dimension (2) must match to dimension (3)." );
618
619 // (4) cross-check array dimensions: faceNormals (C,P,D) vs. worksetJacobians (C,P,D,D)
620 for (auto i=0;i<3;++i) {
621 INTREPID2_TEST_FOR_EXCEPTION( faceNormals.extent(i) != worksetJacobians.extent(i), std::invalid_argument,
622 ">>> ERROR (Intrepid2::CellTools::getPhysicalFaceNormals): faceNormals dimension (i) must match to worksetJacobians dimension (i)." );
623 }
624#endif
625 constexpr bool are_accessible =
626 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
627 typename decltype(faceNormals)::memory_space>::accessible &&
628 Kokkos::Impl::MemorySpaceAccess<MemSpaceType,
629 typename decltype(worksetJacobians)::memory_space>::accessible;
630 static_assert(are_accessible, "CellTools<DeviceType>::getPhysicalFaceNormals(..): input/output views' memory spaces are not compatible with DeviceType");
631
632
633 // this should be provided from users
634 // Storage for physical face tangents: rank-3 (C,P,D) arrays
635 const auto worksetSize = worksetJacobians.extent(0);
636 const auto facePtCount = worksetJacobians.extent(1);
637 const auto dim = parentCell.getDimension();
638
639 auto vcprop = Kokkos::common_view_alloc_prop(faceNormals);
640 using common_value_type = typename decltype(vcprop)::value_type;
641 Kokkos::DynRankView< common_value_type, DeviceType > faceTanU ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanU", vcprop), worksetSize, facePtCount, dim);
642 Kokkos::DynRankView< common_value_type, DeviceType > faceTanV ( Kokkos::view_alloc("CellTools::getPhysicalFaceNormals::faceTanV", vcprop), worksetSize, facePtCount, dim);
643
644 getPhysicalFaceTangents(faceTanU, faceTanV,
645 worksetJacobians,
646 worksetFaceOrd,
647 parentCell);
648
649 RealSpaceTools<DeviceType>::vecprod(faceNormals, faceTanU, faceTanV);
650 }
651}
652
653#endif
static void getReferenceSubcellNodes(Kokkos::DynRankView< subcellNodeValueType, subcellNodeProperties... > subcellNodes, const ordinal_type subcellDim, const ordinal_type subcellOrd, const shards::CellTopology parentCell)
Retrieves the Cartesian coordinates of all nodes of a reference subcell.
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.
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...
static void getReferenceEdgeTangent(Kokkos::DynRankView< refEdgeTangentValueType, refEdgeTangentProperties... > refEdgeTangent, const ordinal_type edgeOrd, const shards::CellTopology parentCell)
Computes constant tangent vectors to edges of 2D or 3D reference cells.
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...
static void getReferenceFaceNormal(Kokkos::DynRankView< refFaceNormalValueType, refFaceNormalProperties... > refFaceNormal, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to faces of 3D reference cell.
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 ...
static void getReferenceCellCenter(Kokkos::DynRankView< cellCenterValueType, cellCenterProperties... > cellCenter, const shards::CellTopology cell)
Computes the Cartesian coordinates of reference cell barycenter.
static void getReferenceFaceTangents(Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanU, Kokkos::DynRankView< refFaceTanValueType, refFaceTanProperties... > refFaceTanV, const ordinal_type faceOrd, const shards::CellTopology parentCell)
Computes pairs of constant tangent vectors to faces of a 3D reference cells.
static void getReferenceSideNormal(Kokkos::DynRankView< refSideNormalValueType, refSideNormalProperties... > refSideNormal, const ordinal_type sideOrd, const shards::CellTopology parentCell)
Computes constant normal vectors to sides of 2D or 3D reference cells.
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.
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.
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 .
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.
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
static ConstViewType get(const unsigned cellTopoKey)
Retrieves the Cartesian coordinates of a reference cell barycenter.
static ConstViewType get(const unsigned cellTopoKey)
Retrieves the Cartesian coordinates of reference cell nodes.
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...