query based meshing on top of GMSH
from meshql.mesh import *
from meshql.transactions import *
from .ql import *
from .preprocessing import *
meshql/__init__.py
OpenOrion/meshql
from enum import Enum
import cadquery as cq
from cadquery.cq import CQObject
from dataclasses import dataclass
from typing import Iterable, Optional, OrderedDict, Sequence, Union, cast
from meshql.utils.cq import CQLinq, CQType, CQ_TYPE_STR_MAPPING
from meshql.utils.types import OrderedSet


ENTITY_DIM_MAPPING: dict[CQType, int] = {
    "vertex": 0,
    "edge": 1,
    "face": 2,
    "solid": 3,
}

@dataclass
class Entity:
    type: CQType
    "dimension type of the entity."
    
    tag: int = -1
    "tag of the entity."

    name: Optional[str] = None
    "name of the entity."

    def __eq__(self, __value: object) -> bool:
        if isinstance(__value, Entity):
            return self.type == __value.type and self.tag == __value.tag
        return False

    @property
    def dim(self):
        if self.type not in ENTITY_DIM_MAPPING:
            raise ValueError(f"Entity type {self.type} not supported, only {ENTITY_DIM_MAPPING.keys()}")

        return ENTITY_DIM_MAPPING[self.type]

    def __hash__(self) -> int:
        return hash((self.type, self.tag))


class CQEntityContext:
    "Maps OCC objects to gmsh entity tags"
    def __init__(self, workplane: cq.Workplane) -> None:
        self.dimension = 3 if len(workplane.solids().vals()) else 2

        self.entity_registries: dict[CQType, OrderedDict[CQObject, Entity]] = {
            "solid": OrderedDict[CQObject, Entity](),
            "shell": OrderedDict[CQObject, Entity](),
            "face":  OrderedDict[CQObject, Entity](),
            "wire": OrderedDict[CQObject, Entity](),
            "edge": OrderedDict[CQObject, Entity](),
            "vertex": OrderedDict[CQObject, Entity](),
        }

        if self.dimension == 3:
            self._init_3d_objs(workplane)
        else:
            self._init_2d_objs(workplane)

    def add(self, obj: CQObject):
        entity_type = CQ_TYPE_STR_MAPPING[type(obj)]
        registry = self.entity_registries[entity_type]
        if obj not in registry:
            tag = len(registry) + 1
            registry[obj] = Entity(entity_type, tag)


    def select(self, obj: CQObject):
        entity_type = CQ_TYPE_STR_MAPPING[type(obj)]
        registry = self.entity_registries[entity_type]
        return registry[obj]
    
    def select_many(self, target: Union[cq.Workplane, Iterable[CQObject]], type: Optional[CQType] = None):
        entities = OrderedSet[Entity]()
        objs = target.vals() if isinstance(target, cq.Workplane) else target
        selected_objs = objs if type is None else CQLinq.select(objs, type)
        for obj in selected_objs:
            try:
                selected_entity = self.select(obj)
                entities.add(selected_entity)
            except:
                ...

        return entities
    
    def select_batch(self, target: Union[cq.Workplane, Iterable[CQObject]], parent_type: CQType, child_type: CQType):
        objs = target.vals() if isinstance(target, cq.Workplane) else target
        selected_batches = CQLinq.select_batch(objs, parent_type, child_type)
        for selected_batch in selected_batches:
            yield self.select_many(selected_batch)

    def _init_3d_objs(self, target: Union[cq.Workplane, Sequence[CQObject]]):
        objs = CQLinq.select(target, "solid")
        for compound in cast(Sequence[cq.Solid], objs):
            for solid in compound.Solids():
                for shell in solid.Shells():
                    self._init_2d_objs(shell.Faces())
                    self.add(shell)
                self.add(solid)

    def _init_2d_objs(self, target: Union[cq.Workplane, Sequence[CQObject]]):
        objs = CQLinq.select(target, "face")
        for face in cast(Sequence[cq.Face], objs):
            for wire in [face.outerWire(), *face.innerWires()]:
                for edge in wire.Edges():
                    for vertex in edge.Vertices():
                        self.add(vertex)
                    self.add(edge)
                self.add(wire)
            self.add(face)





meshql/entity.py
OpenOrion/meshql
from . importers import *
from .exporters import *
from .mesh import *
meshql/mesh/__init__.py
OpenOrion/meshql
from typing import Sequence, Union
from .mesh import Mesh


def export_to_su2(meshes: Union[Mesh, Sequence[Mesh]], file_path: str):
    """Export a mesh from SU2 format"""
    from su2fmt import Mesh as Su2Mesh, export_mesh
    from su2fmt.mesh import ElementType as Su2ElementType, Zone
    if not isinstance(meshes, Sequence):
        meshes = [meshes]
    zones = []
    for izone, mesh in enumerate(meshes):
        element_types: Sequence[Su2ElementType] = []
        for element_type in mesh.element_types:
            try:
                su2_element_type = Su2ElementType[element_type.name]
                element_types.append(su2_element_type)
            except:
                raise ValueError("Warning: Element type not supported: ", element_type.name)
        
        marker_types = {}
        for marker_label, marker_element_types in mesh.marker_types.items():
            marker_types[marker_label] = [Su2ElementType[marker_element_type.name] for marker_element_type in marker_element_types]
        
        zone = Zone(
            izone=izone+1,
            ndime=mesh.dim,
            elements=mesh.elements,
            element_types=element_types,
            points=mesh.points,
            markers=mesh.markers,
            marker_types=marker_types
        )
        zones.append(zone)
    su2_mesh = Su2Mesh(len(zones), zones)
    return export_mesh(su2_mesh, file_path)
meshql/mesh/exporters.py
OpenOrion/meshql

from meshql.entity import ENTITY_DIM_MAPPING
from meshql.utils.types import NumpyFloat
from .mesh import ElementType, Mesh
import numpy.typing as npt
import numpy as np
import gmsh

def import_from_su2(file_path: str):
    """Import a mesh from SU2 format"""
    from su2fmt import parse_mesh
    su2_mesh = parse_mesh(file_path)
    meshes = []
    for zone in su2_mesh.zones:
        element_types = []
        for su2_element_type in zone.element_types:
            element_type = ElementType[su2_element_type.name]
            element_types.append(element_type)

        marker_types = {}
        for marker_label, su2_marker_element_types in zone.marker_types.items():
            marker_types[marker_label] = [ElementType[su2_marker_element_type.name] for su2_marker_element_type in su2_marker_element_types]
        
        mesh = Mesh(
            zone.ndime,
            zone.elements,
            element_types,
            zone.points,
            zone.markers,
            marker_types
        )
        meshes.append(mesh)
    if len(meshes) == 1:
        return meshes[0]
    return meshes

def import_from_msh(file_path: str):
    """Import a mesh from Gmsh format"""
    import gmsh
    gmsh.initialize()
    gmsh.open(file_path)
    mesh = import_from_gmsh()
    gmsh.finalize()
    return mesh

def import_from_file(file_path: str):
    """Import a mesh from a file"""
    if file_path.endswith('.su2'):
        return import_from_su2(file_path)
    elif file_path.endswith('.msh'):
        return import_from_msh(file_path)
    raise ValueError(f"File extension not supported: {file_path}")

def import_from_gmsh() -> Mesh:
    dim = gmsh.model.getDimension()
    elements: list[npt.NDArray[np.uint16]] = []
    element_types: list[ElementType] = []

    node_tags, points_concatted, _ = gmsh.model.mesh.getNodes()
    node_indices = np.argsort(node_tags-1)  # type: ignore
    points = np.array(points_concatted, dtype=NumpyFloat).reshape((-1, 3))[node_indices]

    grouped_concatted_elements = gmsh.model.mesh.getElements()
    for element_type_value, grouped_element_tags, grouped_node_tags_concatted in zip(*grouped_concatted_elements):
        if element_type_value == ElementType.POINT.value or element_type_value == ElementType.LINE.value:
            continue
        num_nodes = gmsh.model.mesh.getElementProperties(element_type_value)[3]
        group_elements = np.array(grouped_node_tags_concatted, dtype=np.uint16).reshape((-1, num_nodes)) - 1
        elements += list(group_elements)
        element_types += [ElementType(element_type_value)] * len(group_elements)

    # get physical groups
    markers = {}
    marker_types = {}
    physical_groups = gmsh.model.getPhysicalGroups()
    for group_dim, group_tag in physical_groups:
        marker_name = gmsh.model.getPhysicalName(group_dim, group_tag)
        if len(marker_name) == 0 or group_dim == ENTITY_DIM_MAPPING["solid"]:
            continue
        entities = gmsh.model.getEntitiesForPhysicalGroup(group_dim, group_tag)
        for entity in entities:
            marker_grouped_concatted_elements = gmsh.model.mesh.getElements(group_dim, tag=entity)
            assert len(marker_grouped_concatted_elements[0]) == 1, "There should only be one group"
            marker_element_type, marker_node_tags_concatted = marker_grouped_concatted_elements[0][0], marker_grouped_concatted_elements[2][0]
            if marker_element_type == ElementType.LINE.value:
                node_count = 2
            elif marker_element_type == ElementType.TRIANGLE.value:
                node_count = 3
            elif marker_element_type == ElementType.QUADRILATERAL.value:
                node_count = 4
            else:
                raise ValueError("Only expecting plane element types")

            marker_elements = np.array(marker_node_tags_concatted, dtype=np.uint16).reshape((-1, node_count)) - 1
            if marker_name not in markers:
                markers[marker_name] = []
                marker_types[marker_name] = []
            
            markers[marker_name] += list(marker_elements)
            marker_types[marker_name] += [ElementType(marker_element_type)] * len(marker_elements)


    return Mesh(
        dim,
        elements,
        element_types,
        points,
        markers,
        marker_types
    )

meshql/mesh/importers.py
OpenOrion/meshql

from dataclasses import dataclass, field
from enum import Enum
from typing import Sequence
import numpy.typing as npt
import numpy as np
from meshql.utils.types import NumpyFloat

class ElementType(Enum):
    LINE = 1
    TRIANGLE = 2
    QUADRILATERAL = 3
    TETRAHEDRON = 4
    HEXAHEDRON = 5
    PRISM = 6
    PYRAMID = 7
    POINT = 15


@dataclass
class Mesh:
    dim: int
    elements: Sequence[npt.NDArray[np.uint16]]
    element_types: Sequence[ElementType]
    points: npt.NDArray[NumpyFloat]
    markers: dict[str, Sequence[npt.NDArray[np.uint16]]]
    marker_types: dict[str, Sequence[ElementType]]
    target_points: dict[str, dict[np.uint16, str]] = field(default_factory=dict)

meshql/mesh/mesh.py
OpenOrion/meshql
import numpy as np
import cadquery as cq
from cadquery.cq import VectorLike
from typing import Literal, Optional, Sequence, Union
from meshql.utils.cq import CQCache, CQExtensions, CQLinq
from meshql.utils.types import LineTuple, VectorTuple
from jupyter_cadquery import show
Axis = Union[Literal["X", "Y", "Z"], VectorTuple, cq.Vector]
def get_normal_from_axis(axis: Axis):
    if isinstance(axis, str):
        return cq.Vector([1 if axis == "X" else 0, 1 if axis == "Y" else 0, 1 if axis == "Z" else 0])        
    elif isinstance(axis, tuple):
        return cq.Vector(axis)
    return axis

def norm_line_tuple(line: LineTuple):
    pnt1 = tuple(float(v) for v in ((*line[0], 0) if len(line[0]) == 2 else line[0]))
    pnt2 = tuple(float(v) for v in ((*line[1], 0) if len(line[1]) == 2 else line[1]))
    return (pnt1, pnt2)

class Split:
    @staticmethod
    def from_plane(
        base_pnt: VectorLike = (0,0,0), 
        angle: VectorTuple = (0,0,0), 
    ):
        return cq.Face.makePlane(None, None, base_pnt, tuple(np.radians(angle)))

    @staticmethod
    def from_faces(
        workplane: cq.Workplane, 
        face_type: Literal['interior', 'exterior'],
        snap_tolerance: Optional[float] = None,
        angle_offset: VectorTuple = (0,0,0),
    ):
        offset = cq.Vector(tuple(np.radians(angle_offset)))
        type_groups = CQLinq.groupByTypes(workplane, only_faces=True, check_splits=False)
        if False:
            yield None
        dir = "away" if face_type == "interior" else "towards"
        face_edge_groups: dict[cq.Edge, set[cq.Face]] = {}

        for face in type_groups[face_type]:
            assert isinstance(face, cq.Face)
            edges = CQLinq.select(face, "edge")
            for edge in edges:
                assert isinstance(edge, cq.Edge)
                if edge not in face_edge_groups:
                    face_edge_groups[edge] = set()
                face_edge_groups[edge].add(face)

        for edge, faces in face_edge_groups.items():
            average_normal = np.average([face.normalAt().toTuple() for face in faces], axis=0)
            average_normal_vec = cq.Vector(tuple(average_normal)) + offset
            edge_vec = cq.Vector(edge.endPoint().toTuple()) - cq.Vector(edge.startPoint().toTuple())
            is_parallel = edge_vec.dot(average_normal_vec) == 0
            if is_parallel:
                yield Split.from_edge(workplane, edge, average_normal_vec, dir, snap_tolerance)


    @staticmethod
    def from_anchor(
        workplane: cq.Workplane, 
        anchor: Union[list[VectorTuple], VectorTuple] = (0,0,0), 
        angle: Union[list[VectorTuple], VectorTuple] = (0,0,0),
        snap_tolerance: Optional[float] = None,
        until: Literal["next", "all"] = "next",
    ):
        anchors = [anchor] if isinstance(anchor, tuple) else anchor
        angles = [angle] if isinstance(angle, tuple) else angle
        assert len(anchors) == len(angles), "anchors and dirs must be the same length"

        edges = []
        for anchor, angle in zip(anchors, angles):
            split_face = Split.from_plane(anchor, angle)
            if until == "next":
                intersect_vertex = CQExtensions.split_intersect(workplane, anchor, split_face, snap_tolerance)
                edges.append((anchor, intersect_vertex.toTuple())) # type: ignore
            else:
                edges.append(split_face)
        return Split.from_lines(workplane, edges)
    
    @staticmethod
    def from_pnts(pnts: Sequence[VectorLike]):
        return cq.Face.makeFromWires(cq.Wire.makePolygon(pnts))

    @staticmethod
    def from_edge(
        workplane: cq.Workplane,
        edge: cq.Edge, 
        axis: Union[Literal["X", "Y", "Z"], VectorTuple, cq.Vector] = "Z",
        dir: Literal["away", "towards", "both"] = "both",
        snap_tolerance: Optional[float] = None,
    ):
        scaled_edge = CQExtensions.scale(edge, z=10)
        maxDim = workplane.findSolid().BoundingBox().DiagonalLength * 10.0
        normal_vector = get_normal_from_axis(axis)        

        max_dim_edge = scaled_edge.translate(normal_vector * maxDim) if dir in ("both", "towards") else scaled_edge
        min_dim_edge = scaled_edge.translate(-normal_vector * maxDim) if dir in ("both", "away") else scaled_edge

        
        split_face = cq.Face.makeFromWires(cq.Wire.assembleEdges([
            min_dim_edge,
            cq.Edge.makeLine(min_dim_edge.endPoint(), max_dim_edge.endPoint()),
            max_dim_edge,
            cq.Edge.makeLine(max_dim_edge.startPoint(), min_dim_edge.startPoint()),
        ]))

        # TODO: for some reason this snap tolerance is not working
        if snap_tolerance:
            intersect_vertex = CQExtensions.split_intersect(workplane, edge.startPoint(), split_face, snap_tolerance)
            assert intersect_vertex, "No intersection found"
            intersect_vec = cq.Vector(intersect_vertex.toTuple()) - cq.Vector(edge.startPoint().toTuple())
            intersect_vec_norm = intersect_vec/intersect_vec.Length
            return Split.from_edge(workplane, edge, intersect_vec_norm, "towards")
        return split_face

    @staticmethod
    def from_lines(
        workplane: cq.Workplane, 
        lines: Union[list[LineTuple], LineTuple], 
        axis: Union[Literal["X", "Y", "Z"], VectorTuple] = "Z",
        dir: Literal["away", "towards", "both"] = "both",
    ):
        if isinstance(lines, tuple):
            edges_pnts = np.array([norm_line_tuple(lines), norm_line_tuple(lines)])
        elif isinstance(lines, list) and len(lines) == 1:
            edges_pnts = np.array([norm_line_tuple(lines[0]), norm_line_tuple(lines[0])])

        else: 
            edges_pnts = [norm_line_tuple(line) for line in lines]
        maxDim = workplane.findSolid().BoundingBox().DiagonalLength * 10.0
        normal_vector = np.array(get_normal_from_axis(axis).toTuple())

        if dir in ("both", "towards"):
            edges_pnts[0] += maxDim * normal_vector
        if dir in ("both", "away"):        
            edges_pnts[-1] -= maxDim * normal_vector

        side1 = edges_pnts[:, 0].tolist()
        side2 = edges_pnts[:, 1].tolist()
        wire_pnts = [side1[0], *side2, *side1[1:][::-1]] 
        return Split.from_pnts(wire_pnts)


def split_workplane(workplane: cq.Workplane, splits: Sequence[cq.Face], use_cache: bool = True):
    shape_combo = [*workplane.vals(), *splits]
    cache_exists = CQCache.get_cache_exists(shape_combo) if use_cache else False
    cache_file_name = CQCache.get_file_name(shape_combo) if use_cache else ""
    if use_cache and cache_exists:
        shape = CQCache.import_brep(cache_file_name)
    else:
        for split in splits:      
            workplane = workplane.split(split)
        shape = CQExtensions.fuse_shapes(workplane.vals())
        if use_cache:
            CQCache.export_brep(shape, cache_file_name)
    return cq.Workplane(shape)

meshql/preprocessing/split.py
OpenOrion/meshql
import gmsh
import numpy as np
import cadquery as cq
from cadquery.cq import CQObject
from typing import Callable, Iterable, Literal, Optional, Sequence, Union, cast
from cadquery.selectors import Selector
from meshql.entity import CQEntityContext, Entity
from meshql.preprocessing.split import split_workplane
from meshql.transaction import Transaction, TransactionContext
from meshql.transactions.algorithm import MeshAlgorithm2DType, MeshAlgorithm3DType, MeshSubdivisionType, SetMeshAlgorithm2D, SetMeshAlgorithm3D, SetSubdivisionAlgorithm
from meshql.transactions.boundary_layer import UnstructuredBoundaryLayer, UnstructuredBoundaryLayer2D, get_boundary_ratio
from meshql.transactions.physical_group import SetPhysicalGroup
from meshql.transactions.refinement import Recombine, Refine, SetMeshSize, SetSmoothing
from meshql.transactions.transfinite import SetTransfiniteEdge, SetTransfiniteFace, SetTransfiniteSolid, TransfiniteArrangementType, TransfiniteMeshType
from meshql.mesh.exporters import export_to_su2
from meshql.utils.cq import CQ_TYPE_RANKING, CQ_TYPE_STR_MAPPING, CQExtensions, CQGroupTypeString, CQLinq, CQType
from meshql.utils.types import OrderedSet
from meshql.visualizer import visualize_mesh
from jupyter_cadquery import show

class GeometryQL:
    _workplane: cq.Workplane
    _initial_workplane: cq.Workplane
    def __init__(self) -> None:
        self._initial_workplane = self._workplane = None # type: ignore
        self._ctx = TransactionContext()
        self.is_structured = False
        self._transfinite_edge_groups = list[set[cq.Edge]]()
    def __enter__(self):
        gmsh.initialize()
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        gmsh.finalize()

    def end(self, num: Optional[int] = None):
        if num is None:
            self._workplane = self._initial_workplane
        else:
            self._workplane = self._workplane.end(num)
        return self

    def load(self, target: Union[cq.Workplane, str, Iterable[CQObject]], splits: Optional[Callable[[cq.Workplane], Sequence[cq.Face]]] = None, use_cache: bool = False):
        assert self._workplane is None, "Workplane is already loaded."

        workplane = self._initial_workplane = CQExtensions.import_workplane(target)
        is_2d = CQExtensions.get_dimension(workplane) == 2

        # extrudes 2D shapes to 3D
        if is_2d:
            workplane = workplane.extrude(-1)

        if splits:
            workplane = split_workplane(workplane, splits(workplane), use_cache)


        if is_2d:
            # fuses top faces to appear as one Compound in GMSH
            faces = cast(Sequence[cq.Face], workplane.faces(">Z").vals())
            fused_face = CQExtensions.fuse_shapes(faces)
            self._initial_workplane = cq.Workplane(fused_face)
        else:
            self._initial_workplane = workplane

        self._type_groups = CQLinq.groupByTypes(self._initial_workplane, exclude_split=is_2d)


        self._workplane = self._initial_workplane
        topods = self._workplane.toOCC()
        gmsh.model.occ.importShapesNativePointer(topods._address())
        gmsh.model.occ.synchronize()

        self._entity_ctx = CQEntityContext(self._workplane)

        self._tag_workplane()

        return self    
    

    def _tag_workplane(self):
        "Tag all gmsh entity tags to workplane"
        for cq_type, registry in self._entity_ctx.entity_registries.items():
            for occ_obj in registry.keys():
                tag = f"{cq_type}/{registry[occ_obj].tag}"
                self._workplane.newObject([occ_obj]).tag(tag)

    def solids(self, selector: Union[Selector, str, None] = None, tag: Union[str, None] = None, type: Optional[CQGroupTypeString] = None, indices: Optional[Sequence[int]] = None):
        obj_type = type and self._type_groups[type]
        selector = CQExtensions.get_selector(selector, obj_type, indices)
        self._workplane = self._workplane.solids(selector, tag)
        return self

    def faces(self, selector: Union[Selector, str, None] = None, tag: Union[str, None] = None, type: Optional[CQGroupTypeString] = None, indices: Optional[Sequence[int]] = None):
        obj_type = type and self._type_groups[type]
        selector = CQExtensions.get_selector(selector, obj_type, indices)
        self._workplane = self._workplane.faces(selector, tag)
        return self
    
    def edges(self, selector: Union[Selector, str, None] = None, tag: Union[str, None] = None, type: Optional[CQGroupTypeString] = None, indices: Optional[Sequence[int]] = None):
        obj_type = type and self._type_groups[type]
        selector = CQExtensions.get_selector(selector, obj_type, indices)
        self._workplane = self._workplane.edges(selector, tag)
        return self

    def wires(self, selector: Union[Selector, str, None] = None, tag: Union[str, None] = None, type: Optional[CQGroupTypeString] = None, indices: Optional[Sequence[int]] = None):
        obj_type = type and self._type_groups[type]
        selector = CQExtensions.get_selector(selector, obj_type, indices)
        self._workplane = self._workplane.wires(selector, tag)
        return self

    def vertices(self, selector: Union[Selector, str, None] = None, tag: Union[str, None] = None, type: Optional[CQGroupTypeString] = None, indices: Optional[Sequence[int]] = None):
        obj_type = type and self._type_groups[type]
        selector = CQExtensions.get_selector(selector, obj_type, indices)
        self._workplane = self._workplane.vertices(selector, tag)
        return self

    def vals(self):
        return self._entity_ctx.select_many(self._workplane)

    def val(self):
        return self._entity_ctx.select(self._workplane.val())


    def tag(self, names: Union[str, Sequence[str]]):
        if isinstance(names, str):
            self._workplane.tag(names)
        else:
            for i, cq_obj in enumerate(self._workplane.vals()):
                self._workplane.newObject([cq_obj]).tag(names[i])
        return self

    def fromTagged(self, tags: Union[str, Iterable[str]], resolve_type: Optional[CQType] = None, invert: bool = True):        
        if isinstance(tags, str) and resolve_type is None:
            self._workplane = self._workplane._getTagged(tags)
        else:
            tagged_objs = list(CQLinq.select_tagged(self._workplane, tags, resolve_type))
            tagged_cq_type = CQ_TYPE_STR_MAPPING[type(tagged_objs[0])]
            workplane_objs = CQLinq.select(self._workplane, tagged_cq_type)
            filtered_objs = CQLinq.filter(workplane_objs, tagged_objs, invert)
            self._workplane = self._workplane.newObject(filtered_objs)
        return self

    def addPhysicalGroup(self, group: Union[str, Sequence[str]]):
        if  isinstance(group, str):
            set_physical_group = SetPhysicalGroup(self.vals(), group)
            self._ctx.add_transaction(set_physical_group)
        else:
            objs = list(self.vals())
            group_entities: dict[str, OrderedSet[Entity]] = {}

            for i, group_name in enumerate(group):
                new_group_entity = objs[i]
                if group_name not in group_entities:
                    group_entities[group_name] = OrderedSet()
                group_entities[group_name].add(new_group_entity)
            
            for group_name, group_objs in group_entities.items():
                set_physical_group = SetPhysicalGroup(group_objs, group_name)
                self._ctx.add_transaction(set_physical_group)

        return self

    def recombine(self, angle: float = 45):
        faces = self._entity_ctx.select_many(self._workplane, "face")
        recombines = [Recombine(face, angle) for face in faces]
        self._ctx.add_transactions(recombines)
        return self

    def setMeshSize(self, size: Union[float, Callable[[float,float,float], float]]):
        points = self._entity_ctx.select_many(self._workplane, "vertex")
        set_size = SetMeshSize(points, size)
        self._ctx.add_transaction(set_size)
        return self

    def setMeshAlgorithm(self, type: MeshAlgorithm2DType, per_face: bool = False):
        faces = self._entity_ctx.select_many(self._workplane, "face")
        if per_face:
            set_algorithms = [SetMeshAlgorithm2D(type, face) for face in faces]
            self._ctx.add_transactions(set_algorithms)
        else:
            set_algorithm = SetMeshAlgorithm2D(type)
            self._ctx.add_transaction(set_algorithm)
        
        return self

    def setMeshAlgorithm3D(self, type: MeshAlgorithm3DType):
        set_algorithm3D = SetMeshAlgorithm3D(type)
        self._ctx.add_transaction(set_algorithm3D)
        return self

    def setSubdivisionAlgorithm(self, type: MeshSubdivisionType):
        set_subdivision_algorithm = SetSubdivisionAlgorithm(type)
        self._ctx.add_transaction(set_subdivision_algorithm)
        return self

    def smooth(self, num_smooths = 1):
        faces = self._entity_ctx.select_many(self._workplane)
        set_smoothings = [SetSmoothing(face, num_smooths) for face in faces]
        self._ctx.add_transactions(set_smoothings)
        return self

    def refine(self, num_refines = 1):
        refine = Refine(num_refines)
        self._ctx.add_transaction(refine)
        return self

    def setTransfiniteEdge(self, num_nodes: Optional[Union[Sequence[int], int]] = None, mesh_type: Optional[Union[TransfiniteMeshType, Sequence[TransfiniteMeshType]]] = None, coef: Optional[Union[float, Sequence[float]]] = None):
        edge_batch = self._entity_ctx.select_batch(self._workplane, "face", "edge")
        for edges in edge_batch:
            for i, edge in enumerate(edges):
                transaction = cast(SetTransfiniteEdge, self._ctx.get_transaction(SetTransfiniteEdge, edge))
                if transaction is not None:
                    if num_nodes is not None:
                        transaction.num_elems = num_nodes if isinstance(num_nodes, int) else num_nodes[i]
                    if mesh_type is not None:
                        transaction.mesh_type = mesh_type if isinstance(mesh_type, str) else mesh_type[i]
                    if coef is not None:
                        transaction.coef = coef if isinstance(coef, (int, float)) else coef[i]
                else:
                    assert num_nodes is not None, "num_nodes must be specified"
                    mesh_type = mesh_type or "Progression"
                    coef = coef or 1.0
                    set_transfinite_edge = SetTransfiniteEdge(
                        edge, 
                        num_nodes if isinstance(num_nodes, int) else num_nodes[i], 
                        mesh_type if isinstance(mesh_type, str) else mesh_type[i], 
                        coef if isinstance(coef, (int, float)) else coef[i]
                    )
                    self._ctx.add_transaction(set_transfinite_edge)

        return self

    def setTransfiniteFace(self, arrangement: TransfiniteArrangementType = "Left"):
        self.is_structured = True    
        cq_face_batch = CQLinq.select_batch(self._workplane, "solid", "face")
        for i, cq_faces in enumerate(cq_face_batch):
            faces = self._entity_ctx.select_many(cq_faces)
            set_transfinite_faces = [SetTransfiniteFace(face, arrangement) for face in faces]
            self._ctx.add_transactions(set_transfinite_faces)

        return self

    def setTransfiniteSolid(self):
        self.is_structured = True    
        solids = self._entity_ctx.select_many(self._workplane, "solid")
        set_transfinite_solids = [SetTransfiniteSolid(solid) for solid in solids]
        self._ctx.add_transactions(set_transfinite_solids)
        return self

    def _getTransfiniteEdgeGroups(self, cq_faces: Sequence[cq.Face]):
        transfinite_edge_groups: list[set[cq.Edge]] = []
        for cq_face in cq_faces:
            sorted_edges = CQLinq.sort(cq_face.Edges())
            for i, path in enumerate(sorted_edges):
                cq_edge = path.edge
                parllel_edge_index = (i+2 if i+2 < len(sorted_edges) else (i+2) - len(sorted_edges))
                cq_parllel_edge = sorted_edges[parllel_edge_index].edge
                found_group: Optional[set] = None
                for i, group in enumerate(transfinite_edge_groups):
                    if not found_group:
                        if cq_edge in group:
                            group.add(cq_parllel_edge)
                            found_group = group
                        elif cq_parllel_edge in group:
                            group.add(cq_edge)
                            found_group = group
                    else: 
                        if cq_edge in group or cq_parllel_edge in group:
                            found_group.update(group)
                            transfinite_edge_groups.remove(group)

                if found_group is None:
                    transfinite_edge_groups.append(set([path.edge, cq_parllel_edge]))
        return transfinite_edge_groups

    def _setTransfiniteFaceAuto(
        self, 
        cq_faces: Sequence[cq.Face], 
        max_nodes: int, 
        min_nodes: int = 1,
        arrangement: TransfiniteArrangementType = "Left"
    ):
        self.is_structured = True    
        for cq_face in cq_faces:
            face = self._entity_ctx.select(cq_face)
            set_transfinite_face = SetTransfiniteFace(face, arrangement)
            self._ctx.add_transaction(set_transfinite_face)
        self._transfinite_edge_groups = self._getTransfiniteEdgeGroups(cq_faces)
        
        for transfinite_group in self._transfinite_edge_groups:
            total_length = sum([cq_edge.Length() for cq_edge in transfinite_group])
            group_max_num_nodes = 0
            for cq_edge in transfinite_group:
                edge_num_nodes = int(np.ceil((cq_edge.Length()/total_length)*max_nodes))
                if edge_num_nodes < min_nodes:
                    edge_num_nodes = min_nodes
                if edge_num_nodes > group_max_num_nodes:
                    group_max_num_nodes = edge_num_nodes
            
            assert group_max_num_nodes > 0, "group_max_num_nodes must be greater than 0, make num_nodes higher"
            group_edges = self._entity_ctx.select_many(transfinite_group)
            set_transfinite_edges = [SetTransfiniteEdge(edge, group_max_num_nodes) for edge in group_edges]
            self._ctx.add_transactions(set_transfinite_edges)


    def addTransaction(self, toTransaction: Callable[["GeometryQL"], Transaction]):
        self._ctx.add_transaction(toTransaction(self))
        return self

    def setTransfiniteAuto(
        self,
        max_nodes: int,
        min_nodes: int = 1,
        auto_recombine: bool = True
    ):
        self.is_structured = True
        if CQExtensions.get_dimension(self._workplane) == 2:
            cq_faces = cast(Sequence[cq.Face], list(CQLinq.select(self._workplane, "face")))
            self._setTransfiniteFaceAuto(cq_faces, max_nodes, min_nodes)

        else:
            for cq_solid in cast(Sequence[cq.Solid], CQLinq.select(self._workplane, "solid")):
                solid = self._entity_ctx.select(cq_solid)
                set_transfinite_solid = SetTransfiniteSolid(solid)
                self._ctx.add_transaction(set_transfinite_solid)
            cq_faces = cast(Sequence[cq.Face], list(CQLinq.select(self._workplane, "face")))
            self._setTransfiniteFaceAuto(cq_faces, max_nodes, min_nodes)

        # transfinite_auto = SetTransfiniteAuto()
        # self._ctx.add_transaction(transfinite_auto)


        if auto_recombine:
            self.recombine()

        return self

    def _addStructuredBoundaryLayer(
            self, 
            cq_objs: Sequence[CQObject], 
            size: Optional[float] = None,
            ratio: Optional[float] = None,
        ):
        assert self.is_structured, "Structured boundary layer can only be applied after setTransfiniteAuto"
        assert (size is None) != (ratio is None), "Either size or ratio must be specified, not both"

        boundary_vertices =  list(CQLinq.select(cq_objs, "vertex"))

        for (cq_edge, edge) in self._entity_ctx.entity_registries["edge"].items():
            transaction = cast(SetTransfiniteEdge, self._ctx.get_transaction(SetTransfiniteEdge, edge))
            assert edge.type == "edge", "StructuredBoundaryLayer only accepts edges"
            if size:
                edge_ratio = get_boundary_ratio(cq_edge.Length(), size, transaction.num_elems) # type: ignore
            elif ratio:
                edge_ratio = ratio
            else:
                raise ValueError("Either size or ratio must be specified, not both")
            cq_curr_edge_vertices =  cq_edge.Vertices() # type: ignore
            if cq_curr_edge_vertices[0] in boundary_vertices and cq_curr_edge_vertices[-1] not in boundary_vertices:
                transaction.coef = edge_ratio

            elif cq_curr_edge_vertices[-1] in boundary_vertices and cq_curr_edge_vertices[0] not in boundary_vertices:
                transaction.coef = -edge_ratio

    def addBoundaryLayer(self, size: float, ratio: Optional[float] = None, num_layers: Optional[int] = None, auto_recombine: bool = True):
        if self.is_structured:
            self._addStructuredBoundaryLayer(self._workplane.vals(), size, ratio)
        else:
            ratio = ratio or 1.0
            assert num_layers is not None and size is not None and ratio is not None, "num_layers, hwall_n and ratio must be specified for unstructured boundary layer"
            if CQ_TYPE_RANKING[type(self._workplane.val())] < CQ_TYPE_RANKING[cq.Face]:
                boundary_layer = UnstructuredBoundaryLayer2D(self.vals(), ratio, size, num_layers)
            else:
                boundary_layer = UnstructuredBoundaryLayer(self.vals(), ratio, size, num_layers)
                if auto_recombine:
                    self.recombine()
            self._ctx.add_transaction(boundary_layer)
        return self

    def generate(self, dim: int = 3):
        self._ctx.generate(dim)
        return self

    def write(self, filename: str, dim: int = 3):
        if self._ctx.mesh is None:
            self.generate(dim)
        assert self._ctx.mesh is not None
        if filename.endswith(".su2"):
            export_to_su2(self._ctx.mesh, filename)
        else:
            gmsh.write(filename)
        return self

    def showTransfiniteGroup(self, group_index: int):
        assert self.is_structured, "Structured boundary layer can only be applied after setTransfiniteAuto"
        assert group_index < len(self._transfinite_edge_groups), f"Group index {group_index} is out of range"
        group = self._transfinite_edge_groups[group_index]
        show(self._workplane.newObject(group), theme="dark")
        return self

    def show(self, type: Literal["gmsh", "mesh", "cq", "plot"] = "cq", only_markers: bool = False):
        if type == "gmsh":
            gmsh.fltk.run()
        elif type == "mesh":
            assert self._ctx.mesh is not None, "Mesh is not generated yet."
            visualize_mesh(self._ctx.mesh, only_markers=only_markers)
        elif type == "plot":
            CQExtensions.plot_cq(self._workplane, ctx=self._entity_ctx)
        elif type == "cq":
            show(self._workplane, theme="dark")
        else:
            raise NotImplementedError(f"Unknown show type {type}")
        return self


    def close(self):
        gmsh.finalize()

meshql/ql.py
OpenOrion/meshql
import uuid
import gmsh
from dataclasses import dataclass
from typing import Optional, OrderedDict, Sequence, Union
from meshql.entity import Entity
from meshql.mesh.importers import import_from_gmsh
from meshql.mesh.mesh import Mesh
from meshql.utils.types import OrderedSet

@dataclass
class Transaction:
    def __post_init__(self):
        self.is_commited: bool = False
        self.is_generated: bool = False
        self.id = uuid.uuid4()
    
    def before_gen(self):
        "completes transaction before mesh generation."
        ...

    def after_gen(self):
        "completes transaction after mesh generation."
        ...

    def __hash__(self) -> int:
        return self.id.__hash__()

    def __eq__(self, __value: object) -> bool:
        return self.id.__eq__(__value)
    


@dataclass(eq=False)
class SingleEntityTransaction(Transaction):
    entity: Entity
    "The entity that transaction will be applied towards"

@dataclass(eq=False)
class MultiEntityTransaction(Transaction):
    entities: OrderedSet[Entity]
    "The entities that transaction will be applied towards"


class TransactionContext:
    def __init__(self):
        self.entity_transactions = OrderedDict[tuple[type[Transaction], Entity], Transaction]()
        self.system_transactions = OrderedDict[type[Transaction], Transaction]()

        self.mesh: Optional[Mesh] = None
    
    def get_transaction(self, transaction_type: type[Transaction], entity: Optional[Entity] = None) -> Optional[Transaction]:
        if entity is None:
            return self.system_transactions.get(transaction_type)
        else:
            return self.entity_transactions.get((transaction_type, entity))


    def add_transaction(self, transaction: Transaction, ignore_duplicates: bool = False):
        """
        transaction: Transaction - transaction to be added
        ignore_duplicates: bool = False - if True, will ignore duplicate entity transactions
        """
        if isinstance(transaction, (SingleEntityTransaction, MultiEntityTransaction) ):
            entities = transaction.entities if isinstance(transaction, MultiEntityTransaction) else OrderedSet([transaction.entity])
            for entity in entities:
                transaction_id = (type(transaction), entity)
                if transaction_id not in self.entity_transactions:
                    self.entity_transactions[transaction_id] = transaction
                else:
                    if not ignore_duplicates:
                        old_transaction = self.entity_transactions[transaction_id]
                        if isinstance(old_transaction, MultiEntityTransaction):
                            old_transaction.entities.remove(entity)
                        self.entity_transactions[transaction_id] = transaction
        else:
            self.system_transactions[type(transaction)] = transaction

    def add_transactions(self, transactions: Sequence[Transaction], ignore_duplicates: bool = False):
        """
        transactions: Sequence[Transaction] - transactions to be added
        ignore_duplicates: bool = False - if True, will ignore duplicate entity transactions
        """
        for transaction in transactions:
            self.add_transaction(transaction, ignore_duplicates)

    def generate(self, dim: int = 3):
        gmsh.model.occ.synchronize()
        transactions = OrderedSet([
            *self.entity_transactions.values(),
            * self.system_transactions.values()
        ])
        for transaction in transactions:
            transaction.before_gen()

        gmsh.model.mesh.generate(dim)

        for transaction in transactions:
            transaction.after_gen()
        
        self.entity_transactions = OrderedDict()
        self.system_transactions = OrderedDict()
        self.mesh = import_from_gmsh()

meshql/transaction.py
OpenOrion/meshql
from .algorithm import *
from .boundary_layer import *
from .physical_group import *
from .refinement import *
from .transfinite import *
meshql/transactions/__init__.py
OpenOrion/meshql
import gmsh
from enum import Enum
from dataclasses import dataclass
from typing import Literal, Optional
from meshql.entity import Entity
from meshql.transaction import SingleEntityTransaction, Transaction


MESH_ALGORITHM_2D_MAPPING = {
    "MeshAdapt": 1,
    "Automatic": 2,
    "InitialMeshOnly": 3,
    "Delaunay": 5,
    "FrontalDelaunay": 6,
    "BAMG": 7,
    "FrontalDelaunayQuads": 8,
    "PackingOfParallelograms": 9,
    "QuasiStructuredQuad": 11,
}

MESH_ALGORITHM_3D_MAPPING = {
   "Delaunay": 1,
   "InitialMeshOnly": 3,
   "Frontal": 4,
   "MMG3D": 7,
   "RTree": 9,
   "HXT": 10,
}

MESH_SUBDIVISION_ALGORITHM_MAPPING = {
    "None": 0,
    "AllQuadrangles": 1,
    "AllHexahedra": 2,
    "Barycentric": 3,
}

MeshAlgorithm2DType =  Literal[
    "MeshAdapt",
    "Automatic",
    "InitialMeshOnly",
    "Delaunay",
    "FrontalDelaunay",
    "BAMG",
    "FrontalDelaunayQuads",
    "PackingOfParallelograms",
    "QuasiStructuredQuad"
]

MeshAlgorithm3DType =  Literal[
    "Delaunay",
    "InitialMeshOnly",
    "Frontal",
    "MMG3D", 
    "RTree", 
    "HXT"
]

MeshSubdivisionType = Literal[
    "None",
    "AllQuadrangles",
    "AllHexahedra",
    "Barycentric",
]

@dataclass(eq=False)
class SetMeshAlgorithm2D(Transaction):
    type: MeshAlgorithm2DType
    "algorithm to use"

    entity: Optional[Entity] = None
    "Entity to set algorithm for"


    def before_gen(self):
        if self.entity:
            assert self.entity.type == "face", "Can only set per face for edges"
            gmsh.model.mesh.setAlgorithm(self.entity.dim, self.entity.tag, MESH_ALGORITHM_2D_MAPPING[self.type])
        else:
            algo_val = MESH_ALGORITHM_2D_MAPPING[self.type]
            gmsh.option.setNumber("Mesh.Algorithm", algo_val)

@dataclass(eq=False)
class SetMeshAlgorithm3D(Transaction):
    type: MeshAlgorithm3DType
    "algorithm to use"

    def before_gen(self):
        algo_val = MESH_ALGORITHM_3D_MAPPING[self.type]
        gmsh.option.setNumber("Mesh.Algorithm3D", algo_val)

@dataclass(eq=False)
class SetSubdivisionAlgorithm(Transaction):
    type: MeshSubdivisionType
    "algorithm to use"

    def before_gen(self):
        algo_val = MESH_SUBDIVISION_ALGORITHM_MAPPING[self.type]
        gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", algo_val)


meshql/transactions/algorithm.py
OpenOrion/meshql
import gmsh
from dataclasses import dataclass
from typing import Optional
import numpy as np
from meshql.entity import ENTITY_DIM_MAPPING, Entity
from meshql.transaction import MultiEntityTransaction
from meshql.utils.types import OrderedSet
def get_boundary_sizes(ratio: float, size: float, num_layers: int):
    sizes = [size]
    for i in range(1, num_layers): 
        # geometric series for boundary layer heights
        sizes.append(sizes[-1] + sizes[0] * ratio**i)
    return sizes

def get_boundary_ratio(length: float, normal_height: float, num_layers: int):
    # sympy solve for an = a1*(r)**(n - 1)
    # solve([(normal_height*(ratio)**(num_layers - 1)) - length], [r], dict=True)
    return (length/normal_height)**(1/(num_layers - 1))
    
@dataclass(eq=False)
class UnstructuredBoundaryLayer(MultiEntityTransaction):
    "geometric series boundary layer for unstructured meshes"

    entities: OrderedSet[Entity]
    "faces to be added to the boundary layer"

    ratio: float
    "size Ratio Between Two Successive Layers"

    size: float
    "mesh size normal to the the wall"

    num_layers: int
    "number of layers"

    def before_gen(self):
        sizes = get_boundary_sizes(self.ratio, self.size, self.num_layers)
        
        dim_tags = []
        for face in self.entities:
            assert face.dim == ENTITY_DIM_MAPPING["face"], "boundary layer can only be applied to faces"
            dim_tags.append((face.dim, face.tag))

        extruded_bnd_layer = gmsh.model.geo.extrudeBoundaryLayer(dim_tags, [1] * self.num_layers, sizes, True)

        top = []
        for i in range(1, len(extruded_bnd_layer)):
            if extruded_bnd_layer[i][0] == 3:
                top.append(extruded_bnd_layer[i-1])
        gmsh.model.geo.synchronize()
        bnd = gmsh.model.getBoundary(top)
        gmsh.model.geo.addCurveLoop([c[1] for c in bnd])



@dataclass(eq=False)
class UnstructuredBoundaryLayer2D(MultiEntityTransaction):
    "boundary layer for unstructured 2D meshes"

    entities: OrderedSet[Entity]
    "edges to be added to the boundary layer"

    ratio: float
    "size Ratio Between Two Successive Layers"

    size: float
    "mesh size Normal to the wall"

    num_layers: int
    "number of layers"

    def before_gen(self):
        edge_tags = []
        for edge in self.entities:
            assert edge.dim == ENTITY_DIM_MAPPING["edge"], "boundary layer can only be applied to edges"
            edge_tags.append(edge.tag)

        self.tag = gmsh.model.mesh.field.add('BoundaryLayer')
        gmsh.model.mesh.field.setAsBoundaryLayer(self.tag)

        gmsh.model.mesh.field.setNumbers(self.tag, 'CurvesList', edge_tags)
        gmsh.model.mesh.field.setNumber(self.tag, "Quads", 1)
        gmsh.model.mesh.field.setNumber(self.tag, "Size", self.size)
        gmsh.model.mesh.field.setNumber(self.tag, "NbLayers", self.num_layers)
        gmsh.model.mesh.field.setNumber(self.tag, "BetaLaw", 1)
meshql/transactions/boundary_layer.py
OpenOrion/meshql
import gmsh
from dataclasses import dataclass
from meshql.entity import Entity
from meshql.transaction import MultiEntityTransaction
from meshql.utils.types import OrderedSet


@dataclass(eq=False)
class SetPhysicalGroup(MultiEntityTransaction):
    entities: OrderedSet[Entity]
    "The entities that will be added to the physical group."

    name: str
    "The name of the physical group."

    def __post_init__(self):
        super().__post_init__()
        self.entity_type = self.entities.first.dim

    def before_gen(self):
        entity_tags = []
        for entity in self.entities:
            assert entity.type == self.entities.first.type, "all entities must be of the same type"
            entity.name = self.name
            entity_tags.append(entity.tag)

        physical_group_tag = gmsh.model.addPhysicalGroup(self.entities.first.dim, entity_tags)
        gmsh.model.set_physical_name(self.entity_type, physical_group_tag, self.name)
meshql/transactions/physical_group.py
OpenOrion/meshql
import gmsh
from typing import Callable, Union
from dataclasses import dataclass
from meshql.entity import Entity
from meshql.transaction import Transaction, SingleEntityTransaction, MultiEntityTransaction
from meshql.utils.types import OrderedSet

@dataclass(eq=False)
class Recombine(SingleEntityTransaction):
    entity: Entity
    "Entity to recombine for"

    angle: float = 45
    "Angle to recombine with"

    def before_gen(self):
        gmsh.model.mesh.setRecombine(self.entity.dim, self.entity.tag, self.angle)

@dataclass(eq=False)
class SetSmoothing(SingleEntityTransaction):
    entity: Entity
    "Entity to smooth for"
    
    num_smooths: int = 1
    "Number of times to smooth the mesh"

    def after_gen(self):
        gmsh.model.mesh.setSmoothing(self.entity.dim, self.entity.tag, self.num_smooths)


@dataclass(eq=False)
class Refine(Transaction):
    num_refines: int = 1
    "Number of times to refine the mesh"
    
    def after_gen(self):
        for _ in range(self.num_refines):
            gmsh.model.mesh.refine()

@dataclass(eq=False)
class SetMeshSize(MultiEntityTransaction):
    entities: OrderedSet[Entity]
    "Entities to set mesh sizes for"

    size: Union[float, Callable[[float, float, float], float]]
    "Size to set points"

    def before_gen(self):
        point_tags = [(entity.dim, entity.tag) for entity in self.entities]
        if isinstance(self.size, float):
            gmsh.model.mesh.setSize(point_tags, self.size)
        else:
            gmsh.model.mesh.setSizeCallback(lambda dim, tag, x, y, z, lc: lc if tag in point_tags else self.size(x, y, z)) # type: ignore

meshql/transactions/refinement.py
OpenOrion/meshql
import gmsh
import numpy as np
from dataclasses import dataclass
from typing import Literal, Optional, Sequence
from meshql.entity import Entity, ENTITY_DIM_MAPPING
from meshql.transaction import SingleEntityTransaction, MultiEntityTransaction, Transaction
from meshql.utils.types import OrderedSet

TransfiniteArrangementType = Literal["Left", "Right", "AlternateLeft", "AlternateRight"]
TransfiniteMeshType = Literal["Progression", "Bump", "Beta"]

def get_num_nodes_for_ratios(total_num_nodes: int, ratios: Sequence[float]):
    assert np.round(np.sum(ratios), 5) == 1, "Ratios must sum to 1"
    assert total_num_nodes > len(ratios), f"Number of nodes must be greater than number of ratios {len(ratios)}"
    allocated_nodes = []
    for ratio in ratios:
        num_nodes = int(np.round(ratio * total_num_nodes)) or 1
        allocated_nodes.append(num_nodes)

    # remove nodes uniformly from highest ratios
    descending_ratio_indexes = sorted(range(len(ratios)), key=lambda i: -ratios[i])
    num_allocated_nodes = sum(allocated_nodes)
    if num_allocated_nodes > total_num_nodes:
        total_node_diff = num_allocated_nodes - total_num_nodes
        for i in descending_ratio_indexes:
            if allocated_nodes[i] > 1:
                allocated_node_diff = int(np.ceil(total_node_diff*ratios[i]))
                allocated_nodes[i] -= allocated_node_diff
                num_allocated_nodes -= allocated_node_diff
            if num_allocated_nodes == total_num_nodes:
                break
    assert sum(allocated_nodes) == total_num_nodes, (
        f"Number of allocated nodes must equal num_nodes, {num_allocated_nodes} != {total_num_nodes}"
    )

    return allocated_nodes


@dataclass(eq=False)
class SetTransfiniteEdge(SingleEntityTransaction):
    entity: Entity
    "edge to be added to the boundary layer"

    num_elems: int
    "number of elems for edge"

    mesh_type: TransfiniteMeshType = "Progression"
    "mesh type for edge"

    coef: float = 1.0
    "coefficients for edge"

    def __post_init__(self):
        super().__post_init__()

    def before_gen(self):
        assert self.entity.type == "edge", "SetTransfiniteEdge only accepts edges"
        num_nodes = self.num_elems + 1
        gmsh.model.mesh.setTransfiniteCurve(self.entity.tag, num_nodes, self.mesh_type, self.coef)



@dataclass(eq=False)
class SetTransfiniteFace(SingleEntityTransaction):
    entity: Entity
    "face to apply field"

    arrangement: TransfiniteArrangementType = "Left"
    "arrangement of transfinite face"

    corners: Optional[OrderedSet[Entity]] = None
    "corner point tags for transfinite face"

    def before_gen(self):
        assert self.entity.type == "face", "SetTransfiniteFace only accepts faces"
        corner_tags = [corner.tag for corner in self.corners] if self.corners else []
        gmsh.model.mesh.setTransfiniteSurface(self.entity.tag, self.arrangement, corner_tags)

@dataclass(eq=False)
class SetTransfiniteSolid(SingleEntityTransaction):
    entity: Entity
    "face to apply field"

    corners: Optional[OrderedSet[Entity]] = None
    "corner point tags for transfinite face"

    def before_gen(self):
        assert self.entity.type == "solid", "SetTransfiniteSolid only accepts solids"
        corner_tags = [corner.tag for corner in self.corners] if self.corners else []
        gmsh.model.mesh.setTransfiniteVolume(self.entity.tag, corner_tags)



@dataclass(eq=False)
class SetCompound(MultiEntityTransaction):
    entities: OrderedSet[Entity]
    "face to apply field"

    def before_gen(self):
        entity_tags = [entity.tag for entity in self.entities]
        gmsh.model.mesh.setCompound(ENTITY_DIM_MAPPING["edge"], entity_tags)


@dataclass(eq=False)
class SetTransfiniteAuto(Transaction):

    def before_gen(self):
        gmsh.option.setNumber('Mesh.MeshSizeMin', 0.5)
        gmsh.option.setNumber('Mesh.MeshSizeMax', 0.5)
        gmsh.model.mesh.setTransfiniteAutomatic()

meshql/transactions/transfinite.py
OpenOrion/meshql
from dataclasses import dataclass, field
import os
import tempfile
import hashlib
import base64
from plotly import graph_objects as go
import cadquery as cq
from cadquery.cq import CQObject, VectorLike
from typing import Callable, Iterable, Literal, Optional, Sequence, TypeVar, Union, cast
import numpy as np
from meshql.utils.plot import add_plot
from meshql.utils.shapes import get_sampling
from meshql.utils.types import OrderedSet, NumpyFloat
from OCP.BRepTools import BRepTools
from OCP.BRep import BRep_Builder
from OCP.TopoDS import TopoDS_Shape
from jupyter_cadquery import show

CQType = Literal["compound", "solid", "shell", "face", "wire", "edge", "vertex"]
CQGroupTypeString = Literal["split", "interior", "exterior"]
CQEdgeOrFace = Union[cq.Edge, cq.Face]
TShape = TypeVar('TShape', bound=cq.Shape)


TEMPDIR_PATH = tempfile.gettempdir()
CACHE_DIR_NAME = "meshql_geom_cache"
CACHE_DIR_PATH = os.path.join(TEMPDIR_PATH, CACHE_DIR_NAME)

@dataclass
class DirectedPath:
    edge: cq.Edge
    "edge or face of path"

    direction: int = 1
    "direction of path"

    def __post_init__(self):
        assert isinstance(self.edge, cq.Edge), "edge must be an edge"
        assert self.direction in [-1, 1], "direction must be -1 or 1"
        self.vertices = self.edge.Vertices()[::self.direction]
        self.start = self.vertices[0]
        self.end = self.vertices[-1]

    def __eq__(self, __value: object) -> bool:
        return self.edge == __value

    def __hash__(self) -> int:
        return self.edge.__hash__()

@dataclass
class Group:
    paths: list[DirectedPath] = field(default_factory=list)
    "elements in group"

    prev_group: Optional["Group"] = None
    "previous group"

    next_group: Optional["Group"] = None
    "next group"

    @property
    def start(self):
        return self.paths[0].start

    @property
    def end(self):
        return self.paths[-1].end

CQ_TYPE_STR_MAPPING: dict[type[CQObject], CQType] = {
    cq.Compound: "compound",
    cq.Solid: "solid",
    cq.Shell: "shell",
    cq.Face: "face",
    cq.Wire: "wire",
    cq.Edge: "edge",
    cq.Vertex: "vertex",
}

CQ_TYPE_RANKING = dict(zip(CQ_TYPE_STR_MAPPING.keys(), range(len(CQ_TYPE_STR_MAPPING)+1)[::-1]))

class CQLinq:    
    @staticmethod
    def select_tagged(workplane: cq.Workplane, tags: Union[str, Iterable[str]], type: Optional[CQType] = None):
        for tag in ([tags] if isinstance(tags, str) else tags):
            if type is None:
                yield from workplane._getTagged(tag).vals()
            else:
               yield from CQLinq.select(workplane._getTagged(tag).vals(), type)    

    @staticmethod
    def select(target: Union[cq.Workplane, Iterable[CQObject], CQObject], select_type: Optional[CQType] = None):
        cq_objs = target.vals() if isinstance(target, cq.Workplane) else (list(target) if isinstance(target, Iterable) else [target])
        
        if type is None or len(cq_objs) > 0 and CQ_TYPE_STR_MAPPING[type(cq_objs[0])] == select_type:
            yield from cq_objs
        
        for cq_obj in cq_objs:
            assert isinstance(cq_obj, cq.Shape), "target must be a shape"
            if select_type is None:
                yield cq_obj
            elif select_type == "compound":
                yield from cq_obj.Compounds()
            elif select_type == "solid":
                yield from cq_obj.Solids()
            elif select_type == "shell":
                yield from cq_obj.Shells()
            elif select_type == "face":
                yield from cq_obj.Faces()
            elif select_type == "wire":
                yield from cq_obj.Wires()
            elif select_type == "edge":
                yield from cq_obj.Edges()
            elif select_type == "vertex":
                yield from cq_obj.Vertices()

    @staticmethod
    def select_batch(target: Union[cq.Workplane, Iterable[CQObject]], parent_type: CQType, child_type: CQType):
        cq_objs = list(target.vals() if isinstance(target, cq.Workplane) else target)
        if CQ_TYPE_STR_MAPPING[type(cq_objs[0])] == child_type:
            yield cq_objs
        else:
            parent_cq_objs = CQLinq.select(cq_objs, parent_type)
            for parent_occ_obj in parent_cq_objs:
                yield CQLinq.select([parent_occ_obj], child_type)

    # TODO: take a look later if this function is even needed
    @staticmethod
    def filter(objs: Iterable[CQObject], filter_objs: Iterable[CQObject], invert: bool):
        filtered = []
        filter_objs = OrderedSet(filter_objs)
        for cq_obj in objs:
            if (invert and cq_obj in filter_objs) or (not invert and cq_obj not in filter_objs):
                filtered.append(cq_obj)

        # sort filtered in the same order as filtered_objs
        return [cq_obj for cq_obj in filter_objs if cq_obj in filtered]
        

    @staticmethod
    def sort(target: Union[cq.Edge, Sequence[cq.Edge]]):
        unsorted_cq_edges = [target] if isinstance(target, cq.Edge) else target
        cq_edges = list(unsorted_cq_edges[1:])
        sorted_paths = [DirectedPath(unsorted_cq_edges[0])]
        while cq_edges:
            for i, cq_edge in enumerate(cq_edges):
                vertices = cq_edge.Vertices()
                if vertices[0].toTuple() == sorted_paths[-1].end.toTuple():
                    sorted_paths.append(DirectedPath(cq_edge))
                    cq_edges.pop(i)
                    break
                elif vertices[-1].toTuple() == sorted_paths[-1].end.toTuple():
                    sorted_paths.append(DirectedPath(cq_edge, direction=-1))
                    cq_edges.pop(i)
                    break
                elif vertices[0].toTuple() == sorted_paths[0].start.toTuple():
                    sorted_paths.insert(0, DirectedPath(cq_edge, direction=-1))
                    cq_edges.pop(i)
                    break
    
            else:
                raise ValueError("Edges do not form a closed loop")
        
        assert sorted_paths[-1].end == sorted_paths[0].start, "Edges do not form a closed loop"
        return sorted_paths
    

    @staticmethod
    def groupByTypes(target: Union[cq.Workplane, Sequence[CQObject]], only_faces=False, check_splits: bool = True, exclude_split: bool = False): 
        workplane = target if isinstance(target, cq.Workplane) else cq.Workplane().add(target)

        add_wire_to_group = lambda wires, group: group.update([
            *wires,
            *CQLinq.select(wires, "edge"),
            *CQLinq.select(wires, "vertex"),
        ])
        
        groups: dict[CQGroupTypeString, OrderedSet[CQObject]] = {
            "split": OrderedSet[CQObject](),
            "interior": OrderedSet[CQObject](),
            "exterior": OrderedSet[CQObject](),
        }
        faces = list(CQLinq.select(target, "face"))
        is_2d = CQExtensions.get_dimension(workplane) == 2
        if is_2d:
            for face in faces:
                add_wire_to_group(face.innerWires(), groups["interior"])
                add_wire_to_group([face.outerWire()], groups["exterior"])
        else:
            for face in faces:
                assert isinstance(face, cq.Face), "object must be a face"
                split_intersect = CQExtensions.split_intersect(workplane, face.Center(), cq.Edge.makeLine(face.Center(), face.Center() + (face.normalAt()*1E-5))) if check_splits else False
                is_interior = CQExtensions.is_interior_face(face) 
                if split_intersect:
                    face_registry = groups["split"]
                else:
                    face_registry = groups["interior" if is_interior else "exterior"]
                face_registry.add(face)
                if not only_faces:
                    add_wire_to_group(face.Wires(), face_registry)

        if exclude_split:
            new_exterior = OrderedSet[CQObject]()
            for obj in groups["exterior"]:
                if obj not in groups["split"]:
                    new_exterior.add(obj)
            groups["exterior"] = new_exterior

            new_interior = OrderedSet[CQObject]()
            for obj in groups["interior"]:
                if obj not in groups["split"]:
                    new_interior.add(obj)
            groups["interior"] = new_interior

        return groups


    @staticmethod
    def groupBy(target: Union[cq.Workplane, Sequence[CQObject]], parent_type: CQType, child_type: CQType): 
        groups: dict[CQObject, set[CQObject]] = {}
        
        cq_objs = target.vals() if isinstance(target, cq.Workplane) else (target if isinstance(target, Iterable) else [target])
        for cq_obj in cq_objs:
            parents = CQLinq.select(cq_obj, parent_type)
            for parent in parents:
                children = CQLinq.select(parent, child_type)
                for child in children:
                    if child not in groups:
                        groups[child] = set[CQObject]()
                    groups[child].add(parent)
        return groups




    @staticmethod
    def find(target: Union[cq.Workplane, Iterable[CQObject]], where: Callable[[CQObject], bool]):
        cq_objs = target.vals() if isinstance(target, cq.Workplane) else (target if isinstance(target, Iterable) else [target])
        for cq_obj in cq_objs:
            if where(cq_obj):
                yield cq_obj

class CQExtensions:
    @staticmethod
    def is_interior_face(face: CQObject, invert: bool = False):
        assert isinstance(face, cq.Face), "object must be a face"
        face_normal = face.normalAt()
        face_centroid = face.Center()
        interior_dot_product = face_normal.dot(face_centroid)
        return (not invert and interior_dot_product < 0) or (invert and interior_dot_product >= 0)

    @staticmethod
    def find_nearest_point(workplane: cq.Workplane, near_point: cq.Vertex, tolerance: float = 1e-2) -> cq.Vertex:
        min_dist_vertex, min_dist = None, float("inf")
        for vertex in workplane.vertices().vals():
            dist = cast(cq.Vertex, vertex).distance(near_point)
            if dist < min_dist and dist <= tolerance:
                min_dist_vertex, min_dist = vertex, dist
        return cast(cq.Vertex, min_dist_vertex)

    @staticmethod
    def split_intersect(
        workplane: cq.Workplane, 
        anchor: VectorLike, 
        splitter: CQObject, 
        snap_tolerance: Optional[float] = None,
    ) -> Optional[cq.Vertex]:

        intersected_edges = workplane.intersect(cq.Workplane(splitter)).edges().vals()
        if len(intersected_edges) == 0:
            return None

        intersected_vertices = OrderedSet(CQLinq.select(intersected_edges, "vertex"))
        min_dist_vertex, min_dist = None, float("inf") 
        for vertex in intersected_vertices:
                intersect_dist = (cq.Vector(vertex.X, vertex.Y, vertex.Z) - cq.Vector(anchor)).Length
                if intersect_dist !=0 and intersect_dist < min_dist:
                    min_dist_vertex, min_dist = cast(cq.Vertex, vertex), intersect_dist
            
        if snap_tolerance and isinstance(min_dist_vertex, cq.Vertex):
            nearest_point = CQExtensions.find_nearest_point(workplane, min_dist_vertex, snap_tolerance)
            if nearest_point:
                return nearest_point
        return min_dist_vertex

    @staticmethod
    def get_angle_between(prev: CQEdgeOrFace, curr: CQEdgeOrFace):
        if isinstance(prev, cq.Edge) and isinstance(curr, cq.Edge):
            prev_tangent_vec = prev.tangentAt(0.5) # type: ignore
            tangent_vec = curr.tangentAt(0.5)      # type: ignore
        else:
            prev_tangent_vec = prev.normalAt() # type: ignore
            tangent_vec = curr.normalAt()      # type: ignore
        angle = prev_tangent_vec.getAngle(tangent_vec)
        assert not np.isnan(angle), "angle should not be NaN"
        return angle

    @staticmethod
    def fuse_shapes(shapes: Sequence[CQObject]) -> cq.Shape:
        fused_shape: Optional[cq.Shape] = None
        for shape in shapes:
            assert isinstance(shape, cq.Shape), "shape must be a shape"
            if fused_shape:
                fused_shape = fused_shape.fuse(shape)
            else:
                fused_shape = shape
        assert fused_shape is not None, "No shapes to fuse"
        return fused_shape

    @staticmethod
    def plot_cq(
        target: Union[cq.Workplane, CQObject, Sequence[CQObject], Sequence[Group], Sequence[Sequence[CQObject]]], 
        title: str = "Plot", 
        samples_per_spline: int = 50,
        ctx = None
    ):
        from meshql.entity import CQEntityContext
        ctx = cast(CQEntityContext, ctx)

        fig = go.Figure(
            layout=go.Layout(title=go.layout.Title(text=title))
        )
        if isinstance(target, cq.Workplane):
            edge_groups = [[edge] for edge in CQLinq.select(target, "edge")]
        elif isinstance(target, CQObject):
            edge_groups = [[edge] for edge in CQLinq.select(target, "edge")]
        elif isinstance(target, Sequence) and isinstance(target[0], CQObject):
            edge_groups = [cast(Sequence[CQObject], target)]
        elif isinstance(target, Sequence) and isinstance(target[0], Group):
            edge_groups = [[path.edge for path in cast(Group, group).paths]for group in target]
        else:
            target = cast(Sequence[Sequence], target)
            edge_groups = cast(Sequence[Sequence[CQObject]], target)

        for i, edges in enumerate(edge_groups):

            edge_name = f"Edge{ctx.select(edges[0]).tag}" if ctx else f"Edge{i}"
            sampling = get_sampling(0, 1, samples_per_spline, False)
            coords = np.concatenate([np.array([vec.toTuple() for vec in edge.positions(sampling)], dtype=NumpyFloat) for edge in edges]) # type: ignore
            add_plot(coords, fig, edge_name)

        fig.layout.yaxis.scaleanchor = "x"  # type: ignore
        fig.show()

    @staticmethod
    def get_dimension(workplane: cq.Workplane):
        return 2 if len(workplane.solids().vals()) == 0 else 3

    @staticmethod
    def import_workplane(target: Union[cq.Workplane, str, Iterable[CQObject]]):
        if isinstance(target, str):
            if target.lower().endswith(".step"):
                workplane = cq.importers.importStep(target)
            elif target.lower().endswith(".dxf"):
                workplane = cq.importers.importDXF(target)
            else:
                raise ValueError(f"Unsupported file type: {target}")
        elif isinstance(target, Iterable):
            workplane = cq.Workplane().newObject(target)
        else:
            workplane = target

        return workplane


    @staticmethod
    def get_selector(selector: Union[cq.Selector, str, None], group: Optional[OrderedSet[CQObject]], indices: Optional[Sequence[int]] = None):
        selectors = []
        if isinstance(selector, str):
            selector = selectors.append(cq.StringSyntaxSelector(selector))
        elif isinstance(selector, cq.Selector):
            selectors.append(selector)

        if group is not None:
            selectors.append(GroupSelector(group))
        if indices is not None:
            selectors.append(IndexSelector(indices))

        if len(selectors) > 0:
            prev_selector = selectors[0]
            for selector in selectors[1:]:
                prev_selector = cq.selectors.AndSelector(prev_selector, selector)
            return prev_selector

    @staticmethod
    def scale(shape: TShape, x: float = 1, y: float = 1, z: float = 1) -> TShape:
        t = cq.Matrix([
            [x, 0, 0, 0],
            [0, y, 0, 0],
            [0, 0, z, 0],
            [0, 0, 0, 1]
        ])
        return shape.transformGeometry(t)

class CQCache:
    @staticmethod
    def import_brep(file_path: str):
        """
        Import a boundary representation model
        Returns a TopoDS_Shape object
        """
        builder = BRep_Builder()
        shape = TopoDS_Shape()
        return_code = BRepTools.Read_s(shape, file_path, builder)
        if return_code is False:
            raise ValueError("Import failed, check file name")
        return cq.Compound(shape)

    @staticmethod
    def get_cache_exists(obj: Union[Sequence[CQObject], CQObject]):
        cache_file_name = CQCache.get_file_name(obj)
        return os.path.isfile(cache_file_name)

    @staticmethod
    def get_file_name(shape: Union[Sequence[CQObject], CQObject]):
        prev_vector = cq.Vector(0,0,0)
        for vertex in cast(Sequence[cq.Vertex], CQLinq.select(shape, "vertex")):
            prev_vector += vertex.Center()
        
        hasher = hashlib.md5()
        hasher.update(bytes(str(tuple(np.round(prev_vector.toTuple(), 4))), "utf-8"))
        # encode the hash as a filesystem safe string
        shape_id = base64.urlsafe_b64encode(hasher.digest()).decode("utf-8")
        return f"{CACHE_DIR_PATH}/{shape_id}.brep"

    @staticmethod
    def export_brep(shape: cq.Shape, file_path: str):
        if CACHE_DIR_NAME not in os.listdir(TEMPDIR_PATH):
            os.mkdir(CACHE_DIR_PATH)
        shape.exportBrep(file_path)

    @staticmethod
    def clear_cache():
        if CACHE_DIR_NAME in os.listdir(TEMPDIR_PATH):
            for file in os.listdir(CACHE_DIR_PATH):
                os.remove(os.path.join(CACHE_DIR_PATH, file))



class IndexSelector(cq.Selector):
    def __init__(self, indices: Sequence[int]):
        self.indices = indices
    def filter(self, objectList):
        return [objectList[i] for i in self.indices]

class GroupSelector(cq.Selector):
    def __init__(self, allow: OrderedSet[CQObject], is_interior: bool = True):
        self.allow = allow
        self.is_interior = is_interior
    def filter(self, objectList):
        filtered_objs = []
        for obj in objectList:
            if obj in self.allow:
                filtered_objs.append(obj)
        return filtered_objs


meshql/utils/cq.py
OpenOrion/meshql
import numpy as np
import numpy.typing as npt
from typing import Sequence, Union
from plotly import graph_objects as go
from meshql.utils.types import NumpyFloat

def add_plot(coords: npt.NDArray[NumpyFloat], fig: go.Figure=go.Figure(), label: str="Plot"):
    dim = 3 if np.all(coords[:,2]) else 2
    if dim == 3:
        fig.add_scatter3d(
            x=coords[:,0],
            y=coords[:,1],
            z=coords[:,2],
            name=label,
            mode="lines+markers",
            marker=dict(size=1)       
        )
    else:
        fig.add_scatter(
            x=coords[:,0],
            y=coords[:,1],
            name=label,
            fill="toself",
            mode="lines+markers",
            marker=dict(size=1)       
        )

    return fig

meshql/utils/plot.py
OpenOrion/meshql
from typing import Sequence
from meshql.utils.types import Number
import numpy as np
import numpy.typing as npt

def get_sampling(start: Number, end: Number, num_samples: int, is_cosine_sampling: bool):
    if is_cosine_sampling:
        beta = np.linspace(0.0,np.pi, num_samples, endpoint=True)
        return 0.5*(1.0-np.cos(beta))*(end-start) + start
    else:
        return np.linspace(start, end, num_samples, endpoint=True)


def generate_circle(r, num_points=100):
    theta = np.linspace(0, 2*np.pi, num_points)
    x = r * np.cos(theta)
    y = r * np.sin(theta)
    return np.column_stack((x, y))

def generate_naca4_airfoil(naca_string: str, num_points: int = 100, use_cosine_sampling=True) -> npt.NDArray[np.float64]:

    """generates NACA4 coordinates

    Parameters
    ==========

    naca_string: str
        NACA4 string

    num_points: int
        number of points to generate
    """
    M = int(naca_string[0]) / 100
    P = int(naca_string[1]) / 10
    XX = int(naca_string[2:]) / 100

    # camber line
    if use_cosine_sampling:
        beta = np.linspace(0.0,np.pi, num_points)
        xc = 0.5*(1.0-np.cos(beta))
    else:
        xc = np.linspace(0.0, 1.0, num_points)
    # thickness distribution from camber line
    a0 = 0.2969
    a1 = -0.1260
    a2 = -0.3516
    a3 = 0.2843
    a4 = -0.1036
    yt = 5.0*XX*(np.sqrt(xc)*a0 + xc**4*a4 + xc**3*a3 + xc**2*a2 + xc*a1)
    
    # camber line slope
    if P == 0:
        xl = xu = xc
        yu = yt
        yl = -yt

    else:
        yc1 = M*(-xc**2 + 2*xc*P)/P**2
        yc2 = M*(-xc**2 + 2*xc*P - 2*P + 1)/(1 - P)**2
        yc = (np.select([np.logical_and.reduce((np.greater_equal(xc, 0),np.less(xc, P))),np.logical_and.reduce((np.greater_equal(xc, P),np.less_equal(xc, 1))),True], [yc1,yc2,1], default=np.nan))

        dyc1dx = M*(-2*xc + 2*P)/P**2
        dyc2dx = M*(-2*xc + 2*P)/(1 - P)**2
        dycdx = (np.select([np.logical_and.reduce((np.greater_equal(xc, 0),np.less(xc, P))),np.logical_and.reduce((np.greater_equal(xc, P),np.less_equal(xc, 1))),True], [dyc1dx,dyc2dx,1], default=np.nan))
        theta = np.arctan(dycdx)
        
        xu = xc - yt*np.sin(theta)
        yu = yc + yt*np.cos(theta)
        xl = xc + yt*np.sin(theta)
        yl = yc - yt*np.cos(theta)

    # thickness lines
    x = np.concatenate((xu[1:-1], xl[::-1]))
    y = np.concatenate((yu[1:-1], yl[::-1]))

    return np.column_stack((x, y)) # type: ignore

meshql/utils/shapes.py
OpenOrion/meshql
from typing import Iterable, Iterator, MutableSet, TypeVar, Union
import numpy as np

NumpyFloat = np.float32
Number = Union[int, float, NumpyFloat]
VectorTuple = Union[tuple[float, float, float], tuple[float, float]]
LineTuple = tuple[VectorTuple, VectorTuple]

T = TypeVar("T")
class OrderedSet(MutableSet[T]):
    """A set that preserves insertion order by internally using a dict."""

    def __init__(self, iterable: Iterable[T] = []):
        self._d = dict.fromkeys(iterable)

    def add(self, x: T) -> None:
        self._d[x] = None

    def discard(self, x: T) -> None:
        self._d.pop(x, None)

    def update(self, iterable: Iterable[T]) -> None:
        self._d.update(dict.fromkeys(iterable))

    @property
    def first(self) -> T:
        return next(iter(self))

    @property
    def last(self) -> T:
        *_, last = iter(self)
        return last

    def __contains__(self, x: object) -> bool:
        return self._d.__contains__(x)

    def __len__(self) -> int:
        return self._d.__len__()

    def __iter__(self) -> Iterator[T]:
        return self._d.__iter__()

    def __str__(self):
        return f"{{{', '.join(str(i) for i in self)}}}"

    def __repr__(self):
        return f"<OrderedSet {self}>"

    # def __getitem__(self, index):
    #     return list(self)[index]

meshql/utils/types.py
OpenOrion/meshql
from typing import Any, Sequence, Union, cast
from .mesh.mesh import Mesh
import pythreejs
from IPython.display import display
from IPython.core.display import HTML
import ipywidgets as widgets
import numpy as np
import colorsys

def generate_color_legend_html(title: str, color_labels: dict[str, list[int]]):
    title = f"<h2>{title}</h2>"
    legend = '<table>'
    for label, color in color_labels.items():
        assert len(color) == 3, "Color must be a list of 3 integers"
        legend += f'<tr><td style="background-color: {to_rgb_str(color)}" width="20"></td><td>{label}</td></tr>'
    legend += '</table>'
    return f'<div style="float: left; padding-right: 50px">{title+legend}</div>'


def generate_rgb_values(n_colors, is_grayscale=False):
    if n_colors == 0:
        return []
    colors=[]
    for i in np.arange(0., 360., 360. / n_colors):
        hue = i/360.
        if is_grayscale:
            min_rgb = 0.5
            rgb = (1 - min_rgb)*hue + min_rgb
            rgb_values = [rgb,rgb,rgb]
        else:
            lightness = (50 + np.random.rand() * 10)/100.
            saturation = (90 + np.random.rand() * 10)/100.
            rgb_values = list(colorsys.hls_to_rgb(hue, lightness, saturation))

        colors.append(rgb_values)


    return colors

def to_rgb_str(color: Sequence[int]):
    return f"rgb({int(color[0]*255)},{int(color[1]*255)},{int(color[2]*255)})"


def visualize_mesh(
    meshes: Union[Mesh, list[Mesh]], 
    view_width=800, 
    view_height=600,
    only_markers=False
):
    coord_html = widgets.HTML("Coords: ()")

    def on_surf_mousemove(change):
        # write coordinates to html container
        if change.new is None:
            coord_html.value = "Coords: ()"
        else:
            coord_html.value = "Coords: (%f, %f, %f)" % change.owner.point

    if not isinstance(meshes, list):
        meshes = [meshes]

    # Legend Colors
    mesh_colors = generate_rgb_values(len(meshes), is_grayscale=True)
    marker_colors = generate_rgb_values(sum([len(mesh.markers) for mesh in meshes]))

    # Legend Color Labels
    marker_color_labels = {}
    mesh_color_labels = {}

    marker_line_segments = []
    buffer_meshes = []
    target_point_spheres = []
    for i, mesh in enumerate(meshes):
        mesh_color = mesh_colors[i]
        mesh_color_labels[f"Zone {i}"] = mesh_color
        max_point = mesh.points.min(axis=0)
        min_point = mesh.points.max(axis=0)
        point_size = max(max_point[0] - min_point[0], max_point[1] - min_point[1])*0.03
        # Marker line segment points and colors
        marker_line_points = []
        marker_segment_colors = []
        marker_elements_to_name = {}
        for marker_name, marker_elements in mesh.markers.items():
            is_line = len(marker_elements) == 2
            for elements in marker_elements:
                # iterate all marker element line combinations
                for i in range(len(elements)):
                    line_from_to = (elements[i], elements[i+1]) if i+1 < len(elements) else (elements[-1], elements[0])
                    
                    # only allow mesh line label override if line is appart of edge if already labeled for edge mesh line
                    if (line_from_to in marker_elements_to_name and is_line) or line_from_to not in marker_elements_to_name:
                        marker_elements_to_name[line_from_to] = marker_name


                if marker_name in mesh.target_points and elements[1] in mesh.target_points[marker_name]:
                    target_point_sphere = pythreejs.Mesh(
                        geometry=pythreejs.SphereGeometry(radius=point_size),
                        material=pythreejs.MeshLambertMaterial(color='red', side='DoubleSide'),
                    )
                    target_point_sphere.position = mesh.points[elements[1]].tolist()
                    target_point_spheres.append(target_point_sphere)

        # Non-marker line segment points
        non_marker_line_points = []
        for point_tags in mesh.elements:
            for i in range(len(point_tags)):
                if i + 1 < len(point_tags):
                    line_point_tags = (point_tags[i+1], point_tags[i])
                else:
                    line_point_tags = (point_tags[0], point_tags[i])
                line_points = [mesh.points[line_point_tags[0]].tolist(), mesh.points[line_point_tags[1]].tolist()]

                marker_point_tags = line_point_tags if line_point_tags in marker_elements_to_name else line_point_tags[::-1]
                if marker_point_tags in marker_elements_to_name:
                    marker_name = marker_elements_to_name[marker_point_tags]
                    if marker_name not in marker_color_labels:
                        marker_color_labels[marker_name] = marker_colors[len(marker_color_labels)]
                    marker_color = marker_color_labels[marker_name]
                    marker_segment_colors.append([marker_color, marker_color])
                    marker_line_points.append(line_points)
                else:
                    non_marker_line_points.append(line_points)

        if not only_markers and len(non_marker_line_points) > 0:
            non_marker_lines = pythreejs.LineSegments2(
                cast(Any, pythreejs.LineSegmentsGeometry(positions=non_marker_line_points)),
                cast(Any, pythreejs.LineMaterial(linewidth=1, color=to_rgb_str(mesh_color)))
            )
            marker_line_segments.append(non_marker_lines)

        if len(marker_line_points) > 0:
            marker_lines = pythreejs.LineSegments2(
                cast(Any, pythreejs.LineSegmentsGeometry(positions=marker_line_points, colors=marker_segment_colors)),
                cast(Any, pythreejs.LineMaterial(linewidth=2, vertexColors='VertexColors'))
            )
            marker_line_segments.append(marker_lines)

        if len(mesh.elements) > 0:
            buffer_geom = pythreejs.BufferGeometry(attributes=dict(
                position=pythreejs.BufferAttribute(mesh.points, normalized=False),
                index=pythreejs.BufferAttribute(np.concatenate(mesh.elements), normalized=False),
            ))

            buffer_mesh = pythreejs.Mesh(
                geometry=buffer_geom,
                material=pythreejs.MeshLambertMaterial(color='white', side='DoubleSide'),
            )
            buffer_meshes.append(buffer_mesh)

    camera = pythreejs.PerspectiveCamera(position=[0, 0, 1], far=100000, near=0.001, aspect=cast(Any, view_width/view_height))
    scene = pythreejs.Scene(children=[*marker_line_segments, *buffer_meshes, *target_point_spheres, pythreejs.AmbientLight(intensity=cast(int, 0.8))], background="black")
    orbit_controls = pythreejs.OrbitControls(controlling=camera)

    pickable_objects = pythreejs.Group()
    for buffer_mesh in buffer_meshes:
        pickable_objects.add(buffer_mesh)

    mousemove_picker = pythreejs.Picker(
        controlling=pickable_objects,
        event='mousemove'
    )
    mousemove_picker.observe(on_surf_mousemove, names=cast(Any, ['faceIndex']))

    renderer = pythreejs.Renderer(
        camera=camera,
        scene=scene,
        controls=[orbit_controls, mousemove_picker],
        width=view_width,
        height=view_height
    )

    # Plot renderer
    display(coord_html, renderer)

    # Plot legend
    marker_legend_html = generate_color_legend_html("Markers", marker_color_labels)
    mesh_legend_html = generate_color_legend_html("Zones", mesh_color_labels)
    display(HTML(marker_legend_html+mesh_legend_html))

meshql/visualizer.py
OpenOrion/meshql


from setuptools import setup

setup(
   name='meshql',
   version='1.0.0',
   description='query based meshing on top of GMSH',
   author='Afshawn Lotfi',
   author_email='',
   packages=['meshql', 'meshql.mesh', 'meshql.preprocessing', 'meshql.transactions', 'meshql.utils'],
   install_requires=[
    "numpy",
    "gmsh",
    "ipywidgets==7.6",
    "ipython_genutils",
    "pythreejs==2.4.2",
    "su2fmt @ git+https://github.com/Turbodesigner/su2fmt.git",
    "shapely",
    "scipy",
    "jupyter_cadquery"
   ]
)
setup.py
OpenOrion/meshql