Intrepid2
Intrepid2_CellToolsDefParametrization.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_PARAMETRIZATION_HPP__
50 #define __INTREPID2_CELLTOOLS_DEF_PARAMETRIZATION_HPP__
51 
52 // disable clang warnings
53 #if defined (__clang__) && !defined (__INTEL_COMPILER)
54 #pragma clang system_header
55 #endif
56 
57 namespace Intrepid2 {
58 
59  //============================================================================================//
60  // //
61  // Parametrization and hgrad objects caching //
62  // //
63  //============================================================================================//
64 
65 
66  template<typename SpT>
67  void
68  CellTools<SpT>::
69  setSubcellParametrization() {
70  if(isSubcellParametrizationSet_)
71  return;
72 
73  {
74  const auto tet = shards::CellTopology(shards::getCellTopologyData<shards::Tetrahedron<4> >());
75  setSubcellParametrization( subcellParamData_.tetFaces, 2, tet );
76  setSubcellParametrization( subcellParamData_.tetEdges, 1, tet );
77  }
78  {
79  const auto hex = shards::CellTopology(shards::getCellTopologyData<shards::Hexahedron<8> >());
80  setSubcellParametrization( subcellParamData_.hexFaces, 2, hex );
81  setSubcellParametrization( subcellParamData_.hexEdges, 1, hex );
82  }
83  {
84  const auto pyr = shards::CellTopology(shards::getCellTopologyData<shards::Pyramid<5> >());
85  setSubcellParametrization( subcellParamData_.pyrFaces, 2, pyr );
86  setSubcellParametrization( subcellParamData_.pyrEdges, 1, pyr );
87  }
88  {
89  const auto wedge = shards::CellTopology(shards::getCellTopologyData<shards::Wedge<6> >());
90  setSubcellParametrization( subcellParamData_.wedgeFaces, 2, wedge );
91  setSubcellParametrization( subcellParamData_.wedgeEdges, 1, wedge );
92  }
93  {
94  const auto tri = shards::CellTopology(shards::getCellTopologyData<shards::Triangle<3> >());
95  setSubcellParametrization( subcellParamData_.triEdges, 1, tri );
96  }
97  {
98  const auto quad = shards::CellTopology(shards::getCellTopologyData<shards::Quadrilateral<4> >());
99  setSubcellParametrization( subcellParamData_.quadEdges, 1, quad );
100  }
101  {
102  const auto line = shards::CellTopology(shards::getCellTopologyData<shards::ShellLine<2> >());
103  setSubcellParametrization( subcellParamData_.lineEdges, 1, line );
104  }
105 
106  Kokkos::push_finalize_hook( [=] {
107  subcellParamData_.dummy = subcellParamViewType();
108  subcellParamData_.lineEdges = subcellParamViewType();
109  subcellParamData_.triEdges = subcellParamViewType();
110  subcellParamData_.quadEdges = subcellParamViewType();
111  subcellParamData_.shellTriEdges = subcellParamViewType();
112  subcellParamData_.shellQuadEdges = subcellParamViewType();
113  subcellParamData_.tetEdges = subcellParamViewType();
114  subcellParamData_.hexEdges = subcellParamViewType();
115  subcellParamData_.pyrEdges = subcellParamViewType();
116  subcellParamData_.wedgeEdges = subcellParamViewType();
117  subcellParamData_.shellTriFaces = subcellParamViewType();
118  subcellParamData_.shellQuadFaces = subcellParamViewType();
119  subcellParamData_.tetFaces = subcellParamViewType();
120  subcellParamData_.hexFaces = subcellParamViewType();
121  subcellParamData_.pyrFaces = subcellParamViewType();
122  subcellParamData_.wedgeFaces = subcellParamViewType();
123  });
124 
125  isSubcellParametrizationSet_= true;
126  }
127 
128  // template<typename SpT>
129  // const Basis*
130  // CellTools<SpT>::
131  // getHgradBasis( const shards::CellTopology cellTopo ) {
132  // Basis<SpT>* basis = NULL;
133  // switch( cellTopo.getKey() ){
134  // case shards::Line<2>::key: basis = &CachedHgradBasis::c1::line; break;
135  // case shards::Triangle<3>::key: basis = &CachedHgradBasis::c1::tri; break;
136  // case shards::Quadrilateral<4>::key: basis = &CachedHgradBasis::c1::quad; break;
137  // case shards::Tetrahedron<4>::key: basis = &CachedHgradBasis::c1::tet; break;
138  // case shards::Hexahedron<8>::key: basis = &CachedHgradBasis::c1::hex; break;
139  // case shards::Wedge<6>::key: basis = &CachedHgradBasis::c1::wedge; break;
140  // case shards::Pyramid<5>::key: basis = &CachedHgradBasis::c1::pyr; break;
141  // case shards::Triangle<6>::key: basis = &CachedHgradBasis::c2::tri; break;
142  // case shards::Quadrilateral<9>::key: basis = &CachedHgradBasis::c2::quad; break;
143  // case shards::Tetrahedron<10>::key: basis = &CachedHgradBasis::c2::tet; break;
144  // case shards::Hexahedron<27>::key: basis = &CachedHgradBasis::c2::hex; break;
145  // case shards::Wedge<18>::key: basis = &CachedHgradBasis::c2::wedge; break;
146  // default: {
147  // INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
148  // ">>> ERROR (Intrepid2::CellTools::getHgradBasis): Cell topology not supported.");
149  // }
150  // }
151  // }
152 
153  template<typename SpT>
154  void
155  CellTools<SpT>::
156  getSubcellParametrization( subcellParamViewType &subcellParam,
157  const ordinal_type subcellDim,
158  const shards::CellTopology parentCell ) {
159 #ifdef HAVE_INTREPID2_DEBUG
160  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
161  ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): the specified cell topology does not have a reference cell.");
162 #endif
163 
164  if (!isSubcellParametrizationSet_)
165  setSubcellParametrization();
166 
167  // Select subcell parametrization according to its parent cell type
168  const auto pcd = parentCell.getDimension(); // parent cell dim
169  INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 || subcellDim > static_cast<ordinal_type>(pcd-1), std::invalid_argument,
170  ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): Parametrizations defined in a range between 1 and (dim-1)");
171 
172  switch (parentCell.getKey() ) {
173  case shards::Tetrahedron<4>::key:
174  case shards::Tetrahedron<8>::key:
175  case shards::Tetrahedron<10>::key:
176  case shards::Tetrahedron<11>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.tetFaces : subcellParamData_.tetEdges ); break;
177 
178  case shards::Hexahedron<8>::key:
179  case shards::Hexahedron<20>::key:
180  case shards::Hexahedron<27>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.hexFaces : subcellParamData_.hexEdges ); break;
181 
182  case shards::Pyramid<5>::key:
183  case shards::Pyramid<13>::key:
184  case shards::Pyramid<14>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.pyrFaces : subcellParamData_.pyrEdges ); break;
185 
186  case shards::Wedge<6>::key:
187  case shards::Wedge<15>::key:
188  case shards::Wedge<18>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.wedgeFaces : subcellParamData_.wedgeEdges ); break;
189 
190  case shards::Triangle<3>::key:
191  case shards::Triangle<4>::key:
192  case shards::Triangle<6>::key: subcellParam = subcellParamData_.triEdges; break;
193 
194  case shards::Quadrilateral<4>::key:
195  case shards::Quadrilateral<8>::key:
196  case shards::Quadrilateral<9>::key: subcellParam = subcellParamData_.quadEdges; break;
197 
198  // case shards::ShellTriangle<3>::key:
199  // case shards::ShellTriangle<6>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.shellTriFaces : subcellParamData_.shellTriEdges ); break;
200 
201  // case shards::ShellQuadrilateral<4>::key:
202  // case shards::ShellQuadrilateral<8>::key:
203  // case shards::ShellQuadrilateral<9>::key: subcellParam = ( subcellDim == 2 ? subcellParamData_.shellQuadFaces : subcellParamData_.shellQuadEdges ); break;
204 
205  case shards::ShellLine<2>::key:
206  case shards::ShellLine<3>::key:
207  case shards::Beam<2>::key:
208  case shards::Beam<3>::key: subcellParam = subcellParamData_.lineEdges; break;
209  default: {
210  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
211  ">>> ERROR (Intrepid2::CellTools::getSubcellParametrization): invalid cell topology.");
212  }
213  }
214  }
215 
216  template<typename SpT>
217  void
218  CellTools<SpT>::
219  setSubcellParametrization( subcellParamViewType &subcellParam,
220  const ordinal_type subcellDim,
221  const shards::CellTopology parentCell ) {
222 #ifdef HAVE_INTREPID2_DEBUG
223  INTREPID2_TEST_FOR_EXCEPTION( !hasReferenceCell(parentCell), std::invalid_argument,
224  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): the specified cell topology does not have a reference cell.");
225 #endif
226  // subcellParametrization is rank-3 FieldContainer with dimensions (SC, PCD, COEF) where:
227  // - SC is the subcell count of subcells with the specified dimension in the parent cell
228  // - PCD is Parent Cell Dimension, which gives the number of coordinate functions in the map
229  // PCD = 2 for standard 2D cells and non-standard 2D cells: shell line and beam
230  // PCD = 3 for standard 3D cells and non-standard 3D cells: shell Tri and Quad
231  // - COEF is number of coefficients needed to specify a coordinate function:
232  // COEFF = 2 for edge parametrizations
233  // COEFF = 3 for both Quad and Tri face parametrizations. Because all Quad reference faces
234  // are affine, the coefficient of the bilinear term u*v is zero and is not stored, i.e.,
235  // 3 coefficients are sufficient to store Quad face parameterization maps.
236  //
237  // Edge parametrization maps [-1,1] to edge defined by (v0, v1)
238  // Face parametrization maps [-1,1]^2 to quadrilateral face (v0, v1, v2, v3), or
239  // standard 2-simplex {(0,0),(1,0),(0,1)} to traingle face (v0, v1, v2).
240  // This defines orientation-preserving parametrizations with respect to reference edge and
241  // face orientations induced by their vertex order.
242 
243  // get subcellParametrization dimensions: (sc, pcd, coeff)
244  const auto sc = parentCell.getSubcellCount(subcellDim);
245  const auto pcd = parentCell.getDimension();
246  const auto coeff = (subcellDim == 1) ? 2 : 3;
247 
248  INTREPID2_TEST_FOR_EXCEPTION( subcellDim < 1 || subcellDim > static_cast<ordinal_type>(pcd-1), std::invalid_argument,
249  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): Parametrizations defined in a range between 1 and (dim-1)");
250 
251 
252  // create a view
253  subcellParam = subcellParamViewType("CellTools::setSubcellParametrization",
254  sc, pcd, coeff);
255 
256  referenceNodeDataViewType
257  v0("CellTools::setSubcellParametrization::v0", Parameters::MaxDimension),
258  v1("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension),
259  v2("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension),
260  v3("CellTools::setSubcellParametrization::v1", Parameters::MaxDimension);
261 
262  if (subcellDim == 1) {
263  // Edge parametrizations of 2D and 3D cells (shell lines and beams are 2D cells with edges)
264  for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
265  // vertexK[0] = x_k; vertexK[1] = y_k; vertexK[2] = z_k; z_k = 0 for 2D cells
266  // Note that ShellLine and Beam are 2D cells!
267  const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
268  const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
269 
270  getReferenceVertex(v0, parentCell, v0ord);
271  getReferenceVertex(v1, parentCell, v1ord);
272 
273  // x(t) = (x0 + x1)/2 + t*(x1 - x0)/2
274  subcellParam(subcellOrd, 0, 0) = (v0[0] + v1[0])/2.0;
275  subcellParam(subcellOrd, 0, 1) = (v1[0] - v0[0])/2.0;
276 
277  // y(t) = (y0 + y1)/2 + t*(y1 - y0)/2
278  subcellParam(subcellOrd, 1, 0) = (v0[1] + v1[1])/2.0;
279  subcellParam(subcellOrd, 1, 1) = (v1[1] - v0[1])/2.0;
280 
281  if( pcd == 3 ) {
282  // z(t) = (z0 + z1)/2 + t*(z1 - z0)/2
283  subcellParam(subcellOrd, 2, 0) = (v0[2] + v1[2])/2.0;
284  subcellParam(subcellOrd, 2, 1) = (v1[2] - v0[2])/2.0;
285  }
286  }
287  }
288  else if (subcellDim == 2) {
289  // Face parametrizations of 3D cells: (shell Tri and Quad are 3D cells with faces)
290  // A 3D cell can have both Tri and Quad faces, but because they are affine images of the
291  // parametrization domain, 3 coefficients are enough to store them in both cases.
292  for (size_type subcellOrd=0;subcellOrd<sc;++subcellOrd) {
293 
294  switch (parentCell.getKey(subcellDim,subcellOrd)) {
295 
296  case shards::Triangle<3>::key:
297  case shards::Triangle<4>::key:
298  case shards::Triangle<6>::key: {
299  const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
300  const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
301  const auto v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
302 
303  getReferenceVertex(v0, parentCell, v0ord);
304  getReferenceVertex(v1, parentCell, v1ord);
305  getReferenceVertex(v2, parentCell, v2ord);
306 
307  // x(u,v) = x0 + (x1 - x0)*u + (x2 - x0)*v
308  subcellParam(subcellOrd, 0, 0) = v0[0];
309  subcellParam(subcellOrd, 0, 1) = v1[0] - v0[0];
310  subcellParam(subcellOrd, 0, 2) = v2[0] - v0[0];
311 
312  // y(u,v) = y0 + (y1 - y0)*u + (y2 - y0)*v
313  subcellParam(subcellOrd, 1, 0) = v0[1];
314  subcellParam(subcellOrd, 1, 1) = v1[1] - v0[1];
315  subcellParam(subcellOrd, 1, 2) = v2[1] - v0[1];
316 
317  // z(u,v) = z0 + (z1 - z0)*u + (z2 - z0)*v
318  subcellParam(subcellOrd, 2, 0) = v0[2];
319  subcellParam(subcellOrd, 2, 1) = v1[2] - v0[2];
320  subcellParam(subcellOrd, 2, 2) = v2[2] - v0[2];
321  break;
322  }
323  case shards::Quadrilateral<4>::key:
324  case shards::Quadrilateral<8>::key:
325  case shards::Quadrilateral<9>::key: {
326  const auto v0ord = parentCell.getNodeMap(subcellDim, subcellOrd, 0);
327  const auto v1ord = parentCell.getNodeMap(subcellDim, subcellOrd, 1);
328  const auto v2ord = parentCell.getNodeMap(subcellDim, subcellOrd, 2);
329  const auto v3ord = parentCell.getNodeMap(subcellDim, subcellOrd, 3);
330 
331  getReferenceVertex(v0, parentCell, v0ord);
332  getReferenceVertex(v1, parentCell, v1ord);
333  getReferenceVertex(v2, parentCell, v2ord);
334  getReferenceVertex(v3, parentCell, v3ord);
335 
336  // 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
337  subcellParam(subcellOrd, 0, 0) = ( v0[0] + v1[0] + v2[0] + v3[0])/4.0;
338  subcellParam(subcellOrd, 0, 1) = (-v0[0] + v1[0] + v2[0] - v3[0])/4.0;
339  subcellParam(subcellOrd, 0, 2) = (-v0[0] - v1[0] + v2[0] + v3[0])/4.0;
340 
341  // 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
342  subcellParam(subcellOrd, 1, 0) = ( v0[1] + v1[1] + v2[1] + v3[1])/4.0;
343  subcellParam(subcellOrd, 1, 1) = (-v0[1] + v1[1] + v2[1] - v3[1])/4.0;
344  subcellParam(subcellOrd, 1, 2) = (-v0[1] - v1[1] + v2[1] + v3[1])/4.0;
345 
346  // 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
347  subcellParam(subcellOrd, 2, 0) = ( v0[2] + v1[2] + v2[2] + v3[2])/4.0;
348  subcellParam(subcellOrd, 2, 1) = (-v0[2] + v1[2] + v2[2] - v3[2])/4.0;
349  subcellParam(subcellOrd, 2, 2) = (-v0[2] - v1[2] + v2[2] + v3[2])/4.0;
350  break;
351  }
352  default: {
353  INTREPID2_TEST_FOR_EXCEPTION( true, std::invalid_argument,
354  ">>> ERROR (Intrepid2::CellTools::setSubcellParametrization): parametrization not defined for the specified face topology.");
355  }
356  }
357  }
358  }
359  }
360 
361 
362 
363  template<typename SpT>
364  bool
365  CellTools<SpT>::
366  isSubcellParametrizationSet_ = false;
367 
368  template<typename SpT>
369  typename CellTools<SpT>::SubcellParamData
370  CellTools<SpT>::
371  subcellParamData_ = typename CellTools<SpT>::SubcellParamData();
372 
373 
374 }
375 
376 #endif
Intrepid2::Parameters::MaxDimension
static constexpr ordinal_type MaxDimension
The maximum ambient space dimension.
Definition: Intrepid2_Types.hpp:135
Intrepid2::CellTools
A stateless class for operations on cell data. Provides methods for:
Definition: Intrepid2_CellTools.hpp:104
Intrepid2::CellTools::SubcellParamData
Parametrization coefficients of edges and faces of reference cells.
Definition: Intrepid2_CellTools.hpp:255
Intrepid2::CellTools::setSubcellParametrization
static void setSubcellParametrization()
Defines orientation-preserving parametrizations of reference edges and faces of cell topologies with ...
Definition: Intrepid2_CellToolsDefParametrization.hpp:69
Intrepid2::CellTools::getSubcellParametrization
static void getSubcellParametrization(subcellParamViewType &subcellParam, const ordinal_type subcellDim, const shards::CellTopology parentCell)
Returns array with the coefficients of the parametrization maps for the edges or faces of a reference...
Definition: Intrepid2_CellToolsDefParametrization.hpp:156
Generated by   doxygen 1.8.5