OGS 6.2.0-97-g4a610c866
NaturalCoordinatesMapping.cpp
Go to the documentation of this file.
1 
11 
12 #include <cassert>
13 #ifndef NDEBUG
14 #include <iostream>
15 #endif // NDEBUG
16 
17 #include "BaseLib/Error.h"
18 
37 
55 
56 #include "ShapeMatrices.h"
57 
58 namespace NumLib
59 {
60 
61 namespace detail
62 {
63 
64 template <ShapeMatrixType FIELD_TYPE> struct FieldType {};
65 
66 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
68  const T_MESH_ELEMENT &/*ele*/,
69  const double* natural_pt,
70  const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
71  T_SHAPE_MATRICES &shapemat,
73 {
74  T_SHAPE_FUNC::computeShapeFunction(natural_pt, shapemat.N);
75 }
76 
77 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
78 inline
79 typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
81  const T_MESH_ELEMENT &/*ele*/,
82  const double* natural_pt,
83  const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
84  T_SHAPE_MATRICES &shapemat,
86 {
87  double* const dNdr = shapemat.dNdr.data();
88  T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
89 }
90 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
91 inline
92 typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
94  const T_MESH_ELEMENT &/*ele*/,
95  const double* /*natural_pt*/,
96  const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
97  T_SHAPE_MATRICES &/*shapemat*/,
99 {
100 }
101 
102 static void checkJacobianDeterminant(const double detJ,
103  MeshLib::Element const& element)
104 {
105  if (detJ > 0)
106  { // The usual case
107  return;
108  }
109 
110  if (detJ < 0)
111  {
112  ERR("det J = %g is negative for element %d.", detJ, element.getID());
113 #ifndef NDEBUG
114  std::cerr << element << "\n";
115 #endif // NDEBUG
116  OGS_FATAL("Please check the node numbering of the element.");
117  }
118 
119  if (detJ == 0)
120  {
121  ERR("det J is zero for element %d.", element.getID());
122 #ifndef NDEBUG
123  std::cerr << element << "\n";
124 #endif // NDEBUG
125  OGS_FATAL(
126  "Please check whether:\n"
127  "\t the element nodes may have the same coordinates,\n"
128  "\t or the coordinates of all nodes are not given on the x-axis "
129  "for a 1D problem or in the xy-plane in the 2D case.");
130  }
131 }
132 
133 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
134 inline
135 typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
137  const T_MESH_ELEMENT &ele,
138  const double* natural_pt,
139  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
140  T_SHAPE_MATRICES &shapemat,
142 {
143  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
144  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR>());
145 
146  auto const dim = T_MESH_ELEMENT::dimension;
147  auto const nnodes = T_MESH_ELEMENT::n_all_nodes;
148 
149  //jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
150  for (auto k = decltype(nnodes){0}; k<nnodes; k++) {
151  const MathLib::Point3d& mapped_pt = ele_local_coord.getMappedCoordinates(k);
152  // outer product of dNdr and mapped_pt for a particular node
153  for (auto i_r = decltype(dim){0}; i_r<dim; i_r++) {
154  for (auto j_x = decltype(dim){0}; j_x<dim; j_x++) {
155  shapemat.J(i_r,j_x) += shapemat.dNdr(i_r,k) * mapped_pt[j_x];
156  }
157  }
158  }
159 
160  shapemat.detJ = shapemat.J.determinant();
161  checkJacobianDeterminant(shapemat.detJ, ele);
162 }
163 
164 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
165 inline
166 typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
168  const T_MESH_ELEMENT &/*ele*/,
169  const double* /*natural_pt*/,
170  const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
171  T_SHAPE_MATRICES &shapemat,
173 {
174  shapemat.detJ = 1.0;
175 }
176 
177 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
179  const T_MESH_ELEMENT &ele,
180  const double* natural_pt,
181  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
182  T_SHAPE_MATRICES &shapemat,
184 {
185  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
186  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
187  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
188  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
189 }
190 
191 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
192 inline
193 typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
195  const T_MESH_ELEMENT &ele,
196  const double* natural_pt,
197  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
198  T_SHAPE_MATRICES &shapemat,
200 {
201  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
202  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
203 
204  checkJacobianDeterminant(shapemat.detJ, ele);
205 
206  //J^-1, dshape/dx
207  shapemat.invJ.noalias() = shapemat.J.inverse();
208 
209  auto const nnodes(shapemat.dNdr.cols());
210  auto const ele_dim(shapemat.dNdr.rows());
211  assert(shapemat.dNdr.rows()==ele.getDimension());
212  const unsigned global_dim = ele_local_coord.getGlobalDimension();
213  if (global_dim==ele_dim) {
214  shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() = shapemat.invJ * shapemat.dNdr;
215  } else {
216  auto const& matR = ele_local_coord.getRotationMatrixToGlobal(); // 3 x 3
217  auto invJ_dNdr = shapemat.invJ * shapemat.dNdr;
218  auto dshape_global = matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; //3 x nnodes
219  shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);;
220  }
221 }
222 
223 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
224 inline
225 typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
227  const T_MESH_ELEMENT &ele,
228  const double* natural_pt,
229  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
230  T_SHAPE_MATRICES &shapemat,
232 {
233  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
234  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
235 }
236 
237 
238 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
240  const T_MESH_ELEMENT &ele,
241  const double* natural_pt,
242  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
243  T_SHAPE_MATRICES &shapemat,
245 {
246  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
247  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
248  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
249  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDX>());
250 }
251 
252 template <class T_MESH_ELEMENT,
253  class T_SHAPE_FUNC,
254  class T_SHAPE_MATRICES,
255  ShapeMatrixType T_SHAPE_MATRIX_TYPE>
256 void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
257  const double* natural_pt,
258  T_SHAPE_MATRICES& shapemat,
259  const unsigned global_dim)
260 {
261  const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, global_dim);
262 
264  T_MESH_ELEMENT,
265  T_SHAPE_FUNC,
266  T_SHAPE_MATRICES>
267  (ele,
268  natural_pt,
269  ele_local_coord,
270  shapemat,
272 }
273 
274 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
275  RULE, SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
276  template void naturalCoordinatesMappingComputeShapeMatrices< \
277  MeshLib::TemplateElement<MeshLib::RULE>, \
278  NumLib::SHAPE, \
279  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
280  ShapeMatrixType::WHICHPART>( \
281  MeshLib::TemplateElement<MeshLib::RULE> const&, \
282  double const*, \
283  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
284  const unsigned global_dim)
285 
286 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(RULE, SHAPE) \
287  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
288  RULE, SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
289  /* Those instantiations are needed in unit tests only */ \
290  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
291  RULE, SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
292  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
293  RULE, SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
294  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
295  RULE, SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
296  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
297  RULE, SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
298  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
299  RULE, SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
300 
301 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(RULE, SHAPE, DIM) \
302  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
303  RULE, SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
304  /* Those instantiations are needed in unit tests only */ \
305  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
306  RULE, SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
307  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
308  RULE, SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
309  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
310  RULE, SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
311  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
312  RULE, SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
313  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
314  RULE, SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
315 
320 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PointRule1, ShapePoint1);
321 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule15, ShapePrism15);
322 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule6, ShapePrism6);
323 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule13, ShapePyra13);
324 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule5, ShapePyra5);
332 
333 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 1);
335 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 1);
336 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 1);
337 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 1);
338 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 1);
339 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 1);
340 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 1);
341 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 1);
342 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 1);
343 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 1);
344 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 1);
345 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 1);
349 
350 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 2);
352 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 2);
353 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 2);
354 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 2);
355 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 2);
356 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 2);
357 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 2);
358 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 2);
359 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 2);
360 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 2);
361 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 2);
362 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 2);
366 
367 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 3);
369 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 3);
370 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 3);
371 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 3);
372 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 3);
373 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 3);
374 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 3);
375 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 3);
376 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 3);
377 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 3);
378 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 3);
379 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 3);
383 
384 } // namespace detail
385 
386 } // namespace NumLib
static void checkJacobianDeterminant(const double detJ, MeshLib::Element const &element)
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 1)
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(HexRule20, ShapeHex20)
void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT &ele, const double *natural_pt, T_SHAPE_MATRICES &shapemat, const unsigned global_dim)
Used to explicitly instantiate the NaturalCoordinatesMapping class template.
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:90
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
void computeMappingMatrices(const T_MESH_ELEMENT &, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::N >)
ShapeMatrixType
Shape matrix type to be calculated.
Definition: ShapeMatrices.h:26