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 constexpr nnodes = T_SHAPE_FUNC::NPOINTS;
148
149 // jacobian: J=[dx/dr dy/dr // dx/ds dy/ds]
150 for (auto k = decltype(nnodes){0}; k < nnodes; k++)
151 {
152 const MathLib::Point3d& mapped_pt =
153 ele_local_coord.getMappedCoordinates(k);
154 shapemat.J +=
155 shapemat.dNdr.col(k) *
156 mapped_pt.asEigenVector3d().head(T_SHAPE_FUNC::DIM).transpose();
157 }
158
159 shapemat.detJ = shapemat.J.determinant();
160 checkJacobianDeterminant(shapemat.detJ, ele);
161 }
162 else
163 {
164 shapemat.detJ = 1.0;
165 }
166}
167
168template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
169inline void computeMappingMatrices(
170 const MeshLib::Element& ele,
171 const double* natural_pt,
172 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
173 T_SHAPE_MATRICES& shapemat,
174 FieldType<ShapeMatrixType::N_J> /*unused*/)
175{
176 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
177 ele,
178 natural_pt,
179 ele_local_coord,
180 shapemat,
181 FieldType<ShapeMatrixType::N>());
182 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
183 ele,
184 natural_pt,
185 ele_local_coord,
186 shapemat,
187 FieldType<ShapeMatrixType::DNDR_J>());
188}
189
190template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
191inline void computeMappingMatrices(
192 const MeshLib::Element& ele,
193 const double* natural_pt,
194 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
195 T_SHAPE_MATRICES& shapemat,
196 FieldType<ShapeMatrixType::DNDX> /*unused*/)
197{
198 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
199 ele,
200 natural_pt,
201 ele_local_coord,
202 shapemat,
203 FieldType<ShapeMatrixType::DNDR_J>());
204 if constexpr (T_SHAPE_FUNC::DIM != 0)
205 {
206 checkJacobianDeterminant(shapemat.detJ, ele);
207
208 // J^-1, dshape/dx
209 shapemat.invJ.noalias() = shapemat.J.inverse();
210
211 assert(shapemat.dNdr.rows() == ele.getDimension());
212 const unsigned global_dim = ele_local_coord.getGlobalDimension();
213 if (global_dim == T_SHAPE_FUNC::DIM)
214 {
215 shapemat.dNdx
216 .template topLeftCorner<T_SHAPE_FUNC::DIM,
217 T_SHAPE_FUNC::NPOINTS>()
218 .noalias() = shapemat.invJ * shapemat.dNdr;
219 }
220 else
221 {
222 auto const& matR =
223 (ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
224 global_dim, T_SHAPE_FUNC::DIM))
225 .eval();
226
227 auto const invJ_dNdr = shapemat.invJ * shapemat.dNdr;
228 auto const dshape_global = matR * invJ_dNdr;
229 shapemat.dNdx =
230 dshape_global.topLeftCorner(global_dim, T_SHAPE_FUNC::NPOINTS);
231 }
232 }
233}
234
235template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
236inline void computeMappingMatrices(
237 const MeshLib::Element& ele,
238 const double* natural_pt,
239 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
240 T_SHAPE_MATRICES& shapemat,
241 FieldType<ShapeMatrixType::ALL> /*unused*/)
242{
243 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
244 ele,
245 natural_pt,
246 ele_local_coord,
247 shapemat,
248 FieldType<ShapeMatrixType::N>());
249 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
250 ele,
251 natural_pt,
252 ele_local_coord,
253 shapemat,
254 FieldType<ShapeMatrixType::DNDX>());
255}
256
257template <class T_SHAPE_FUNC,
258 class T_SHAPE_MATRICES,
259 ShapeMatrixType T_SHAPE_MATRIX_TYPE>
260void naturalCoordinatesMappingComputeShapeMatrices(const MeshLib::Element& ele,
261 const double* natural_pt,
262 T_SHAPE_MATRICES& shapemat,
263 const unsigned global_dim)
264{
265 const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele,
266 global_dim);
267
268 detail::computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
269 ele,
270 natural_pt,
271 ele_local_coord,
272 shapemat,
273 detail::FieldType<T_SHAPE_MATRIX_TYPE>());
274}
275
276#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
277 SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
278 template void naturalCoordinatesMappingComputeShapeMatrices< \
279 NumLib::SHAPE, \
280 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
281 ShapeMatrixType::WHICHPART>( \
282 MeshLib::Element const&, \
283 double const*, \
284 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
285 const unsigned global_dim)
286
287#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE) \
288 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
289 SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
290 /* Those instantiations are needed in unit tests only */ \
291 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
292 SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
293 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
294 SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
295 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
296 SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
297 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
298 SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
299 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
300 SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
301
302#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM) \
303 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
304 SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
305 /* Those instantiations are needed in unit tests only */ \
306 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
307 SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
308 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
309 SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
310 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
311 SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
312 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
313 SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
314 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
315 SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
316
317OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeHex20);
318OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeHex8);
319OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeLine2);
320OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeLine3);
321OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePoint1);
322OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePrism15);
323OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePrism6);
324OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePyra13);
325OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapePyra5);
326OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad4);
327OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad8);
328OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeQuad9);
329OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTet10);
330OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTet4);
331OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri3);
332OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeTri6);
333
334OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 0);
335
336OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 1);
337OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 1);
338OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 1);
339
340OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 2);
341OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 2);
342OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 2);
343OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 2);
344OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 2);
345OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 2);
346OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 2);
347OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 2);
348
349OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex20, 3);
350OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex8, 3);
351OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine2, 3);
352OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeLine3, 3);
353OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePoint1, 3);
354OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism15, 3);
355OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePrism6, 3);
356OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra13, 3);
357OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapePyra5, 3);
358OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad4, 3);
359OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad8, 3);
360OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeQuad9, 3);
361OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet10, 3);
362OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTet4, 3);
363OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri3, 3);
364OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeTri6, 3);
365
366} // namespace detail
367
368} // namespace NumLib
OGS_FATAL
#define OGS_FATAL(...)
Definition Error.h:26
ERR
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM)
Definition NaturalCoordinatesMapping.cpp:302
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE)
Definition NaturalCoordinatesMapping.cpp:287
MathLib::Point3d::asEigenVector3d
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
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:91
NumLib::ShapePoint1
Shape function for a point element in natural coordinates.
Definition ShapePoint1.h:20
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:260
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