OGS 6.2.1-97-g73d1aeda3
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 typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
80  const T_MESH_ELEMENT& /*ele*/,
81  const double* natural_pt,
82  const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
83  T_SHAPE_MATRICES& shapemat,
85 {
86  double* const dNdr = shapemat.dNdr.data();
87  T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
88 }
89 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
90 inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
92  const T_MESH_ELEMENT& /*ele*/,
93  const double* /*natural_pt*/,
94  const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
95  T_SHAPE_MATRICES& /*shapemat*/,
97 {
98 }
99 
100 static void checkJacobianDeterminant(const double detJ,
101  MeshLib::Element const& element)
102 {
103  if (detJ > 0)
104  { // The usual case
105  return;
106  }
107 
108  if (detJ < 0)
109  {
110  ERR("det J = %g is negative for element %d.", detJ, element.getID());
111 #ifndef NDEBUG
112  std::cerr << element << "\n";
113 #endif // NDEBUG
114  OGS_FATAL("Please check the node numbering of the element.");
115  }
116 
117  if (detJ == 0)
118  {
119  ERR("det J is zero for element %d.", element.getID());
120 #ifndef NDEBUG
121  std::cerr << element << "\n";
122 #endif // NDEBUG
123  OGS_FATAL(
124  "Please check whether:\n"
125  "\t the element nodes may have the same coordinates,\n"
126  "\t or the coordinates of all nodes are not given on the x-axis "
127  "for a 1D problem or in the xy-plane in the 2D case.");
128  }
129 }
130 
131 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
132 inline typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
134  const T_MESH_ELEMENT& ele,
135  const double* natural_pt,
136  const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
137  T_SHAPE_MATRICES& shapemat,
139 {
140  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
141  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR>());
142 
143  auto const dim = T_MESH_ELEMENT::dimension;
144  auto const nnodes = T_MESH_ELEMENT::n_all_nodes;
145 
146  //jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
147  for (auto k = decltype(nnodes){0}; k<nnodes; k++) {
148  const MathLib::Point3d& mapped_pt = ele_local_coord.getMappedCoordinates(k);
149  // outer product of dNdr and mapped_pt for a particular node
150  for (auto i_r = decltype(dim){0}; i_r<dim; i_r++) {
151  for (auto j_x = decltype(dim){0}; j_x<dim; j_x++) {
152  shapemat.J(i_r,j_x) += shapemat.dNdr(i_r,k) * mapped_pt[j_x];
153  }
154  }
155  }
156 
157  shapemat.detJ = shapemat.J.determinant();
158  checkJacobianDeterminant(shapemat.detJ, ele);
159 }
160 
161 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
162 inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
164  const T_MESH_ELEMENT& /*ele*/,
165  const double* /*natural_pt*/,
166  const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
167  T_SHAPE_MATRICES& shapemat,
169 {
170  shapemat.detJ = 1.0;
171 }
172 
173 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
175  const T_MESH_ELEMENT& ele,
176  const double* natural_pt,
177  const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
178  T_SHAPE_MATRICES& shapemat,
180 {
181  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
182  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
183  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
184  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
185 }
186 
187 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
188 inline typename std::enable_if<T_SHAPE_FUNC::DIM != 0>::type
190  const T_MESH_ELEMENT& ele,
191  const double* natural_pt,
192  const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
193  T_SHAPE_MATRICES& shapemat,
195 {
196  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
197  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
198 
199  checkJacobianDeterminant(shapemat.detJ, ele);
200 
201  //J^-1, dshape/dx
202  shapemat.invJ.noalias() = shapemat.J.inverse();
203 
204  auto const nnodes(shapemat.dNdr.cols());
205  auto const ele_dim(shapemat.dNdr.rows());
206  assert(shapemat.dNdr.rows()==ele.getDimension());
207  const unsigned global_dim = ele_local_coord.getGlobalDimension();
208  if (global_dim==ele_dim) {
209  shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() = shapemat.invJ * shapemat.dNdr;
210  } else {
211  auto const& matR = ele_local_coord.getRotationMatrixToGlobal(); // 3 x 3
212  auto invJ_dNdr = shapemat.invJ * shapemat.dNdr;
213  auto dshape_global = matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; //3 x nnodes
214  shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);;
215  }
216 }
217 
218 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
219 inline typename std::enable_if<T_SHAPE_FUNC::DIM == 0>::type
221  const T_MESH_ELEMENT& ele,
222  const double* natural_pt,
223  const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
224  T_SHAPE_MATRICES& shapemat,
226 {
227  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
228  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
229 }
230 
231 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
233  const T_MESH_ELEMENT& ele,
234  const double* natural_pt,
235  const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
236  T_SHAPE_MATRICES& shapemat,
238 {
239  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
240  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
241  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
242  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDX>());
243 }
244 
245 template <class T_MESH_ELEMENT,
246  class T_SHAPE_FUNC,
247  class T_SHAPE_MATRICES,
248  ShapeMatrixType T_SHAPE_MATRIX_TYPE>
249 void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
250  const double* natural_pt,
251  T_SHAPE_MATRICES& shapemat,
252  const unsigned global_dim)
253 {
254  const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, global_dim);
255 
257  T_MESH_ELEMENT,
258  T_SHAPE_FUNC,
259  T_SHAPE_MATRICES>
260  (ele,
261  natural_pt,
262  ele_local_coord,
263  shapemat,
265 }
266 
267 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
268  RULE, SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
269  template void naturalCoordinatesMappingComputeShapeMatrices< \
270  MeshLib::TemplateElement<MeshLib::RULE>, \
271  NumLib::SHAPE, \
272  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
273  ShapeMatrixType::WHICHPART>( \
274  MeshLib::TemplateElement<MeshLib::RULE> const&, \
275  double const*, \
276  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
277  const unsigned global_dim)
278 
279 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(RULE, SHAPE) \
280  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
281  RULE, SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
282  /* Those instantiations are needed in unit tests only */ \
283  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
284  RULE, SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
285  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
286  RULE, SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
287  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
288  RULE, SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
289  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
290  RULE, SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
291  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
292  RULE, SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
293 
294 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(RULE, SHAPE, DIM) \
295  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
296  RULE, SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
297  /* Those instantiations are needed in unit tests only */ \
298  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
299  RULE, SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
300  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
301  RULE, SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
302  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
303  RULE, SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
304  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
305  RULE, SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
306  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
307  RULE, SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
308 
313 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PointRule1, ShapePoint1);
314 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule15, ShapePrism15);
315 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule6, ShapePrism6);
316 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule13, ShapePyra13);
317 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule5, ShapePyra5);
325 
326 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 1);
328 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 1);
329 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 1);
330 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 1);
331 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 1);
332 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 1);
333 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 1);
334 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 1);
335 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 1);
336 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 1);
337 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 1);
338 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 1);
342 
343 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 2);
345 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 2);
346 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 2);
347 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 2);
348 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 2);
349 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 2);
350 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 2);
351 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 2);
352 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 2);
353 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 2);
354 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 2);
355 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 2);
359 
360 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 3);
362 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 3);
363 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 3);
364 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 3);
365 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 3);
366 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 3);
367 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 3);
368 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 3);
369 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 3);
370 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 3);
371 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 3);
372 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 3);
376 
377 } // namespace detail
378 
379 } // 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