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