OGS
NaturalCoordinatesMapping.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14#ifndef NDEBUG
15#include <iostream>
16#endif // NDEBUG
17
18#include "BaseLib/Error.h"
54#include "ShapeMatrices.h"
55
56namespace NumLib
57{
58namespace detail
59{
60template <ShapeMatrixType FIELD_TYPE>
62{
63};
64
65template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
67 const MeshLib::Element& /*ele*/,
68 const double* natural_pt,
69 const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
70 T_SHAPE_MATRICES& shapemat,
72{
73 T_SHAPE_FUNC::computeShapeFunction(natural_pt, shapemat.N);
74}
75
76template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
78 const MeshLib::Element& /*ele*/,
79 [[maybe_unused]] const double* natural_pt,
80 const MeshLib::ElementCoordinatesMappingLocal& /*ele_local_coord*/,
81 [[maybe_unused]] T_SHAPE_MATRICES& shapemat,
83{
84 if constexpr (T_SHAPE_FUNC::DIM != 0)
85 {
86 double* const dNdr = shapemat.dNdr.data();
87 T_SHAPE_FUNC::computeGradShapeFunction(natural_pt, dNdr);
88 }
89}
90
91static void checkJacobianDeterminant(const double detJ,
92 MeshLib::Element const& element)
93{
94 if (detJ > 0)
95 { // The usual case
96 return;
97 }
98
99 if (detJ < 0)
100 {
101 ERR("det J = {:g} is negative for element {:d}.",
102 detJ,
103 element.getID());
104#ifndef NDEBUG
105 std::cerr << element << "\n";
106#endif // NDEBUG
107 OGS_FATAL(
108 "Please check whether the node numbering of the element is correct,"
109 "or additional elements (like boundary elements) are still present "
110 "in the mesh.");
111 }
112
113 if (detJ == 0)
114 {
115 ERR("det J is zero for element {:d}.", element.getID());
116#ifndef NDEBUG
117 std::cerr << element << "\n";
118#endif // NDEBUG
119 OGS_FATAL(
120 "Please check whether:\n"
121 "\t the element nodes may have the same coordinates,\n"
122 "\t or the coordinates of all nodes are not given on the x-axis "
123 "for a 1D problem or in the xy-plane in the 2D case.\n"
124 "The first case can occur, if not all boundary elements"
125 "were removed from the bulk mesh.");
126 }
127}
128
129template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
131 const MeshLib::Element& ele,
132 [[maybe_unused]] const double* natural_pt,
133 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
134 T_SHAPE_MATRICES& shapemat,
136{
137 if constexpr (T_SHAPE_FUNC::DIM != 0)
138 {
139 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
140 ele,
141 natural_pt,
142 ele_local_coord,
143 shapemat,
145
146 auto const dim = T_SHAPE_FUNC::DIM;
147 auto const 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 // outer product of dNdr and mapped_pt for a particular node
155 for (auto i_r = decltype(dim){0}; i_r < dim; i_r++)
156 {
157 for (auto j_x = decltype(dim){0}; j_x < dim; j_x++)
158 {
159 shapemat.J(i_r, j_x) +=
160 shapemat.dNdr(i_r, k) * mapped_pt[j_x];
161 }
162 }
163 }
164
165 shapemat.detJ = shapemat.J.determinant();
166 checkJacobianDeterminant(shapemat.detJ, ele);
167 }
168 else
169 {
170 shapemat.detJ = 1.0;
171 }
172}
173
174template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
176 const MeshLib::Element& ele,
177 const double* natural_pt,
178 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
179 T_SHAPE_MATRICES& shapemat,
181{
182 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
183 ele,
184 natural_pt,
185 ele_local_coord,
186 shapemat,
188 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
189 ele,
190 natural_pt,
191 ele_local_coord,
192 shapemat,
194}
195
196template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
198 const MeshLib::Element& ele,
199 const double* natural_pt,
200 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
201 T_SHAPE_MATRICES& shapemat,
203{
204 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
205 ele,
206 natural_pt,
207 ele_local_coord,
208 shapemat,
210 if constexpr (T_SHAPE_FUNC::DIM != 0)
211 {
212 checkJacobianDeterminant(shapemat.detJ, ele);
213
214 // J^-1, dshape/dx
215 shapemat.invJ.noalias() = shapemat.J.inverse();
216
217 assert(shapemat.dNdr.rows() == ele.getDimension());
218 const unsigned global_dim = ele_local_coord.getGlobalDimension();
219 if (global_dim == T_SHAPE_FUNC::DIM)
220 {
221 shapemat.dNdx
222 .template topLeftCorner<T_SHAPE_FUNC::DIM,
223 T_SHAPE_FUNC::NPOINTS>()
224 .noalias() = shapemat.invJ * shapemat.dNdr;
225 }
226 else
227 {
228 auto const& matR =
229 (ele_local_coord.getRotationMatrixToGlobal().topLeftCorner(
230 global_dim, T_SHAPE_FUNC::DIM))
231 .eval();
232
233 auto const invJ_dNdr = shapemat.invJ * shapemat.dNdr;
234 auto const dshape_global = matR * invJ_dNdr;
235 shapemat.dNdx =
236 dshape_global.topLeftCorner(global_dim, T_SHAPE_FUNC::NPOINTS);
237 }
238 }
239}
240
241template <class T_SHAPE_FUNC, class T_SHAPE_MATRICES>
243 const MeshLib::Element& ele,
244 const double* natural_pt,
245 const MeshLib::ElementCoordinatesMappingLocal& ele_local_coord,
246 T_SHAPE_MATRICES& shapemat,
248{
249 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
250 ele,
251 natural_pt,
252 ele_local_coord,
253 shapemat,
255 computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
256 ele,
257 natural_pt,
258 ele_local_coord,
259 shapemat,
261}
262
263template <class T_SHAPE_FUNC,
264 class T_SHAPE_MATRICES,
265 ShapeMatrixType T_SHAPE_MATRIX_TYPE>
267 const double* natural_pt,
268 T_SHAPE_MATRICES& shapemat,
269 const unsigned global_dim)
270{
271 const MeshLib::ElementCoordinatesMappingLocal ele_local_coord(ele,
272 global_dim);
273
274 detail::computeMappingMatrices<T_SHAPE_FUNC, T_SHAPE_MATRICES>(
275 ele,
276 natural_pt,
277 ele_local_coord,
278 shapemat,
280}
281
282#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
283 SHAPE, DIM, WHICHPART, SHAPEMATRIXPOLICY) \
284 template void naturalCoordinatesMappingComputeShapeMatrices< \
285 NumLib::SHAPE, \
286 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices, \
287 ShapeMatrixType::WHICHPART>( \
288 MeshLib::Element const&, \
289 double const*, \
290 SHAPEMATRIXPOLICY<NumLib::SHAPE, DIM>::ShapeMatrices&, \
291 const unsigned global_dim)
292
293#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(SHAPE) \
294 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
295 SHAPE, 0, ALL, EigenDynamicShapeMatrixPolicy); \
296 /* Those instantiations are needed in unit tests only */ \
297 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
298 SHAPE, 0, N, EigenDynamicShapeMatrixPolicy); \
299 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
300 SHAPE, 0, DNDR, EigenDynamicShapeMatrixPolicy); \
301 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
302 SHAPE, 0, N_J, EigenDynamicShapeMatrixPolicy); \
303 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
304 SHAPE, 0, DNDR_J, EigenDynamicShapeMatrixPolicy); \
305 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
306 SHAPE, 0, DNDX, EigenDynamicShapeMatrixPolicy)
307
308#define OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(SHAPE, DIM) \
309 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
310 SHAPE, DIM, ALL, EigenFixedShapeMatrixPolicy); \
311 /* Those instantiations are needed in unit tests only */ \
312 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
313 SHAPE, DIM, N, EigenFixedShapeMatrixPolicy); \
314 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
315 SHAPE, DIM, DNDR, EigenFixedShapeMatrixPolicy); \
316 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
317 SHAPE, DIM, N_J, EigenFixedShapeMatrixPolicy); \
318 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
319 SHAPE, DIM, DNDR_J, EigenFixedShapeMatrixPolicy); \
320 OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_PART( \
321 SHAPE, DIM, DNDX, EigenFixedShapeMatrixPolicy)
322
339
356
373
390
391} // namespace detail
392
393} // namespace NumLib
#define OGS_FATAL(...)
Definition: Error.h:26
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
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.
virtual std::size_t getID() const final
Returns the ID of the element.
Definition: Element.h:89
Shape function for a point element in natural coordinates.
Definition: ShapePoint1.h:19
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_FIX(ShapeHex20, 1)
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 >)
OGS_INSTANTIATE_NATURAL_COORDINATES_MAPPING_DYN(ShapeHex20)
static void checkJacobianDeterminant(const double detJ, MeshLib::Element const &element)
ShapeMatrixType
Shape matrix type to be calculated.
Definition: ShapeMatrices.h:24