feat: Add export of mean curvature for dolfin

Signed-off-by: Christopher T. Lee <[email protected]>

include/gamer/TetMesh.h

@@ -79,7 +79,7 @@ struct TMVertexProperties {
 struct TMVertex : Vertex, TMVertexProperties {
   TMVertex() : TMVertex(Vertex(), TMVertexProperties()) {}
   template <typename... Args>
-  TMVertex(Args &&... args) : TMVertex(Vertex(std::forward<Args>(args)...)) {}
+  TMVertex(Args &&...args) : TMVertex(Vertex(std::forward<Args>(args)...)) {}
   TMVertex(Vertex v) : TMVertex(v, TMVertexProperties(-1)) {}
   TMVertex(Vertex v, TMVertexProperties p) : Vertex(v), TMVertexProperties(p) {}
 
@@ -323,6 +323,16 @@ makeTetMesh(const std::vector<SurfaceMesh const *> &surfmeshes,
  */
 std::unique_ptr<SurfaceMesh> extractSurface(const TetMesh &mesh);
 
+/**
+ * @brief      Extracts the boundary surfaces of a tetrahedral mesh from
+ * boundary markers
+ *
+ * @param[in]  mesh  The mesh
+ *
+ * @return     Bounding surface meshes
+ */
+std::unique_ptr<SurfaceMesh> extractSurfaceFromBoundary(const TetMesh &mesh);
+
 /**
  * @brief      Laplacian smoothing of tetrahedral mesh
  *
@@ -382,4 +392,15 @@ void writeTriangle(const std::string &filename, const TetMesh &mesh);
  * @return     Tetrahedral mesh
  */
 std::unique_ptr<TetMesh> readDolfin(const std::string &filename);
+
+/**
+ * @brief Compute curvatures and write them to dolfin file
+ *
+ * @param filename The filename
+ * @param mesh     The mesh
+ * @param tetmesh  Tetrahedral mesh
+ */
+void curvatureMDSBtoDolfin(const std::string &filename, const SurfaceMesh &mesh,
+                           const TetMesh &tetmesh);
+
 } // end namespace gamer

pygamer/src/TetMesh.cpp

@@ -517,6 +517,19 @@ void init_TetMesh(py::module& mod){
         )delim"
     );
 
+      TetMeshCls.def("extractSurfaceFromBoundary",
+        &extractSurfaceFromBoundary,
+        R"delim(
+            Extract the surface of the TetMesh from boundary markers
+
+            Args:
+                tetmesh (TetMesh): Tetrahedral mesh to extract from.
+
+            Returns:
+                :py:class:`SurfaceMesh`: Surface meshes.
+        )delim"
+    );
+
     /************************************
      *  ITERATORS
      ************************************/

pygamer/src/pygamer.cpp

@@ -325,6 +325,18 @@ PYBIND11_MODULE(pygamer, pygamer) {
         )delim"
     );
 
+    pygamer.def("curvatureMDSBtoDolfin", &curvatureMDSBtoDolfin,
+            py::arg("filename"), py::arg("mesh"), py::arg("tetmesh"),
+        R"delim(
+            Write curvature to dolfin
+
+            Args:
+                filename (:py:class:`str`): Filename to write to
+                mesh (:py:class:`surfacemesh.SurfaceMesh`): list of surface meshes
+                tetmesh (:py:class:`tetmesh.TetMesh`): Tet mesh
+        )delim"
+    );
+
     pygamer.def("__version__",
         [](){
             extern const std::string gVERSION;

src/TetMesh.cpp

@@ -385,6 +385,32 @@ std::unique_ptr<SurfaceMesh> extractSurface(const TetMesh &tetmesh) {
   return surfmesh;
 }
 
+std::unique_ptr<SurfaceMesh>
+extractSurfaceFromBoundary(const TetMesh &tetmesh) {
+  std::unique_ptr<SurfaceMesh> mesh(new SurfaceMesh);
+  std::vector<int> keys;
+
+  for (auto faceID : tetmesh.get_level_id<3>()) {
+    auto name = tetmesh.get_name(faceID);
+    auto data = *faceID;
+    if (data.marker != 0) {
+      mesh->insert(name, SMFace(data.marker, false));
+    }
+  }
+
+  for (auto vertexID : mesh->get_level_id<1>()) {
+    auto &data = *vertexID;
+    auto name = mesh->get_name(vertexID); // Same as in tetmesh
+    data.position = (*tetmesh.get_simplex_up(name)).position;
+  }
+  casc::compute_orientation(*mesh);
+
+  if (getVolume(*mesh) < 0)
+    flipNormals(*mesh);
+
+  return mesh;
+}
+
 void writeVTK(const std::string &filename, const TetMesh &mesh) {
   std::ofstream fout(filename);
   if (!fout.is_open()) {
@@ -884,4 +910,66 @@ std::unique_ptr<TetMesh> readDolfin(const std::string &filename) {
   return mesh;
 }
 
+void curvatureMDSBtoDolfin(const std::string &filename, const SurfaceMesh &mesh,
+                           const TetMesh &tetmesh) {
+  double *kh;
+  double *kg;
+  double *k1;
+  double *k2;
+  std::map<typename SurfaceMesh::KeyType, typename SurfaceMesh::KeyType> sigma;
+  std::tie(kh, kg, k1, k2, sigma) = curvatureViaMDSB(mesh);
+
+  std::ofstream fout(filename + "kh.xml");
+  if (!fout.is_open()) {
+    std::stringstream ss;
+    ss << "File '" << filename << "' could not be written to.";
+    gamer_runtime_error(ss.str());
+  }
+
+  fout << "<?xml version=\"1.0\"?>\n"
+       << "<dolfin xmlns:dolfin=\"http://fenicsproject.org\">\n";
+  fout << "  <mesh_function>\n";
+  fout << "    <mesh_value_collection name=\"kh\" type=\"double\" dim=\"0\" "
+          "size=\""
+       << mesh.size<1>() << "\">\n";
+
+  // hacky way to regenerate cell ids...
+  std::map<TetMesh::SimplexID<4>, std::size_t> simplex_map;
+  std::size_t cnt = 0;
+  for (const auto tetID : tetmesh.get_level_id<4>()) {
+    simplex_map.emplace(tetID, cnt++);
+  }
+
+  for (const auto vertexID : mesh.get_level_id<1>()) {
+    auto idx = mesh.get_name(vertexID)[0];
+
+    auto tet =
+        *(tetmesh.up(tetmesh.up(tetmesh.up(tetmesh.get_simplex_up({idx}))))
+              .begin());
+    auto tetname = tetmesh.get_name(tet);
+    std::size_t local_entity = 0;
+    for (; local_entity < 4; ++local_entity) {
+      if (tetname[local_entity] == idx)
+        break;
+    }
+
+    // std::cout << "vid:" << idx << " cid:" << simplex_map[tet] << " "
+    //           << casc::to_string(tetname) << " " << local_entity <<
+    //           std::endl;
+
+    fout << "        <value cell_index=\"" << simplex_map[tet] << "\" "
+         << " local_entity=\"" << local_entity << "\" "
+         << " value=\"" << kh[sigma[idx]] << "\" />\n";
+  }
+
+  fout << "    </mesh_value_collection>\n";
+  fout << "  </mesh_function>\n";
+  fout << "</dolfin>\n";
+
+  delete[] kh;
+  delete[] kg;
+  delete[] k2;
+  delete[] k1;
+}
+
 } // end namespace gamer

tools/blendgamer/src/tetrahedralization.py

@@ -299,6 +299,12 @@ class GAMerTetrahedralizationPropertyGroup(bpy.types.PropertyGroup):
         description="(Untested) Generate Comsol mphtxt output",
     )
 
+    export_mean_curvature = BoolProperty(
+        name="Dolfin kH",
+        default=False,
+        description="(Untested) Write out mean curvatures in dolfin format",
+    )
+
     status = StringProperty(name="status", default="")
 
     def associate_region_point(self, report, context):
@@ -528,9 +534,11 @@ def tetrahedralize(self, context, report):
                     except Exception as ex:
                         report({"ERROR"}, str(ex))
 
-                sm = tetmesh.extractSurface()
-
-                g.writeOFF("test.off", sm)
+                if self.export_mean_curvature:
+                    surfaces = tetmesh.extractSurfaceFromBoundary()
+                    g.curvatureMDSBtoDolfin(
+                        f"{filename}_mean_curvature", surfaces, tetmesh
+                    )
 
         print("######################## End Tetrahedralize ########################")
 

tools/blendgamer/src/ui.py

@@ -566,6 +566,8 @@ def draw(self, context):
                 col.prop(tetprops, "paraview")
                 col = row.column()
                 col.prop(tetprops, "comsol")
+                col = row.column()
+                col.prop(tetprops, "export_mean_curvature")
 
             row = layout.row()
             icon = "PROP_OFF"