From 1225e7be17297cd034df4c89070494caf9886ada Mon Sep 17 00:00:00 2001
From: marioalexis <mario.passaglia@gmail.com>
Date: Sat, 9 Mar 2024 15:05:53 -0300
Subject: [PATCH] Fem: Export/import mesh elements to Vtk using proper order

---
 src/Mod/Fem/App/FemVTKTools.cpp               | 392 +++++++-----------
 .../Fem/femtest/data/mesh/tetra10_mesh.vtk    |   2 +-
 2 files changed, 157 insertions(+), 237 deletions(-)

diff --git a/src/Mod/Fem/App/FemVTKTools.cpp b/src/Mod/Fem/App/FemVTKTools.cpp
index ac88e59518..fd5df42bc9 100644
--- a/src/Mod/Fem/App/FemVTKTools.cpp
+++ b/src/Mod/Fem/App/FemVTKTools.cpp
@@ -97,6 +97,52 @@ void writeVTKFile(const char* filename, vtkSmartPointer<vtkUnstructuredGrid> dat
     writer->Write();
 }
 
+namespace
+{
+
+// Helper function to fill vtkCellArray from SMDS_Mesh using vtk cell order
+template<typename T, typename E>
+void fillVtkArray(vtkSmartPointer<vtkCellArray>& elemArray, std::vector<int>& types, const E* elem)
+{
+    vtkSmartPointer<T> cell = vtkSmartPointer<T>::New();
+    const std::vector<int>& order = SMDS_MeshCell::toVtkOrder(elem->GetEntityType());
+    if (!order.empty()) {
+        for (int i = 0; i < elem->NbNodes(); ++i) {
+            cell->GetPointIds()->SetId(i, elem->GetNode(order[i])->GetID() - 1);
+        }
+    }
+    else {
+        for (int i = 0; i < elem->NbNodes(); ++i) {
+            cell->GetPointIds()->SetId(i, elem->GetNode(i)->GetID() - 1);
+        }
+    }
+    elemArray->InsertNextCell(cell);
+    types.push_back(SMDS_MeshCell::toVtkType(elem->GetEntityType()));
+}
+
+// Helper function to fill SMDS_Mesh elements ID from vtk cell
+void fillMeshElementIds(vtkCell* cell, std::vector<int>& ids)
+{
+    VTKCellType cellType = static_cast<VTKCellType>(cell->GetCellType());
+    const std::vector<int>& order = SMDS_MeshCell::fromVtkOrder(cellType);
+    vtkIdType* vtkIds = cell->GetPointIds()->GetPointer(0);
+    ids.clear();
+    int nbPoints = cell->GetNumberOfPoints();
+    ids.resize(nbPoints);
+    if (!order.empty()) {
+        for (int i = 0; i < nbPoints; ++i) {
+            ids[i] = vtkIds[order[i]] + 1;
+        }
+    }
+    else {
+        for (int i = 0; i < nbPoints; ++i) {
+            ids[i] = vtkIds[i] + 1;
+        }
+    }
+}
+
+}  // namespace
+
 
 void FemVTKTools::importVTKMesh(vtkSmartPointer<vtkDataSet> dataset, FemMesh* mesh, float scale)
 {
@@ -105,10 +151,6 @@ void FemVTKTools::importVTKMesh(vtkSmartPointer<vtkDataSet> dataset, FemMesh* me
     Base::Console().Log("%d nodes/points and %d cells/elements found!\n", nPoints, nCells);
     Base::Console().Log("Build SMESH mesh out of the vtk mesh data.\n", nPoints, nCells);
 
-    // vtkSmartPointer<vtkCellArray> cells = dataset->GetCells();
-    //  works only for vtkUnstructuredGrid
-    vtkSmartPointer<vtkIdList> idlist = vtkSmartPointer<vtkIdList>::New();
-
     // Now fill the SMESH datastructure
     SMESH_Mesh* smesh = mesh->getSMesh();
     SMESHDS_Mesh* meshds = smesh->GetMeshDS();
@@ -120,138 +162,120 @@ void FemVTKTools::importVTKMesh(vtkSmartPointer<vtkDataSet> dataset, FemMesh* me
     }
 
     for (vtkIdType iCell = 0; iCell < nCells; iCell++) {
-        idlist->Reset();
-        idlist = dataset->GetCell(iCell)->GetPointIds();
-        vtkIdType* ids = idlist->GetPointer(0);
-        switch (dataset->GetCellType(iCell)) {
+        vtkCell* cell = dataset->GetCell(iCell);
+        std::vector<int> ids;
+        fillMeshElementIds(cell, ids);
+        switch (cell->GetCellType()) {
             // 2D faces
             case VTK_TRIANGLE:  // tria3
-                meshds->AddFaceWithID(ids[0] + 1, ids[1] + 1, ids[2] + 1, iCell + 1);
+                meshds->AddFaceWithID(ids[0], ids[1], ids[2], iCell + 1);
                 break;
             case VTK_QUADRATIC_TRIANGLE:  // tria6
-                meshds->AddFaceWithID(ids[0] + 1,
-                                      ids[1] + 1,
-                                      ids[2] + 1,
-                                      ids[3] + 1,
-                                      ids[4] + 1,
-                                      ids[5] + 1,
-                                      iCell + 1);
+                meshds->AddFaceWithID(ids[0], ids[1], ids[2], ids[3], ids[4], ids[5], iCell + 1);
                 break;
             case VTK_QUAD:  // quad4
-                meshds->AddFaceWithID(ids[0] + 1, ids[1] + 1, ids[2] + 1, ids[3] + 1, iCell + 1);
+                meshds->AddFaceWithID(ids[0], ids[1], ids[2], ids[3], iCell + 1);
                 break;
             case VTK_QUADRATIC_QUAD:  // quad8
-                meshds->AddFaceWithID(ids[0] + 1,
-                                      ids[1] + 1,
-                                      ids[2] + 1,
-                                      ids[3] + 1,
-                                      ids[4] + 1,
-                                      ids[5] + 1,
-                                      ids[6] + 1,
-                                      ids[7] + 1,
+                meshds->AddFaceWithID(ids[0],
+                                      ids[1],
+                                      ids[2],
+                                      ids[3],
+                                      ids[4],
+                                      ids[5],
+                                      ids[6],
+                                      ids[7],
                                       iCell + 1);
                 break;
-
             // 3D volumes
             case VTK_TETRA:  // tetra4
-                meshds->AddVolumeWithID(ids[0] + 1, ids[1] + 1, ids[2] + 1, ids[3] + 1, iCell + 1);
+                meshds->AddVolumeWithID(ids[0], ids[1], ids[2], ids[3], iCell + 1);
                 break;
             case VTK_QUADRATIC_TETRA:  // tetra10
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        ids[6] + 1,
-                                        ids[7] + 1,
-                                        ids[8] + 1,
-                                        ids[9] + 1,
+                meshds->AddVolumeWithID(ids[0],
+                                        ids[1],
+                                        ids[2],
+                                        ids[3],
+                                        ids[4],
+                                        ids[5],
+                                        ids[6],
+                                        ids[7],
+                                        ids[8],
+                                        ids[9],
                                         iCell + 1);
                 break;
             case VTK_HEXAHEDRON:  // hexa8
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        ids[6] + 1,
-                                        ids[7] + 1,
+                meshds->AddVolumeWithID(ids[0],
+                                        ids[1],
+                                        ids[2],
+                                        ids[3],
+                                        ids[4],
+                                        ids[5],
+                                        ids[6],
+                                        ids[7],
                                         iCell + 1);
                 break;
             case VTK_QUADRATIC_HEXAHEDRON:  // hexa20
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        ids[6] + 1,
-                                        ids[7] + 1,
-                                        ids[8] + 1,
-                                        ids[9] + 1,
-                                        ids[10] + 1,
-                                        ids[11] + 1,
-                                        ids[12] + 1,
-                                        ids[13] + 1,
-                                        ids[14] + 1,
-                                        ids[15] + 1,
-                                        ids[16] + 1,
-                                        ids[17] + 1,
-                                        ids[18] + 1,
-                                        ids[19] + 1,
+                meshds->AddVolumeWithID(ids[0],
+                                        ids[1],
+                                        ids[2],
+                                        ids[3],
+                                        ids[4],
+                                        ids[5],
+                                        ids[6],
+                                        ids[7],
+                                        ids[8],
+                                        ids[9],
+                                        ids[10],
+                                        ids[11],
+                                        ids[12],
+                                        ids[13],
+                                        ids[14],
+                                        ids[15],
+                                        ids[16],
+                                        ids[17],
+                                        ids[18],
+                                        ids[19],
                                         iCell + 1);
                 break;
             case VTK_WEDGE:  // penta6
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        iCell + 1);
+                meshds->AddVolumeWithID(ids[0], ids[1], ids[2], ids[3], ids[4], ids[5], iCell + 1);
                 break;
             case VTK_QUADRATIC_WEDGE:  // penta15
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        ids[6] + 1,
-                                        ids[7] + 1,
-                                        ids[8] + 1,
-                                        ids[9] + 1,
-                                        ids[10] + 1,
-                                        ids[11] + 1,
-                                        ids[12] + 1,
-                                        ids[13] + 1,
-                                        ids[14] + 1,
+                meshds->AddVolumeWithID(ids[0],
+                                        ids[1],
+                                        ids[2],
+                                        ids[3],
+                                        ids[4],
+                                        ids[5],
+                                        ids[6],
+                                        ids[7],
+                                        ids[8],
+                                        ids[9],
+                                        ids[10],
+                                        ids[11],
+                                        ids[12],
+                                        ids[13],
+                                        ids[14],
                                         iCell + 1);
                 break;
             case VTK_PYRAMID:  // pyra5
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        iCell + 1);
+                meshds->AddVolumeWithID(ids[0], ids[1], ids[2], ids[3], ids[4], iCell + 1);
                 break;
             case VTK_QUADRATIC_PYRAMID:  // pyra13
-                meshds->AddVolumeWithID(ids[0] + 1,
-                                        ids[1] + 1,
-                                        ids[2] + 1,
-                                        ids[3] + 1,
-                                        ids[4] + 1,
-                                        ids[5] + 1,
-                                        ids[6] + 1,
-                                        ids[7] + 1,
-                                        ids[8] + 1,
-                                        ids[9] + 1,
-                                        ids[10] + 1,
-                                        ids[11] + 1,
-                                        ids[12] + 1,
+                meshds->AddVolumeWithID(ids[0],
+                                        ids[1],
+                                        ids[2],
+                                        ids[3],
+                                        ids[4],
+                                        ids[5],
+                                        ids[6],
+                                        ids[7],
+                                        ids[8],
+                                        ids[9],
+                                        ids[10],
+                                        ids[11],
+                                        ids[12],
                                         iCell + 1);
                 break;
 
@@ -314,58 +338,23 @@ void exportFemMeshFaces(vtkSmartPointer<vtkUnstructuredGrid> grid,
     vtkSmartPointer<vtkCellArray> elemArray = vtkSmartPointer<vtkCellArray>::New();
     std::vector<int> types;
 
-    for (; aFaceIter->more();) {
+    while (aFaceIter->more()) {
         const SMDS_MeshFace* aFace = aFaceIter->next();
-
         // triangle
-        if (aFace->NbNodes() == 3) {
-            vtkSmartPointer<vtkTriangle> tria = vtkSmartPointer<vtkTriangle>::New();
-            tria->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID() - 1);
-            tria->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID() - 1);
-            tria->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID() - 1);
-
-            elemArray->InsertNextCell(tria);
-            types.push_back(VTK_TRIANGLE);
+        if (aFace->GetEntityType() == SMDSEntity_Triangle) {
+            fillVtkArray<vtkTriangle>(elemArray, types, aFace);
         }
         // quad
-        else if (aFace->NbNodes() == 4) {
-            vtkSmartPointer<vtkQuad> quad = vtkSmartPointer<vtkQuad>::New();
-            quad->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID() - 1);
-            quad->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID() - 1);
-            quad->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID() - 1);
-            quad->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID() - 1);
-
-            elemArray->InsertNextCell(quad);
-            types.push_back(VTK_QUAD);
+        else if (aFace->GetEntityType() == SMDSEntity_Quadrangle) {
+            fillVtkArray<vtkQuad>(elemArray, types, aFace);
         }
         // quadratic triangle
-        else if (aFace->NbNodes() == 6) {
-            vtkSmartPointer<vtkQuadraticTriangle> tria =
-                vtkSmartPointer<vtkQuadraticTriangle>::New();
-            tria->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID() - 1);
-            tria->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID() - 1);
-            tria->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID() - 1);
-            tria->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID() - 1);
-            tria->GetPointIds()->SetId(4, aFace->GetNode(4)->GetID() - 1);
-            tria->GetPointIds()->SetId(5, aFace->GetNode(5)->GetID() - 1);
-
-            elemArray->InsertNextCell(tria);
-            types.push_back(VTK_QUADRATIC_TRIANGLE);
+        else if (aFace->GetEntityType() == SMDSEntity_Quad_Triangle) {
+            fillVtkArray<vtkQuadraticTriangle>(elemArray, types, aFace);
         }
         // quadratic quad
-        else if (aFace->NbNodes() == 8) {
-            vtkSmartPointer<vtkQuadraticQuad> quad = vtkSmartPointer<vtkQuadraticQuad>::New();
-            quad->GetPointIds()->SetId(0, aFace->GetNode(0)->GetID() - 1);
-            quad->GetPointIds()->SetId(1, aFace->GetNode(1)->GetID() - 1);
-            quad->GetPointIds()->SetId(2, aFace->GetNode(2)->GetID() - 1);
-            quad->GetPointIds()->SetId(3, aFace->GetNode(3)->GetID() - 1);
-            quad->GetPointIds()->SetId(4, aFace->GetNode(4)->GetID() - 1);
-            quad->GetPointIds()->SetId(5, aFace->GetNode(5)->GetID() - 1);
-            quad->GetPointIds()->SetId(6, aFace->GetNode(6)->GetID() - 1);
-            quad->GetPointIds()->SetId(7, aFace->GetNode(7)->GetID() - 1);
-
-            elemArray->InsertNextCell(quad);
-            types.push_back(VTK_QUADRATIC_QUAD);
+        else if (aFace->GetEntityType() == SMDSEntity_Quad_Quadrangle) {
+            fillVtkArray<vtkQuadraticQuad>(elemArray, types, aFace);
         }
         else {
             throw Base::TypeError("Face not yet supported by FreeCAD's VTK mesh builder\n");
@@ -387,101 +376,32 @@ void exportFemMeshCells(vtkSmartPointer<vtkUnstructuredGrid> grid,
     vtkSmartPointer<vtkCellArray> elemArray = vtkSmartPointer<vtkCellArray>::New();
     std::vector<int> types;
 
-    for (; aVolIter->more();) {
+    while (aVolIter->more()) {
         const SMDS_MeshVolume* aVol = aVolIter->next();
 
-        if (aVol->NbNodes() == 4) {  // tetra4
-            Base::Console().Log("    Volume tetra4\n");
-            vtkSmartPointer<vtkTetra> cell = vtkSmartPointer<vtkTetra>::New();
-            cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID() - 1);
-            cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID() - 1);
-            cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID() - 1);
-            cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID() - 1);
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_TETRA);
+        if (aVol->GetEntityType() == SMDSEntity_Tetra) {  // tetra4
+            fillVtkArray<vtkTetra>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 5) {  // pyra5
-            Base::Console().Log("    Volume pyra5\n");
-            vtkSmartPointer<vtkPyramid> cell = vtkSmartPointer<vtkPyramid>::New();
-            cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID() - 1);
-            cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID() - 1);
-            cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID() - 1);
-            cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID() - 1);
-            cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID() - 1);
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_PYRAMID);
+        else if (aVol->GetEntityType() == SMDSEntity_Pyramid) {  // pyra5
+            fillVtkArray<vtkPyramid>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 6) {  // penta6
-            Base::Console().Log("    Volume penta6\n");
-            vtkSmartPointer<vtkWedge> cell = vtkSmartPointer<vtkWedge>::New();
-            cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID() - 1);
-            cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID() - 1);
-            cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID() - 1);
-            cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID() - 1);
-            cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID() - 1);
-            cell->GetPointIds()->SetId(5, aVol->GetNode(5)->GetID() - 1);
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_WEDGE);
+        else if (aVol->GetEntityType() == SMDSEntity_Penta) {  // penta6
+            fillVtkArray<vtkWedge>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 8) {  // hexa8
-            Base::Console().Log("    Volume hexa8\n");
-            vtkSmartPointer<vtkHexahedron> cell = vtkSmartPointer<vtkHexahedron>::New();
-            cell->GetPointIds()->SetId(0, aVol->GetNode(0)->GetID() - 1);
-            cell->GetPointIds()->SetId(1, aVol->GetNode(1)->GetID() - 1);
-            cell->GetPointIds()->SetId(2, aVol->GetNode(2)->GetID() - 1);
-            cell->GetPointIds()->SetId(3, aVol->GetNode(3)->GetID() - 1);
-            cell->GetPointIds()->SetId(4, aVol->GetNode(4)->GetID() - 1);
-            cell->GetPointIds()->SetId(5, aVol->GetNode(5)->GetID() - 1);
-            cell->GetPointIds()->SetId(6, aVol->GetNode(6)->GetID() - 1);
-            cell->GetPointIds()->SetId(7, aVol->GetNode(7)->GetID() - 1);
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_HEXAHEDRON);
+        else if (aVol->GetEntityType() == SMDSEntity_Hexa) {  // hexa8
+            fillVtkArray<vtkHexahedron>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 10) {  // tetra10
-            Base::Console().Log("    Volume tetra10\n");
-            vtkSmartPointer<vtkQuadraticTetra> cell = vtkSmartPointer<vtkQuadraticTetra>::New();
-            for (int i = 0; i < 10; i++) {
-                cell->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID() - 1);
-            }
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_QUADRATIC_TETRA);
+        else if (aVol->GetEntityType() == SMDSEntity_Quad_Tetra) {  // tetra10
+            fillVtkArray<vtkQuadraticTetra>(elemArray, types, aVol);
         }
-
-        else if (aVol->NbNodes() == 13) {  // pyra13
-            Base::Console().Log("    Volume pyra13\n");
-            vtkSmartPointer<vtkQuadraticPyramid> cell = vtkSmartPointer<vtkQuadraticPyramid>::New();
-            for (int i = 0; i < 13; i++) {
-                cell->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID() - 1);
-            }
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_QUADRATIC_PYRAMID);
+        else if (aVol->GetEntityType() == SMDSEntity_Quad_Pyramid) {  // pyra13
+            fillVtkArray<vtkQuadraticPyramid>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 15) {  // penta15
-            Base::Console().Log("    Volume penta15\n");
-            vtkSmartPointer<vtkQuadraticWedge> cell = vtkSmartPointer<vtkQuadraticWedge>::New();
-            for (int i = 0; i < 15; i++) {
-                cell->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID() - 1);
-            }
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_QUADRATIC_WEDGE);
+        else if (aVol->GetEntityType() == SMDSEntity_Quad_Penta) {  // penta15
+            fillVtkArray<vtkQuadraticWedge>(elemArray, types, aVol);
         }
-        else if (aVol->NbNodes() == 20) {  // hexa20
-            Base::Console().Log("    Volume hexa20\n");
-            vtkSmartPointer<vtkQuadraticHexahedron> cell =
-                vtkSmartPointer<vtkQuadraticHexahedron>::New();
-            for (int i = 0; i < 20; i++) {
-                cell->GetPointIds()->SetId(i, aVol->GetNode(i)->GetID() - 1);
-            }
-
-            elemArray->InsertNextCell(cell);
-            types.push_back(VTK_QUADRATIC_HEXAHEDRON);
+        else if (aVol->GetEntityType() == SMDSEntity_Quad_Hexa) {  // hexa20
+            fillVtkArray<vtkQuadraticHexahedron>(elemArray, types, aVol);
         }
         else {
             throw Base::TypeError("Volume not yet supported by FreeCAD's VTK mesh builder\n");
diff --git a/src/Mod/Fem/femtest/data/mesh/tetra10_mesh.vtk b/src/Mod/Fem/femtest/data/mesh/tetra10_mesh.vtk
index 453413b5b1..bc946d3b70 100644
--- a/src/Mod/Fem/femtest/data/mesh/tetra10_mesh.vtk
+++ b/src/Mod/Fem/femtest/data/mesh/tetra10_mesh.vtk
@@ -8,7 +8,7 @@ POINTS 10 float
 9 6 18 6 9 9 3 3 9 
 9 3 9 
 CELLS 1 11
-10 0 1 2 3 4 5 6 7 8 9 
+10 0 2 1 3 6 5 4 7 9 8
 
 CELL_TYPES 1
 24