further work

examples/scripts/meshing/gmsh_basic.py

@@ -0,0 +1,41 @@
+import gmsh
+
+# Initialize Gmsh
+gmsh.initialize()
+gmsh.model.add("split_line_occ_example")
+
+# Define the OCC geometry
+# Add two points to define the original line
+p1 = gmsh.model.occ.addPoint(0, 0, 0)
+p2 = gmsh.model.occ.addPoint(1, 0, 0)
+p3 = gmsh.model.occ.addPoint(0.5, -2, 0)
+p4 = gmsh.model.occ.addPoint(0.5, 1, 0)
+
+# Create a line segment between these points using OCC
+line = gmsh.model.occ.addLine(p1, p2)
+line2 = gmsh.model.occ.addLine(p3, p4)
+
+# Add the splitting point (vertex) along the line
+split_point = gmsh.model.occ.addPoint(0.5, 0, 0)
+gmsh.model.occ.synchronize()
+gmsh.fltk.run()
+
+# Perform the fragmentation using occ.fragment
+# The first argument is the list of objects to fragment (the line)
+# The second argument is the list of objects that will be used to fragment the first list (the vertex)
+result = gmsh.model.occ.fragment([(1, line2)], [(0, split_point)])
+
+# The result of the fragment operation contains the new entities, which need to be extracted
+fragmented_lines = [entity for entity in result[0] if entity[0] == 1]
+
+# Synchronize to apply the changes
+gmsh.model.occ.synchronize()
+
+# Optionally, generate a mesh to see the result
+gmsh.model.mesh.generate(1)
+
+# Launch the Gmsh GUI to visualize the result
+gmsh.fltk.run()
+
+# Finalize Gmsh
+gmsh.finalize()

examples/scripts/meshing/partitioning_edge_intersect.py

@@ -22,7 +22,7 @@ def edges_intersect(use_xact=False):
     Config().update_config_globally(
         "meshing_open_viewer_breakpoint_names",
         [
-            # "partition_bm_split_pre",
+            #"partition_isect_bm_loop",
             # "partition_isect_bm_pre",
             # "partition_bm_split_cut_1"
         ],

examples/scripts/meshing/partitioning_perp_plates_and_beams.py

@@ -22,10 +22,13 @@ def main():
     pl3 = ada.Plate("pl3", p1x1, 0.01, orientation=ada.Placement(xdir=(1, 0, 0), zdir=(0, -1, 0)))
 
     p = ada.Part("MyFem") / [pl1_5, pl3, bm_1, bm_2]
+
     fem = p.to_fem_obj(1, bm_repr="line", pl_repr="shell")
+    fem.show(solid_beams=True)
     p.fem = fem
+
     a = ada.Assembly("Test") / p
-    a.to_fem("ADA_pl_mesh", "usfos", overwrite=True)
+    # a.to_fem("ADA_pl_mesh", "usfos", overwrite=True)
 
     assert len(fem.nodes) == 10
     assert len(fem.elements) == 13

src/ada/api/spatial/part.py

@@ -673,8 +673,9 @@ def beam_clash_check(self, margins=5e-5):
 
         all_parts = self.get_all_subparts() + [self]
         all_beams = [bm for p in all_parts for bm in p.beams]
+        all_bm_containers = [p.beams for p in all_parts]
 
-        return filter(None, [basic_intersect(bm, margins, all_parts) for bm in all_beams])
+        return filter(None, [basic_intersect(bm, margins, all_bm_containers) for bm in all_beams])
 
     def move_all_mats_and_sec_here_from_subparts(self):
         for p in self.get_all_subparts():
@@ -732,6 +733,7 @@ def to_fem_obj(
 
         options = GmshOptions(Mesh_Algorithm=8) if options is None else options
         masses: list[Shape] = []
+
         with GmshSession(silent=silent, options=options, debug_mode=debug_mode) as gs:
             for obj in self.get_all_physical_objects(sub_elements_only=False):
                 if isinstance(obj, Beam):
@@ -748,8 +750,8 @@ def to_fem_obj(
             if interactive is True:
                 gs.open_gui()
 
-            gs.partition_plates()
 
+            gs.partition_plates()
             gs.partition_beams()
 
             if interactive is True:

src/ada/core/clash_check.py

@@ -3,20 +3,24 @@
 import traceback
 from dataclasses import dataclass
 from itertools import chain
-from typing import Iterable, List
+from typing import Iterable, List, TYPE_CHECKING
 
 import numpy as np
 
 import ada
-from ada import Assembly, Beam, Node, Part, Pipe, PipeSegStraight, Plate, PrimCyl
+
 from ada.config import logger
 
 from ..api.transforms import EquationOfPlane
 from .utils import Counter
 from .vector_utils import intersect_calc, is_parallel, vector_length
 
+if TYPE_CHECKING:
+    from ada import Assembly, Beam, Node, Part, Pipe, PipeSegStraight, Plate, PrimCyl
+    from ada.api.containers import Beams
+
 
-def basic_intersect(bm: Beam, margins, all_parts: [Part]):
+def basic_intersect(bm: Beam, margins, all_beam_containers: [Beams]):
     if bm.section.type == "gensec":
         return bm, []
     try:
@@ -27,7 +31,7 @@ def basic_intersect(bm: Beam, margins, all_parts: [Part]):
     vol_in = [x for x in zip(vol[0], vol[1])]
     beams = filter(
         lambda x: x != bm,
-        chain.from_iterable([p.beams.get_beams_within_volume(vol_in, margins=margins) for p in all_parts]),
+        chain.from_iterable([beams.get_beams_within_volume(vol_in, margins=margins) for beams in all_beam_containers]),
     )
     return bm, beams
 
@@ -63,6 +67,8 @@ def beam_cross_check(bm1: Beam, bm2: Beam, outofplane_tol=0.1):
 
 def are_beams_connected(bm1: Beam, beams: List[Beam], out_of_plane_tol, point_tol, nodes, nmap) -> None:
     # TODO: Function should be renamed, or return boolean. Unclear what the function does at the moment
+    from ada import Node
+
     for bm2 in beams:
         if bm1 == bm2:
             continue
@@ -140,6 +146,7 @@ def filter_beams_along_plate_edges(pl: Plate, beams: Iterable[Beam]):
 def find_beams_connected_to_plate(pl: Plate, beams: list[Beam]) -> list[Beam]:
     """Return all beams with their midpoints inside a specified plate for a given list of beams"""
     from ada.api.containers import Nodes
+    from ada import Node
 
     nid = Counter(1)
     nodes = Nodes(

src/ada/fem/concept.py

@@ -392,7 +392,7 @@ def show(
         unique_id=None,
         purpose: FilePurposeDC = FilePurposeDC.ANALYSIS,
         params_override: RenderParams = None,
-        solid_beams=False,
+        solid_beams=True,
         ping_timeout=1,
     ) -> None:
 

src/ada/fem/meshing/concepts.py

@@ -177,11 +177,13 @@ def make_cuts(self, keep_cutting_planes=True):
         self.model.geo.synchronize()
 
     def partition_beams(self):
-        from ada.fem.meshing.partitioning.partition_intersecting_beams import split_crossing_beams
+        from ada.fem.meshing.partitioning.partition_beams import split_crossing_beams, split_intersecting_beams
         plates = [obj for obj in self.model_map.keys() if type(obj) is Plate]
 
-        if len(plates) == 0: # For some reason this
+        if len(plates) == 0: # For some reason this breaks the model whenever plates are in the model
             split_crossing_beams(self)
+        else:
+            split_intersecting_beams(self)
 
     def partition_plates(self):
         """Split plates that are intersecting each other"""

src/ada/fem/meshing/partitioning/partition_beam_interiors.py

@@ -0,0 +1,67 @@
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+from ada import CurvePoly2d, Placement
+
+if TYPE_CHECKING:
+    from ada import Beam
+
+    from ..concepts import GmshData, GmshSession
+
+
+def ibeam(model: "GmshData", gmsh_session: "GmshSession"):
+    # pass
+    bm_obj: "Beam" = model.obj
+    for cut in make_ig_cutplanes(bm_obj):
+        gmsh_session.add_cutting_plane(cut, [model])
+
+    gmsh_session.make_cuts()
+
+    for dim, tag in gmsh_session.model.get_entities():
+        if dim == 2:
+            gmsh_session.model.mesh.set_transfinite_surface(tag)
+            gmsh_session.model.mesh.setRecombine(dim, tag)
+
+    for dim, tag in gmsh_session.model.get_entities():
+        if dim == 3:
+            gmsh_session.model.mesh.set_transfinite_volume(tag)
+            gmsh_session.model.mesh.setRecombine(dim, tag)
+
+    # gmsh_session.open_gui()
+
+    # raise NotImplementedError()
+
+
+def make_ig_cutplanes(bm: "Beam"):
+    from ..concepts import CutPlane
+
+    points2d = bm.section.get_section_profile().outer_curve.points2d
+    sec_place = Placement(bm.n1.p, bm.yvec, bm.up)
+    points3d = sec_place.transform_local_points_back_to_global(points2d)
+
+    minz = min([x[2] for x in points3d])
+    maxz = max([x[2] for x in points3d])
+    pmin, pmax = bm.bbox().p1, bm.bbox().p2
+    dx, dy, dz = (np.array(pmax) - np.array(pmin)) * 1.0
+    x, y, _ = pmin
+
+    sec = bm.section
+
+    cut1 = CutPlane((x, y, minz + sec.t_fbtn), dx=dx, dy=dy)
+    cut2 = CutPlane((x, y, maxz - sec.t_fbtn), dx=dx, dy=dy)
+
+    web_left = bm.n1.p - (sec.t_w / 2) * bm.yvec - (sec.h / 2) * bm.up
+    web_right = bm.n1.p + (sec.t_w / 2) * bm.yvec - (sec.h / 2) * bm.up
+    dy = sec.h
+    cut3 = CutPlane(web_left, dx=dx, dy=dy, plane="XZ")
+    cut4 = CutPlane(web_right, dx=dx, dy=dy, plane="XZ")
+
+    return [cut1, cut2, cut3, cut4]
+
+
+def get_bm_section_curve(bm: "Beam", origin=None) -> CurvePoly2d:
+    origin = origin if origin is not None else bm.n1.p
+    section_profile = bm.section.get_section_profile(True)
+    points2d = section_profile.outer_curve.points2d
+    return CurvePoly2d(points2d=points2d, origin=origin, xdir=bm.yvec, normal=bm.xvec, parent=bm.parent)

src/ada/fem/meshing/partitioning/partition_beams.py

@@ -1,67 +1,68 @@
-from typing import TYPE_CHECKING
+from ada.config import Config
+from ada.fem.meshing import GmshSession
 
-import numpy as np
 
-from ada import CurvePoly2d, Placement
-
-if TYPE_CHECKING:
+def split_crossing_beams(gmsh_session: GmshSession):
     from ada import Beam
+    br_names = Config().meshing_open_viewer_breakpoint_names
 
-    from ..concepts import GmshData, GmshSession
-
-
-def ibeam(model: "GmshData", gmsh_session: "GmshSession"):
-    # pass
-    bm_obj: "Beam" = model.obj
-    for cut in make_ig_cutplanes(bm_obj):
-        gmsh_session.add_cutting_plane(cut, [model])
-
-    gmsh_session.make_cuts()
+    beams = [obj for obj in gmsh_session.model_map.keys() if type(obj) is Beam]
+    if len(beams) == 1:
+        return None
 
-    for dim, tag in gmsh_session.model.get_entities():
-        if dim == 2:
-            gmsh_session.model.mesh.set_transfinite_surface(tag)
-            gmsh_session.model.mesh.setRecombine(dim, tag)
+    if br_names is not None and "partition_isect_bm_pre" in br_names:
+        gmsh_session.open_gui()
 
-    for dim, tag in gmsh_session.model.get_entities():
-        if dim == 3:
-            gmsh_session.model.mesh.set_transfinite_volume(tag)
-            gmsh_session.model.mesh.setRecombine(dim, tag)
+    intersecting_beams = []
+    int_bm_map = dict()
+    for bm in beams:
+        bm_gmsh_obj = gmsh_session.model_map[bm]
+        for li_dim, li_ent in bm_gmsh_obj.entities:
+            intersecting_beams.append((li_dim, li_ent))
+            int_bm_map[(li_dim, li_ent)] = bm_gmsh_obj
 
-    # gmsh_session.open_gui()
+    res, res_map = gmsh_session.model.occ.fragment(
+        intersecting_beams, intersecting_beams, removeTool=True, removeObject=True
+    )
 
-    # raise NotImplementedError()
+    for i, int_bm in enumerate(intersecting_beams):
+        bm_gmsh_obj = int_bm_map[int_bm]
+        new_ents = res_map[i]
+        bm_gmsh_obj.entities = new_ents
 
+    gmsh_session.model.occ.synchronize()
 
-def make_ig_cutplanes(bm: "Beam"):
-    from ..concepts import CutPlane
 
-    points2d = bm.section.get_section_profile().outer_curve.points2d
-    sec_place = Placement(bm.n1.p, bm.yvec, bm.up)
-    points3d = sec_place.transform_local_points_back_to_global(points2d)
+def split_intersecting_beams(gmsh_session: GmshSession, margins=5e-5, out_of_plane_tol=0.1, point_tol=Config().general_point_tol):
+    from ada import Beam, Node
+    from ada.api.containers import Beams, Nodes
+    from ada.core.clash_check import basic_intersect, are_beams_connected
 
-    minz = min([x[2] for x in points3d])
-    maxz = max([x[2] for x in points3d])
-    pmin, pmax = bm.bbox().p1, bm.bbox().p2
-    dx, dy, dz = (np.array(pmax) - np.array(pmin)) * 1.0
-    x, y, _ = pmin
+    br_names = Config().meshing_open_viewer_breakpoint_names
 
-    sec = bm.section
+    all_beams = [obj for obj in gmsh_session.model_map.keys() if type(obj) is Beam]
+    bm_cont = Beams(all_beams)
+    if len(all_beams) == 1:
+        return None
 
-    cut1 = CutPlane((x, y, minz + sec.t_fbtn), dx=dx, dy=dy)
-    cut2 = CutPlane((x, y, maxz - sec.t_fbtn), dx=dx, dy=dy)
+    nodes = Nodes()
+    nmap: dict[Node, list[Beam]] = dict()
+    for bm, cbeams in filter(None, [basic_intersect(bm, margins, [bm_cont]) for bm in all_beams]):
+        are_beams_connected(bm, cbeams, out_of_plane_tol, point_tol, nodes, nmap)
 
-    web_left = bm.n1.p - (sec.t_w / 2) * bm.yvec - (sec.h / 2) * bm.up
-    web_right = bm.n1.p + (sec.t_w / 2) * bm.yvec - (sec.h / 2) * bm.up
-    dy = sec.h
-    cut3 = CutPlane(web_left, dx=dx, dy=dy, plane="XZ")
-    cut4 = CutPlane(web_right, dx=dx, dy=dy, plane="XZ")
+    for n, beams in nmap.items():
+        split_point = gmsh_session.model.occ.addPoint(n.x, n.y, n.z)
+        for bm in beams:
+            if n == bm.n1 or n == bm.n2:
+                continue
+            bm_gmsh_obj = gmsh_session.model_map[bm]
+            # entities_1 = gmsh_session.model.occ.get_entities(1)
 
-    return [cut1, cut2, cut3, cut4]
+            res, res_map = gmsh_session.model.occ.fragment(bm_gmsh_obj.entities, [(0, split_point)])
 
+            bm_gmsh_obj.entities = [x for x in res_map[0] if x[0] == 1]
+            gmsh_session.model.occ.synchronize()
+            if br_names is not None and "partition_isect_bm_loop" in br_names:
+                gmsh_session.open_gui()
 
-def get_bm_section_curve(bm: "Beam", origin=None) -> CurvePoly2d:
-    origin = origin if origin is not None else bm.n1.p
-    section_profile = bm.section.get_section_profile(True)
-    points2d = section_profile.outer_curve.points2d
-    return CurvePoly2d(points2d=points2d, origin=origin, xdir=bm.yvec, normal=bm.xvec, parent=bm.parent)
+    gmsh_session.model.occ.synchronize()