OGS
NaturalCoordinatesMapping.cpp
Go to the documentation of this file.
1
11#include "NaturalCoordinatesMapping.h"
12
13#include <Eigen/LU>
14#include <cassert>
15#ifndef NDEBUG
16#include <iostream>
17#endif // NDEBUG
18
19#include "BaseLib/Error.h"
20#include "MeshLib/ElementCoordinatesMappingLocal.h"
21#include "MeshLib/Elements/HexRule20.h"
22#include "MeshLib/Elements/HexRule8.h"
23#include "MeshLib/Elements/LineRule2.h"
24#include "MeshLib/Elements/LineRule3.h"
25#include "MeshLib/Elements/PointRule1.h"
26#include "MeshLib/Elements/PrismRule15.h"
27#include "MeshLib/Elements/PrismRule6.h"
28#include "MeshLib/Elements/PyramidRule13.h"
29#include "MeshLib/Elements/PyramidRule5.h"
30#include "MeshLib/Elements/QuadRule4.h"
31#include "MeshLib/Elements/QuadRule8.h"
32#include "MeshLib/Elements/QuadRule9.h"
33#include "MeshLib/Elements/TemplateElement.h"
34#include "MeshLib/Elements/TetRule10.h"
35#include "MeshLib/Elements/TetRule4.h"
36#include "MeshLib/Elements/TriRule3.h"
37#include "MeshLib/Elements/TriRule6.h"
38#include "NumLib/Fem/ShapeFunction/ShapeHex20.h"
39#include "NumLib/Fem/ShapeFunction/ShapeHex8.h"
40#include "NumLib/Fem/ShapeFunction/ShapeLine2.h"
41#include "NumLib/Fem/ShapeFunction/ShapeLine3.h"
42#include "NumLib/Fem/ShapeFunction/ShapePoint1.h"
43#include "NumLib/Fem/ShapeFunction/ShapePrism15.h"
44#include "NumLib/Fem/ShapeFunction/ShapePrism6.h"
45#include "NumLib/Fem/ShapeFunction/ShapePyra13.h"
46#include "NumLib/Fem/ShapeFunction/ShapePyra5.h"
47#include "NumLib/Fem/ShapeFunction/ShapeQuad4.h"
48#include "NumLib/Fem/ShapeFunction/ShapeQuad8.h"
49#include "NumLib/Fem/ShapeFunction/ShapeQuad9.h"
50#include "NumLib/Fem/ShapeFunction/ShapeTet10.h"
51#include "NumLib/Fem/ShapeFunction/ShapeTet4.h"
52#include "NumLib/Fem/ShapeFunction/ShapeTri3.h"
53#include "NumLib/Fem/ShapeFunction/ShapeTri6.h"
54#include "NumLib/Fem/ShapeMatrixPolicy.h"
55#include "ShapeMatrices.h"
56
57namespace NumLib
58{
59namespace detail
60{
61template <ShapeMatrixType FIELD_TYPE>
62struct FieldType
63{
64};
65
66template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
67inline void computeMappingMatrices(
68 const MeshLib::Element& /*ele*/,
69 const double* natural_pt,
70 const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
71 T_SHAPE_MATRICES& shapemat,
72 FieldType<ShapeMatrixType::N> /*unused*/)
73{
74 T_SHAPE_FUNC::computeShapeFunction(natural_pt, shapemat.N);
75}
76
77template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
78inline void computeMappingMatrices(
79 const MeshLib::Element& /*ele*/,
80 [[maybe_unused]] const double* natural_pt,
81 const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
82 [[maybe_unused]] T_SHAPE_MATRICES& shapemat,
83 FieldType<ShapeMatrixType::DNDR> /*unused*/)
84{
85 if constexpr (T_SHAPE_FUNC::DIM != 0)
86 {
87 double* const dNdr = shapemat.dNdr.data();
88 T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
89 }
90}
91
92static void checkJacobianDeterminant(const double detJ,
93 MeshLib::Element const& element)
94{
95 if (detJ > 0)
96 { // The usual case
97 return;
98 }
99
100 if (detJ < 0)
101 {
102 ERR("det J = {:g} is negative for element {:d}.",
103 detJ,
104 element.getID());
105#ifndef NDEBUG
106 std::cerr << element << "\n";
107#endif // NDEBUG
108 OGS_FATAL(
109 "Please check whether the node numbering of the element is correct,"
110 "or additional elements (like boundary elements) are still present "
111 "in the mesh.");
112 }
113
114 if (detJ == 0)
115 {
116 ERR("det J is zero for element {:d}.", element.getID());
117#ifndef NDEBUG
118 std::cerr << element << "\n";
119#endif // NDEBUG
120 OGS_FATAL(
121 "Please check whether:\n"
122 "\t the element nodes may have the same coordinates,\n"
123 "\t or the coordinates of all nodes are not given on the x-axis "
124 "for a 1D problem or in the xy-plane in the 2D case.\n"
125 "The first case can occur, if not all boundary elements"
126 "were removed from the bulk mesh.");
127 }
128}
129
130template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
131inline void computeMappingMatrices(
132 const MeshLib::Element& ele,
133 [[maybe_unused]] const double* natural_pt,
134 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
135 T_SHAPE_MATRICES& shapemat,
136 FieldType<ShapeMatrixType::DNDR_J> /*unused*/)
137{
138 if constexpr (T_SHAPE_FUNC::DIM != 0)
139 {
140 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
141 ele,
142 natural_pt,
143 ele_local_coord,
144 shapemat,
145 FieldType<ShapeMatrixType::DNDR>());
146
147 auto const dim = T_SHAPE_FUNC::DIM;
148 auto const nnodes = T_SHAPE_FUNC::NPOINTS;
149
150 // jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
151 for (auto k = decltype(nnodes){0}; k < nnodes; k++)
152 {
153 const MathLib::Point3d& mapped_pt =
154 ele_local_coord.getMappedCoordinates(k);
155 // outer product of dNdr and mapped_pt for a particular node
156 for (auto i_r = decltype(dim){0}; i_r < dim; i_r++)
157 {
158 for (auto j_x = decltype(dim){0}; j_x < dim; j_x++)
159 {
160 shapemat.J(i_r, j_x) +=
161 shapemat.dNdr(i_r, k) * mapped_pt[j_x];
162 }
163 }
164 }
165
166 shapemat.detJ = shapemat.J.determinant();
167 checkJacobianDeterminant(shapemat.detJ, ele);
168 }
169 else
170 {
171 shapemat.detJ = 1.0;
172 }
173}
174
175template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
176inline void computeMappingMatrices(
177 const MeshLib::Element& ele,
178 const double* natural_pt,
179 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
180 T_SHAPE_MATRICES& shapemat,
181 FieldType<ShapeMatrixType::N_J> /*unused*/)
182{
183 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
184 ele,
185 natural_pt,
186 ele_local_coord,
187 shapemat,
188 FieldType<ShapeMatrixType::N>());
189 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
190 ele,
191 natural_pt,
192 ele_local_coord,
193 shapemat,
194 FieldType<ShapeMatrixType::DNDR_J>());
195}
196
197template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
198inline void computeMappingMatrices(
199 const MeshLib::Element& ele,
200 const double* natural_pt,
201 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
202 T_SHAPE_MATRICES& shapemat,
203 FieldType<ShapeMatrixType::DNDX> /*unused*/)
204{
205 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
206 ele,
207 natural_pt,
208 ele_local_coord,
209 shapemat,
210 FieldType<ShapeMatrixType::DNDR_J>());
211 if constexpr (T_SHAPE_FUNC::DIM != 0)
212 {
213 checkJacobianDeterminant(shapemat.detJ, ele);
214
215 // J^-1, dshape/dx
216 shapemat.invJ.noalias() = shapemat.J.inverse();
217
218 assert(shapemat.dNdr.rows() == ele.getDimension());
219 const unsigned global_dim = ele_local_coord.getGlobalDimension();
220 if (global_dim == T_SHAPE_FUNC::DIM)
221 {
222 shapemat.dNdx
223 .template topLeftCorner<T_SHAPE_FUNC::DIM,
224 T_SHAPE_FUNC::NPOINTS>()
225 .noalias() = shapemat.invJ * shapemat.dNdr;
226 }
227 else
228 {
229 auto const& matR =
230 (ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
231 global_dim, T_SHAPE_FUNC::DIM))
232 .eval();
233
234 auto const invJ_dNdr = shapemat.invJ * shapemat.dNdr;
235 auto const dshape_global = matR * invJ_dNdr;
236 shapemat.dNdx =
237 dshape_global.topLeftCorner(global_dim, T_SHAPE_FUNC::NPOINTS);
238 }
239 }
240}
241
242template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
243inline void computeMappingMatrices(
244 const MeshLib::Element& ele,
245 const double* natural_pt,
246 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
247 T_SHAPE_MATRICES& shapemat,
248 FieldType<ShapeMatrixType::ALL> /*unused*/)
249{
250 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
251 ele,
252 natural_pt,
253 ele_local_coord,
254 shapemat,
255 FieldType<ShapeMatrixType::N>());
256 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
257 ele,
258 natural_pt,
259 ele_local_coord,
260 shapemat,
261 FieldType<ShapeMatrixType::DNDX>());
262}
263
264template <class T_SHAPE_FUNC,
265 class T_SHAPE_MATRICES,
266 ShapeMatrixType T_SHAPE_MATRIX_TYPE>
267void naturalCoordinatesMappingComputeShapeMatrices(const MeshLib::Element& ele,
268 const double* natural_pt,
269 T_SHAPE_MATRICES& shapemat,
270 const unsigned global_dim)
271{
272 const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele,
273 global_dim);
274
275 detail::computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
276 ele,
277 natural_pt,
278 ele_local_coord,
279 shapemat,
280 detail::FieldType<T_SHAPE_MATRIX_TYPE>());
281}
282
283#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
284 SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
285 template void naturalCoordinatesMappingComputeShapeMatrices< \
286 NumLib::SHAPE, \
287 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
288 ShapeMatrixType::WHICHPART>( \
289 MeshLib::Element const&, \
290 double const*, \
291 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
292 const unsigned global_dim)
293
294#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE) \
295 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
296 SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
297 /* Those instantiations are needed in unit tests only */ \
298 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
299 SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
300 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
301 SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
302 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
303 SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
304 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
305 SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
306 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
307 SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
308
309#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM) \
310 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
311 SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
312 /* Those instantiations are needed in unit tests only */ \
313 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
314 SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
315 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
316 SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
317 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
318 SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
319 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
320 SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
321 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
322 SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
323
324OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeHex20);
325OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeHex8);
326OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeLine2);
327OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeLine3);
328OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePoint1);
329OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePrism15);
330OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePrism6);
331OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePyra13);
332OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePyra5);
333OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad4);
334OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad8);
335OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad9);
336OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTet10);
337OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTet4);
338OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri3);
339OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri6);
340
341OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 0);
342
343OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 1);
344OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 1);
345OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 1);
346
347OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 2);
348OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 2);
349OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 2);
350OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 2);
351OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 2);
352OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 2);
353OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 2);
354OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 2);
355
356OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex20, 3);
357OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex8, 3);
358OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 3);
359OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 3);
360OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 3);
361OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism15, 3);
362OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism6, 3);
363OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra13, 3);
364OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra5, 3);
365OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 3);
366OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 3);
367OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 3);
368OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet10, 3);
369OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet4, 3);
370OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 3);
371OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 3);
372
373} // namespace detail
374
375} // namespace NumLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
ERR
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM)
Definition NaturalCoordinatesMapping.cpp:309
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE)
Definition NaturalCoordinatesMapping.cpp:294
MeshLib::ElementCoordinatesMappingLocal::getRotationMatrixToGlobal
const RotationMatrix & getRotationMatrixToGlobal() const
return a rotation matrix converting to global coordinates
Definition ElementCoordinatesMappingLocal.h:49
MeshLib::ElementCoordinatesMappingLocal::getMappedCoordinates
MathLib::Point3d const & getMappedCoordinates(std::size_t node_id) const
return mapped coordinates of the node
Definition ElementCoordinatesMappingLocal.h:43
MeshLib::ElementCoordinatesMappingLocal::getGlobalDimension
unsigned getGlobalDimension() const
return the global dimension
Definition ElementCoordinatesMappingLocal.h:40
MeshLib::Element::getDimension
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
MeshLib::Element::getID
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
NumLib::ShapePoint1
Shape function for a point element in natural coordinates.
Definition ShapePoint1.h:19
NumLib::detail::naturalCoordinatesMappingComputeShapeMatrices
void naturalCoordinatesMappingComputeShapeMatrices(const MeshLib::Element &ele, const double *natural_pt, T_SHAPE_MATRICES &shapemat, const unsigned global_dim)
Used to explicitly instantiate the NaturalCoordinatesMapping class template.
Definition NaturalCoordinatesMapping.cpp:267
NumLib::detail::computeMappingMatrices
void computeMappingMatrices(const MeshLib::Element &, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::N >)
Definition NaturalCoordinatesMapping.cpp:67
NumLib::detail::checkJacobianDeterminant
static void checkJacobianDeterminant(const double detJ, MeshLib::Element const &element)
Definition NaturalCoordinatesMapping.cpp:92
NumLib::ShapeMatrixType
ShapeMatrixType
Shape matrix type to be calculated.
Definition ShapeMatrices.h:24