48 #ifndef __INTREPID2_CELLTOOLS_SERIAL_HPP__
49 #define __INTREPID2_CELLTOOLS_SERIAL_HPP__
51 #include "Intrepid2_ConfigDefs.hpp"
53 #include "Shards_CellTopology.hpp"
68 typedef Kokkos::DynRankView<double,Kokkos::HostSpace> NodeDataHostView;
69 typedef Kokkos::DynRankView<const double,Kokkos::HostSpace,Kokkos::MemoryUnmanaged> ConstUnmanagedNodeDataHostView;
72 double line[2][3], line_3[3][3];
73 double triangle[3][3], triangle_4[4][3], triangle_6[6][3];
74 double quadrilateral[4][3], quadrilateral_8[8][3], quadrilateral_9[9][3];
75 double tetrahedron[4][3], tetrahedron_8[8][3], tetrahedron_10[10][3], tetrahedron_11[10][3];
76 double hexahedron[8][3], hexahedron_20[20][3], hexahedron_27[27][3];
77 double pyramid[5][3], pyramid_13[13][3], pyramid_14[14][3];
78 double wedge[6][3], wedge_15[15][3], wedge_18[18][3];
87 {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}
90 {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 0.0, 0.0}
94 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}
97 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 1/3, 1/3, 0.0}
100 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
101 { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
105 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}
108 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
109 { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}
112 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
113 { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0}, { 0.0, 0.0, 0.0}
117 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
120 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
121 { 1/3, 0.0, 1/3}, { 1/3, 1/3, 1/3}, { 1/3, 1/3, 0.0}, { 0.0, 1/3, 1/3}
124 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
125 { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
128 { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
129 { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}, { 0.0, 0.0, 0.5}, { 0.5, 0.0, 0.5}, { 0.0, 0.5, 0.5}
133 {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
134 {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0}
137 {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
138 {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
139 { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0},
140 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
141 { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0}
144 {-1.0,-1.0,-1.0}, { 1.0,-1.0,-1.0}, { 1.0, 1.0,-1.0}, {-1.0, 1.0,-1.0},
145 {-1.0,-1.0, 1.0}, { 1.0,-1.0, 1.0}, { 1.0, 1.0, 1.0}, {-1.0, 1.0, 1.0},
146 { 0.0,-1.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, {-1.0, 0.0,-1.0},
147 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0},
148 { 0.0,-1.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}, {-1.0, 0.0, 1.0},
150 { 0.0, 0.0,-1.0}, { 0.0, 0.0, 1.0}, {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, {0.0,-1.0, 0.0}, {0.0, 1.0, 0.0}
154 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0}
157 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
158 { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
159 {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}
162 {-1.0,-1.0, 0.0}, { 1.0,-1.0, 0.0}, { 1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}, { 0.0, 0.0, 1.0},
163 { 0.0,-1.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}, {-1.0, 0.0, 0.0},
164 {-0.5,-0.5, 0.5}, { 0.5,-0.5, 0.5}, { 0.5, 0.5, 0.5}, {-0.5, 0.5, 0.5}, { 0.0, 0.0, 0.0}
168 { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0}
171 { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
172 { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
173 { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0}
176 { 0.0, 0.0,-1.0}, { 1.0, 0.0,-1.0}, { 0.0, 1.0,-1.0}, { 0.0, 0.0, 1.0}, { 1.0, 0.0, 1.0}, { 0.0, 1.0, 1.0},
177 { 0.5, 0.0,-1.0}, { 0.5, 0.5,-1.0}, { 0.0, 0.5,-1.0}, { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0},
178 { 0.5, 0.0, 1.0}, { 0.5, 0.5, 1.0}, { 0.0, 0.5, 1.0},
179 { 0.5, 0.0, 0.0}, { 0.5, 0.5, 0.0}, { 0.0, 0.5, 0.0}
185 template<
typename refNodeViewType>
188 getReferenceNode(
const refNodeViewType &nodes,
189 const shards::CellTopology &cell,
190 const ordinal_type nodeOrd ) {
191 ConstUnmanagedNodeDataHostView ref;
192 switch (cell.getKey() ) {
193 case shards::Line<2>::key:
194 case shards::ShellLine<2>::key:
195 case shards::Beam<2>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().line, 2, 3);
break;
196 case shards::Line<3>::key:
197 case shards::ShellLine<3>::key:
198 case shards::Beam<3>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().line_3, 3, 3);
break;
200 case shards::Triangle<3>::key:
201 case shards::ShellTriangle<3>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().triangle, 3, 3);
break;
202 case shards::Triangle<4>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().triangle_4, 4, 3);
break;
203 case shards::Triangle<6>::key:
204 case shards::ShellTriangle<6>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().triangle_6, 6, 3);
break;
206 case shards::Quadrilateral<4>::key:
207 case shards::ShellQuadrilateral<4>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().quadrilateral, 4, 3);
break;
208 case shards::Quadrilateral<8>::key:
209 case shards::ShellQuadrilateral<8>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().quadrilateral_8, 8, 3);
break;
210 case shards::Quadrilateral<9>::key:
211 case shards::ShellQuadrilateral<9>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().quadrilateral_9, 9, 3);
break;
213 case shards::Tetrahedron<4>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().tetrahedron, 4, 3);
break;
214 case shards::Tetrahedron<8>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().tetrahedron_8, 8, 3);
break;
215 case shards::Tetrahedron<10>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().tetrahedron_10, 10, 3);
break;
216 case shards::Tetrahedron<11>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().tetrahedron_11, 11, 3);
break;
218 case shards::Hexahedron<8>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().hexahedron, 8, 3);
break;
219 case shards::Hexahedron<20>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().hexahedron_20, 20, 3);
break;
220 case shards::Hexahedron<27>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().hexahedron_27, 27, 3);
break;
222 case shards::Pyramid<5>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().pyramid, 5, 3);
break;
223 case shards::Pyramid<13>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().pyramid_13, 13, 3);
break;
224 case shards::Pyramid<14>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().pyramid_14, 14, 3);
break;
226 case shards::Wedge<6>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().wedge, 6, 3);
break;
227 case shards::Wedge<15>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().wedge_15, 15, 3);
break;
228 case shards::Wedge<18>::key: ref = ConstUnmanagedNodeDataHostView((
const double*)getRefNodeData().wedge_18, 18, 3);
break;
231 INTREPID2_TEST_FOR_ABORT(
true,
"invalid input cell topology.");
236 const ordinal_type D = cell.getDimension();
237 for (ordinal_type i=0;i<D;++i)
238 nodes(i) = ref(nodeOrd, i);
242 NodeDataHostView dummy;
243 NodeDataHostView lineEdges;
244 NodeDataHostView triEdges, quadEdges;
245 NodeDataHostView shellTriEdges, shellQuadEdges;
246 NodeDataHostView tetEdges, hexEdges, pyrEdges, wedgeEdges;
247 NodeDataHostView shellTriFaces, shellQuadFaces;
248 NodeDataHostView tetFaces, hexFaces, pyrFaces, wedgeFaces;
254 Kokkos::push_finalize_hook( [=] {
255 subcellParamData.dummy = NodeDataHostView();
256 subcellParamData.lineEdges = NodeDataHostView();
257 subcellParamData.triEdges = NodeDataHostView();
258 subcellParamData.quadEdges = NodeDataHostView();
259 subcellParamData.shellTriEdges = NodeDataHostView();
260 subcellParamData.shellQuadEdges = NodeDataHostView();
261 subcellParamData.tetEdges = NodeDataHostView();
262 subcellParamData.hexEdges = NodeDataHostView();
263 subcellParamData.pyrEdges = NodeDataHostView();
264 subcellParamData.wedgeEdges = NodeDataHostView();
265 subcellParamData.shellTriFaces = NodeDataHostView();
266 subcellParamData.shellQuadFaces = NodeDataHostView();
267 subcellParamData.tetFaces = NodeDataHostView();
268 subcellParamData.hexFaces = NodeDataHostView();
269 subcellParamData.pyrFaces = NodeDataHostView();
270 subcellParamData.wedgeFaces = NodeDataHostView();
272 return subcellParamData;
277 setSubcellParametrization() {
278 static bool isSubcellParametrizationSet =
false;
279 if (!isSubcellParametrizationSet) {
281 const auto tet = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >());
282 setSubcellParametrization( getSubcellParamData().tetFaces, 2, tet );
283 setSubcellParametrization( getSubcellParamData().tetEdges, 1, tet );
286 const auto hex = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
287 setSubcellParametrization( getSubcellParamData().hexFaces, 2, hex );
288 setSubcellParametrization( getSubcellParamData().hexEdges, 1, hex );
291 const auto pyr = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >());
292 setSubcellParametrization( getSubcellParamData().pyrFaces, 2, pyr );
293 setSubcellParametrization( getSubcellParamData().pyrEdges, 1, pyr );
296 const auto wedge = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >());
297 setSubcellParametrization( getSubcellParamData().wedgeFaces, 2, wedge );
298 setSubcellParametrization( getSubcellParamData().wedgeEdges, 1, wedge );
301 const auto tri = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >());
302 setSubcellParametrization( getSubcellParamData().triEdges, 1, tri );
305 const auto quad = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
306 setSubcellParametrization( getSubcellParamData().quadEdges, 1, quad );
309 const auto line = shards::CellTopology(shards::getCellTopologyData<shards::ShellLine<2> >());
310 setSubcellParametrization( getSubcellParamData().lineEdges, 1, line );
313 isSubcellParametrizationSet =
true;
318 setSubcellParametrization( NodeDataHostView &subcellParam,
319 const ordinal_type subcellDim,
320 const shards::CellTopology &parentCell ) {
322 const auto sc = parentCell.getSubcellCount(subcellDim);
323 const auto pcd = parentCell.getDimension();
324 const auto coeff = (subcellDim == 1) ? 2 : 3;
327 subcellParam = NodeDataHostView(
"subcellParam",
330 NodeDataHostView v0(
"v0", 3), v1(
"v1", 3), v2(
"v1", 3), v3(
"v1", 3);
331 if (subcellDim == 1) {
333 for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
336 const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
337 const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
339 getReferenceNode(v0, parentCell, v0ord);
340 getReferenceNode(v1, parentCell, v1ord);
343 subcellParam(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
344 subcellParam(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
347 subcellParam(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
348 subcellParam(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
352 subcellParam(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
353 subcellParam(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
357 else if (subcellDim == 2) {
361 for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
363 switch (parentCell.getKey(subcellDim,subcellOrd)) {
365 case shards::Triangle<3>::key:
366 case shards::Triangle<4>::key:
367 case shards::Triangle<6>::key: {
368 const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
369 const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
370 const auto v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
372 getReferenceNode(v0, parentCell, v0ord);
373 getReferenceNode(v1, parentCell, v1ord);
374 getReferenceNode(v2, parentCell, v2ord);
377 subcellParam(subcellOrd, 0, 0) = v0[0];
378 subcellParam(subcellOrd, 0, 1) = v1[0] - v0[0];
379 subcellParam(subcellOrd, 0, 2) = v2[0] - v0[0];
382 subcellParam(subcellOrd, 1, 0) = v0[1];
383 subcellParam(subcellOrd, 1, 1) = v1[1] - v0[1];
384 subcellParam(subcellOrd, 1, 2) = v2[1] - v0[1];
387 subcellParam(subcellOrd, 2, 0) = v0[2];
388 subcellParam(subcellOrd, 2, 1) = v1[2] - v0[2];
389 subcellParam(subcellOrd, 2, 2) = v2[2] - v0[2];
392 case shards::Quadrilateral<4>::key:
393 case shards::Quadrilateral<8>::key:
394 case shards::Quadrilateral<9>::key: {
395 const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
396 const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
397 const auto v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
398 const auto v3ord = parentCell.getNodeMap(subcellDim, subcellOrd, 3);
400 getReferenceNode(v0, parentCell, v0ord);
401 getReferenceNode(v1, parentCell, v1ord);
402 getReferenceNode(v2, parentCell, v2ord);
403 getReferenceNode(v3, parentCell, v3ord);
406 subcellParam(subcellOrd, 0, 0) = ( v0[0] + v1[0] + v2[0] + v3[0])/4.0;
407 subcellParam(subcellOrd, 0, 1) = (-v0[0] + v1[0] + v2[0] - v3[0])/4.0;
408 subcellParam(subcellOrd, 0, 2) = (-v0[0] - v1[0] + v2[0] + v3[0])/4.0;
410 subcellParam(subcellOrd, 1, 0) = ( v0[1] + v1[1] + v2[1] + v3[1])/4.0;
411 subcellParam(subcellOrd, 1, 1) = (-v0[1] + v1[1] + v2[1] - v3[1])/4.0;
412 subcellParam(subcellOrd, 1, 2) = (-v0[1] - v1[1] + v2[1] + v3[1])/4.0;
415 subcellParam(subcellOrd, 2, 0) = ( v0[2] + v1[2] + v2[2] + v3[2])/4.0;
416 subcellParam(subcellOrd, 2, 1) = (-v0[2] + v1[2] + v2[2] - v3[2])/4.0;
417 subcellParam(subcellOrd, 2, 2) = (-v0[2] - v1[2] + v2[2] + v3[2])/4.0;
421 INTREPID2_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
422 ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): parametrization not defined for the specified face topology.");
436 template<
typename jacobianViewType,
437 typename basisGradViewType,
438 typename nodeViewType>
439 KOKKOS_INLINE_FUNCTION
441 computeJacobian(
const jacobianViewType &jacobian,
442 const basisGradViewType &grads,
443 const nodeViewType &nodes) {
444 const auto N = nodes.extent(0);
446 const auto D = jacobian.extent(0);
447 const auto sD = jacobian.extent(1);
449 INTREPID2_TEST_FOR_ABORT( N != grads.extent(0),
"grad dimension_0 does not match to cardinality.");
450 INTREPID2_TEST_FOR_ABORT(sD != grads.extent(1),
"grad dimension_1 does not match to space dim.");
451 INTREPID2_TEST_FOR_ABORT( D != nodes.extent(1),
"node dimension_1 does not match to space dim.");
453 Kernels::Serial::gemm_trans_notrans(1.0, nodes, grads, 0.0, jacobian);
461 template<
typename pointViewType,
462 typename basisValViewType,
463 typename nodeViewType>
464 KOKKOS_INLINE_FUNCTION
466 mapToPhysicalFrame(
const pointViewType &point,
467 const basisValViewType &vals,
468 const nodeViewType &nodes) {
469 const auto N = vals.extent(0);
470 const auto D = point.extent(0);
472 INTREPID2_TEST_FOR_ABORT(N != nodes.extent(0),
"nodes dimension_0 does not match to vals dimension_0.");
473 INTREPID2_TEST_FOR_ABORT(D != nodes.extent(1),
"node dimension_1 does not match to space dim.");
475 Kernels::Serial::gemv_trans(1.0, nodes, vals, 0.0, point);
485 template<
typename implBasisType,
486 typename refPointViewType,
487 typename phyPointViewType,
488 typename nodeViewType>
489 KOKKOS_INLINE_FUNCTION
491 mapToReferenceFrame(
const refPointViewType &xref,
492 const phyPointViewType &xphy,
493 const nodeViewType &nodes) {
494 const ordinal_type sD = xref.extent(0);
495 const ordinal_type D = xphy.extent(0);
496 const ordinal_type N = nodes.extent(0);
498 INTREPID2_TEST_FOR_ABORT(sD > D,
"subcell dimension is greater than physical cell dimension.");
499 INTREPID2_TEST_FOR_ABORT(D != static_cast<ordinal_type>(nodes.extent(1)),
"xphy dimension_0 does not match to space dim.");
501 typedef typename refPointViewType::non_const_value_type value_type;
505 value_type buf[27*3 + 27 + 9 + 9 + 9 + 9 + 3 + 3] = {}, *ptr = &buf[0];
506 Kokkos::DynRankView<value_type,
507 Kokkos::Impl::ActiveExecutionMemorySpace,
508 Kokkos::MemoryUnmanaged> grads(ptr, N, sD); ptr += N*sD;
510 Kokkos::DynRankView<value_type,
511 Kokkos::Impl::ActiveExecutionMemorySpace,
512 Kokkos::MemoryUnmanaged> vals(ptr, N); ptr += N;
514 Kokkos::DynRankView<value_type,
515 Kokkos::Impl::ActiveExecutionMemorySpace,
516 Kokkos::MemoryUnmanaged> jac(ptr, D, sD); ptr += D*sD;
518 Kokkos::DynRankView<value_type,
519 Kokkos::Impl::ActiveExecutionMemorySpace,
520 Kokkos::MemoryUnmanaged> metric(ptr, sD, sD); ptr += sD*sD;
522 Kokkos::DynRankView<value_type,
523 Kokkos::Impl::ActiveExecutionMemorySpace,
524 Kokkos::MemoryUnmanaged> invmetric(ptr, sD, sD); ptr += sD*sD;
526 Kokkos::DynRankView<value_type,
527 Kokkos::Impl::ActiveExecutionMemorySpace,
528 Kokkos::MemoryUnmanaged> invdf(ptr, sD, D); ptr += sD*D;
530 Kokkos::DynRankView<value_type,
531 Kokkos::Impl::ActiveExecutionMemorySpace,
532 Kokkos::MemoryUnmanaged> xtmp(ptr, sD); ptr += sD;
534 Kokkos::DynRankView<value_type,
535 Kokkos::Impl::ActiveExecutionMemorySpace,
536 Kokkos::MemoryUnmanaged> xold(ptr, sD); ptr += sD;
539 for (ordinal_type j=0;j<D;++j) xold(j) = 0;
541 const double tol = tolerence();
545 mapToPhysicalFrame(xtmp, vals, nodes);
549 CellTools::Serial::computeJacobian(jac, grads, nodes);
551 Kernels::Serial::gemm_trans_notrans(1.0, jac, jac, 0.0, metric);
552 Kernels::Serial::inverse(invmetric, metric);
553 Kernels::Serial::gemm_notrans_trans(1.0, invmetric, jac, 0.0, invdf);
556 Kernels::Serial::z_is_axby(xtmp, 1.0, xphy, -1.0, xtmp);
557 Kernels::Serial::gemv_notrans(1.0, invdf, xtmp, 0.0, xref);
558 Kernels::Serial::z_is_axby(xref, 1.0, xold, 1.0, xref);
561 Kernels::Serial::z_is_axby(xtmp, 1.0, xold, -1.0, xref);
563 double err = Kernels::Serial::norm(xtmp, NORM_ONE);
568 Kernels::Serial::copy(xold, xref);
577 static ConstUnmanagedNodeDataHostView
579 const shards::CellTopology parent_cell) {
583 Kokkos::DynRankView<const double,Kokkos::HostSpace,Kokkos::MemoryUnmanaged> r_val;
584 if (subcell_dim == 2) {
585 switch (parent_cell.getBaseKey()) {
586 case shards::Tetrahedron<>::key: r_val = getSubcellParamData().tetFaces;
break;
587 case shards::Hexahedron<>::key: r_val = getSubcellParamData().hexFaces;
break;
588 case shards::Pyramid<>::key: r_val = getSubcellParamData().pyrFaces;
break;
589 case shards::Wedge<18>::key: r_val = getSubcellParamData().wedgeFaces;
break;
592 else if (subcell_dim == 1) {
593 switch (parent_cell.getBaseKey()) {
594 case shards::Tetrahedron<>::key: r_val = getSubcellParamData().tetEdges;
break;
595 case shards::Hexahedron<>::key: r_val = getSubcellParamData().hexEdges;
break;
596 case shards::Pyramid<>::key: r_val = getSubcellParamData().pyrEdges;
break;
597 case shards::Wedge<>::key: r_val = getSubcellParamData().wedgeEdges;
break;
598 case shards::Triangle<>::key: r_val = getSubcellParamData().triEdges;
break;
599 case shards::Quadrilateral<>::key: r_val = getSubcellParamData().quadEdges;
break;
600 case shards::Line<>::key: r_val = getSubcellParamData().lineEdges;
break;
603 INTREPID2_TEST_FOR_ABORT(r_val.rank() == 0,
"subcell param is not set up before.");
607 template<
typename refEdgeTanViewType>
610 getReferenceEdgeTangent(
const refEdgeTanViewType &ref_edge_tangent,
611 const ordinal_type edge_ordinal,
612 const shards::CellTopology &parent_cell ) {
615 const ordinal_type D = parent_cell.getDimension();
616 for (ordinal_type i=0;i<D;++i)
617 ref_edge_tangent(i) = edge_map(edge_ordinal, i, 1);
620 template<
typename refFaceTanViewType>
622 getReferenceFaceTangent(
const refFaceTanViewType &ref_face_tan_u,
623 const refFaceTanViewType &ref_face_tan_v,
624 const ordinal_type face_ordinal,
625 const shards::CellTopology &parent_cell) {
630 const ordinal_type D = parent_cell.getDimension();
631 for (ordinal_type i=0;i<D;++i) {
632 ref_face_tan_u(i) = face_map(face_ordinal, i, 1);
633 ref_face_tan_v(i) = face_map(face_ordinal, i, 2);
637 template<
typename edgeTangentViewType,
638 typename jacobianViewType>
641 getPhysicalEdgeTangent(
const edgeTangentViewType &edge_tangent,
642 const jacobianViewType &jacobian,
643 const ordinal_type edge_ordinal,
644 const shards::CellTopology &parent_cell) {
645 typedef typename edgeTangentViewType::non_const_value_type value_type;
646 const ordinal_type D = parent_cell.getDimension();
648 Kokkos::DynRankView<value_type,
649 Kokkos::Impl::ActiveExecutionMemorySpace,
650 Kokkos::MemoryUnmanaged> ref_edge_tangent(&buf[0], D);
652 getReferenceEdgeTangent(ref_edge_tangent, edge_ordinal, parent_cell);
653 Kernels::Serial::matvec_product(edge_tangent, jacobian, ref_edge_tangent);
656 template<
typename faceTanViewType,
657 typename jacobianViewType>
660 getPhysicalFaceTangent(
const faceTanViewType &face_tan_u,
661 const faceTanViewType &face_tan_v,
662 const jacobianViewType &jacobian,
663 const ordinal_type face_ordinal,
664 const shards::CellTopology &parent_cell) {
665 typedef typename faceTanViewType::non_const_value_type value_type;
666 const ordinal_type D = parent_cell.getDimension();
667 INTREPID2_TEST_FOR_ABORT(D != 3,
"computing face normal requires dimension 3.");
669 Kokkos::DynRankView<value_type,
670 Kokkos::Impl::ActiveExecutionMemorySpace,
671 Kokkos::MemoryUnmanaged> ref_face_tan_u(&buf[0], D), ref_face_tan_v(&buf[3], D);
673 getReferenceFaceTangent(ref_face_tan_u,
678 Kernels::Serial::matvec_product_d3(face_tan_u, jacobian, ref_face_tan_u);
679 Kernels::Serial::matvec_product_d3(face_tan_v, jacobian, ref_face_tan_v);
683 template<
typename faceNormalViewType,
684 typename jacobianViewType>
687 getPhysicalFaceNormal(
const faceNormalViewType &face_normal,
688 const jacobianViewType &jacobian,
689 const ordinal_type face_ordinal,
690 const shards::CellTopology &parent_cell) {
691 typedef typename faceNormalViewType::non_const_value_type value_type;
692 const ordinal_type D = parent_cell.getDimension();
693 INTREPID2_TEST_FOR_ABORT(D != 3,
"computing face normal requires dimension 3.");
695 Kokkos::DynRankView<value_type,
696 Kokkos::Impl::ActiveExecutionMemorySpace,
697 Kokkos::MemoryUnmanaged> face_tan_u(&buf[0], D), face_tan_v(&buf[3], D);
699 getPhysicalFaceTangent(face_tan_u, face_tan_v,
703 Kernels::Serial::vector_product_d3(face_normal, face_tan_u, face_tan_v);
706 template<
typename sideNormalViewType,
707 typename jacobianViewType>
710 getPhysicalSideNormal(
const sideNormalViewType &side_normal,
711 const jacobianViewType &jacobian,
712 const ordinal_type side_ordinal,
713 const shards::CellTopology &parent_cell ) {
714 const ordinal_type D = parent_cell.getDimension();
715 typedef typename sideNormalViewType::non_const_value_type value_type;
719 Kokkos::DynRankView<value_type,
720 Kokkos::Impl::ActiveExecutionMemorySpace,
721 Kokkos::MemoryUnmanaged> edge_tangent(&buf[0], D);
722 getPhysicalEdgeTangent(edge_tangent, jacobian, side_ordinal, parent_cell);
723 side_normal(0) = edge_tangent(1);
724 side_normal(1) = -edge_tangent(0);
728 getPhysicalFaceNormal(side_normal, jacobian, side_ordinal, parent_cell);
732 INTREPID2_TEST_FOR_ABORT(
true,
"cell dimension is out of range.");
Header function for Intrepid2::Util class and other utility functions.
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
Contains definitions of custom data types in Intrepid2.
Header file for small functions used in Intrepid2.