Generate Minimal Surfaces with Intersection of an Input Geometry

gen_minimal_surfaces.py

''' 
Initial code from [MaSMaker](https://github.com/CONMAD-CIDESIMX/MaSMaker)
Reference: Tenorio-Suárez, M. I., Gómez-Ortega, A., Canales, H., Piedra, S., & Pérez-Barrera, J. (2022). MaSMaker: An open-source, portable software to create and integrate maze-like surfaces into arbitrary geometries. SoftwareX, 19, 101203.

Adapted by [Yang Liu](yangliu.mit.edu) in HTM(A)A, Fall 2022.

'''
import os
import trimesh
import pycork
import numpy as np
from pathlib import Path
from numpy import sin, cos, pi
from skimage.measure import marching_cubes

# Triply periodic minimal surface (TPMS) Algebraic Equations
def tpms_equation(x, y, z, a, t, type="gyroid"):
    cox = cos(2.0*pi*x/a)
    siy = sin(2.0*pi*y/a)
    coy = cos(2.0*pi*y/a)
    siz = sin(2.0*pi*z/a)
    coz = cos(2.0*pi*z/a)
    six = sin(2.0*pi*x/a)

    co2x = cos(4.0*pi*x/a)
    co2y = cos(4.0*pi*y/a)
    co2z = cos(4.0*pi*z/a)
    si2x = sin(4.0*pi*x/a)
    si2y = sin(4.0*pi*y/a)
    si2z = sin(4.0*pi*z/a)

    if type == "gyroid":
        # Gyroid surface [Schoen-Gyroid]: sin X cos Y + sin Y cos Z + sin Z cos X = c
        return (six*coy + siy*coz + siz*cox)**2 - t**2
    elif type == "diamond":
        # Diamond surface [Schwarz-Diamond]: cos X cos Y cos Z −sin X sin Y sin Z = c
        return (cox*coy*coz - six*siy*siz)**2 - t**2
    elif type == "primitive":
        # Primitive surface [Schwarz-Primitive]: cos X + cos Y + cos Z = c
        return (cox + coy + coz)**2 - t**2
    elif type == "iwp":
        # I-WP surface [Schwarz-IWP]: 2*(cos X cos Y + cos Y cos Z + cos Z cos X) - (cos 2X + cos 2Y + cos 2Z) = c
        return (2*(cox*coy + coy*coz + coz*cox) - (co2x + co2y + co2z))**2 - t**2
    elif type == "neovius":
        # Neovius surface [Neovius]: 3*(cos X + cos Y + cos Z) + 4*cos X cos Y cos Z = c
        return (3*(cox + coy + coz) - 4*cox*coy*coz)**2 - t**2
    elif type == "s":
        # S surface [Fischer Koch S]: cos 2X sin Y cos Z + cos X cos 2Y sin Z + sin X cos Y cos 2Z = c
        return (co2x*siy*coz + cox*co2y*siz + six*coy*co2z)**2 - t**2
    elif type == "frd":
        # F-RD surface [Schoen-FRD]: 4(cos X cos Y cos Z) − (cos 2 X cos 2Y + cos 2Y cos 2Z + cos 2Z cos 2X ) = c
        return (4*cox*coy*coz - (co2x*co2y + co2y*co2z + co2z*co2x))**2 - t**2
    elif type == "pmy":
        # PMY surface [PMY]: 2 cos X cos Y cos Z + sin 2X sin Y + sin X sin 2Z + sin 2Y sin Z = c
        return (2*cox*coy*coz + si2x*siy + six*si2z + si2y*siz)**2 - t**2

path = "./common-3d-test-models/data/" 
name = "spot.obj"
# name = "bimba.obj"

path_out = path + "output/"
Path(path_out).mkdir(parents=True, exist_ok=True)

# [0] Minimal Surfaces Parameters
type = "gyroid"  # Minimal surfaces type - "gyroid" / 
# type = "s"  # Minimal surfaces type
cell_size = 2e-4  # [m] Unit cell size
sheet_thickness = 0.7  # "t" value for the gyroid equation. Level-set value for the gyroid surface   
resolution = 20 # Resolution of a single unit cell. Number of voxels [initially 20]                    
res = resolution*1j 
# Size correction and number of units
cell_size = 1000 * cell_size  # Factor to get the correct size

# [1] Input Geometry_input
meshg=trimesh.load_mesh(path+name)

xmax=meshg.bounds[:,0][1]
xmin=meshg.bounds[:,0][0]
ymax=meshg.bounds[:,1][1]
ymin=meshg.bounds[:,1][0]
zmax=meshg.bounds[:,2][1]
zmin=meshg.bounds[:,2][0]

cenx = (xmax + xmin) / 2
ceny = (ymax + ymin) / 2
cenz = (zmax + zmin) / 2
centro = np.array([cenx, ceny, cenz])

meshg.vertices -= centro

total_l = (xmax-xmin)+1.2
total_w = (ymax-ymin)+1.2
total_h = (zmax-zmin)+1.2

geometry_input=trimesh.exchange.off.export_off(meshg, digits=3)
with open(path+"geometry_input.off", "w") as text_file:
    text_file.write("%s" % geometry_input)

nx = total_l / cell_size  
ny = total_w / cell_size 
nz = total_h / cell_size 
  
# [2] Marching Cubes to Create Minimal surfaces
X, Y, Z = np.mgrid[0:total_l:res * nx, 0:total_w:res * ny, 0:total_h:res * nz]
vol = tpms_equation(X, Y, Z, cell_size, sheet_thickness, type)
verts, faces, normals, values = marching_cubes(vol, 0, spacing=(total_l / (int(nx * resolution)-1), total_w / (int(ny * resolution)-1), total_h / (int(nz * resolution)-1)))  #
meshf = trimesh.Trimesh(vertices=verts[:], faces=faces[:], vertex_normals=normals[:])
meshf.vertices -= meshf.centroid        
meshoff=trimesh.exchange.off.export_off(meshf, digits=3)
with open(path+"minimal_surfaces.off", "w") as text_file:
    text_file.write("%s" % meshoff)


# [3] Boolean Intersection using Pycork
meshA = trimesh.load_mesh(path+"geometry_input.off", process=False)
meshB = trimesh.load_mesh(path+"minimal_surfaces.off", process=False)
vertsA = meshA.vertices
trisA = meshA.faces
vertsB = meshB.vertices
trisB = meshB.faces
vertsD, trisD = pycork.intersection(vertsA, trisA, vertsB, trisB)
meshC = trimesh.Trimesh(vertices=vertsD, faces=trisD, process=False)
meshC.export(path_out + Path(name).stem + "_tpms.stl")
os.remove(path + "geometry_input.off")
os.remove(path + "minimal_surfaces.off")

# [4] Output Minimal Surfaces Properties
area = round((meshC.area),4) # [mm^2] Total surface 
vol  = round((meshC.volume),4) # [mm^2] Total volume 
RD   = round(meshC.volume/meshA.volume*100, 4) # [%] Volume fraction 

print("Total surface [mm^2]: ", area)
print("Total volume  [mm^2]: ", vol)
print("Volume fraction  (%): ", RD)