OGS 6.1.0-1699-ge946d4c5f
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) // The usual case
106  return;
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
133 typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
135  const T_MESH_ELEMENT &ele,
136  const double* natural_pt,
137  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
138  T_SHAPE_MATRICES &shapemat,
140 {
141  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
142  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR>());
143 
144  auto const dim = T_MESH_ELEMENT::dimension;
145  auto const nnodes = T_MESH_ELEMENT::n_all_nodes;
146 
147  //jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
148  for (auto k = decltype(nnodes){0}; k<nnodes; k++) {
149  const MathLib::Point3d& mapped_pt = ele_local_coord.getMappedCoordinates(k);
150  // outer product of dNdr and mapped_pt for a particular node
151  for (auto i_r = decltype(dim){0}; i_r<dim; i_r++) {
152  for (auto j_x = decltype(dim){0}; j_x<dim; j_x++) {
153  shapemat.J(i_r,j_x) += shapemat.dNdr(i_r,k) * mapped_pt[j_x];
154  }
155  }
156  }
157 
158  shapemat.detJ = shapemat.J.determinant();
159  checkJacobianDeterminant(shapemat.detJ, ele);
160 }
161 
162 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
163 inline
164 typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
166  const T_MESH_ELEMENT &/*ele*/,
167  const double* /*natural_pt*/,
168  const MeshLib::ElementCoordinatesMappingLocal &/*ele_local_coord*/,
169  T_SHAPE_MATRICES &shapemat,
171 {
172  shapemat.detJ = 1.0;
173 }
174 
175 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
177  const T_MESH_ELEMENT &ele,
178  const double* natural_pt,
179  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
180  T_SHAPE_MATRICES &shapemat,
182 {
183  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
184  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
185  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
186  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
187 }
188 
189 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
190 inline
191 typename std::enable_if<T_SHAPE_FUNC::DIM!=0>::type
193  const T_MESH_ELEMENT &ele,
194  const double* natural_pt,
195  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
196  T_SHAPE_MATRICES &shapemat,
198 {
199  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
200  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
201 
202  checkJacobianDeterminant(shapemat.detJ, ele);
203 
204  //J^-1, dshape/dx
205  shapemat.invJ.noalias() = shapemat.J.inverse();
206 
207  auto const nnodes(shapemat.dNdr.cols());
208  auto const ele_dim(shapemat.dNdr.rows());
209  assert(shapemat.dNdr.rows()==ele.getDimension());
210  const unsigned global_dim = ele_local_coord.getGlobalDimension();
211  if (global_dim==ele_dim) {
212  shapemat.dNdx.topLeftCorner(ele_dim, nnodes).noalias() = shapemat.invJ * shapemat.dNdr;
213  } else {
214  auto const& matR = ele_local_coord.getRotationMatrixToGlobal(); // 3 x 3
215  auto invJ_dNdr = shapemat.invJ * shapemat.dNdr;
216  auto dshape_global = matR.topLeftCorner(3u, ele_dim) * invJ_dNdr; //3 x nnodes
217  shapemat.dNdx = dshape_global.topLeftCorner(global_dim, nnodes);;
218  }
219 }
220 
221 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
222 inline
223 typename std::enable_if<T_SHAPE_FUNC::DIM==0>::type
225  const T_MESH_ELEMENT &ele,
226  const double* natural_pt,
227  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
228  T_SHAPE_MATRICES &shapemat,
230 {
231  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
232  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDR_J>());
233 }
234 
235 
236 template <class T_MESH_ELEMENT, class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
238  const T_MESH_ELEMENT &ele,
239  const double* natural_pt,
240  const MeshLib::ElementCoordinatesMappingLocal &ele_local_coord,
241  T_SHAPE_MATRICES &shapemat,
243 {
244  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
245  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::N>());
246  computeMappingMatrices<T_MESH_ELEMENT, T_SHAPE_FUNC, T_SHAPE_MATRICES>
247  (ele, natural_pt, ele_local_coord, shapemat, FieldType<ShapeMatrixType::DNDX>());
248 }
249 
250 template <class T_MESH_ELEMENT,
251  class T_SHAPE_FUNC,
252  class T_SHAPE_MATRICES,
253  ShapeMatrixType T_SHAPE_MATRIX_TYPE>
254 void naturalCoordinatesMappingComputeShapeMatrices(const T_MESH_ELEMENT& ele,
255  const double* natural_pt,
256  T_SHAPE_MATRICES& shapemat,
257  const unsigned global_dim)
258 {
259  const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele, global_dim);
260 
262  T_MESH_ELEMENT,
263  T_SHAPE_FUNC,
264  T_SHAPE_MATRICES>
265  (ele,
266  natural_pt,
267  ele_local_coord,
268  shapemat,
270 }
271 
272 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
273  RULE, SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
274  template void naturalCoordinatesMappingComputeShapeMatrices< \
275  MeshLib::TemplateElement<MeshLib::RULE>, \
276  NumLib::SHAPE, \
277  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
278  ShapeMatrixType::WHICHPART>( \
279  MeshLib::TemplateElement<MeshLib::RULE> const&, \
280  double const*, \
281  SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
282  const unsigned global_dim)
283 
284 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(RULE, SHAPE) \
285  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
286  RULE, SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
287  /* Those instantiations are needed in unit tests only */ \
288  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
289  RULE, SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
290  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
291  RULE, SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
292  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
293  RULE, SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
294  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
295  RULE, SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
296  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
297  RULE, SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
298 
299 #define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(RULE, SHAPE, DIM) \
300  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
301  RULE, SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
302  /* Those instantiations are needed in unit tests only */ \
303  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
304  RULE, SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
305  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
306  RULE, SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
307  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
308  RULE, SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
309  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
310  RULE, SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
311  OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
312  RULE, SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
313 
318 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PointRule1, ShapePoint1);
319 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule15, ShapePrism15);
320 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PrismRule6, ShapePrism6);
321 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule13, ShapePyra13);
322 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(PyramidRule5, ShapePyra5);
330 
331 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 1);
333 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 1);
334 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 1);
335 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 1);
336 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 1);
337 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 1);
338 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 1);
339 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 1);
340 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 1);
341 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 1);
342 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 1);
343 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 1);
347 
348 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 2);
350 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 2);
351 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 2);
352 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 2);
353 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 2);
354 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 2);
355 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 2);
356 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 2);
357 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 2);
358 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 2);
359 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 2);
360 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 2);
364 
365 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(HexRule20, ShapeHex20, 3);
367 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule2, ShapeLine2, 3);
368 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(LineRule3, ShapeLine3, 3);
369 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PointRule1, ShapePoint1, 3);
370 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule15, ShapePrism15, 3);
371 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PrismRule6, ShapePrism6, 3);
372 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule13, ShapePyra13, 3);
373 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(PyramidRule5, ShapePyra5, 3);
374 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule4, ShapeQuad4, 3);
375 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule8, ShapeQuad8, 3);
376 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(QuadRule9, ShapeQuad9, 3);
377 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(TetRule10, ShapeTet10, 3);
381 
382 } // detail
383 
384 } // 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:71
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