Support for COMSOL (.STL) compatible mesh output file

[maxnyf]

I'm wondering if anyone has ideas for how to create a mesh file that is compatible with COMSOL, which accepts .ply, .stl and some others. COMSOL only accepts triangular mesh's for .ply, so the polymesh output to .ply does not work. Trying to convert .vtk from the tri-mesh to .ply or .stl with meshio does not work either.

Ideally I want a polymesh converted to an STL file, not sure if this has been attempted.

[kip-hart]

Hi Max, thanks for reaching out. To make an STL for a polymesh, you could split each 3D polygon facet into triangles. If a facet has 5 vertices, you could form triangles from vertices (0, 1, 2), vertices (0, 2, 3), and vertices (0, 3, 4). The code to write the STL would look something like:

with open('mySTL.stl', 'w') as file:
    file.write('solid mySTL\n')
    for vertices in pmesh.facets:
        pt0 = pmesh.points[vertices[0]]
        for idx in range(2, len(vertices)):
            pt1 = pmesh.points[vertices[idx-1]]
            pt2 = pmesh.points[vertices[idx]]
            write_tri_in_STL_format(file, pt0, pt1, pt2)
    file.write('endsolid mySTL\n')

The write_tri_in_STL_format function would write the triangle in STL format

There's a risk there of generating thin triangles with this method. Instead, you may want to create a center point for each facet and define triangles that emanate from it. This would look something like:

with open('mySTL.stl', 'w') as file:
    file.write('solid mySTL\n')
    for vertices in pmesh.facets:
        pts = pmesh.points[vertices]
        center = np.mean(pts, axis=2)
        n = len(vertices)
        for idx_pt1 in range(n):
            idx_pt2 = idx_pt1 + 1
            if idx_pt2 == n:
                idx_pt2 = 0
            pt1 = pts[idx_pt1]
            pt2 = pts[idx_pt2]
            write_tri_in_STL_format(file, center, pt1, pt2)
    file.write('endsolid mySTL\n')

You could also include an if ... else that catches cases where the original polygon is a triangle. This second option will write 2 additional triangles to the file for each facet, but it will prevent numerical issues caused by thin triangles.

If you'd like to contribute to the project and add an STL option to the polymesh.write function, I would be happy to review and merge it in.

[maxnyf]

Thanks for the quick response, I will look into it further!

[maxnyf]

This worked well! Both methods worked and imported into COMSOL! I just attached a rough script to demonstrate it in case anyone else ever needs. Assigning normal direction does not seem to be necessary for COMSOL. Not sure if more complicated geometry will ever cause issues, but this basic example works well. The below image is from COMSOL.

import matplotlib.pyplot as plt
import microstructpy as msp
import scipy.stats
import numpy as np

def write_stl_file(filename: str,pmesh: msp.meshing.polymesh.PolyMesh):
    solid_name = filename.split('.')[0]
    with open(filename, 'w') as file:
        file.write(f'solid {solid_name}\n')
        for vertices in pmesh.facets:
            pts = [pmesh.points[i] for i in vertices]
            center = np.mean(pts, axis=0)
            n = len(vertices)
            if n == 3:
                write_tri_in_STL_format(file, pts[0], pts[1], pts[2])
            else:
                for idx_pt1 in range(n):
                    idx_pt2 = idx_pt1 + 1
                    if idx_pt2 == n:
                        idx_pt2 = 0
                    pt1 = pts[idx_pt1]
                    pt2 = pts[idx_pt2]
                    write_tri_in_STL_format(file, center, pt1, pt2)
        file.write(f'endsolid {solid_name}\n')

def write_stl_file2(filename: str, pmesh: msp.meshing.polymesh.PolyMesh):
    solid_name = filename.split('.')[0]
    with open(filename, 'w') as file:
        file.write(f'solid {solid_name}\n')
        for vertices in pmesh.facets:
            pt0 = pmesh.points[vertices[0]]
            for idx in range(2, len(vertices)):
                pt1 = pmesh.points[vertices[idx - 1]]
                pt2 = pmesh.points[vertices[idx]]
                write_tri_in_STL_format(file, pt0, pt1, pt2)
        file.write(f'endsolid {solid_name}\n')

def write_tri_in_STL_format(file, pt1, pt2, pt3):
    ni = 0
    nj = 0
    nk = 0

    file.write("facet normal %.1f %.1f %.1f\n" % (ni,nj,nk))
    file.write(" outer loop\n")
    file.write("  vertex %.5f %.5f %.5f\n" % tuple(pt1))
    file.write("  vertex %.5f %.5f %.5f\n" % tuple(pt2))
    file.write("  vertex %.5f %.5f %.5f\n" % tuple(pt3))
    file.write(" endloop\n")
    file.write("endfacet\n")

material_1 = {
    'name': 'Matrix',
    'material_type': 'solid',
    'fraction': 2,
    'shape': 'sphere',
    'size': scipy.stats.uniform(loc=0, scale=1.5)
}

materials = [material_1]
domain = msp.geometry.Cube(side_length=5, center=(0, 0,0))

seed_area = domain.volume
rng_seeds = {'size': np.random.randint(100)}
seeds = msp.seeding.SeedList.from_info(materials,
                                       seed_area,
                                       rng_seeds)

seeds.position(domain)
pmesh = msp.meshing.PolyMesh.from_seeds(seeds, domain)

write_stl_file("test_stl1.stl",pmesh)
write_stl_file2("test_stl2.stl",pmesh)