updated geometry tolerance

src/pyhexspline/gmsh_mesh_builder.py

@@ -957,18 +957,19 @@ def mesh_generate(self, dim: int, element_order: int):
         self.option.setNumber("Mesh.Recombine3DLevel", 1)
         self.option.setNumber("Mesh.ElementOrder", element_order)
         self.option.setNumber("Mesh.Smoothing", 100000)
+        self.option.setNumber("Geometry.Tolerance", 1e-3)
 
         self.model.mesh.generate(dim)
         self.model.mesh.removeDuplicateNodes()
         self.model.mesh.removeDuplicateElements()
         self.model.occ.synchronize()
 
+        self.logger.info("Optimising mesh")
         if element_order == 1:
-            pass
-            # self.logger.info("Optimising mesh")
-            # self.model.mesh.optimize(method="HighOrderElastic", force=True)
+            self.model.mesh.optimize(
+                method="UntangleMeshGeometry", niter=20, force=True
+            )
         elif element_order > 1:
-            self.logger.info("Optimising mesh")
             self.model.mesh.optimize(method="HighOrderFastCurving", force=False)
         return None
 
@@ -981,6 +982,7 @@ def analyse_mesh_quality(self, hiding_thresh: float) -> None:
         return None
 
     def gmsh_get_nodes(self) -> Dict[uint64, ndarray]:
+        self.factory.synchronize()
         nTagsCoord_dict = {}
         for physical_group in self.model.getPhysicalGroups():
             nTags_s, nCoord_s = self.model.mesh.getNodesForPhysicalGroup(
@@ -994,6 +996,7 @@ def gmsh_get_nodes(self) -> Dict[uint64, ndarray]:
         return nTagsCoord_dict
 
     def gmsh_get_elms(self) -> Dict[uint64, ndarray]:
+        self.factory.synchronize()
         elms_dict = {}
         nodeTags = None
         for vol in self.model.getEntities(3):

src/pyhexspline/spline_mesher.py

@@ -46,13 +46,24 @@ def __init__(
 
     def make_compound(self, dim: int, tags: list):
         gmsh.model.occ.synchronize()
-        # get a list of the tags of the entities to be combined
-        compound = []
-        for entity in tags:
-            compound.append(entity[0][1])
+        # make sure that tags is a list of integers or numpy array of integers
+        if all(isinstance(item, int) for item in tags) or all(
+            isinstance(item, (int, np.integer)) for item in tags
+        ):
+            compound = tags
+        else:
+            # get a list of the tags of the entities to be combined
+            compound = []
+            for entity in tags:
+                compound.append(entity[0][1])
 
         # make compound of size 'dim' and tags 'compound'
-        gmsh.model.mesh.setCompound(dim, compound)
+        if dim == 1:
+            subcompounds = np.array_split(compound, 10)  # TODO: generalise this
+            for subcomp in subcompounds:
+                gmsh.model.mesh.setCompound(dim, subcomp)
+        else:
+            gmsh.model.mesh.setCompound(dim, compound)
         gmsh.model.occ.synchronize()
         return None
 
@@ -176,9 +187,9 @@ def split_thrusection_entities(self, thrusections_tags):
     def add_thrusection(
         self,
         surf_tags,
-        MAX_DEGREE=2,
+        MAX_DEGREE=1,
         MAKE_SOLID=True,
-        CONTINUITY="G2",
+        CONTINUITY="G0",
     ):
         gmsh.model.occ.synchronize()
         thrusections_volumes = []
@@ -315,9 +326,11 @@ def mesher(
         )
 
         cortical_v.plot_mhd_slice()
-        print(f'Image size before entering binary_threshold(): {cortical_v.sitk_image.GetSize()}')
+        print(
+            f"Image size before entering binary_threshold(): {cortical_v.sitk_image.GetSize()}"
+        )
         image_pad = cortical_v.binary_threshold(img_path=cortical_v.img_path)
-        print(f'Image size after entering binary_threshold(): {image_pad.GetSize()}')
+        print(f"Image size after entering binary_threshold(): {image_pad.GetSize()}")
         cortical_ext_split, cortical_int_split, cortical_int_sanity = (
             cortical_v.volume_spline_fast_implementation(image_pad)
         )
@@ -407,11 +420,10 @@ def mesher(
         # *ThruSections for Cortical Slices
         cort_slice_thrusections_volumes = self.add_thrusection(
             slices_tags,
-            MAX_DEGREE=2,
+            MAX_DEGREE=1,
             MAKE_SOLID=True,
-            CONTINUITY="G2",
+            CONTINUITY="G0",
         )
-
         cort_lines, cort_splines_radial, cort_splines_longitudinal = (
             self.split_thrusection_entities(cort_slice_thrusections_volumes)
         )
@@ -448,9 +460,9 @@ def mesher(
 
         trab_inner_vol_thrusection = self.add_thrusection(
             trab_curveloop_h_list,
-            MAX_DEGREE=2,
+            MAX_DEGREE=1,
             MAKE_SOLID=True,
-            CONTINUITY="G2",
+            CONTINUITY="G0",
         )
 
         trab_inner_lines, trab_inner_splines_radial, trab_inner_splines_longitudinal = (
@@ -475,9 +487,9 @@ def mesher(
         # add thrusection
         trab_slice_vol_thrusections = self.add_thrusection(
             trabecular_slice_curve_loop_tags,
-            MAX_DEGREE=2,
+            MAX_DEGREE=1,
             MAKE_SOLID=True,
-            CONTINUITY="G2",
+            CONTINUITY="G0",
         )
         # split entities
         trab_lines, trab_splines_radial, trab_splines_longitudinal = (
@@ -488,21 +500,21 @@ def mesher(
         radial_lines.extend(trab_splines_radial)
         longitudinal_lines.extend(trab_splines_longitudinal)
 
-        # mesher.factory.synchronize()
-        # all_line_entities = gmsh.model.getEntities(1)
-        # all_line_entities = [line[1] for line in all_line_entities]
-        # lines_to_delete = [
-        #     line
-        #     for line in all_line_entities
-        #     if line not in transverse_lines
-        #     or line not in radial_lines
-        #     or line not in longitudinal_lines
-        #     or line not in cort_lines
-        # ]
-        # print(f"Deleting {len(lines_to_delete)} entities")
-        # [gmsh.model.occ.remove([(1, line)]) for line in lines_to_delete]
-        # # [gmsh.model.occ.remove([(1, line)]) for line in all_line_entities]
-        # mesher.factory.synchronize()
+        mesher.factory.synchronize()
+        all_line_entities = gmsh.model.getEntities(1)
+        all_line_entities = [line[1] for line in all_line_entities]
+        lines_to_delete = [
+            line
+            for line in all_line_entities
+            if line not in transverse_lines
+            or line not in radial_lines
+            or line not in longitudinal_lines
+            or line not in cort_lines
+        ]
+        print(f"Deleting {len(lines_to_delete)} entities")
+        [gmsh.model.occ.remove([(1, line)]) for line in lines_to_delete]
+        # [gmsh.model.occ.remove([(1, line)]) for line in all_line_entities]
+        mesher.factory.synchronize()
 
         trab_refinement = None
         quadref_vols = None
@@ -543,11 +555,34 @@ def mesher(
         )
 
         # * Make ThruSection Compound
-        self.make_compound(1, cortical_ext_bspline)
-        self.make_compound(1, cortical_int_bspline)
-        self.make_compound(2, trab_surfs)
-        self.make_compound(3, cort_slice_thrusections_volumes)
-        # self.make_compound(3, trab_slice_vol_thrusections)
+        cort_thrusection_surfs = []
+        for subvols in cort_slice_thrusections_volumes:
+            _, surfs = mesher.model.getAdjacencies(3, subvols[0][1])
+            cort_thrusection_surfs.extend(surfs)
+        # self.make_compound(1, cortical_ext_bspline)
+        # self.make_compound(1, cortical_int_bspline)
+        # self.make_compound(2, [2, 8, 14, 20])
+        # self.make_compound(2, [4, 10, 16, 22])
+        # self.make_compound(2, [5, 11, 17, 23])
+        # self.make_compound(2, [6, 12, 18, 24])
+        # self.make_compound(3, cort_slice_thrusections_volumes)
+        # self.make_compound(3, [1, 2, 3, 4])
+        # self.make_compound(3, [5, 6, 7, 8])
+
+        # gmsh.model.mesh.setCompound(1, [73, 116])
+        # gmsh.model.mesh.setCompound(1, [75, 118])
+        # gmsh.model.mesh.setCompound(1, [74, 114])
+        # gmsh.model.mesh.setCompound(1, [76, 117])
+        # gmsh.model.mesh.setCompound(2, [1, 21])
+        # gmsh.model.mesh.setCompound(2, [4, 22])
+        # gmsh.model.mesh.setCompound(2, [2, 20])
+
+        # gmsh.model.mesh.setCompound(1, [83, 119])
+        # gmsh.model.mesh.setCompound(1, [84, 120])
+        # gmsh.model.mesh.setCompound(1, [77, 113])
+        # gmsh.model.mesh.setCompound(1, [79, 115])
+        # gmsh.model.mesh.setCompound(2, [5, 23])
+        # gmsh.model.mesh.setCompound(2, [6, 24])
 
         # * make transfinite
         mesher.make_thrusections_transfinite(
@@ -571,7 +606,6 @@ def mesher(
         )
 
         mesher.factory.synchronize()
-        gmsh.model.geo.removeAllDuplicates()
         mesher.mesh_generate(dim=3, element_order=ELM_ORDER)
 
         if MESH_ANALYSIS:

src/pyhexspline/spline_volume.py

@@ -903,7 +903,7 @@ def evaluate_bspline(self, contour_coords):
         tckp, _ = splprep(
             [contour_coords[:, 0], contour_coords[:, 1]],
             s=self.S,
-            k=3,
+            k=self.K,
             per=1,
         )
         x_new, y_new = splev(np.linspace(0, 1, self.INTERP_POINTS_S), tckp)

standalone.py

@@ -41,21 +41,21 @@ def main():
         "aspect": 100,  # aspect ratio of the plots
         "_slice": 1,  # slice of the image to be plotted
         "undersampling": 1,  # undersampling factor of the image
-        "slicing_coefficient": 5,  # using every nth slice of the image for the spline reconstruction
+        "slicing_coefficient": 6,  # using every nth slice of the image for the spline reconstruction
         "inside_val": int(0),  # threshold value for the inside of the mask
         "outside_val": int(1),  # threshold value for the outside of the mask
         "lower_thresh": float(0),  # lower threshold for the mask
         "upper_thresh": float(0.9),  # upper threshold for the mask
-        "s": 1000,  # smoothing factor of the spline
+        "s": 700,  # smoothing factor of the spline
         "k": 3,  # degree of the spline
         "interp_points": 500,  # number of points to interpolate the spline
-        "thickness_tol": 5e-1,  # minimum cortical thickness tolerance: 3 * XCTII voxel size
+        "thickness_tol": 1,  # minimum cortical thickness tolerance: 3 * XCTII voxel size
         "phases": 2,  # 1: only external contour, 2: external and internal contour
         "center_square_length_factor": 0.4,  # size ratio of the refinement square: 0 < l_f < 1
         "mesh_order": 1,  # set element order (1: linear, 2: quadratic, >2: higher order, not tested)
-        "n_elms_longitudinal": 60,  # number of elements in the longitudinal direction
-        "n_elms_transverse_trab": 15,  # number of elements in the transverse direction for the trabecular compartment
-        "n_elms_transverse_cort": 4,  # number of elements in the transverse direction for the cortical compartment
+        "n_elms_longitudinal": 30,  # number of elements in the longitudinal direction
+        "n_elms_transverse_trab": 20,  # number of elements in the transverse direction for the trabecular compartment
+        "n_elms_transverse_cort": 3,  # number of elements in the transverse direction for the cortical compartment
         "n_elms_radial": 20,  # number of elements in the radial direction # ! Should be 10 if trab_refinement is True
         "ellipsoid_fitting": True,  # True: perform ellipsoid fitting
         "show_plots": False,  # show plots during construction
@@ -69,14 +69,14 @@ def main():
     #     path_np_s="/Users/msb/Documents/01_PHD/03_Methods/Meshing/01_AIM/C0003094/C0003094_CORT_MASK_UNCOMP.npy"
     # )
 
-    # sitk_image_s = sitk.ReadImage(
-    #     "99_testing_prototyping/repro_sim_errors/C0001643_CORTMASK.mhd"
-    # )
-
     sitk_image_s = sitk.ReadImage(
-        "/home/simoneponcioni/Desktop/radius-validation-ubelix/C0002209_CORTMASK.mhd"
+        "99_testing_prototyping/repro_sim_errors/C0001643_CORTMASK.mhd"
     )
 
+    # sitk_image_s = sitk.ReadImage(
+    #     "/home/simoneponcioni/Desktop/radius-validation-ubelix/C0002209_CORTMASK.mhd"
+    # )
+
     sitk_image_s = sitk_image_s[:, :, 5:-5]
     print(sitk_image_s.GetSize())