OGS
NaturalCoordinatesMapping.cpp
Go to the documentation of this file.
1
12
13#include <Eigen/LU>
14#include <cassert>
15#ifndef NDEBUG
16#include <iostream>
17#endif // NDEBUG
18
19#include "BaseLib/Error.h"
55#include "ShapeMatrices.h"
56
57namespace NumLib
58{
59namespace detail
60{
61template <ShapeMatrixType FIELD_TYPE>
63{
64};
65
66template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
68 const MeshLib::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
77template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
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,
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>
132 const MeshLib::Element& ele,
133 [[maybe_unused]] const double* natural_pt,
134 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
135 T_SHAPE_MATRICES& shapemat,
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,
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>
177 const MeshLib::Element& ele,
178 const double* natural_pt,
179 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
180 T_SHAPE_MATRICES& shapemat,
182{
183 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
184 ele,
185 natural_pt,
186 ele_local_coord,
187 shapemat,
189 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
190 ele,
191 natural_pt,
192 ele_local_coord,
193 shapemat,
195}
196
197template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
199 const MeshLib::Element& ele,
200 const double* natural_pt,
201 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
202 T_SHAPE_MATRICES& shapemat,
204{
205 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
206 ele,
207 natural_pt,
208 ele_local_coord,
209 shapemat,
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>
244 const MeshLib::Element& ele,
245 const double* natural_pt,
246 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
247 T_SHAPE_MATRICES& shapemat,
249{
250 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
251 ele,
252 natural_pt,
253 ele_local_coord,
254 shapemat,
256 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
257 ele,
258 natural_pt,
259 ele_local_coord,
260 shapemat,
262}
263
264template <class T_SHAPE_FUNC,
265 class T_SHAPE_MATRICES,
266 ShapeMatrixType T_SHAPE_MATRIX_TYPE>
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,
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
340
357
374
391
392} // namespace detail
393
394} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM)
#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE)
const RotationMatrix & getRotationMatrixToGlobal() const
return a rotation matrix converting to global coordinates
MathLib::Point3d const & getMappedCoordinates(std::size_t node_id) const
return mapped coordinates of the node
unsigned getGlobalDimension() const
return the global dimension
virtual constexpr unsigned getDimension() const =0
Get dimension of the mesh element.
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:89
Shape function for a point element in natural coordinates.
Definition ShapePoint1.h:19
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.
void computeMappingMatrices(const MeshLib::Element &, const double *natural_pt, const MeshLib::ElementCoordinatesMappingLocal &, T_SHAPE_MATRICES &shapemat, FieldType< ShapeMatrixType::N >)
static void checkJacobianDeterminant(const double detJ, MeshLib::Element const &element)
ShapeMatrixType
Shape matrix type to be calculated.