From 932f56f081d4d7d9d6c9a3bb3b086861a626ab70 Mon Sep 17 00:00:00 2001 From: Simone Poncioni Date: Sat, 22 Jun 2024 11:50:23 +0200 Subject: [PATCH] updated geometry tolerance --- src/pyhexspline/gmsh_mesh_builder.py | 11 ++- src/pyhexspline/spline_mesher.py | 108 ++++++++++++++++++--------- src/pyhexspline/spline_volume.py | 2 +- standalone.py | 22 +++--- 4 files changed, 90 insertions(+), 53 deletions(-) diff --git a/src/pyhexspline/gmsh_mesh_builder.py b/src/pyhexspline/gmsh_mesh_builder.py index 718b2b9..ce75308 100644 --- a/src/pyhexspline/gmsh_mesh_builder.py +++ b/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): diff --git a/src/pyhexspline/spline_mesher.py b/src/pyhexspline/spline_mesher.py index 295b65d..5ddbe35 100644 --- a/src/pyhexspline/spline_mesher.py +++ b/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: diff --git a/src/pyhexspline/spline_volume.py b/src/pyhexspline/spline_volume.py index 2397aa7..32f99f9 100644 --- a/src/pyhexspline/spline_volume.py +++ b/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) diff --git a/standalone.py b/standalone.py index c5aa299..72e9811 100644 --- a/standalone.py +++ b/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())