OGS
NodeReordering.cpp
Go to the documentation of this file.
1
11
12#include <tclap/CmdLine.h>
13
14#include <Eigen/Dense>
15#include <algorithm>
16#include <array>
17#include <concepts>
18#include <functional>
19#include <vector>
20
21#include "BaseLib/Algorithm.h"
22#include "BaseLib/Logging.h"
23#include "BaseLib/MPI.h"
25#include "InfoLib/GitInfo.h"
30#include "MeshLib/Mesh.h"
31#include "MeshLib/MeshEnums.h"
32#include "MeshLib/Node.h"
50
51template <typename GradShapeFunction, int Dim>
53 int const mesh_space_dimension,
54 std::array<double, Dim> const& xi,
55 bool const check_reordered = false)
56{
57 Eigen::Matrix<double, GradShapeFunction::DIM, GradShapeFunction::NPOINTS,
58 Eigen::RowMajor>
59 dNdxi;
60 Eigen::Map<Eigen::VectorXd> dN_vec(dNdxi.data(), dNdxi.size());
61 GradShapeFunction::computeGradShapeFunction(xi.data(), dN_vec);
63 e, mesh_space_dimension);
64
65 Eigen::MatrixXd J = Eigen::MatrixXd::Zero(Dim, Dim);
66 assert(e.getNumberOfNodes() == GradShapeFunction::NPOINTS);
67 for (unsigned k = 0; k < GradShapeFunction::NPOINTS; k++)
68 {
69 const MathLib::Point3d& mapped_pt =
70 ele_local_coord.getMappedCoordinates(k);
71 J += dNdxi.col(k) * mapped_pt.asEigenVector3d().head(Dim).transpose();
72 }
73
74 if (!check_reordered)
75 {
76 return !(J.determinant() < 0);
77 }
78
79 if (J.determinant() < 0)
80 {
82 "Element {:d} has negative Jacobian determinant {:g}. "
83 "NodeReordering fails.",
84 e.getID(), J.determinant());
85 }
86
87 return true;
88}
89
90template <typename T>
91concept ShapeFunction = requires(double* xi, Eigen::Map<Eigen::VectorXd> dN) {
92 { T::DIM } -> std::convertible_to<int>;
93 { T::NPOINTS } -> std::convertible_to<int>;
94 T::computeGradShapeFunction(xi, dN);
95};
96
98{
99 std::function<bool(MeshLib::Element&, int, bool)> is_node_ordering_correct;
100 std::function<void(MeshLib::Element&, const std::vector<MeshLib::Node*>&)>
102};
103
104template <ShapeFunction ShapeFunc, int Dim>
106 std::array<double, Dim> xi,
107 std::function<void(MeshLib::Element&, const std::vector<MeshLib::Node*>&)>
108 reorder)
109{
111 [xi](MeshLib::Element& e, int mesh_space_dimension,
112 bool check_reordered)
113 {
115 e, mesh_space_dimension, xi, check_reordered);
116 },
117 reorder};
118}
119
120static const std::array<ElementReorderConfigBase,
121 static_cast<int>(MeshLib::CellType::enum_length)>
123{
124 auto swap_nodes_i_j = [](MeshLib::Element& element,
125 const std::vector<MeshLib::Node*>& nodes,
126 unsigned const i, unsigned const j)
127 {
128 element.setNode(i, nodes[j]);
129 element.setNode(j, nodes[i]);
130 };
131
132 auto order_nodes_quadratic_quad =
133 [swap_nodes_i_j](MeshLib::Element& element,
134 const std::vector<MeshLib::Node*>& nodes)
135 {
136 swap_nodes_i_j(element, nodes, 1, 3);
137 swap_nodes_i_j(element, nodes, 5, 6);
138 swap_nodes_i_j(element, nodes, 4, 7);
139 };
140
141 std::array<ElementReorderConfigBase,
142 static_cast<int>(MeshLib::CellType::enum_length)>
143 arr{};
144
145 arr[static_cast<int>(MeshLib::CellType::LINE2)] =
147 {0.5},
148 [&swap_nodes_i_j](MeshLib::Element& element,
149 const std::vector<MeshLib::Node*>& nodes)
150 { swap_nodes_i_j(element, nodes, 0, 1); });
151
152 arr[static_cast<int>(MeshLib::CellType::LINE3)] =
154 {0.5},
155 [&swap_nodes_i_j](MeshLib::Element& element,
156 const std::vector<MeshLib::Node*>& nodes)
157 { swap_nodes_i_j(element, nodes, 0, 1); });
158
159 arr[static_cast<int>(MeshLib::CellType::TRI3)] =
161 {1.0 / 3.0, 1.0 / 3.0},
162 [&swap_nodes_i_j](MeshLib::Element& element,
163 const std::vector<MeshLib::Node*>& nodes)
164 { swap_nodes_i_j(element, nodes, 1, 2); });
165
166 arr[static_cast<int>(MeshLib::CellType::TRI6)] =
168 {1.0 / 3.0, 1.0 / 3.0},
169 [&swap_nodes_i_j](MeshLib::Element& element,
170 const std::vector<MeshLib::Node*>& nodes)
171 {
172 swap_nodes_i_j(element, nodes, 1, 2);
173 swap_nodes_i_j(element, nodes, 3, 5);
174 });
175
176 arr[static_cast<int>(MeshLib::CellType::QUAD4)] =
178 {0.0, 0.0},
179 [&swap_nodes_i_j](MeshLib::Element& element,
180 const std::vector<MeshLib::Node*>& nodes)
181 { swap_nodes_i_j(element, nodes, 0, 2); });
182
183 arr[static_cast<int>(MeshLib::CellType::QUAD8)] =
185 {0.0, 0.0},
186 [&order_nodes_quadratic_quad](
187 MeshLib::Element& element,
188 const std::vector<MeshLib::Node*>& nodes)
189 { order_nodes_quadratic_quad(element, nodes); });
190 arr[static_cast<int>(MeshLib::CellType::QUAD9)] =
192 {0.0, 0.0},
193 [&order_nodes_quadratic_quad](
194 MeshLib::Element& element,
195 const std::vector<MeshLib::Node*>& nodes)
196 { order_nodes_quadratic_quad(element, nodes); });
197 arr[static_cast<int>(MeshLib::CellType::TET4)] =
199 {0.25, 0.25, 0.25},
200 [&swap_nodes_i_j](MeshLib::Element& element,
201 const std::vector<MeshLib::Node*>& nodes)
202 { swap_nodes_i_j(element, nodes, 1, 2); });
203
204 arr[static_cast<int>(MeshLib::CellType::TET10)] =
206 {0.25, 0.25, 0.25},
207 [&swap_nodes_i_j](MeshLib::Element& element,
208 const std::vector<MeshLib::Node*>& nodes)
209 {
210 swap_nodes_i_j(element, nodes, 1, 2);
211 swap_nodes_i_j(element, nodes, 4, 6);
212 swap_nodes_i_j(element, nodes, 8, 9);
213 });
214 arr[static_cast<int>(MeshLib::CellType::PRISM6)] =
216 {1.0 / 3.0, 1.0 / 3.0, 0.5},
217 [&swap_nodes_i_j](MeshLib::Element& element,
218 const std::vector<MeshLib::Node*>& nodes)
219 {
220 swap_nodes_i_j(element, nodes, 1, 2);
221 swap_nodes_i_j(element, nodes, 4, 5);
222 });
223
224 arr[static_cast<int>(MeshLib::CellType::PRISM15)] =
226 {1.0 / 3.0, 1.0 / 3., 0.5},
227 [&swap_nodes_i_j](MeshLib::Element& element,
228 const std::vector<MeshLib::Node*>& nodes)
229 {
230 swap_nodes_i_j(element, nodes, 0, 3);
231 swap_nodes_i_j(element, nodes, 1, 4);
232 swap_nodes_i_j(element, nodes, 2, 5);
233 swap_nodes_i_j(element, nodes, 6, 9);
234 swap_nodes_i_j(element, nodes, 7, 10);
235 swap_nodes_i_j(element, nodes, 8, 11);
236 });
237
238 arr[static_cast<int>(MeshLib::CellType::PYRAMID5)] =
240 {0.25, 0.25, 0.5},
241 [&swap_nodes_i_j](MeshLib::Element& element,
242 const std::vector<MeshLib::Node*>& nodes)
243 { swap_nodes_i_j(element, nodes, 0, 2); });
244
245 arr[static_cast<int>(MeshLib::CellType::PYRAMID13)] =
247 {0.25, 0.25, 0.5},
248 [&swap_nodes_i_j](MeshLib::Element& element,
249 const std::vector<MeshLib::Node*>& nodes)
250 {
251 swap_nodes_i_j(element, nodes, 0, 2);
252 swap_nodes_i_j(element, nodes, 9, 11);
253 swap_nodes_i_j(element, nodes, 5, 6);
254 swap_nodes_i_j(element, nodes, 7, 8);
255 });
256 arr[static_cast<int>(MeshLib::CellType::HEX8)] =
258 {0.5, 0.5, 0.5},
259 [&swap_nodes_i_j](MeshLib::Element& element,
260 const std::vector<MeshLib::Node*>& nodes)
261 {
262 swap_nodes_i_j(element, nodes, 0, 2);
263 swap_nodes_i_j(element, nodes, 4, 6);
264 });
265
266 arr[static_cast<int>(MeshLib::CellType::HEX20)] =
268 {0.5, 0.5, 0.5},
269 [&swap_nodes_i_j](MeshLib::Element& element,
270 const std::vector<MeshLib::Node*>& nodes)
271 {
272 swap_nodes_i_j(element, nodes, 0, 2);
273 swap_nodes_i_j(element, nodes, 4, 6);
274 swap_nodes_i_j(element, nodes, 16, 18);
275 swap_nodes_i_j(element, nodes, 8, 9);
276 swap_nodes_i_j(element, nodes, 10, 11);
277 swap_nodes_i_j(element, nodes, 12, 13);
278 swap_nodes_i_j(element, nodes, 14, 15);
279 });
280 return arr;
281}();
282
292void reverseNodeOrder(std::vector<MeshLib::Element*>& elements,
293 int const mesh_space_dimension, bool const forced)
294{
295 std::size_t n_corrected_elements = 0;
296
297 for (auto* element : elements)
298 {
299 auto const cell_type = element->getCellType();
300
301 if (cell_type == MeshLib::CellType::INVALID)
302 {
303 OGS_FATAL("Element type `{:s}' does not exist.",
304 MeshLib::CellType2String(cell_type));
305 }
306 // This check is already done in MeshLib::readMeshFromFile(). Therefore,
307 // this is just a safeguard.
308 if (cell_type == MeshLib::CellType::HEX27 ||
309 cell_type == MeshLib::CellType::PRISM18)
310 {
311 OGS_FATAL("Element type `{:s}' is not supported in OGS.",
312 MeshLib::CellType2String(cell_type));
313 }
314
315 if (cell_type == MeshLib::CellType::POINT1)
316 {
317 continue;
318 }
319
320 const auto& element_config =
321 element_configs_array[static_cast<int>(cell_type)];
322
323 if (element_config.is_node_ordering_correct(
324 *element, mesh_space_dimension, false /*check_reordered*/))
325 {
326 continue;
327 }
328
329 // Save nodes before reordering
330 const unsigned nElemNodes = element->getNumberOfNodes();
331 std::vector<MeshLib::Node*> nodes(element->getNodes(),
332 element->getNodes() + nElemNodes);
333
334 double const element_volume_origin = element->computeVolume();
335
336 element_config.reorder_element_nodes(*element, nodes);
337
338 // Ensure that the element volume did not change.
339 double const element_volume = element->computeVolume();
340 // We use a fixed tolerance here because for very small elements the
341 // machine epsilon might be too small.
342 if (std::fabs(element_volume - element_volume_origin) > 1e-14)
343 {
344 OGS_FATAL(
345 "Reordering of element {:d} nodes failed as its volume changed "
346 "from {:.12e} to {:.12e}.",
347 element->getID(), element_volume_origin, element_volume);
348 }
349 // Element::computeVolume() uses the element vertecies to compute
350 // the element volume, so the change of edge nodes are not
351 // considered. Therefore, we need to additionally check the Jacobian
352 // determinant here if forced is true.
353 if (forced)
354 {
355 element_config.is_node_ordering_correct(
356 *element, mesh_space_dimension, true /*check_reordered*/);
357 }
358
359 ++n_corrected_elements;
360 }
361
362 INFO("Corrected {:d} elements.", n_corrected_elements);
363}
364
369void fixVtkInconsistencies(std::vector<MeshLib::Element*>& elements)
370{
371 for (auto* const element : elements)
372 {
373 const unsigned nElemNodes(element->getNumberOfBaseNodes());
374 std::vector<MeshLib::Node*> nodes(element->getNodes(),
375 element->getNodes() + nElemNodes);
376
377 for (std::size_t j = 0; j < nElemNodes; ++j)
378 {
379 if (element->getGeomType() == MeshLib::MeshElemType::PRISM)
380 {
381 for (std::size_t k = 0; k < 3; ++k)
382 {
383 element->setNode(k, nodes[k + 3]);
384 element->setNode(k + 3, nodes[k]);
385 }
386 break;
387 }
388 }
389 }
390}
391
394{
395 std::vector<MeshLib::Node*> base_nodes;
396 std::vector<MeshLib::Node*> nonlinear_nodes;
397 for (MeshLib::Element const* e : mesh.getElements())
398 {
399 for (unsigned i = 0; i < e->getNumberOfBaseNodes(); i++)
400 {
401 base_nodes.push_back(const_cast<MeshLib::Node*>(e->getNode(i)));
402 }
403 for (unsigned i = e->getNumberOfBaseNodes(); i < e->getNumberOfNodes();
404 i++)
405 {
406 nonlinear_nodes.push_back(
407 const_cast<MeshLib::Node*>(e->getNode(i)));
408 }
409 }
410
411 BaseLib::makeVectorUnique(base_nodes,
413 BaseLib::makeVectorUnique(nonlinear_nodes,
415
416 std::vector<MeshLib::Node*>& allnodes =
417 const_cast<std::vector<MeshLib::Node*>&>(mesh.getNodes());
418 allnodes.clear();
419
420 allnodes.insert(allnodes.end(), base_nodes.begin(), base_nodes.end());
421 allnodes.insert(allnodes.end(), nonlinear_nodes.begin(),
422 nonlinear_nodes.end());
423
424 mesh.resetNodeIDs();
425}
426
427int main(int argc, char* argv[])
428{
429 enum class ExpectedCellType
430 {
431 INVALID = 0,
432 POINT1 = 1,
433 LINE2 = 2,
434 LINE3 = 3,
435 TRI3 = 4,
436 TRI6 = 5,
437 QUAD4 = 6,
438 QUAD8 = 7,
439 QUAD9 = 8,
440 TET4 = 9,
441 TET10 = 10,
442 HEX8 = 11,
443 HEX20 = 12,
444 HEX27 = 13,
445 PRISM6 = 14,
446 PRISM15 = 15,
447 PRISM18 = 16,
448 PYRAMID5 = 17,
449 PYRAMID13 = 18,
451 };
452
453 constexpr bool are_expected_cell_types =
454 static_cast<int>(ExpectedCellType::enum_length) ==
455 static_cast<int>(MeshLib::CellType::enum_length) &&
456 static_cast<int>(ExpectedCellType::INVALID) ==
457 static_cast<int>(MeshLib::CellType::INVALID) &&
458 static_cast<int>(ExpectedCellType::POINT1) ==
459 static_cast<int>(MeshLib::CellType::POINT1) &&
460 static_cast<int>(ExpectedCellType::LINE2) ==
461 static_cast<int>(MeshLib::CellType::LINE2) &&
462 static_cast<int>(ExpectedCellType::LINE3) ==
463 static_cast<int>(MeshLib::CellType::LINE3) &&
464 static_cast<int>(ExpectedCellType::TRI3) ==
465 static_cast<int>(MeshLib::CellType::TRI3) &&
466 static_cast<int>(ExpectedCellType::TRI6) ==
467 static_cast<int>(MeshLib::CellType::TRI6) &&
468 static_cast<int>(ExpectedCellType::QUAD4) ==
469 static_cast<int>(MeshLib::CellType::QUAD4) &&
470 static_cast<int>(ExpectedCellType::QUAD8) ==
471 static_cast<int>(MeshLib::CellType::QUAD8) &&
472 static_cast<int>(ExpectedCellType::QUAD9) ==
473 static_cast<int>(MeshLib::CellType::QUAD9) &&
474 static_cast<int>(ExpectedCellType::TET4) ==
475 static_cast<int>(MeshLib::CellType::TET4) &&
476 static_cast<int>(ExpectedCellType::TET10) ==
477 static_cast<int>(MeshLib::CellType::TET10) &&
478 static_cast<int>(ExpectedCellType::HEX8) ==
479 static_cast<int>(MeshLib::CellType::HEX8) &&
480 static_cast<int>(ExpectedCellType::HEX20) ==
481 static_cast<int>(MeshLib::CellType::HEX20) &&
482 static_cast<int>(ExpectedCellType::HEX27) ==
483 static_cast<int>(MeshLib::CellType::HEX27) &&
484 static_cast<int>(ExpectedCellType::PRISM6) ==
485 static_cast<int>(MeshLib::CellType::PRISM6) &&
486 static_cast<int>(ExpectedCellType::PRISM15) ==
487 static_cast<int>(MeshLib::CellType::PRISM15) &&
488 static_cast<int>(ExpectedCellType::PRISM18) ==
489 static_cast<int>(MeshLib::CellType::PRISM18) &&
490 static_cast<int>(ExpectedCellType::PYRAMID5) ==
491 static_cast<int>(MeshLib::CellType::PYRAMID5) &&
492 static_cast<int>(ExpectedCellType::PYRAMID13) ==
493 static_cast<int>(MeshLib::CellType::PYRAMID13);
494 // This error should never occur unless someone has changed the
495 // MeshLib::CellType enum class. These assertions ensure that the
496 // array 'element_configs_array' is up to date.
497 if (!are_expected_cell_types)
498 {
499 OGS_FATAL(
500 "The enum class MeshLib::CellType has changed. Please adapt the "
501 "array 'element_configs_array' in NodeReordering.cpp accordingly.");
502 }
503
504 TCLAP::CmdLine cmd(
505 "Reorders mesh nodes in elements to make old or incorrectly ordered "
506 "meshes compatible with OGS6.\n"
507 "Three options are available:\n"
508 "Method 0: Reversing order of nodes and checking again for all "
509 "elements.\n"
510 "Method 1: Reversing order of nodes unless it's perceived correct by "
511 "OGS6 standards. This is the default selection.\n"
512 "Method 2: Fixing node ordering issues between VTK and OGS6 (only "
513 "applies to prism-elements)\n"
514 "Method 3: Re-ordering of mesh node vector such that all base nodes "
515 "are sorted before all nonlinear nodes.\n\n"
516 "OpenGeoSys-6 software, version " +
518 ".\n"
519 "Copyright (c) 2012-2025, OpenGeoSys Community "
520 "(http://www.opengeosys.org)",
522
523 std::vector<int> method_ids{0, 1, 2, 3};
524 TCLAP::ValuesConstraint<int> allowed_values(method_ids);
525 TCLAP::ValueArg<int> method_arg("m", "method",
526 "reordering method selection", false, 1,
527 &allowed_values);
528 cmd.add(method_arg);
529 TCLAP::ValueArg<std::string> output_mesh_arg(
530 "o", "output_mesh", "Output (.vtu). The name of the output mesh file",
531 true, "", "OUTPUT_FILE");
532 cmd.add(output_mesh_arg);
533 TCLAP::ValueArg<std::string> input_mesh_arg(
534 "i", "input_mesh",
535 "Input (.vtu | .vtk | .msh). The name of the input mesh file", true, "",
536 "INPUT_FILE");
537 cmd.add(input_mesh_arg);
538 auto log_level_arg = BaseLib::makeLogLevelArg();
539 cmd.add(log_level_arg);
540 cmd.parse(argc, argv);
541
542 BaseLib::MPI::Setup mpi_setup(argc, argv);
543 BaseLib::initOGSLogger(log_level_arg.getValue());
544
545 std::unique_ptr<MeshLib::Mesh> mesh(
546 MeshLib::IO::readMeshFromFile(input_mesh_arg.getValue()));
547
548 if (!mesh)
549 {
550 return EXIT_FAILURE;
551 }
552
553 INFO("Reordering nodes... ");
554 if (!method_arg.isSet() || method_arg.getValue() < 2)
555 {
556 bool const forced = (method_arg.getValue() == 0);
557
558 if (forced)
559 {
560 INFO("Method 0: Reversing order of nodes will be checked again.");
561 }
562 INFO(
563 "Method: Reversing order of nodes unless it is considered correct "
564 "by the OGS6 standard, i.e. such that det(J) > 0, where J is the "
565 "Jacobian of the global-to-local coordinate transformation.");
566 int const mesh_space_dimension =
567 MeshLib::getSpaceDimension(mesh->getNodes());
569 const_cast<std::vector<MeshLib::Element*>&>(mesh->getElements()),
570 mesh_space_dimension, forced);
571 }
572 else if (method_arg.getValue() == 2)
573 {
575 const_cast<std::vector<MeshLib::Element*>&>(mesh->getElements()));
576 }
577 else if (method_arg.getValue() == 3)
578 {
580 }
581
582 MeshLib::IO::writeMeshToFile(*mesh, output_mesh_arg.getValue());
583
584 INFO("VTU file written.");
585
586 return EXIT_SUCCESS;
587}
Definition of the Element class.
#define OGS_FATAL(...)
Definition Error.h:26
Git information.
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
Definition of mesh-related Enumerations.
Definition of the Mesh class.
int main(int argc, char *argv[])
bool checkJacobianDeterminant(MeshLib::Element const &e, int const mesh_space_dimension, std::array< double, Dim > const &xi, bool const check_reordered=false)
static const std::array< ElementReorderConfigBase, static_cast< int >(MeshLib::CellType::enum_length)> element_configs_array
void reverseNodeOrder(std::vector< MeshLib::Element * > &elements, int const mesh_space_dimension, bool const forced)
Reverses order of nodes. In particular, this fixes issues between (Gmsh or OGS5) and OGS6 meshes.
void fixVtkInconsistencies(std::vector< MeshLib::Element * > &elements)
void reorderNonlinearNodes(MeshLib::Mesh &mesh)
Orders the base nodes of each elements before its non-linear nodes.
ElementReorderConfigBase makeElementConfig(std::array< double, Dim > xi, std::function< void(MeshLib::Element &, const std::vector< MeshLib::Node * > &)> reorder)
Definition of the Node class.
Eigen::Vector3d const & asEigenVector3d() const
Definition Point3d.h:64
MathLib::Point3d const & getMappedCoordinates(std::size_t node_id) const
return mapped coordinates of the node
virtual unsigned getNumberOfNodes() const =0
virtual void setNode(unsigned idx, Node *node)=0
std::size_t getID() const
Returns the ID of the element.
Definition Element.h:91
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition Mesh.h:108
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition Mesh.h:111
void resetNodeIDs()
Resets the IDs of all mesh-nodes to their position in the node vector.
Definition Mesh.cpp:160
TCLAP::ValueArg< std::string > makeLogLevelArg()
void initOGSLogger(std::string const &log_level)
Definition Logging.cpp:64
void makeVectorUnique(std::vector< T > &v)
Definition Algorithm.h:180
GITINFOLIB_EXPORT const std::string ogs_version
MeshLib::Mesh * readMeshFromFile(const std::string &file_name, bool const compute_element_neighbors)
int writeMeshToFile(const MeshLib::Mesh &mesh, std::filesystem::path const &file_path, std::set< std::string > variable_output_names)
bool idsComparator(T const a, T const b)
Definition Mesh.h:208
std::string CellType2String(const CellType t)
Given a MeshElemType this returns the appropriate string.
int getSpaceDimension(std::vector< Node * > const &nodes)
Computes dimension of the embedding space containing the set of given points.
Definition of readMeshFromFile function.
std::function< bool(MeshLib::Element &, int, bool)> is_node_ordering_correct
std::function< void(MeshLib::Element &, const std::vector< MeshLib::Node * > &)> reorder_element_nodes