Merge pull request #5 from ctlee/development

Fixing tetrahedralization in Blender

.travis.yml

@@ -13,7 +13,7 @@ branches:
 matrix:
   include:
     - os: osx
-      osx_image: xcode9.3beta
+      osx_image: xcode
     - os: linux
       addons:
         apt:

examples/MeshProtein.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+
+import pygamer as g
+
+# Generate the solvent excluded  molecular surface
+mesh = g.ReadPDB_molsurf('2jho.pdb')
+# mesh = g.ReadPDB_gauss('2jho.pdb', -0.2, 2.5) # filename, blobbiness, isovalue
+
+# Set outward facing normals
+g.compute_orientation(mesh)
+g.correctNormals(mesh)
+
+# Get the root node and set some metadata
+gInfo = mesh.getRoot()
+gInfo.closed = True
+gInfo.ishole = True # Don't tetrahedralize this domain
+gInfo.marker = -1
+
+# Perform some coarsening and smoothing operations
+g.coarse(mesh, 0.016, 1, 0)
+g.coarse(mesh, 2.5, 0, 10)
+g.smoothMesh(mesh, 15, 165, 5, True)
+
+
+# Set boundary markers of the mesh to 23
+for faceID in mesh.faceIDs():
+    faceID.data().marker = 23
+
+# for face in mesh.faces():
+#     print(face.marker)
+
+# Generate a surrounding box
+box = g.meshCube(5)
+
+# Set box boundary markers to 50
+for faceID in box.faceIDs():
+    faceID.data().marker = 50
+
+# Get and set box metadata
+gInfo = box.getRoot()
+gInfo.closed = True
+gInfo.ishole = False
+gInfo.marker = 5
+
+# Get the radius of the protein mesh and scale the box to be big enough.
+radius = g.getRadius(mesh)
+g.scale(box, radius*2)
+
+
+# Generate list of meshes
+meshes = g.VectorSM(); # Holder for list of meshes
+meshes.push_back(mesh)
+meshes.push_back(box)
+
+# Call tetgen
+tetmesh = g.MakeTetMesh(meshes, "pq1.4zYAQ")
+
+g.writeVTK('outputTetMesh.vtk', tetmesh)
+g.writeDolfin('outputTetMesh.xml', tetmesh)

include/SurfaceMesh.h

@@ -111,12 +111,6 @@ struct Global
     float volumeConstraint;
     /// flag that determines if the volume constraint is used.
     bool  useVolumeConstraint;
-    /// Minimum coordinate of vertices
-    float min[3];
-    /// Max coordinate of vertices
-    float max[3];
-    /// Average edge length
-    float avglen;
     /// Flag that determines if the mesh represents a hole or not
     bool  ishole;
 };

include/TetMesh.h

@@ -99,6 +99,14 @@ struct Global
     bool higher_order;
 };
 
+struct TetVertex : Vertex
+{
+    using Vertex::Vertex; // Inherit constructor from Vertex
+    TetVertex() {}
+    TetVertex(Vertex v) : Vertex(v) {}
+    double error = 0;
+};
+
 /**
  * @brief      A helper struct containing the traits/types in the simplicial
  *             complex
@@ -108,7 +116,7 @@ struct complex_traits
     /// The index type
     using KeyType = int;
     /// The types of each node
-    using NodeTypes = util::type_holder<Global, Vertex, Edge, Face, Cell>;
+    using NodeTypes = util::type_holder<Global, TetVertex, Edge, Face, Cell>;
     /// The types of each edge
     using EdgeTypes = util::type_holder<casc::Orientable, casc::Orientable, casc::Orientable, casc::Orientable>;
 };
@@ -131,11 +139,11 @@ std::unique_ptr<TetMesh> makeTetMesh(
 #endif //SWIG
 
 
+void smoothMesh(TetMesh & mesh);
 void writeVTK(const std::string& filename, const TetMesh &mesh);
 void writeOFF(const std::string& filename, const TetMesh &mesh);
-
 void writeDolfin(const std::string &filename, const TetMesh &mesh);
-
+void writeTriangle(const std::string &filename, const TetMesh &mesh);
 
 // void writeMCSF(const std::string &filename, const TetMesh &mesh);
 //void writeDiffPack

include/stringutil.h

@@ -29,6 +29,7 @@
 #include <cctype>
 #include <functional>
 #include <string>
+#include <vector>
 
 // TODO: (2) Do we need this stringutil for gamer?
 namespace stringutil

libraries/CMakeLists.txt

@@ -26,4 +26,4 @@ if( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/casc/.git" )
 )
 endif()
 add_subdirectory(tetgen)
-# add_subdirectory(triangle)
+add_subdirectory(triangle)

src/TetMesh.cpp

@@ -48,7 +48,6 @@ std::unique_ptr<TetMesh> makeTetMesh(
 
     size_t nVertices = 0, nFaces = 0, nRegions = 0, nHoles = 0;
     int i = 0;
-
     for (auto &surfmesh : surfmeshes){
         size_t nverts = surfmesh->template size<1>();
         size_t nfaces = surfmesh->template size<3>();
@@ -157,9 +156,9 @@ std::unique_ptr<TetMesh> makeTetMesh(
         normal /= std::sqrt(normal|normal);
 
         auto fname = surfmesh->get_name(faceID);
-        Vertex a = *surfmesh->get_simplex_up({fname[0]});
-        Vertex b = *surfmesh->get_simplex_up({fname[1]});
-        Vertex c = *surfmesh->get_simplex_up({fname[2]});
+        auto a = *surfmesh->get_simplex_up({fname[0]});
+        auto b = *surfmesh->get_simplex_up({fname[1]});
+        auto c = *surfmesh->get_simplex_up({fname[2]});
 
         Vector d = a-b;
         double weight = std::sqrt(d|d);
@@ -279,7 +278,7 @@ std::unique_ptr<TetMesh> tetgenioToTetMesh(tetgenio &tetio){
         if (vertex != nullptr){
             auto &vdata = *vertex;
             // std::cout << casc::to_string(std::array<double,3>({ptr[0],ptr[1],ptr[2]})) << std::endl;
-            vdata = Vertex(ptr[0], ptr[1], ptr[2], tetio.pointmarkerlist[i], false);
+            vdata = tetmesh::TetVertex(ptr[0], ptr[1], ptr[2], tetio.pointmarkerlist[i], false);
         }
     }
 
@@ -553,5 +552,88 @@ void writeDolfin(const std::string &filename, const TetMesh &mesh){
     fout.close();
 }
 
+void writeTriangle(const std::string &filename, const TetMesh &mesh){
+    std::ofstream fout(filename + ".node");
+    if(!fout.is_open())
+    {
+        std::cerr   << "File '" << filename + ".node"
+                    << "' could not be writen to." << std::endl;
+        exit(1);
+    }
+
+    std::map<typename TetMesh::KeyType,typename TetMesh::KeyType> sigma;
+    size_t cnt = 1;
+
+    // Print out Vertices
+    // std::cout << "Printing Vertices" << std::endl;
+    // nVertices, dimension, nattributes, nmarkers
+    fout << mesh.size<1>() << " 3 " << " 0 " << " 1\n";
+
+    fout.precision(6);
+    for(const auto vertexID : mesh.get_level_id<1>()){
+        size_t idx = cnt;
+        sigma[mesh.get_name(vertexID)[0]] = cnt++;
+        auto vertex = *vertexID;
+
+        fout    << idx << " "
+                << vertex[0] << " "
+                << vertex[1] << " "
+                << vertex[2] << " "
+                << (*vertexID).marker << "\n";
+    }
+    fout.close(); // Close .node file
+
+
+    // Open file filename.ele
+    std::ofstream foutEle(filename + ".ele");
+    if (!foutEle.is_open())
+    {
+        std::cerr   << "File '" << filename + ".ele"
+                    << "' could not be writen to." << std::endl;
+        exit(1);
+    }
+
+    // nTetrahedra, nodes per tet, nAttributes
+    foutEle << mesh.size<4>() << " 4 1\n";
+    cnt = 1;
+    for (const auto tetID :  mesh.get_level_id<4>()){
+        std::size_t idx = cnt++;
+        auto tetName = mesh.get_name(tetID);
+        foutEle << idx << " "
+                << sigma[tetName[0]] << " "
+                << sigma[tetName[1]] << " "
+                << sigma[tetName[2]] << " "
+                << sigma[tetName[3]] << " "
+                << (*tetID).marker << "\n";
+    }
+    foutEle.close(); // Close .ele file
+}
+
+void smoothMesh(TetMesh &mesh){
+    std::set<TetMesh::SimplexID<1>> vertexIDs;
+
+    for (auto faceID : mesh.get_level_id<3>()){
+        if(mesh.up(faceID).size() == 1){
+            auto bdryVertices = mesh.down(mesh.down(faceID)); // Set of vertices
+            vertexIDs.insert(bdryVertices.begin(), bdryVertices.end());
+        }
+    }
+
+    for (auto vID :  mesh.get_level_id<1>()){
+        if (vertexIDs.find(vID) == vertexIDs.end()){
+            std::set<TetMesh::SimplexID<1>> nbhd;
+            neighbors_up(mesh, vID, std::inserter(nbhd, nbhd.end()));
+
+            Vector barycenter;
+            for(auto nbor : nbhd){
+                barycenter += *nbor;
+            }
+            barycenter /= nbhd.size();
+            auto& pos  = *vID;
+            pos = barycenter;
+        }
+    }
+}
+
 
 

swig/TetMesh.i

@@ -34,7 +34,6 @@ wrap_unique_ptr(TMUniquePtr, TetMesh);
 
 %include "TetMesh.h"
 
-
 %inline %{
 TetMesh* MakeTetMesh(const std::vector<SurfaceMesh*> &surfmeshes, std::string tetgen_params){
   return makeTetMesh(surfmeshes, tetgen_params).release();

tools/blender_addon/src/markers.py

@@ -454,7 +454,6 @@ def copyBoundaries(self, fromObject, toObject):
 class GAMerBoundaryMarkersList(bpy.types.PropertyGroup):
     boundary_list = CollectionProperty(type=GAMerBoundaryMarker, name="Boundary List")
     active_bnd_index = IntProperty(name="Active Boundary Index", default=0)
-    include = BoolProperty(name="Include Domain in Model", default=False)
     get_boundary_info = BoolProperty(
         name="Toggle to enable/disable boundary_info", default=False)