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