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