Intrepid2
Intrepid2_CellTools_Serial.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 
48 #ifndef __INTREPID2_CELLTOOLS_SERIAL_HPP__
49 #define __INTREPID2_CELLTOOLS_SERIAL_HPP__
50 
51 #include "Intrepid2_ConfigDefs.hpp"
52 
53 #include "Shards_CellTopology.hpp"
54 
55 #include "Intrepid2_Types.hpp"
56 #include "Intrepid2_Utils.hpp"
57 #include "Intrepid2_Kernels.hpp"
58 
59 namespace Intrepid2 {
60 
61  namespace Impl {
62 
66  class CellTools {
67  public:
68  typedef Kokkos::DynRankView<double,Kokkos::HostSpace> NodeDataHostView;
69  typedef Kokkos::DynRankView<const double,Kokkos::HostSpace,Kokkos::MemoryUnmanaged> ConstUnmanagedNodeDataHostView;
70 
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];
79  };
80 
81  inline
82  static const ReferenceNodeDataType&
83  getRefNodeData() {
84  const static ReferenceNodeDataType refNodeData = {
85  // line
86  { // 2
87  {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}
88  },
89  { // 3
90  {-1.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 0.0, 0.0}
91  },
92  // triangle
93  { // 3
94  { 0.0, 0.0, 0.0}, { 1.0, 0.0, 0.0}, { 0.0, 1.0, 0.0}
95  },
96  { // 4
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}
98  },
99  { // 6
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}
102  },
103  // quad
104  { // 4
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}
106  },
107  { // 8
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}
110  },
111  { // 9
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}
114  },
115  // tet
116  { // 4
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}
118  },
119  { // 8
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}
122  },
123  { // 10
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}
126  },
127  { // 11
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}
130  },
131  // hex
132  { // 8
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}
135  },
136  { // 20
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}
142  },
143  { // 27
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},
149  { 0.0, 0.0, 0.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}
151  },
152  // pyramid
153  { // 5
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}
155  },
156  { // 13
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}
160  },
161  { // 14
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}
165  },
166  // wedge
167  { // 6
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}
169  },
170  { // 15
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}
174  },
175  { // 18
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}
180  }
181  };
182  return refNodeData;
183  }
184 
185  template<typename refNodeViewType>
186  static
187  void
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;
199 
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;
205 
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;
212 
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;
217 
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;
221 
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;
225 
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;
229 
230  default: {
231  INTREPID2_TEST_FOR_ABORT( true, "invalid input cell topology.");
232  break;
233  }
234  }
235 
236  const ordinal_type D = cell.getDimension();
237  for (ordinal_type i=0;i<D;++i)
238  nodes(i) = ref(nodeOrd, i);
239  }
240 
242  NodeDataHostView dummy;
243  NodeDataHostView lineEdges; // edge maps for 2d non-standard cells; shell line and beam
244  NodeDataHostView triEdges, quadEdges; // edge maps for 2d standard cells
245  NodeDataHostView shellTriEdges, shellQuadEdges; // edge maps for 3d non-standard cells; shell tri and quad
246  NodeDataHostView tetEdges, hexEdges, pyrEdges, wedgeEdges; // edge maps for 3d standard cells
247  NodeDataHostView shellTriFaces, shellQuadFaces; // face maps for 3d non-standard cells
248  NodeDataHostView tetFaces, hexFaces, pyrFaces, wedgeFaces; // face maps for 3d standard cells
249  };
250 
251  inline
252  static SubcellParamDataType& getSubcellParamData() {
253  static SubcellParamDataType subcellParamData;
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();
271  });
272  return subcellParamData;
273  }
274 
275  inline
276  static void
277  setSubcellParametrization() {
278  static bool isSubcellParametrizationSet = false;
279  if (!isSubcellParametrizationSet) {
280  {
281  const auto tet = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >());
282  setSubcellParametrization( getSubcellParamData().tetFaces, 2, tet );
283  setSubcellParametrization( getSubcellParamData().tetEdges, 1, tet );
284  }
285  {
286  const auto hex = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
287  setSubcellParametrization( getSubcellParamData().hexFaces, 2, hex );
288  setSubcellParametrization( getSubcellParamData().hexEdges, 1, hex );
289  }
290  {
291  const auto pyr = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >());
292  setSubcellParametrization( getSubcellParamData().pyrFaces, 2, pyr );
293  setSubcellParametrization( getSubcellParamData().pyrEdges, 1, pyr );
294  }
295  {
296  const auto wedge = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >());
297  setSubcellParametrization( getSubcellParamData().wedgeFaces, 2, wedge );
298  setSubcellParametrization( getSubcellParamData().wedgeEdges, 1, wedge );
299  }
300  {
301  const auto tri = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >());
302  setSubcellParametrization( getSubcellParamData().triEdges, 1, tri );
303  }
304  {
305  const auto quad = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
306  setSubcellParametrization( getSubcellParamData().quadEdges, 1, quad );
307  }
308  {
309  const auto line = shards::CellTopology(shards::getCellTopologyData<shards::ShellLine<2> >());
310  setSubcellParametrization( getSubcellParamData().lineEdges, 1, line );
311  }
312  }
313  isSubcellParametrizationSet = true;
314  }
315 
316  inline
317  static void
318  setSubcellParametrization( NodeDataHostView &subcellParam,
319  const ordinal_type subcellDim,
320  const shards::CellTopology &parentCell ) {
321  // get subcellParametrization dimensions: (sc, pcd, coeff)
322  const auto sc = parentCell.getSubcellCount(subcellDim);
323  const auto pcd = parentCell.getDimension();
324  const auto coeff = (subcellDim == 1) ? 2 : 3;
325 
326  // create a view
327  subcellParam = NodeDataHostView("subcellParam",
328  sc, pcd, coeff);
329 
330  NodeDataHostView v0("v0", 3), v1("v1", 3), v2("v1", 3), v3("v1", 3);
331  if (subcellDim == 1) {
332  // Edge parametrizations of 2D and 3D cells (shell lines and beams are 2D cells with edges)
333  for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
334  // vertexK[0] = x_k; vertexK[1] = y_k; vertexK[2] = z_k; z_k = 0 for 2D cells
335  // Note that ShellLine and Beam are 2D cells!
336  const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
337  const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
338 
339  getReferenceNode(v0, parentCell, v0ord);
340  getReferenceNode(v1, parentCell, v1ord);
341 
342  // x(t) = (x0 + x1)/2 + t*(x1 - x0)/2
343  subcellParam(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
344  subcellParam(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
345 
346  // y(t) = (y0 + y1)/2 + t*(y1 - y0)/2
347  subcellParam(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
348  subcellParam(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
349 
350  if( pcd == 3 ) {
351  // z(t) = (z0 + z1)/2 + t*(z1 - z0)/2
352  subcellParam(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
353  subcellParam(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
354  }
355  }
356  }
357  else if (subcellDim == 2) {
358  // Face parametrizations of 3D cells: (shell Tri and Quad are 3D cells with faces)
359  // A 3D cell can have both Tri and Quad faces, but because they are affine images of the
360  // parametrization domain, 3 coefficients are enough to store them in both cases.
361  for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
362 
363  switch (parentCell.getKey(subcellDim,subcellOrd)) {
364 
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);
371 
372  getReferenceNode(v0, parentCell, v0ord);
373  getReferenceNode(v1, parentCell, v1ord);
374  getReferenceNode(v2, parentCell, v2ord);
375 
376  // x(u,v) = x0 + (x1 - x0)*u + (x2 - x0)*v
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];
380 
381  // y(u,v) = y0 + (y1 - y0)*u + (y2 - y0)*v
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];
385 
386  // z(u,v) = z0 + (z1 - z0)*u + (z2 - z0)*v
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];
390  break;
391  }
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);
399 
400  getReferenceNode(v0, parentCell, v0ord);
401  getReferenceNode(v1, parentCell, v1ord);
402  getReferenceNode(v2, parentCell, v2ord);
403  getReferenceNode(v3, parentCell, v3ord);
404 
405  // x(u,v) = (x0+x1+x2+x3)/4+u*(-x0+x1+x2-x3)/4+v*(-x0-x1+x2+x3)/4+uv*(0=x0-x1+x2-x3)/4
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;
409  // y(u,v) = (y0+y1+y2+y3)/4+u*(-y0+y1+y2-y3)/4+v*(-y0-y1+y2+y3)/4+uv*(0=y0-y1+y2-y3)/4
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;
413 
414  // z(u,v) = (z0+z1+z2+z3)/4+u*(-z0+z1+z2-z3)/4+v*(-z0-z1+z2+z3)/4+uv*(0=z0-z1+z2-z3)/4
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;
418  break;
419  }
420  default: {
421  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
422  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): parametrization not defined for the specified face topology.");
423  }
424  }
425  }
426  }
427  }
428 
429  struct Serial {
430 
431  // output:
432  // jacobian (D,sD) - jacobian matrix evaluated at a single point
433  // input:
434  // grads (N,sD) - hgrad basis grad values evaluated at a single point (C1/C2 element only)
435  // nodes (N,D) - cell element-to-node connectivity
436  template<typename jacobianViewType,
437  typename basisGradViewType,
438  typename nodeViewType>
439  KOKKOS_INLINE_FUNCTION
440  static void
441  computeJacobian(const jacobianViewType &jacobian, // D,sD
442  const basisGradViewType &grads, // N,sD
443  const nodeViewType &nodes) { // N,D
444  const auto N = nodes.extent(0);
445 
446  const auto D = jacobian.extent(0);
447  const auto sD = jacobian.extent(1);
448 
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.");
452 
453  Kernels::Serial::gemm_trans_notrans(1.0, nodes, grads, 0.0, jacobian);
454  }
455 
456  // output:
457  // point (D) - mapped physical point
458  // input:
459  // vals (N) - hgrad basis values evaluated at a single point (C1/C2 element only)
460  // nodes (N,D) - cell element-to-node connectivity
461  template<typename pointViewType,
462  typename basisValViewType,
463  typename nodeViewType>
464  KOKKOS_INLINE_FUNCTION
465  static void
466  mapToPhysicalFrame(const pointViewType &point, // D
467  const basisValViewType &vals, // N
468  const nodeViewType &nodes) { // N,D
469  const auto N = vals.extent(0);
470  const auto D = point.extent(0);
471 
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.");
474 
475  Kernels::Serial::gemv_trans(1.0, nodes, vals, 0.0, point);
476  }
477 
478  // template:
479  // implBasisType - impl basis function type e.g., Impl::Basis_HGRAD_QUAD_C1_FEM
480  // output:
481  // xref (sD) - point mapped to reference frame (subcell Dim)
482  // input:
483  // xphy (D) - point in physical frame
484  // nodes (N,D) - cell element-to-node connectivity
485  template<typename implBasisType,
486  typename refPointViewType,
487  typename phyPointViewType,
488  typename nodeViewType>
489  KOKKOS_INLINE_FUNCTION
490  static void
491  mapToReferenceFrame(const refPointViewType &xref, // sD
492  const phyPointViewType &xphy, // D
493  const nodeViewType &nodes) { // N,D
494  const ordinal_type sD = xref.extent(0);
495  const ordinal_type D = xphy.extent(0);
496  const ordinal_type N = nodes.extent(0);
497 
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.");
500 
501  typedef typename refPointViewType::non_const_value_type value_type;
502 
503  // I want to use view instead of dynrankview
504  // NMAX = 28, MAXDIM = 3
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;
509 
510  Kokkos::DynRankView<value_type,
511  Kokkos::Impl::ActiveExecutionMemorySpace,
512  Kokkos::MemoryUnmanaged> vals(ptr, N); ptr += N;
513 
514  Kokkos::DynRankView<value_type,
515  Kokkos::Impl::ActiveExecutionMemorySpace,
516  Kokkos::MemoryUnmanaged> jac(ptr, D, sD); ptr += D*sD;
517 
518  Kokkos::DynRankView<value_type,
519  Kokkos::Impl::ActiveExecutionMemorySpace,
520  Kokkos::MemoryUnmanaged> metric(ptr, sD, sD); ptr += sD*sD;
521 
522  Kokkos::DynRankView<value_type,
523  Kokkos::Impl::ActiveExecutionMemorySpace,
524  Kokkos::MemoryUnmanaged> invmetric(ptr, sD, sD); ptr += sD*sD;
525 
526  Kokkos::DynRankView<value_type,
527  Kokkos::Impl::ActiveExecutionMemorySpace,
528  Kokkos::MemoryUnmanaged> invdf(ptr, sD, D); ptr += sD*D;
529 
530  Kokkos::DynRankView<value_type,
531  Kokkos::Impl::ActiveExecutionMemorySpace,
532  Kokkos::MemoryUnmanaged> xtmp(ptr, sD); ptr += sD;
533 
534  Kokkos::DynRankView<value_type,
535  Kokkos::Impl::ActiveExecutionMemorySpace,
536  Kokkos::MemoryUnmanaged> xold(ptr, sD); ptr += sD;
537 
538  // set initial guess
539  for (ordinal_type j=0;j<D;++j) xold(j) = 0;
540 
541  const double tol = tolerence();
542  for (ordinal_type iter=0;iter<Parameters::MaxNewton;++iter) {
543  // xtmp := F(xold);
544  implBasisType::template Serial<OPERATOR_VALUE>::getValues(vals, xold);
545  mapToPhysicalFrame(xtmp, vals, nodes);
546 
547  // DF^{-1}
548  implBasisType::template Serial<OPERATOR_GRAD>::getValues(grads, xold);
549  CellTools::Serial::computeJacobian(jac, grads, nodes);
550 
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);
554 
555  // Newton
556  Kernels::Serial::z_is_axby(xtmp, 1.0, xphy, -1.0, xtmp); // xtmp := xphy - F(xold);
557  Kernels::Serial::gemv_notrans(1.0, invdf, xtmp, 0.0, xref); // xref := DF^{-1}( xphy - F(xold))
558  Kernels::Serial::z_is_axby(xref, 1.0, xold, 1.0, xref); // xref += xold
559 
560  // l2 error
561  Kernels::Serial::z_is_axby(xtmp, 1.0, xold, -1.0, xref);
562 
563  double err = Kernels::Serial::norm(xtmp, NORM_ONE);
564 
565  if (err < tol)
566  break;
567 
568  Kernels::Serial::copy(xold, xref);
569  }
570  }
571 
575 
576  inline
577  static ConstUnmanagedNodeDataHostView
578  getSubcellParametrization(const ordinal_type subcell_dim,
579  const shards::CellTopology parent_cell) {
580  // all serial interface assumes that they can be called inside parallel for
581  // lazy initilization is not a good idea; init is necessary
582  // setSubcellParametrization();
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;
590  }
591  }
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;
601  }
602  }
603  INTREPID2_TEST_FOR_ABORT(r_val.rank() == 0, "subcell param is not set up before.");
604  return r_val;
605  }
606 
607  template<typename refEdgeTanViewType>
608  inline
609  static void
610  getReferenceEdgeTangent(const refEdgeTanViewType &ref_edge_tangent,
611  const ordinal_type edge_ordinal,
612  const shards::CellTopology &parent_cell ) {
613  const auto edge_map = getSubcellParametrization(1, parent_cell);
614 
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);
618  }
619 
620  template<typename refFaceTanViewType>
621  static void
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) {
626  const auto face_map = getSubcellParametrization(2, parent_cell);
627 
628  // set refFaceTanU -> C_1(*)
629  // set refFaceTanV -> C_2(*)
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);
634  }
635  }
636 
637  template<typename edgeTangentViewType,
638  typename jacobianViewType>
639  inline
640  static void
641  getPhysicalEdgeTangent(const edgeTangentViewType &edge_tangent, // D
642  const jacobianViewType &jacobian, // D, D
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();
647  value_type buf[3];
648  Kokkos::DynRankView<value_type,
649  Kokkos::Impl::ActiveExecutionMemorySpace,
650  Kokkos::MemoryUnmanaged> ref_edge_tangent(&buf[0], D);
651 
652  getReferenceEdgeTangent(ref_edge_tangent, edge_ordinal, parent_cell);
653  Kernels::Serial::matvec_product(edge_tangent, jacobian, ref_edge_tangent);
654  }
655 
656  template<typename faceTanViewType,
657  typename jacobianViewType>
658  inline
659  static void
660  getPhysicalFaceTangent(const faceTanViewType &face_tan_u, // D
661  const faceTanViewType &face_tan_v, // D
662  const jacobianViewType &jacobian, // D, D
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.");
668  value_type buf[6];
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);
672 
673  getReferenceFaceTangent(ref_face_tan_u,
674  ref_face_tan_v,
675  face_ordinal,
676  parent_cell);
677 
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);
680  }
681 
682 
683  template<typename faceNormalViewType,
684  typename jacobianViewType>
685  inline
686  static void
687  getPhysicalFaceNormal(const faceNormalViewType &face_normal, // D
688  const jacobianViewType &jacobian, // D, D
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.");
694  value_type buf[6];
695  Kokkos::DynRankView<value_type,
696  Kokkos::Impl::ActiveExecutionMemorySpace,
697  Kokkos::MemoryUnmanaged> face_tan_u(&buf[0], D), face_tan_v(&buf[3], D);
698 
699  getPhysicalFaceTangent(face_tan_u, face_tan_v,
700  jacobian,
701  face_ordinal,
702  parent_cell);
703  Kernels::Serial::vector_product_d3(face_normal, face_tan_u, face_tan_v);
704  }
705 
706  template<typename sideNormalViewType,
707  typename jacobianViewType>
708  inline
709  static void
710  getPhysicalSideNormal(const sideNormalViewType &side_normal, // D
711  const jacobianViewType &jacobian, // D, D
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;
716  switch (D) {
717  case 2: {
718  value_type buf[3];
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);
725  break;
726  }
727  case 3: {
728  getPhysicalFaceNormal(side_normal, jacobian, side_ordinal, parent_cell);
729  break;
730  }
731  default: {
732  INTREPID2_TEST_FOR_ABORT(true, "cell dimension is out of range.");
733  break;
734  }
735  }
736  }
737  };
738  };
739  }
740 }
741 
742 #endif
743 
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.
static ConstUnmanagedNodeDataHostView getSubcellParametrization(const ordinal_type subcell_dim, const shards::CellTopology parent_cell)
Header file for small functions used in Intrepid2.