Open mdovgialo opened 1 week ago
mdovgialo
commented
1 week ago Managed to load our 4 spheres model, without ground place, experimenting with boundary conditions and poin charges. Right now only boundary conditions seem to work, not sure how:
Making only the outside sphere have 0 potential results in zero potential everywhere, despite a point charge present...
Next: try adding ground plane like in FENICS solver
mdovgialo
commented
1 week ago MFEM seems to have predefined (but flexible) integrators for different classes of differential equations. It might be an issue to write a custom equation to solve in python
mdovgialo
commented
1 week ago Solved a bug which didn't notice yesterday. Now 3 spheres with a point in it can look like this:
COefficients all equall 0.33
coeeficients equal: [0.33, 1.65, 0.0165, 0.33])
and coefficients which are reciprocity of above:
Unity Point source on the side. Outside sphere has a boundary condition of being 0. Can't figure out how to add a grounding plane yet...
mdovgialo
commented
2 days ago Implemented solver, which seems to do the correct thing. Implemented an experimental voxel based resampling of a point cloud of the solution, might not be needed for equally sampled volumetric solutions in spirit of "MRI as a hex mesh" . For now allows to decently sample any mesh at any resolution, smoothed with boxcar kernel of a size of biggest mesh element
We put an electrode on top of the brain, just below the skull/CSF.
First, let's look at analytic solution of assuming
V_kcsd = 1.0 / (4 np.pi sigma_base * distance_to_electrode)
Which we define on vertices of vertices of the mesh we are using, and then sampled using voxel based resampling (kinda sanity check of a MFEM solution postprocessing methods):
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline
sigma = 0.33.
Lets compare it with MFEM simulations of a point charge at the electrode positions, using 3 meshes:
To simulate point charges a DiffusionIntegrator, which solves poisson equation was used. Physical Volumes had assigned coefficients:
conductivities_vector = mfem.Vector(list(1.0 / (4 np.pi np.array([0.33, 1.65, 0.0165, 0.33, 1e-10])))) (conductances, of brain matter, CSF, bone, skin, air (when model had air box))
Additionally a domain integrator was used to place a point charge using DeltaCoefficient.
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline of analytical solution
Contrast set to 1V so we can see the the effect of ground
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline of analytical solution
Contrast set to 1V so we can see the the effect of ground electrode
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline of analytical solution
Contrast set to 400V so we can see the the effect of ground electrode
We can see that without surrounding air, the potential inside the spheres is a lot higher
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline of analytical solution
Contrast set to 400V so we can see the the effect of ground electrode
The overall potential in the spheres is lower compared to electrode in the brain tissue, but still is way higher compared to free space analytical solution.
This image is a difference between FEM solution and free space solution. The mesh is 4 spheres without air, but with a forced infinte grounding plane touching the bottom of the spheres.
Let's try to have a more decent comparisson, by addign the free space solution to the leadfield:
Contrast is set that at 10 Volts it's exactly white, so we can follow 10 Volt isolines. Cross is put on an isoline of analytical solution. We can see that the whole sphere is way above 10 Volt line, except the part, where the sphere meets the grounding plate.
Contrast set at 400 volts, we can see the weird interpolation, slighly triangular in shape, which gets visible even better, if we just look at the difference between FEM solution and free space solution, calculated by the FENICS zoomed in
We can observe triangular, hexagonal artifacts with protrusions on connecting edges, which means, this image is sampled in a way higher resolution than the size of FEM elements (this was calculated using "fine" preset of KESI 4 spheres template)
As we can see We can recreate similar point source potential distribution, on way higher resolution meshes, than FENICS, with calculations taking a fraction of time. Mmesh with maximum element size of 0.005 meters (around 5 milion elements after refining) takes a dozen of seconds to simulate, and 30 to sample using KDTrees, and MFEM is not even compiled with CUDA or multithreading support.
There doesn't seem to be problems with the electrode position singularity. But lack of air volume in legacy KESI models, seems to elevate potential inside the sphere by orders of magnitude. Adding air box, makes the solution similar in shape and orders of magnitude to free space analytical (in the parts of isotropic, even material), and seem to have similar isolines at V=10, 20, but different at low voltages (1, 2 5). The isolines on the FEM solution are deformed by the mesh boundaries, which is what we wanted to see.
Lack of air volume might have introduced elevated potential, and changed the profiles
@danek8317 What do you think? I guess we'll discuss it on monday.
mdovgialo
commented
2 days ago comparison of profiles going vertically through 4 spheres through the electrode
Start by recreating sphere solution from tutorial