Skip to content

Commit

Permalink
Hp refinement (#53)
Browse files Browse the repository at this point in the history
* Add hp refinement

* refactor AdaptiveMeshRefinement

* add examples

---------

Co-authored-by: Kenichi Maeda <[email protected]>
Co-authored-by: JPPayonk <[email protected]>
Co-authored-by: Julius Zimmermann <[email protected]>
  • Loading branch information
4 people authored Aug 30, 2024
1 parent 89a857b commit 640584c
Show file tree
Hide file tree
Showing 10 changed files with 303 additions and 74 deletions.
53 changes: 25 additions & 28 deletions examples/MeshAPI/BrainVerciseMesh.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,50 +5,47 @@
The mesh is generated and then
refined at one surface.
"""
import ossdbs
from ngsolve import Draw, Redraw

settings = \
{"Electrodes":
[{"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
"EncapsulationLayer":
{"Thickness[mm]": 0.0, # indicates that no encapsulation is modelled
},
import ossdbs

},
],
"MaterialDistribution":
{"MRIPath": "../BrainGeometryAPI/segmask.nii.gz"
},
"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":
settings = {
"BrainRegion": {
"Center": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"Dimension": {"x[mm]": 40.0, "y[mm]": 40.0, "z[mm]": 40.0},
"Shape": "Ellipsoid",
},
"Electrodes": [
{
"LoadMesh": False,
"SaveMesh": False
}
}
"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
"EncapsulationLayer": {"Thickness[mm]": 0.0},
},
],
"MaterialDistribution": {"MRIPath": "../../input_files/Butenko_segmask.nii.gz"},
"Mesh": {
"LoadMesh": False,
"SaveMesh": False,
},
"ExportElectrode": False,
}

mesh = ossdbs.generate_mesh(settings)
print(mesh.boundaries)
print(mesh.materials)
print(mesh.ngsolvemesh.ne)
print("Number of elements:", mesh.ngsolvemesh.ne)
Draw(mesh.ngsolvemesh)
# Clip the drawing to observe the refinement
input("Hit enter to refine mesh")
mesh.refine_at_boundaries(["E1C1"])
print(mesh.ngsolvemesh.ne)
print("Number of elements:", mesh.ngsolvemesh.ne)
Redraw()

input("Hit enter to refine mesh")
mesh.refine_at_boundaries(["E1C1"])
print(mesh.ngsolvemesh.ne)
print("Number of elements:", mesh.ngsolvemesh.ne)
Redraw()

input("Hit to curve mesh")
Expand Down
117 changes: 117 additions & 0 deletions examples/MeshAPI/BrainVerciseMesh_AdaptiveMeshRefinement.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
"""
Example of a full model geometry
with one Vercise electrode without
encapsulation layer.
The mesh is generated and adaptive
refinement is used.
"""
import logging

import ngsolve
from ngsolve import Draw

import ossdbs
from ossdbs.api import (
create_bounding_box,
generate_electrodes,
load_images,
prepare_dielectric_properties,
prepare_solver,
prepare_stimulation_signal,
prepare_volume_conductor_model,
run_volume_conductor_model,
set_contact_and_encapsulation_layer_properties,
)
from ossdbs.fem import ConductivityCF
from ossdbs.model_geometry import BrainGeometry, ModelGeometry
from ossdbs.utils.settings import Settings

ossdbs.set_logger(logging.INFO)

input_settings = {
"BrainRegion": {
"Center": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"Dimension": {"x[mm]": 40.0, "y[mm]": 40.0, "z[mm]": 40.0},
"Shape": "Ellipsoid",
},
"Electrodes": [
{
"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"EncapsulationLayer": {"Thickness[mm]": 0.0},
"Contacts": [
{
"Contact_ID": 1,
"Active": True,
"Voltage[V]": 1.0,
},
{
"Contact_ID": 3,
"Active": True,
"Voltage[V]": 0.0,
},
],
},
],
"MaterialDistribution": {"MRIPath": "../../input_files/Butenko_segmask.nii.gz"},
"Mesh": {
"AdaptiveMeshRefinement": {
"Active": True,
"MaxIterations": 2,
"ErrorTolerance": 0.1,
},
"HPRefinement": {
"Active": False,
"Levels": 2,
"Factor": 0.125,
},
},
"StimulationSignal": {
"Type": "Multisine",
"ListOfFrequencies": [1000.0],
},
"ComputeImpedance": True,
}

settings = Settings(input_settings).complete_settings()

mri_image, dti_image = load_images(settings)

electrodes = generate_electrodes(settings)

brain_region = create_bounding_box(settings["BrainRegion"])
brain_model = BrainGeometry(settings["BrainRegion"]["Shape"], brain_region)

geometry = ModelGeometry(brain_model, electrodes)

set_contact_and_encapsulation_layer_properties(settings, geometry)
dielectric_properties = prepare_dielectric_properties(settings)
materials = settings["MaterialDistribution"]["MRIMapping"]

conductivity = ConductivityCF(
mri_image,
brain_region,
dielectric_properties,
materials,
geometry.encapsulation_layers,
complex_data=settings["EQSMode"],
dti_image=dti_image,
)

mesh = ossdbs.generate_mesh(settings)

with ngsolve.TaskManager():
solver = prepare_solver(settings)
volume_conductor = prepare_volume_conductor_model(
settings, geometry, conductivity, solver
)
frequency_domain_signal = prepare_stimulation_signal(settings)

run_volume_conductor_model(settings, volume_conductor, frequency_domain_signal)

print("Number of elements before refinement:", mesh.ngsolvemesh.ne)
print("Number of elements after refinement:", volume_conductor.mesh.ngsolvemesh.ne)
print("Number degrees of freedom:", ngsolve.H1(volume_conductor.mesh.ngsolvemesh).ndof)
Draw(volume_conductor.mesh.ngsolvemesh)
25 changes: 8 additions & 17 deletions examples/MeshAPI/BrainVerciseMesh_SetMeshSizes.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,60 +12,51 @@
import ossdbs

settings = {
"BrainRegion": {
"Center": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"Dimension": {"x[mm]": 40.0, "y[mm]": 40.0, "z[mm]": 40.0},
"Shape": "Ellipsoid",
},
"Electrodes": [
{
"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
"TipPosition": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"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},
"MaxMeshSize": 0.1,
},
],
}
],
"MaterialDistribution": {"MRIPath": "../BrainGeometryAPI/segmask.nii.gz"},
"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",
},
"MaterialDistribution": {"MRIPath": "../../input_files/Butenko_segmask.nii.gz"},
"Mesh": {"LoadMesh": False, "SaveMesh": False},
"ExportElectrode": False,
}

mesh = ossdbs.generate_mesh(settings)
print(mesh.boundaries)
print(mesh.materials)
print(mesh.ngsolvemesh.ne)
print("Number of elements:", mesh.ngsolvemesh.ne)

Draw(mesh.ngsolvemesh)

Expand Down
20 changes: 7 additions & 13 deletions examples/MeshAPI/BrainVerciseMesh_SetMeshSizes_Encap.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,41 +14,35 @@
ossdbs.set_logger()

settings = {
"BrainRegion": {
"Center": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"Dimension": {"x[mm]": 40.0, "y[mm]": 40.0, "z[mm]": 40.0},
"Shape": "Ellipsoid",
},
"Electrodes": [
{
"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 0},
"TipPosition": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"EncapsulationLayer": {"Thickness[mm]": 1.0, "MaxMeshSize": 0.2},
"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,
},
],
},
],
"MaterialDistribution": {"MRIPath": "../BrainGeometryAPI/segmask.nii.gz"},
"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",
},
"MaterialDistribution": {"MRIPath": "../../input_files/Butenko_segmask.nii.gz"},
"Mesh": {"LoadMesh": False, "SaveMesh": False},
"ExportElectrode": False,
}
Expand Down
63 changes: 63 additions & 0 deletions examples/MeshAPI/BrainVerciseMesh_hpRefinement.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
"""
Example of a full model geometry
with one Vercise electrode without
encapsulation layer.
The mesh is generated and then
refined at one surface.
"""
from ngsolve import Draw

import ossdbs

settings = {
"BrainRegion": {
"Center": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"Dimension": {"x[mm]": 40.0, "y[mm]": 40.0, "z[mm]": 40.0},
"Shape": "Ellipsoid",
},
"Electrodes": [
{
"Name": "BostonScientificVercise",
"Rotation[Degrees]": 0,
"Direction": {"x[mm]": 0, "y[mm]": 0, "z[mm]": 1},
"TipPosition": {"x[mm]": -9.48, "y[mm]": 11.61, "z[mm]": 4.68},
"EncapsulationLayer": {"Thickness[mm]": 0.0},
"Contacts": [
{
"Contact_ID": 1,
"Active": True,
"Voltage[V]": 1.0,
},
{
"Contact_ID": 3,
"Active": True,
"Voltage[V]": 0.0,
},
],
},
],
"Mesh": {
"LoadMesh": False,
"SaveMesh": False,
"AdaptiveMeshRefinement": {
"Active": False,
"MaxIterations": 2,
"ErrorTolerance": 0.1,
},
"HPRefinement": {
"Active": True,
"Levels": 2,
"Factor": 0.125,
},
},
"ExportElectrode": False,
}

hp_mesh = ossdbs.generate_mesh(settings)
Draw(hp_mesh.ngsolvemesh)

settings["Mesh"]["HPRefinement"]["Active"] = False
mesh = ossdbs.generate_mesh(settings)

print("Number of elements before refinement:", mesh.ngsolvemesh.ne)
print("Number of elements after refinement:", hp_mesh.ngsolvemesh.ne)
Loading

0 comments on commit 640584c

Please sign in to comment.