fusion-energy / cad_to_dagmc

Convert CAD geometry (STP files) or Cadquery assemblies to DAGMC h5m files
MIT License
20 stars 5 forks source link

use cadquery to facet #28

Open shimwell opened 1 year ago

shimwell commented 1 year ago

instead of installing gmsh we could potential use cadquery to facet the geometry.

here is some code that doesn't quite produce watertight geometry but looks like it is close


this didn't produce non overlapping water tight parts when the geometry was in contact with other surfaces
def tessellate(parts, tolerance: float = 0.1, angularTolerance: float = 0.1):
    """Creates a mesh / faceting / tessellation of the surface"""

    parts.mesh(tolerance, angularTolerance)

    offset = 0

    vertices: List[Vector] = []
    triangles = {}

    for f in parts.Faces():

        loc = TopLoc_Location()
        poly = BRep_Tool.Triangulation_s(f.wrapped, loc)
        Trsf = loc.Transformation()

        reverse = (
            True
            if f.wrapped.Orientation() == TopAbs_Orientation.TopAbs_REVERSED
            else False
        )

        # add vertices
        face_verticles = [
            (v.X(), v.Y(), v.Z()) for v in (v.Transformed(Trsf) for v in poly.Nodes())
        ]
        vertices += face_verticles

        face_triangles = [
            (
                t.Value(1) + offset - 1,
                t.Value(3) + offset - 1,
                t.Value(2) + offset - 1,
            )
            if reverse
            else (
                t.Value(1) + offset - 1,
                t.Value(2) + offset - 1,
                t.Value(3) + offset - 1,
            )
            for t in poly.Triangles()
        ]
        triangles[f.hashCode()] = face_triangles

        offset += poly.NbNodes()

    list_of_triangles_per_solid = []
    for s in parts.Solids():
        triangles_on_solid = []
        for f in s.Faces():
            triangles_on_solid += triangles[f.hashCode()]
        list_of_triangles_per_solid.append(triangles_on_solid)
    for vert in vertices:
    for tri in list_of_triangles_per_solid:
    return vertices, list_of_triangles_per_solid
shimwell commented 9 months ago

this would need to work for an imprinted assembly instead of parts

import cadquery as cq
assembly = cq.Assembly(name="TwoTouchingCuboids")
cuboid1 = cq.Workplane().box(10,10,10)
assembly.add(cuboid1)
cuboid2 = cq.Workplane().moveTo(5,5).box(3,3,3)
assembly.add(cuboid2)
imprinted_assembly, _ = cq.occ_impl.assembly.imprint(assembly)