ENH: Adds surface extraction and oneToOne map

Fixes bugs and tests surface extraction.
Fixes bugs and tests one to one nodal map.

Base/Geometry/Map/imstkGeometryMap.cpp

@@ -67,6 +67,8 @@ std::string GeometryMap::getTypeName() const
             return "Hexa-Triangle map";
         case GeometryMapType::TetraTetra:
             return "Tetra-Tetra map";
+        case GeometryMapType::OneToOne:
+            return "One-to-One nodal map";
         default:
             return "Map type not determined!";
     }

Base/Geometry/Map/imstkOneToOneMap.cpp

@@ -33,6 +33,19 @@ OneToOneMap::compute()
         return;
     }
 
+    // returns the first matching vertex
+    auto findMatchingVertex = [](std::shared_ptr<Mesh> masterMesh, const Vec3d& p) -> int
+    {
+        for (size_t nodeId = 0; nodeId < masterMesh->getNumVertices(); ++nodeId)
+        {
+            if (masterMesh->getInitialVertexPosition(nodeId) == p)
+            {
+                return nodeId;
+            }
+        }
+        return -1;
+    };
+
     auto meshMaster = std::dynamic_pointer_cast<Mesh>(m_master);
     auto meshSlave = std::dynamic_pointer_cast<Mesh>(m_slave);
 
@@ -40,7 +53,7 @@ OneToOneMap::compute()
     for (size_t nodeId = 0; nodeId < meshSlave->getNumVertices(); ++nodeId)
     {
         // Find the enclosing or closest tetrahedron
-        int matchingNodeId = findMatchingVertex(meshMaster, meshSlave->getInitialVertexPosition(nodeId));
+        int matchingNodeId = findMatchingVertex(meshMaster, meshSlave->getVertexPosition(nodeId));
 
         if (matchingNodeId < 0)
         {
@@ -54,20 +67,6 @@ OneToOneMap::compute()
     }
 }
 
-int
-OneToOneMap::findMatchingVertex(std::shared_ptr<Mesh> masterMesh, const Vec3d& p) const
-{
-    // returns the first matching vertex
-    for (size_t nodeId = 0; nodeId < masterMesh->getNumVertices(); ++nodeId)
-    {
-        if (masterMesh->getInitialVertexPosition(nodeId) == p)
-        {
-            return nodeId;
-        }
-    }
-    return -1;
-}
-
 bool
 OneToOneMap::isValid() const
 {
@@ -132,7 +131,7 @@ OneToOneMap::print() const
     GeometryMap::print();
 
     // Print the one-to-one map
-    LOG(INFO) << "[vertexId1, vertexId2]\n";
+    LOG(INFO) << "[slaveVertId, masterVertexId]\n";
     for (auto const& mapValue : this->m_oneToOneMap)
     {
         LOG(INFO) << "[" << mapValue.first << ", " << mapValue.second << "]\n";
@@ -142,9 +141,9 @@ OneToOneMap::print() const
 void
 OneToOneMap::setMaster(std::shared_ptr<Geometry> master)
 {
-    if (master->getType() == GeometryType::HexahedralMesh ||
-        master->getType() == GeometryType::SurfaceMesh ||
-        master->getType() == GeometryType::TetrahedralMesh)
+    if (master->getType() != GeometryType::HexahedralMesh &&
+        master->getType() != GeometryType::SurfaceMesh &&
+        master->getType() != GeometryType::TetrahedralMesh)
     {
         LOG(WARNING) << "The geometry provided is not a mesh!\n";
         return;
@@ -155,9 +154,9 @@ OneToOneMap::setMaster(std::shared_ptr<Geometry> master)
 void
 OneToOneMap::setSlave(std::shared_ptr<Geometry> slave)
 {
-    if (slave->getType() == GeometryType::HexahedralMesh ||
-        slave->getType() == GeometryType::SurfaceMesh ||
-        slave->getType() == GeometryType::TetrahedralMesh)
+    if (slave->getType() != GeometryType::HexahedralMesh &&
+        slave->getType() != GeometryType::SurfaceMesh &&
+        slave->getType() != GeometryType::TetrahedralMesh)
     {
         LOG(WARNING) << "The geometry provided is not a mesh!\n";
         return;

Base/Geometry/Map/imstkOneToOneMap.h

@@ -49,11 +49,6 @@ public:
     ///
     void compute() override;
 
-    ///
-    /// \brief Find the matching node
-    ///
-    int findMatchingVertex(std::shared_ptr<Mesh> tetraMesh, const Vec3d& p) const;
-
     ///
     /// \brief Check the validity of the map
     ///

Base/Geometry/Map/imstkTetraTriangleMap.cpp

@@ -38,7 +38,7 @@ TetraTriangleMap::compute()
 
     m_verticesEnclosingTetraId.clear();
     m_verticesWeights.clear();
-    for (const Vec3d& surfVertPos : triMesh->getInitialVerticesPositions())
+    for (const Vec3d& surfVertPos : triMesh->getVerticesPositions())
     {
         // Find the enclosing or closest tetrahedron
         int closestTetId = findEnclosingTetrahedron(tetMesh, surfVertPos);

Base/Geometry/Mesh/imstkSurfaceMesh.cpp

@@ -34,9 +34,12 @@ namespace imstk {
         setTrianglesVertices(triangles);
         setTextureCoordinates(texCoords);
 
-        computeVerticesNormals();
-        computeTrianglesNormals();
-        computeVerticesTangents();
+        if (computDerivedData)
+        {
+            computeVerticesNormals();
+            computeTrianglesNormals();
+            computeVerticesTangents();
+        }
     }
 
     void
@@ -260,6 +263,12 @@ SurfaceMesh::getVerticesTangents() const
     return m_verticesTangents;
 }
 
+int
+SurfaceMesh::getNumTriangles() const
+{
+    return this->m_trianglesVertices.size();
+}
+
 const Vec4d&
 SurfaceMesh::getVerticeTangent(size_t i) const
 {
@@ -292,4 +301,23 @@ SurfaceMesh::clear()
     Mesh::clear();
 }
 
+void
+SurfaceMesh::print() const
+{
+    LOG(INFO) << "Number of vertices: " << this->getNumVertices() << "\n";
+    LOG(INFO) << "Number of triangles: " << this->getNumTriangles() << "\n";
+
+    LOG(INFO) << "Triangles:\n";
+    for (auto &triVerts : this->getTrianglesVertices())
+    {
+        LOG(INFO) << "(" << triVerts[0] << ", " << triVerts[1] << "," << triVerts[2] << ")\n";
+    }
+
+    LOG(INFO) << "Vertex positions:\n";
+    for (auto &verts : this->getInitialVerticesPositions())
+    {
+        LOG(INFO) << "(" << verts.x() << ", " << verts.y() << "," << verts.z() << ")\n";
+    }
+}
+
 }

Base/Geometry/Mesh/imstkSurfaceMesh.h

@@ -89,6 +89,11 @@ public:
     ///
     void clear();
 
+    ///
+    /// \brief Print the mesh
+    ///
+    void print() const;
+
     // Accessors
 
     ///
@@ -128,6 +133,12 @@ public:
     ///
     const std::vector<Vec4d>& getVerticesTangents() const;
     const Vec4d& getVerticeTangent(size_t i) const;
+
+    ///
+    /// \brief Returns the number of triangles
+    ///
+    int getNumTriangles() const;
+
 protected:
 
     std::vector<TriangleArray> m_trianglesVertices; ///> Triangle connectivity

Base/Geometry/Mesh/imstkTetrahedralMesh.cpp

@@ -127,55 +127,68 @@ TetrahedralMesh::computeTetrahedronBoundingBox(const size_t& tetId, Vec3d& min,
     max[2] = *std::max_element(arrayz.begin(), arrayz.end());
 }
 
-// TODO: Test pending
 bool
 TetrahedralMesh::extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh)
 {
-    if (surfaceMesh)
+    if (!surfaceMesh)
     {
         LOG(WARNING) << "Cannot extract SurfaceMesh: The surface mesh provided is not instantiated!";
         return false;
     }
 
-    const SurfaceMesh::TriangleArray facePattern[4] = { { 0, 1, 2 }, { 0, 1, 3}, {0, 2, 3}, {1, 2, 3} };
+    std::vector<SurfaceMesh::TriangleArray> facePattern;
+    facePattern.push_back({ { 0, 1, 2 } });
+    facePattern.push_back({ { 0, 1, 3 } });
+    facePattern.push_back({ { 0, 2, 3 } });
+    facePattern.push_back({ { 1, 2, 3 } });
 
     // Find number of common vertices
-    auto getNumCommonVerts = [](const TetraArray& array1, const TetraArray& array2, SurfaceMesh::TriangleArray& commonFace) -> int
+    auto getNumCommonVerts = [facePattern](const TetraArray& array1, const TetraArray& array2, SurfaceMesh::TriangleArray& commonFace) -> int
     {
         int numCommonVerts = 0;
-        SurfaceMesh::TriangleArray tmpFace;
+        std::array<bool, 4> tmpFace = {{0,0,0,0}};
         for (size_t i = 0; i < 4; ++i)
         {
             if (array1[i] == array2[0] || array1[i] == array2[1] || array1[i] == array2[2] || array1[i] == array2[3])
             {
-                tmpFace.at(numCommonVerts) = i;
+                tmpFace[i] = true;
                 numCommonVerts++;
             }
         }
         if (numCommonVerts == 3)
         {
-            commonFace = tmpFace;
+            for (size_t i = 0; i < 4; ++i)
+            {
+                if (!tmpFace[i])
+                {
+                    for (size_t j = 0; j < 3; ++j)
+                    {
+                        commonFace[j] = array1[facePattern[3-i][j]];// this is specific to the above pattern
+                    }
+                }
+            }
         }
         return numCommonVerts;
     };
 
-    // Find the common face
+    // Find the common face irrespecive of the order
     auto findCommonFace = [facePattern](const TetraArray& tetVertArray, const SurfaceMesh::TriangleArray& triVertArray) -> int
     {
-        for (size_t i; i < 4; ++i)
+        for (size_t i = 0; i < 4; ++i)
         {
-            auto vertId = tetVertArray[i];
-            if (vertId != triVertArray[0] && vertId != triVertArray[1] && vertId != triVertArray[2])
+            if (tetVertArray[i] != triVertArray[0] && tetVertArray[i] != triVertArray[1] && tetVertArray[i] != triVertArray[2])
             {
-                for (size_t j; j < 4; ++j)
+                for (size_t j = 0; j < 4; ++j)
                 {
-                    if (i != facePattern[j][0] && i != facePattern[j][1] && i != facePattern[j][2])
+                    if (i != facePattern[j].at(0) && i != facePattern[j].at(1) && i != facePattern[j].at(2))
                     {
-                        return j;
+                        return (int)j;
                     }
                 }
             }
         }
+        LOG(WARNING) << "There is no common face!";
+        return -1;// something wrong if you are here
     };
 
     // Find and store the tetrahedral faces that are unique
@@ -187,13 +200,18 @@ TetrahedralMesh::extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh)
     SurfaceMesh::TriangleArray commonFace;
     bool foundFaces[4];
 
-    for (size_t tetId; tetId < this->getNumTetrahedra(); ++tetId)
+    std::cout << this->getNumTetrahedra() << std::endl;
+    for (size_t tetId = 0; tetId < this->getNumTetrahedra(); ++tetId)
     {
         auto tetVertArray = vertArray.at(tetId);
         foundFaces[0] = foundFaces[1] = foundFaces[2] = foundFaces[3] = false;
 
-        for (size_t tetIdInner; tetIdInner < this->getNumTetrahedra(); ++tetIdInner)
+        for (size_t tetIdInner = 0 ; tetIdInner < this->getNumTetrahedra(); ++tetIdInner)
         {
+            if (tetId == tetIdInner)
+            {
+                continue;
+            }
             auto tetVertArrayInner = vertArray.at(tetIdInner);
 
             // check if there is common face
@@ -246,14 +264,14 @@ TetrahedralMesh::extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh)
         centroid = (v0 + v1 + v2) / 3;
 
         normal = ((v0 - v1).cross(v0 - v2));
-        normal.normalize();
+        //normal.normalize();
 
         if (normal.dot(centroid - this->getVertexPosition(tetRemainingVert.at(faceId))) > 0)
         {
             // swap
-            int tmpIndex = surfaceTri[faceId][3];
-            surfaceTri[faceId][3] = surfaceTri[faceId][2];
-            surfaceTri[faceId][3] = tmpIndex;
+            int tmpIndex = surfaceTri[faceId][2];
+            surfaceTri[faceId][2] = surfaceTri[faceId][1];
+            surfaceTri[faceId][2] = tmpIndex;
         }
     }
 
@@ -265,17 +283,18 @@ TetrahedralMesh::extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh)
         uniqueVertIdList.push_back(face[1]);
         uniqueVertIdList.push_back(face[2]);
     }
+    uniqueVertIdList.sort();
     uniqueVertIdList.unique();
 
     int vertId;
     std::list<int>::iterator it;
     std::vector<Vec3d> vertPositions;
-    for (vertId, it = uniqueVertIdList.begin(); it != uniqueVertIdList.end(); ++vertId, it++)
+    for (vertId=0, it = uniqueVertIdList.begin(); it != uniqueVertIdList.end(); ++vertId, it++)
     {
         vertPositions.push_back(this->getVertexPosition(*it));
         for (auto &face : surfaceTri)
         {
-            for (size_t i = 0; i < 4; ++i)
+            for (size_t i = 0; i < 3; ++i)
             {
                 if (face[i] == *it)
                 {

Base/Geometry/Mesh/imstkTetrahedralMesh.h

@@ -24,6 +24,7 @@
 
 #include <set>
 #include <list>
+#include <iostream>
 
 #include "imstkVolumetricMesh.h"
 
@@ -83,7 +84,10 @@ public:
     void computeTetrahedronBoundingBox(const size_t& tetId, Vec3d& min, Vec3d& max) const;
 
     ///
-    /// \brief Extract the confirming triangular mesh from the tetrahedral mesh
+    /// \brief This method
+    /// (a)Extract the confirming triangular mesh from the tetrahedral mesh
+    /// (b) Checks and flips the triangle connectivity order if it is not consistent
+    /// (c) Renumbers the vertices
     ///
     bool extractSurfaceMesh(std::shared_ptr<SurfaceMesh> surfaceMesh);
 

Examples/Sandbox/main.cpp

@@ -18,9 +18,10 @@
 #include "imstkSurfaceMesh.h"
 #include "imstkMeshReader.h"
 
-// Map
+// Maps
 #include "imstkTetraTriangleMap.h"
 #include "imstkIsometricMap.h"
+#include "imstkOneToOneMap.h"
 
 #include "g3log/g3log.hpp"
 
@@ -30,6 +31,8 @@ void testAnalyticalGeometry();
 void testScenesManagement();
 void testIsometricMap();
 void testTetraTriangleMap();
+void testOneToOneNodalMap();
+void testExtractSurfaceMesh();
 
 int main()
 {
@@ -38,11 +41,13 @@ int main()
               << "****************\n";
 
     //testViewer();
-    testReadMesh();
+    //testReadMesh();
     //testAnalyticalGeometry();
     //testScenesManagement();
     //testIsometricMap();
     //testTetraTriangleMap();
+    testExtractSurfaceMesh();
+    //testOneToOneNodalMap();
 
     return 0;
 }
@@ -298,3 +303,95 @@ void testTetraTriangleMap()
 
     getchar();
 }
+
+void testExtractSurfaceMesh()
+{
+    auto sdk = std::make_shared<imstk::SimulationManager>();
+
+    // a. Construct a sample tetrahedral mesh
+
+    // a.1 add vertex positions
+    auto tetMesh = std::make_shared<imstk::TetrahedralMesh>();
+    std::vector<imstk::Vec3d> vertList;
+    vertList.push_back(imstk::Vec3d(0, 0, 0));
+    vertList.push_back(imstk::Vec3d(1.0, 0, 0));
+    vertList.push_back(imstk::Vec3d(0, 1.0, 0));
+    vertList.push_back(imstk::Vec3d(0, 0, 1.0));
+    vertList.push_back(imstk::Vec3d(1.0, 1.0, 1.0));
+    tetMesh->setInitialVerticesPositions(vertList);
+    tetMesh->setVerticesPositions(vertList);
+
+    // a.2 add connectivity
+    std::vector<imstk::TetrahedralMesh::TetraArray> tetConnectivity;
+    imstk::TetrahedralMesh::TetraArray tet1 = { 0, 1, 2, 3 };
+    imstk::TetrahedralMesh::TetraArray tet2 = { 1, 2, 3, 4 };
+    tetConnectivity.push_back(tet1);
+    tetConnectivity.push_back(tet2);
+    tetMesh->setTetrahedraVertices(tetConnectivity);
+
+    // b. Extract the surface mesh
+    auto extractedSurfMesh = std::make_shared<imstk::SurfaceMesh>();
+    tetMesh->extractSurfaceMesh(extractedSurfMesh);
+
+    // c. Print the resulting mesh
+    extractedSurfMesh->print();
+
+    getchar();
+}
+
+void testOneToOneNodalMap()
+{
+    auto sdk = std::make_shared<imstk::SimulationManager>();
+
+    // a. Construct a sample tetrahedral mesh
+
+    // a.1 add vertex positions
+    auto tetMesh = std::make_shared<imstk::TetrahedralMesh>();
+    std::vector<imstk::Vec3d> vertList;
+    vertList.push_back(imstk::Vec3d(0, 0, 0));
+    vertList.push_back(imstk::Vec3d(1.0, 0, 0));
+    vertList.push_back(imstk::Vec3d(0, 1.0, 0));
+    vertList.push_back(imstk::Vec3d(0, 0, 1.0));
+    vertList.push_back(imstk::Vec3d(1.0, 1.0, 1.0));
+    tetMesh->setInitialVerticesPositions(vertList);
+    tetMesh->setVerticesPositions(vertList);
+
+    // b. Construct a surface mesh
+    auto triMesh = std::make_shared<imstk::SurfaceMesh>();
+
+    // b.1 Add vertex positions
+    std::vector<imstk::Vec3d> SurfVertList;
+    SurfVertList.push_back(imstk::Vec3d(0, 0, 0));
+    SurfVertList.push_back(imstk::Vec3d(1.0, 0, 0));
+    SurfVertList.push_back(imstk::Vec3d(0, 1.0, 0));
+    SurfVertList.push_back(imstk::Vec3d(0, 0, 1.0));
+    SurfVertList.push_back(imstk::Vec3d(1.0, 1.0, 1.0));
+    triMesh->setInitialVerticesPositions(SurfVertList);
+    triMesh->setVerticesPositions(SurfVertList);
+
+    // b.2 Add vertex connectivity
+    std::vector<imstk::SurfaceMesh::TriangleArray> triConnectivity;
+    triConnectivity.push_back({ { 0, 1, 2 } });
+    triConnectivity.push_back({ { 0, 1, 3 } });
+    triConnectivity.push_back({ { 0, 2, 3 } });
+    triConnectivity.push_back({ { 1, 2, 4 } });
+    triConnectivity.push_back({ { 1, 3, 4 } });
+    triConnectivity.push_back({ { 2, 3, 4 } });
+    triMesh->setTrianglesVertices(triConnectivity);
+
+    // c. Construct the one to one nodal map based on the above meshes
+    auto oneToOneNodalMap = std::make_shared<imstk::OneToOneMap>();
+    oneToOneNodalMap->setMaster(tetMesh);
+    oneToOneNodalMap->setSlave(triMesh);
+
+    // d. Compute the map
+    oneToOneNodalMap->compute();
+
+    // e. Print the computed nodal map if valid
+    if (oneToOneNodalMap->isValid())
+    {
+        oneToOneNodalMap->print();
+    }
+
+    getchar();
+}
\ No newline at end of file