add mesh scripts

examples/ConvergenceStudy/FiguresReproducibility/MeshRefinement/generate_meshes.py

@@ -0,0 +1,151 @@
+# gold standard:
+# impedance error 0.1% and maximum element size equal
+# to mesh size
+
+import json
+import logging
+import os
+from copy import deepcopy
+
+import ngsolve
+import numpy as np
+
+import ossdbs
+
+
+def remove_file_handler(logger):
+    """Remove file handler from logger instance."""
+    for h in logger.handlers:
+        if isinstance(h, logging.FileHandler):
+            logger.removeHandler(h)
+
+
+def save_input_dict(base_input_dict):
+    """Save input dictionary for convergence run."""
+    input_dict = deepcopy(base_input_dict)
+    input_dict["OutputPath"] = "./"
+
+    with open(
+        os.path.join(base_input_dict["OutputPath"], "input_dict.json"), "w"
+    ) as fp:
+        json.dump(input_dict, fp, indent=2)
+
+
+BASE_CONTACT = {
+    "Contact_ID": 1,
+    "Active": False,
+    "Current[A]": 0.0,
+    "Voltage[V]": 0.0,
+    "Floating": False,
+}
+
+ossdbs.set_logger(logging.INFO)
+_logger = logging.getLogger("ossdbs")
+
+electrode_name = "Medtronic3389"
+
+with open("oss_dbs_parameters.json") as fp:
+    base_input_dict = json.load(fp)
+
+
+# add electrode
+base_input_dict["Electrodes"][0]["Name"] = electrode_name
+
+# initially no adaptive refinement
+base_input_dict["Mesh"]["AdaptiveMeshRefinement"] = {}
+base_input_dict["Mesh"]["AdaptiveMeshRefinement"]["Active"] = False
+
+# mri image
+mri_image, _ = ossdbs.load_images(base_input_dict)
+# pathway info
+pw = ossdbs.point_analysis.Pathway(base_input_dict["PointModel"]["Pathway"]["FileName"])
+pw.write_netgen_meshsize_file(
+    meshsize=min(mri_image.voxel_sizes), filename="meshsizes.txt"
+)
+
+# first refinement level: Default
+base_input_dict["Mesh"]["MeshingHypothesis"]["Type"] = "Default"
+base_input_dict["Mesh"]["MeshingHypothesis"]["MaxMeshSize"] = 1e6
+base_input_dict["OutputPath"] = "Results_PAM_default"
+electrodes = ossdbs.api.generate_electrodes(base_input_dict)
+brain_model = ossdbs.api.build_brain_model(base_input_dict, mri_image)
+try:
+    geometry = ossdbs.ModelGeometry(brain_model, electrodes)
+except RuntimeError:
+    brain_model = ossdbs.api.build_brain_model(
+        base_input_dict, mri_image, rotate_initial_geo=True
+    )
+    geometry = ossdbs.ModelGeometry(brain_model, electrodes)
+
+
+solver = ossdbs.api.prepare_solver(base_input_dict)
+dielectric_properties = ossdbs.api.prepare_dielectric_properties(base_input_dict)
+materials = base_input_dict["MaterialDistribution"]["MRIMapping"]
+conductivity = ossdbs.ConductivityCF(
+    mri_image,
+    brain_model.brain_region,
+    dielectric_properties,
+    materials,
+    geometry.encapsulation_layers,
+    complex_data=base_input_dict["EQSMode"],
+)
+
+
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+
+# refinement level: default and refinement around pathway
+base_input_dict["Mesh"]["MeshingHypothesis"]["MeshSizeFilename"] = "meshsizes.txt"
+base_input_dict["OutputPath"] = "Results_PAM_default_meshsize"
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+# reset mesh size file name
+base_input_dict["Mesh"]["MeshingHypothesis"]["MeshSizeFilename"] = ""
+
+
+# fourth refinement: material refinement
+base_input_dict["OutputPath"] = "Results_PAM_material_refinement"
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    # two refinement steps
+    volume_conductor.refine_mesh_by_material(2)
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()
+
+# sixth refinement: more edge refinement
+base_input_dict["Mesh"]["MeshingHypothesis"]["Type"] = "Default"
+lead_diameter = ossdbs.electrodes.default_electrode_parameters[
+    electrode_name
+].lead_diameter
+perimeter = np.pi * lead_diameter
+edge_size = perimeter / 50.0
+print(edge_size)
+base_input_dict["Electrodes"][0]["Contacts"][0]["MaxMeshSizeEdge"] = edge_size
+base_input_dict["Electrodes"][0]["Contacts"][1]["MaxMeshSizeEdge"] = edge_size
+base_input_dict["OutputPath"] = "Results_PAM_fine_edge_refinement"
+ossdbs.api.set_contact_and_encapsulation_layer_properties(base_input_dict, geometry)
+
+
+with ngsolve.TaskManager():
+    volume_conductor = ossdbs.api.prepare_volume_conductor_model(
+        base_input_dict, geometry, conductivity, solver
+    )
+    geo = volume_conductor.mesh.ngsolvemesh.ngmesh.GetGeometry()
+    for edge in geo.shape.edges:
+        print(edge.name, edge.maxh)
+    print(volume_conductor.mesh.ngsolvemesh.ne)
+    print("Save material")
+    volume_conductor.export_material_distribution_to_vtk()

examples/ConvergenceStudy/FiguresReproducibility/MeshRefinement/oss_dbs_parameters.json

@@ -0,0 +1,150 @@
+{
+    "BrainRegion": {
+        "Center": {"x[mm]": 13.5035, "y[mm]": -9.49908, "z[mm]": -9.87216},
+        "Dimension": {"x[mm]": 60, "y[mm]": 60, "z[mm]": 80},
+        "Shape": "Ellipsoid"
+    },
+    "Electrodes": [
+        {
+            "Name": "",
+            "Rotation[Degrees]": 0.15176099498051523,
+            "Direction": {
+                "x[mm]": 0.0,
+                "y[mm]": 0.0,
+                "z[mm]": 1.0
+            },
+            "TipPosition": {
+                "x[mm]": 12.956271802141353,
+                "y[mm]": -9.901870551007098,
+                "z[mm]": -12.01710780157648
+            },
+            "EncapsulationLayer": {
+                "Thickness[mm]": 0.0,
+                "Material": "Gray matter",
+                "DielectricModel": "ColeCole4",
+                "DielectricParameters": null,
+                "MaxMeshSize": 0.1
+            },
+            "Contacts": [
+                {
+                    "Contact_ID": 1,
+                    "Active": false,
+                    "Current[A]": 0.0,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 2,
+                    "Active": false,
+                    "Current[A]": 0.012,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 3,
+                    "Active": false,
+                    "Current[A]": -0.004,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                },
+                {
+                    "Contact_ID": 4,
+                    "Active": false,
+                    "Current[A]": -0.008,
+                    "Voltage[V]": 0.0,
+                    "Floating": true,
+                    "SurfaceImpedance[Ohmm]": {"real": 0.0, "imag": 0.0}
+                }
+            ]
+        }
+    ],
+    "Surfaces": [
+        {"Name": "BrainSurface", "Active": true, "Current[A]": 0.0, "Voltage[V]": 0.0}
+    ],
+    "MaterialDistribution": {
+        "MRIPath": "../../PAM/segmask.nii.gz",
+        "MRIMapping": {
+            "Blood": 4,
+            "Gray matter": 1,
+            "White matter": 2,
+            "CSF": 3,
+            "Unknown": 0
+        },
+        "DiffusionTensorActive": false,
+        "DTIPath": ""
+    },
+    "DielectricModel": {"Type": "ColeCole4", "CustomParameters": null},
+    "Mesh": {
+        "LoadMesh": false,
+        "LoadPath": "",
+        "MeshingHypothesis": {"Type": "Default"},
+        "HPRefinement": {"Active": false, "Levels": 2, "Factor": 0.125},
+        "AdaptiveMeshRefinement": {
+            "Active": false,
+            "MaxIterations": 10,
+            "ErrorTolerance": 0.1
+        },
+        "MeshSize": {"Edges": {}, "Faces": {}, "Volumes": {}},
+        "SaveMesh": false,
+        "SavePath": "mesh"
+    },
+    "EQSMode": false,
+    "FEMOrder": 2,
+    "StimulationSignal": {
+        "Type": "Rectangle",
+        "ListOfFrequencies": [10000],
+        "Frequency[Hz]": 130.0,
+        "PulseWidth[us]": 30.0,
+        "PulseTopWidth[us]": 0.0,
+        "CounterPulseWidth[us]": 0.0,
+        "InterPulseWidth[us]": 0.0,
+        "SpectrumMode": "OctaveBand",
+        "CounterAmplitude": 1.0,
+        "CutoffFrequency": 250000.0,
+        "CurrentControlled": true
+    },
+    "Solver": {
+        "Type": "CG",
+        "Preconditioner": "local",
+        "PreconditionerKwargs": {},
+        "MaximumSteps": 5000,
+        "Precision": 1e-10
+    },
+    "PointModel": {
+        "Pathway": {
+            "Active": true,
+            "FileName": "../../PAM/Allocated_axons.h5",
+            "ExportField": false
+        },
+        "Lattice": {
+            "Active": false,
+            "Center": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
+            "Shape": {"x": 10, "y": 10, "z": 10},
+            "Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
+            "PointDistance[mm]": 0.1,
+            "CollapseVTA": false,
+            "ExportField": false
+        },
+        "VoxelLattice": {
+            "Active": false,
+            "Shape": {"x": 10, "y": 10, "z": 10},
+            "TimeDomain": false,
+            "ExportField": false
+        }
+    },
+    "OutputPath": "Results_rh",
+    "ComputeImpedance": false,
+    "ExportVTK": true,
+    "ExportFrequency": null,
+    "ExportElectrode": false,
+    "ModelSide": 0,
+    "CalcAxonActivation": true,
+    "ActivationThresholdVTA[V-per-m]": 200.0,
+    "FailFlag": "rh",
+    "OutOfCore": false,
+    "PathwayFile": "../../PAM/Allocated_axons_parameters.json",
+    "StimSets": {"Active": false, "StimSetsFile": null}
+}