feat: Add support for region point definition

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

include/gamer/SurfaceMesh.h

@@ -25,6 +25,7 @@
 #pragma once
 
 #include <iostream>
+#include <limits>
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -52,6 +53,8 @@ struct SMGlobal {
   bool useVolumeConstraint;
   /// Flag that determines if the mesh represents a hole or not
   bool ishole;
+  /// Region point
+  Eigen::Vector3d regionPoint;
 
   /**
    * @brief      Default constructor
@@ -65,7 +68,10 @@ struct SMGlobal {
   SMGlobal(int marker = -1, float volumeConstraint = -1,
            bool useVolumeConstraint = false, bool ishole = false)
       : marker(marker), volumeConstraint(volumeConstraint),
-        useVolumeConstraint(useVolumeConstraint), ishole(ishole) {}
+        useVolumeConstraint(useVolumeConstraint), ishole(ishole),
+        regionPoint(std::numeric_limits<double>::max(),
+                    std::numeric_limits<double>::max(),
+                    std::numeric_limits<double>::max()) {}
 };
 
 struct SMVertex : Vertex {
@@ -545,7 +551,7 @@ void selectFlipEdges(
           // faces. "
           //           << "Returning..." << std::endl;
           gamer_runtime_error("SurfaceMesh is not pseudomanifold. Found "
-                                   "an edge connected to more than 2 faces.");
+                              "an edge connected to more than 2 faces.");
         } else if (up.size() < 2) // Edge is a boundary
         {
           // std::cerr << "This edge participates in fewer than 2
@@ -1102,9 +1108,9 @@ std::vector<std::unique_ptr<SurfaceMesh>> splitSurfaces(SurfaceMesh &mesh);
 
 /**
  * @brief Cache face and vertex normals
- * 
+ *
  * A zero vector normal will be stored instead of raising an error
- * 
+ *
  * @param mesh Surface mesh of interest
  */
 void cacheNormals(SurfaceMesh &mesh);

pygamer/src/SMGlobal.cpp

@@ -19,26 +19,31 @@
 
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <pybind11/eigen.h>
 
 #include "gamer/SurfaceMesh.h"
 
 /// Namespace for all things gamer
-namespace gamer
-{
+namespace gamer {
 
 namespace py = pybind11;
 
-void init_SMGlobal(py::module& mod){
-    py::class_<SMGlobal> global(mod, "Global",
-        R"delim(
+void init_SMGlobal(py::module &mod) {
+  py::class_<SMGlobal> global(mod, "Global",
+                              R"delim(
             Wrapper around a :cpp:class:`gamer::SMGlobal`.
-        )delim"
-    );
+        )delim");
 
-    global.def_readwrite("marker", &SMGlobal::marker, "Domain marker to use when tetrahedralizing.");
-    global.def_readwrite("volumeConstraint", &SMGlobal::volumeConstraint, "Domain volumeConstraint to use when tetrahedralizing.");
-    global.def_readwrite("useVolumeConstraint", &SMGlobal::useVolumeConstraint, "Use a volume constraint when tetrahedralizing.");
-    global.def_readwrite("ishole", &SMGlobal::ishole, "Does this domain represent a hole or not?");
+  global.def_readwrite("marker", &SMGlobal::marker,
+                       "Domain marker to use when tetrahedralizing.");
+  global.def_readwrite("volumeConstraint", &SMGlobal::volumeConstraint,
+                       "Domain volumeConstraint to use when tetrahedralizing.");
+  global.def_readwrite("useVolumeConstraint", &SMGlobal::useVolumeConstraint,
+                       "Use a volume constraint when tetrahedralizing.");
+  global.def_readwrite("ishole", &SMGlobal::ishole,
+                       "Does this domain represent a hole or not?");
+  global.def_readwrite("regionPoint", &SMGlobal::regionPoint,
+                       "Point within the region.");
 }
 
 } // end namespace gamer

src/TetMesh.cpp

@@ -43,10 +43,6 @@ namespace gamer {
 std::unique_ptr<TetMesh>
 makeTetMesh(const std::vector<SurfaceMesh const *> &surfmeshes,
             std::string tetgen_params) {
-
-  // Create new tetmesh object
-  std::unique_ptr<TetMesh> tetmesh(new TetMesh);
-
   size_t nVertices = 0, nFaces = 0, nRegions = 0, nHoles = 0;
   int i = 0;
   for (auto &surfmesh : surfmeshes) {
@@ -151,23 +147,30 @@ makeTetMesh(const std::vector<SurfaceMesh const *> &surfmeshes,
 
     auto metadata = *surfmesh->get_simplex_up();
 
-    // TODO: (25) Improve region point picking strategy
-    // Pick a point inside the region
-    auto faceID = *surfmesh->template get_level_id<3>().begin();
-    Vector normal = getNormal(*surfmesh, faceID);
-    normal /= std::sqrt(normal | normal);
-
-    auto fname = surfmesh->get_name(faceID);
-    Vector a = (*surfmesh->get_simplex_up({fname[0]})).position;
-    Vector b = (*surfmesh->get_simplex_up({fname[1]})).position;
-    Vector c = (*surfmesh->get_simplex_up({fname[2]})).position;
-
-    Vector d = a - b;
-    double weight = std::sqrt(d | d);
-
-    // flip normal and scale by weight
-    normal *= weight;
-    Vector regionPoint((a + b + c) / 3.0 - normal);
+    Vector regionPoint;
+    if (metadata.regionPoint.isApprox(
+            Eigen::Vector3d(std::numeric_limits<double>::max(),
+                            std::numeric_limits<double>::max(),
+                            std::numeric_limits<double>::max()))) {
+      // Pick a point inside the region
+      auto faceID = *surfmesh->template get_level_id<3>().begin();
+      Vector normal = getNormal(*surfmesh, faceID);
+      normal /= std::sqrt(normal | normal);
+
+      auto fname = surfmesh->get_name(faceID);
+      Vector a = (*surfmesh->get_simplex_up({fname[0]})).position;
+      Vector b = (*surfmesh->get_simplex_up({fname[1]})).position;
+      Vector c = (*surfmesh->get_simplex_up({fname[2]})).position;
+
+      Vector d = a - b;
+      double weight = std::sqrt(d | d);
+
+      // flip normal and scale by weight
+      normal *= weight / 2;
+      regionPoint = (a + b + c) / 3.0 - normal;
+    } else {
+      regionPoint = metadata.regionPoint;
+    }
 
     std::cout << "Region point: " << regionPoint << std::endl;
 
@@ -597,6 +600,7 @@ void writeDolfin(const std::string &filename, const TetMesh &mesh) {
   fout << "    </cells>\n";
   fout << "    <domains>\n";
 
+  // Write face markers
   fout << "      <mesh_value_collection name=\"m\" type=\"uint\" dim=\"2\" "
           "size=\""
        << faceMarkerList.size() << "\">\n";
@@ -608,8 +612,9 @@ void writeDolfin(const std::string &filename, const TetMesh &mesh) {
          << " local_entity=\"" << local_entity << "\" "
          << " value=\"" << marker << "\" />\n";
   }
-
   fout << "      </mesh_value_collection>\n";
+
+  // Write cell markers
   fout << "      <mesh_value_collection name=\"m\" type=\"uint\" dim=\"3\" "
           "size=\""
        << mesh.size<4>() << "\">\n";
@@ -623,6 +628,7 @@ void writeDolfin(const std::string &filename, const TetMesh &mesh) {
          << " value=\"" << marker << "\" />\n";
   }
   fout << "      </mesh_value_collection>\n";
+
   fout << "    </domains>\n";
   fout << "  </mesh>\n";
   fout << "</dolfin>\n";

tools/blendgamer/__init__.py.in

@@ -24,7 +24,7 @@ bl_info = {
         "description": "Utilities for generating finite elements simulation compatible meshes",
         "author": "Christopher T. Lee, Justin Laughlin, John B. Moody, Zeyun Yu, Tom Bartol, Johan Hake, and Michael Holst",
         "version": (@PROJECT_VERSION_MAJOR@, @PROJECT_VERSION_MINOR@, @PROJECT_VERSION_PATCH@),
-        "blender": (3, 6, 0),
+        "blender": (3, 6, 5),
         "location": "3D View > Tool Shelf",
         "wiki_url": "https://github.com/ctlee/gamer",
         "tracker_url": "https://github.com/ctlee/gamer/issues",

tools/blendgamer/src/tetrahedralization.py

@@ -41,6 +41,47 @@
 import numpy as np
 
 
+class GAMER_OT_select_domain(bpy.types.Operator):
+    """Select domain"""
+
+    bl_idname = "gamer.select_domain"
+    bl_label = "Select Domain"
+    bl_options = {"INTERNAL"}
+
+    index = StringProperty()
+
+    def execute(self, context):
+        new_active = bpy.data.objects[self.index]
+        context.view_layer.objects.active.select_set(False)
+        context.view_layer.objects.active = new_active
+        new_active.select_set(True)
+        return {"FINISHED"}
+
+
+class GAMER_OT_associate_region_point(bpy.types.Operator):
+    bl_idname = "gamer.associate_region_point"
+    bl_label = "Associate region point"
+    bl_description = "Associate selected region point with currently selected domain"
+    bl_options = {"REGISTER", "UNDO"}
+
+    def execute(self, context):
+        if context.scene.gamer.tet_group.associate_region_point(self.report, context):
+            return {"FINISHED"}
+        else:
+            return {"CANCELLED"}
+
+
+class GAMER_OT_dissociate_region_point(bpy.types.Operator):
+    bl_idname = "gamer.dissociate_region_point"
+    bl_label = "Dissociate region point"
+    bl_description = "Dissociate region point with currently selected domain"
+    bl_options = {"REGISTER", "UNDO"}
+
+    def execute(self, context):
+        context.scene.gamer.tet_group.dissociate_region_point(self.report, context)
+        return {"FINISHED"}
+
+
 class GAMER_OT_tet_domain_add(bpy.types.Operator):
     bl_idname = "gamer.tet_domain_add"
     bl_label = "Add a Tet Domain"
@@ -126,6 +167,13 @@ class GAMerTetDomainPropertyGroup(bpy.types.PropertyGroup):
     object_pointer = PointerProperty(
         type=bpy.types.Object, name="Test", description="Object"
     )
+
+    region_point = PointerProperty(
+        type=bpy.types.Object,
+        name="Region point",
+        description="Empty referencing a point within the region",
+    )
+
     marker = IntProperty(name="Marker", default=-1, description="Domain Marker Integer")
     is_hole = BoolProperty(
         name="Hole", default=False, description="Use this domain as a hole"
@@ -168,7 +216,12 @@ def draw_item_in_row(self, row):
             return
 
         col = row.column()
-        col.label(text=name)
+        # col.label(text=name)
+        col.operator(
+            "gamer.select_domain",
+            text=name,
+            icon="OUTLINER_DATA_MESH",
+        ).index = self.object_pointer.name
         col = row.column()
         col.label(text="Domain ID: " + str(self.domain_id))
         col = row.column()
@@ -177,6 +230,17 @@ def draw_item_in_row(self, row):
         else:
             col.label(text="Domain Marker: " + str(self.marker))
 
+        col = row.column()
+        if self.region_point is not None:
+            # col.label(text=self.region_point.name)
+            col.operator(
+                "gamer.select_domain",
+                text=self.region_point.name,
+                icon="EMPTY_DATA",
+            ).index = self.region_point.name
+        else:
+            col.label(text="Not set")
+
 
 class GAMerTetrahedralizationPropertyGroup(bpy.types.PropertyGroup):
     export_path = StringProperty(
@@ -237,6 +301,24 @@ class GAMerTetrahedralizationPropertyGroup(bpy.types.PropertyGroup):
 
     status = StringProperty(name="status", default="")
 
+    def associate_region_point(self, report, context):
+        """Associate region point with selected domain"""
+        if not context.active_object.type == "EMPTY":
+            report(
+                {"ERROR"}, "Active object denoting region point must be of EMPTY type"
+            )
+            return False
+        # print(self.domain_list[self.active_domain_index].name)
+
+        self.domain_list[self.active_domain_index].region_point = context.active_object
+        report({"INFO"}, "Associated region point with domain")
+        return True
+
+    def dissociate_region_point(self, report, context):
+        """Dissociate region point with selected domain"""
+        self.domain_list[self.active_domain_index].region_point = None
+        report({"INFO"}, "Dissociated region point with domain")
+
     def add_tet_domain(self, report, context):
         """Add a new tet domain to the list of tet domains for each selected object"""
         # From the list of selected objects, only add MESH objects.
@@ -395,6 +477,9 @@ def tetrahedralize(self, context, report):
                     globalInfo.useVolumeConstraint = d.constrain_vol
                     globalInfo.volumeConstraint = d.vol_constraint
 
+                    if d.region_point is not None:
+                        globalInfo.regionPoint = d.region_point.location
+                        # print(d.region_point.location)
                     # Write surface mesh to file for debug
                     # g.writeOFF("surfmesh_%s.off"%(obj_name), gmesh)
 
@@ -403,7 +488,6 @@ def tetrahedralize(self, context, report):
 
             # Tetrahedralize mesh
             if len(gmeshes) > 0:
-
                 quality_str = "q%.4f/%.4fO8/7AYVCa" % (
                     self.max_aspect_ratio,
                     self.min_dihedral,
@@ -444,10 +528,17 @@ def tetrahedralize(self, context, report):
                     except Exception as ex:
                         report({"ERROR"}, str(ex))
 
+                sm = tetmesh.extractSurface()
+
+                g.writeOFF("test.off", sm)
+
         print("######################## End Tetrahedralize ########################")
 
 
 classes = [
+    GAMER_OT_select_domain,
+    GAMER_OT_associate_region_point,
+    GAMER_OT_dissociate_region_point,
     GAMER_OT_tet_domain_add,
     GAMER_OT_tet_domain_remove,
     GAMER_OT_tet_domain_remove_all,

tools/blendgamer/src/ui.py

@@ -522,6 +522,13 @@ def draw(self, context):
             # row = layout.row()
             # row.label ( "Active Index = " + str ( self.active_domain_index ) + ", ID = " + str ( domain.domain_id ) )
 
+            row = layout.row()
+            col = row.column()
+            col.operator("gamer.associate_region_point", icon="EMPTY_DATA")
+            if domain.region_point is not None:
+                col = row.column()
+                col.operator("gamer.dissociate_region_point", icon="X")
+
             domain.draw_layout(layout)
 
             box = layout.box()