Building Voxel-picture out of multiple voxel-pictures - the size of single voxel-pictures should decrease in z-direction while the number of data points in each voxel-picture stays constant

[ghost]

Hello everyone,

with the following code I want to generate a functionally graded gyroid body. The gyroid should be graded in cell size.
To get a very detailed and clean structure an the side where the cells are very small, the data points in this area should increase so that the cleaning and smoothing algorithm can do their job.
Now I've got the problem, that I want to define the number of voxel-pictures in the base of the gyroid, where the struture is not that detailed. The Voxel-size should decrease in z-direction while the number of data points increase.

Does anyone have an idea how to realize that? I already tried in the code with Create_surface() Funktion.
Would be great! Thank you!

'
class Surfaces:

@staticmethod
def gyroid(x,y,z,m,Z_max,Z_min):
    k1 = (m-1)/(Z_max-Z_min)
    C1 = (-Z_min*k1) +1
    C0 = 0.5*k1*(Z_min**2)
    alpha_z = (k1*z) + C1
    beta_z = (k1*z) + C1
    gamma_z = ((k1/2)*z) + C1 + (C0/z)
    return np.sin(beta_z*y)*np.cos(alpha_z*x)+np.sin(gamma_z*z)*np.cos(beta_z*y)+np.sin(alpha_z*x)*np.cos(gamma_z*z)

class Body:

@staticmethod
def create_body(thickness=0.6, dots_per_dimension=200, surface = Surfaces.gyroid, m=1, Z_max=2, Z_min=1, number_slices=6, voxel_size_factor=1.5, unit_block_multiplier=1, target_reduction=0.9, verbose=False):
    '''Creates the surface'''
    gyroid = Body()
    gyroid.set_unit_block_multiplier(unit_block_multiplier)
    gyroid.set_dots_per_dimension(dots_per_dimension)
    gyroid.thickness = thickness
    #gyroid.thickness = gyroid.get_thickness_for_porosity(surface=surface, desired_porosity=porosity, verbose=verbose) --> ändere thickness zu porosity
    gyroid.create_surface(surface=surface, m=m, Z_max=Z_max, Z_min=Z_min, number_slices=number_slices, voxel_size_factor=voxel_size_factor, verbose=verbose)
    gyroid.clean_smooth_surface(verbose=verbose)
    gyroid.reduce_surface(target_reduction=target_reduction, verbose=verbose)
    return gyroid

def __init__(self):
    self.thickness=0.04
    self.set_dots_per_dimension(200)
    self.set_unit_block_multiplier(1)
    self.level_set = None
    self.mask = None

def set_dots_per_dimension(self, n):
    self.nx = self.ny = self.nz = n

def set_unit_block_multiplier(self, a):
    
    self.a = 2*np.pi * a

def create_surface(self, surface=Surfaces.gyroid, m=2, Z_max=2, Z_min=1, number_slices=6, voxel_size_factor=1.5, verbose=False):

    if verbose:
        print("Create Surface.")

    x_start = 0
    x_end = 0
    x_step = self.a/number_slices

    y_start = 0
    y_end = 0
    y_step = self.a/number_slices

    z_start = 0
    z_end = 0
    z_step = self.a/number_slices

    verts_collect = []
    faces_collect = []
    normals_collect = []
    values_collect = []

    while z_end < self.a:

        z_end = z_start + z_step
        if (self.a - z_end) < z_end:
            z_end = self.a

        while y_end < self.a:

            y_end = y_start + y_step
            if (self.a - y_end) < y_end:
                y_end = self.a

            while x_end < self.a:
                
                x_end = x_start + x_step
                if (self.a - x_end) < x_step:
                    x_end = self.a

                self.x, self.y, self.z = np.mgrid[x_start:x_end:self.nx*1j, 
                                                  y_start:y_end:self.ny*1j, 
                                                  z_start:z_end:self.nz*1j]
                
                self.level_set = np.abs(surface(self.x, self.y, self.z, m=m, Z_max=Z_max, Z_min=Z_min))-self.thickness
                self.mask = np.abs(surface(self.x, self.y, self.z, m=m, Z_max=Z_max, Z_min=Z_min))<=self.thickness

                self.level_set[:,:,:]=1
                self.level_set[self.mask]=0
                self.level_set[0,:,:]=1
                self.level_set[:,0,:]=1
                self.level_set[:,:,0]=1
                self.level_set[-1,:,:]=1
                self.level_set[:,-1,:]=1
                self.level_set[:,:,-1]=1

                vertices, faces, normals, values = measure.marching_cubes(self.level_set, 0)
                
                verts_collect.append(vertices)
                faces_collect.append(faces)
                normals_collect.append(normals)
                values_collect.append(values)

                x_start = x_end
                print("Done1")
            
            y_start = y_end
            print("Done2")
        
        z_start = z_end
        z_step /= voxel_size_factor
        
        y_start = 0
        y_step /= voxel_size_factor

        x_start = 0
        x_step /= voxel_size_factor
        print("Done3")

    self.vertices = np.concatenate(verts_collect, axis=0)
    self.faces = np.concatenate(faces_collect, axis=0)
    self.normals = np.concatenate(normals_collect, axis=0)
    self.values = np.concatenate(values_collect, axis=0)
    print("DoneDoneDone")

def volume(self):
    return np.sum(self.mask) / (self.nx*self.ny*self.nz)

def porosity(self):
    return (1-self.volume())*100

def get_thickness_for_porosity(self, surface=Surfaces.gyroid, desired_porosity=70.0, verbose=False):
    max_thickness = np.pi
    min_thickness = 0.01
    
    if verbose:
        print("Determine the thickness for the given porosity.")
    
    self.calculate_level_set(surface=surface, verbose=False)
    current_porosity = self.porosity()
    if verbose:
        print(f"Porosity is {current_porosity}% at thickness {self.thickness}.")

    while abs(current_porosity - desired_porosity) > 0.1:
        mid_thickness = (max_thickness + min_thickness) / 2
        self.thickness = mid_thickness
        self.calculate_level_set(surface=surface, verbose=False)
        current_porosity = self.porosity()

        if current_porosity > desired_porosity:
            min_thickness = mid_thickness
            if verbose:
                print(f"current_porosity: {current_porosity}%, min_thickness: {min_thickness}, max_thickness: {max_thickness}")
        
        elif current_porosity < desired_porosity:
            max_thickness = mid_thickness
            if verbose:
                print(f"current_porosity: {current_porosity}%, min_thickness: {min_thickness}, max_thickness: {max_thickness}")

    if verbose:
        print(f"Thickness of {self.thickness} for a porosity of {desired_porosity}%.")
    
    return self.thickness

def clean_smooth_surface(self, clean=True, smooth=True, verbose=False):    
    if clean:
        if verbose:
            print("Clean.")
        surface_mesh_to_clean = pymeshfix.PyTMesh()
        surface_mesh_to_clean.load_array(self.vertices, self.faces)
        surface_mesh_to_clean.clean(max_iters=10, inner_loops=1)
        if verbose:
            print('There are {:d} boundaries'.format(surface_mesh_to_clean.boundaries()))
        verts_cleaned, faces_cleaned = surface_mesh_to_clean.return_arrays()

    self.vertices = verts_cleaned
    self.faces = faces_cleaned

    if smooth:
        if verbose:
            print("Smoothen.")
        verts_to_smooth = self.vertices
        faces_to_smooth = self.faces
        surface_mesh_to_smooth = trimesh.Trimesh(vertices=verts_to_smooth, faces=faces_to_smooth)
        smoothed_surface_mesh = trimesh.smoothing.filter_mut_dif_laplacian(surface_mesh_to_smooth)
        if verbose:
            print(f"Is watertight: {smoothed_surface_mesh.is_watertight}")
        
    self.surface = smoothed_surface_mesh

def reduce_surface(self, target_reduction=0.9, verbose=False):
    if verbose:
        print("Reduce.")
    if not isinstance(self.surface, pv.PolyData):
        surface_mesh_to_reduce = pv.wrap(self.surface)
    reduced_surface_mesh = surface_mesh_to_reduce.decimate(target_reduction)
    self.surface = reduced_surface_mesh.triangulate()
    if verbose:
        print(f"Mesh has been reduced by {target_reduction * 100}%")

def plot(self):
    self.surface.plot(show_edges=True, color=True, jupyter_backend='ipygany') #'pythreejs', 'ipygany', 'panel'

def save_to_stl(self, filename = "gyroid.stl"):
    self.surface.save(filename)

`