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 
71  struct ReferenceNodeDataType {
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 
241  struct SubcellParamDataType {
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 
Intrepid2::Impl::CellTools::Serial::getSubcellParametrization
static ConstUnmanagedNodeDataHostView getSubcellParametrization(const ordinal_type subcell_dim, const shards::CellTopology parent_cell)
Definition: Intrepid2_CellTools_Serial.hpp:578
Intrepid2::Impl::CellTools
See Intrepid2::CellTools.
Definition: Intrepid2_CellTools_Serial.hpp:66
Intrepid2::Impl::CellTools::SubcellParamDataType
Definition: Intrepid2_CellTools_Serial.hpp:241
Intrepid2::Parameters::MaxNewton
static constexpr ordinal_type MaxNewton
Maximum number of Newton iterations used internally in methods such as computing the action of the in...
Definition: Intrepid2_Types.hpp:137
Intrepid2_Utils.hpp
Header function for Intrepid2::Util class and other utility functions.
Intrepid2_Kernels.hpp
Header file for small functions used in Intrepid2.
Intrepid2::Impl::CellTools::ReferenceNodeDataType
Definition: Intrepid2_CellTools_Serial.hpp:71
Intrepid2_Types.hpp
Contains definitions of custom data types in Intrepid2.
Intrepid2::Impl::CellTools::Serial
Definition: Intrepid2_CellTools_Serial.hpp:429
Generated by   doxygen 1.8.16