Fix edge refinement

examples/MeshAPI/VerciseMeshDirected_SetMeshSizes.py

@@ -0,0 +1,109 @@
+"""
+Example of a full model geometry
+with one Vercise electrode without
+encapsulation layer.
+The mesh is generated and
+local mesh refinements are taken
+into account during the meshing
+process.
+"""
+import ngsolve
+
+import ossdbs
+from ossdbs.utils.vtk_export import FieldSolution
+
+
+def save(cf, filename: str) -> None:
+    """Save solution in VTK format."""
+    names = [f"{cf.label}_real"]
+    if cf.is_complex:
+        names.append(f"{cf.label}_imag")
+
+    coefficients = [cf.solution.real]
+    if cf.is_complex:
+        coefficients.append(cf.solution.imag)
+
+    vtk = ngsolve.VTKOutput(
+        ma=cf.mesh,
+        coefs=coefficients,
+        names=names,
+        filename=filename,
+        subdivision=0,
+    )
+    vtk.Do(vb=ngsolve.BND)
+
+
+settings = {
+    "Electrodes": [
+        {
+            "Name": "BostonScientificVerciseDirected",
+            "Rotation[Degrees]": 0,
+            "Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
+            "TipPosition": {"x[mm]": 0, "y[mm]": 15, "z[mm]": -3},
+            "EncapsulationLayer": {
+                "Thickness[mm]": 0.0,  # indicates that no encapsulation is modelled
+            },
+            "Contacts": [
+                {
+                    "Contact_ID": 1,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 2,
+                    "MaxMeshSizeEdge": 0.01,
+                },
+                {
+                    "Contact_ID": 3,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 4,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 5,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 6,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 7,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 8,
+                    "MaxMeshSizeEdge": 0.05,
+                },
+            ],
+        }
+    ],
+    "BrainRegion": {
+        "Center": {"x[mm]": 5, "y[mm]": 14, "z[mm]": -4.5},
+        "Dimension": {"x[mm]": 50.0, "y[mm]": 50.0, "z[mm]": 50.0},
+        "Shape": "Ellipsoid",
+    },
+    "Mesh": {"LoadMesh": False, "SaveMesh": False},
+    "ExportElectrode": False,
+}
+
+# generate mesh with brain
+mesh = ossdbs.generate_mesh(settings)
+
+print(mesh.boundaries)
+print(mesh.materials)
+
+bnd_dict = {}
+for idx in range(1, 9):
+    bnd_dict[f"E1C{idx}"] = idx
+
+ngsolve.Draw(mesh.boundary_coefficients(bnd_dict), mesh.ngsolvemesh, "bnd")
+
+cf = FieldSolution(
+    solution=mesh.boundary_coefficients(bnd_dict),
+    label="bnd",
+    mesh=mesh.ngsolvemesh,
+    is_complex=False,
+)
+save(cf, "electrode_in_brain")

examples/MeshAPI/VerciseMeshDirected_SetMeshSizes_RefineOne.py

@@ -0,0 +1,76 @@
+"""
+Example of a full model geometry
+with one Vercise electrode without
+encapsulation layer.
+The mesh is generated and
+local mesh refinements are taken
+into account during the meshing
+process.
+"""
+import ossdbs
+
+settings = {
+    "Electrodes": [
+        {
+            "Name": "BostonScientificVerciseDirected",
+            "Rotation[Degrees]": 0,
+            "Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
+            "TipPosition": {"x[mm]": 0, "y[mm]": 15, "z[mm]": -3},
+            "EncapsulationLayer": {
+                "Thickness[mm]": 0.0,  # indicates that no encapsulation is modelled
+            },
+            "Contacts": [
+                {
+                    "Contact_ID": 1,
+                    "Active": True,
+                    "Current[A]": 0.0,
+                    "Voltage[V]": 1.0,
+                    "Floating": False,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0},
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 3,
+                    "Active": True,
+                    "Current[A]": 0.0,
+                    "Voltage[V]": 0.0,
+                    "Floating": False,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0},
+                    "MaxMeshSizeEdge": 0.05,
+                },
+                {
+                    "Contact_ID": 5,
+                    "Active": False,
+                    "Current[A]": 0.0,
+                    "Voltage[V]": 0.0,
+                    "Floating": False,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0},
+                    "MaxMeshSizeEdge": 0.05,
+                },
+            ],
+        }
+    ],
+    "BrainRegion": {
+        "Center": {"x[mm]": 5, "y[mm]": 14, "z[mm]": -4.5},
+        "Dimension": {"x[mm]": 50.0, "y[mm]": 50.0, "z[mm]": 50.0},
+        "Shape": "Ellipsoid",
+    },
+    "ExportElectrode": False,
+}
+
+electrodes = ossdbs.generate_electrodes(settings)
+vercise = electrodes[0]
+electrode_settings = settings["Electrodes"][0]
+# set edge mesh sizes
+for edge in vercise.geometry.edges:
+    if edge.name is not None:
+        contact_number = int(edge.name.replace("Contact_", ""))
+        for contact_setting in electrode_settings["Contacts"]:
+            if contact_setting["Contact_ID"] == contact_number:
+                if "MaxMeshSizeEdge" in contact_setting:
+                    edge.maxh = contact_setting["MaxMeshSizeEdge"]
+                continue
+
+vercise.export_electrode(
+    output_path=".", brain_dict=settings["BrainRegion"], n_electrode=0
+)

ossdbs/electrodes/abbott_stjude.py

@@ -233,7 +233,7 @@ def _contacts(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
             else:
                 # Label all the named contacts appropriately
                 for edge in contact.edges:
-                    if edge.name is not None:
+                    if edge.name == "Rename":
                         edge.name = name
         return netgen.occ.Fuse(contacts)
 
@@ -253,36 +253,22 @@ def _contact_directed(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
         # Centering contact to label edges
         contact = contact.Rotate(axis, angle)
 
-        # TODO refactor / wrap in function
-        # Find  max z, min z, max x, and max y values and label min x and min y edge
-        max_z_val = max_y_val = max_x_val = float("-inf")
-        min_z_val = float("inf")
+        # Label all outer edges
         for edge in contact.edges:
-            if edge.center.z > max_z_val:
-                max_z_val = edge.center.z
-            if edge.center.z < min_z_val:
-                min_z_val = edge.center.z
-            if edge.center.x > max_x_val:
-                max_x_val = edge.center.x
-                max_x_edge = edge
-            if edge.center.y > max_y_val:
-                max_y_val = edge.center.y
-                max_y_edge = edge
-        max_x_edge.name = "max x"
-        max_y_edge.name = "max y"
-        # Label only the outer edges of the contact with min z and max z values
-        for edge in contact.edges:
-            if np.isclose(edge.center.z, max_z_val) and not (
-                np.isclose(edge.center.x, radius / 2)
-                or np.isclose(edge.center.y, radius / 2)
-            ):
-                edge.name = "max z"
-            elif np.isclose(edge.center.z, min_z_val) and not (
-                np.isclose(edge.center.x, radius / 2)
-                or np.isclose(edge.center.y, radius / 2)
-            ):
-                edge.name = "min z"
+            edge_center = np.array([edge.center.x, edge.center.y, edge.center.z])
+
+            # Skip center edge
+            if np.allclose(np.cross(edge_center, self._direction), 0):
+                continue
+
+            new_center = np.dot(edge_center, self._direction) * np.array(
+                self._direction
+            )
+
+            # Mark only outer edges
+            if not np.isclose(np.linalg.norm(edge_center - new_center), radius / 2):
+                edge.name = "Rename"
+
         contact = contact.Rotate(axis, -angle)
-        # TODO check that the starting axis of the contacts
-        # are correct according to the documentation
+
         return contact

ossdbs/electrodes/boston_scientific_cartesia.py

@@ -168,7 +168,7 @@ def _contacts(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
             else:
                 # Label all the named contacts appropriately
                 for edge in contact.edges:
-                    if edge.name is not None:
+                    if edge.name == "Rename":
                         edge.name = name
         return netgen.occ.Fuse(contacts)
 
@@ -202,20 +202,20 @@ def _contact_directed(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
             if edge.center.y > max_y_val:
                 max_y_val = edge.center.y
                 max_y_edge = edge
-        max_x_edge.name = "max x"
-        max_y_edge.name = "max y"
+        max_x_edge.name = "Rename"
+        max_y_edge.name = "Rename"
         # Label only the outer edges of the contact with min z and max z values
         for edge in contact.edges:
             if np.isclose(edge.center.z, max_z_val) and not (
                 np.isclose(edge.center.x, radius / 2)
                 or np.isclose(edge.center.y, radius / 2)
             ):
-                edge.name = "max z"
+                edge.name = "Rename"
             elif np.isclose(edge.center.z, min_z_val) and not (
                 np.isclose(edge.center.x, radius / 2)
                 or np.isclose(edge.center.y, radius / 2)
             ):
-                edge.name = "min z"
+                edge.name = "Rename"
 
         # TODO check that the starting axis of the contacts
         # are correct according to the documentation

ossdbs/electrodes/boston_scientific_vercise.py

@@ -145,10 +145,11 @@ def _contacts(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
             else:
                 # Label all the named contacts appropriately
                 for edge in contact.edges:
-                    if edge.name is not None:
+                    if edge.name == "Rename":
                         edge.name = name
         return netgen.occ.Fuse(contacts)
 
+    # ruff: noqa: C901
     def _contact_directed(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
         point = (0, 0, 0)
         radius = self._parameters.lead_diameter * 0.5
@@ -164,38 +165,23 @@ def _contact_directed(self) -> netgen.libngpy._NgOCC.TopoDS_Shape:
         contact = body - eraser.Rotate(axis, angle) - eraser.Rotate(axis, -angle)
         # Centering contact to label edges
         contact = contact.Rotate(axis, angle)
-        # TODO refactor / wrap in function
-        # Find  max z, min z, max x, and max y values and label min x and min y edge
-        max_z_val = max_y_val = max_x_val = float("-inf")
-        min_z_val = float("inf")
-        for edge in contact.edges:
-            if edge.center.z > max_z_val:
-                max_z_val = edge.center.z
-            if edge.center.z < min_z_val:
-                min_z_val = edge.center.z
-            if edge.center.x > max_x_val:
-                max_x_val = edge.center.x
-                max_x_edge = edge
-            if edge.center.y > max_y_val:
-                max_y_val = edge.center.y
-                max_y_edge = edge
-        max_x_edge.name = "max x"
-        max_y_edge.name = "max y"
-        # Label only the outer edges of the contact with min z and max z values
+
+        # Label all outer edges
         for edge in contact.edges:
-            if np.isclose(edge.center.z, max_z_val) and not (
-                np.isclose(edge.center.x, radius / 2)
-                or np.isclose(edge.center.y, radius / 2)
-            ):
-                edge.name = "max z"
-            elif np.isclose(edge.center.z, min_z_val) and not (
-                np.isclose(edge.center.x, radius / 2)
-                or np.isclose(edge.center.y, radius / 2)
-            ):
-                edge.name = "min z"
-
-        # TODO check that the starting axis of the contacts
-        # are correct according to the documentation
+            edge_center = np.array([edge.center.x, edge.center.y, edge.center.z])
+
+            # Skip center edge
+            if np.allclose(np.cross(edge_center, self._direction), 0):
+                continue
+
+            new_center = np.dot(edge_center, self._direction) * np.array(
+                self._direction
+            )
+
+            # Mark only outer edges
+            if not np.isclose(np.linalg.norm(edge_center - new_center), radius / 2):
+                edge.name = "Rename"
+
         contact = contact.Rotate(axis, -angle)
 
         return contact