I have a simple electrical current problem I’ve been working on where I have two electrodes around a nerve, and I need to apply voltage to each electrode and measure the current produced.
I am getting the correct solution, but my surface integral to calculate current is not working as expected.
The actual solution is symmetric, which I’ve verified, but I get very different current calculations on each electrode. They should be nearly identical.
I tried manually looping through the faces on the exterior of each electrode and computing the integral myself and I get identical answers for the two electrodes, although the actual magnitude of the answer is not accurate.
Here’s the code to compute it both ways.
import dolfin as do
from fenics import *
import numpy as np
import math
# Import mesh and define function space
mesh = do.Mesh("mesh.xml");
markers = MeshFunction('size_t', mesh, 'mesh_physical_region.xml');
dim = mesh.topology().dim();
faces = MeshFunction('size_t', mesh, 'mesh_facet_region.xml');
V = FunctionSpace(mesh, 'P', 1);
W = VectorFunctionSpace(mesh, 'P', 1)
def surface(x, on_boundary):
return (math.sqrt(x[0]*x[0] + x[1]*x[1]) > 1.95) and on_boundary
u_L = Constant(-1000);
u_R = Constant(1000);
bc_L = DirichletBC(V, u_L, faces, 8)
bc_R = DirichletBC(V, u_R, faces, 9)
bc_outer = DirichletBC(V, 0, surface)
bcs = [bc_L, bc_R, bc_outer];
# Assign Material Properties
k_saline = 2
k_epineurium = 0.008
k_perineurium = 0.00336
k_fascicle = 0.5
k_cuff = 2E-10
k_platinum = 9E6
class Kappa(UserExpression):
def __init__(self, markers, **kwargs):
self.markers = markers
super().__init__(**kwargs)
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = k_fascicle
elif self.markers[cell.index] == 2:
values[0] = k_perineurium
elif self.markers[cell.index] == 3:
values[0] = k_epineurium
elif self.markers[cell.index] == 6:
values[0] = k_cuff
elif self.markers[cell.index] == 7:
values[0] = k_saline
else:
values[0] = k_platinum
kappa = Kappa(markers, degree=0)
# Define variational problem
dx = Measure('dx', domain=mesh, subdomain_data=markers)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = kappa * dot(grad(u), grad(v))*dx
L = f*v*dx
# Compute solution
u = Function(V)
solve(a == L, u, bcs, solver_parameters={"linear_solver": "gmres", "preconditioner": "sor"}, form_compiler_parameters={"optimize": True})
ds = Measure('dS', domain=mesh, subdomain_data=faces)
n = FacetNormal(mesh)
grad_u = project(grad(u), W, solver_type='gmres', preconditioner_type='sor')
flux_1 = assemble(inner(grad_u('+'), n('+'))*ds(8))
flux_2 = assemble(inner(grad_u('+'), n('+'))*ds(9))
print("Flux 1: %f" % flux_1)
print("Flux 2: %f" % flux_2)
ons = np.array([1, 1, 1])
volt_diff1 = np.array([])
volt_diff2 = np.array([])
mesh.init(dim - 1, dim)
mesh.init(dim - 3, dim - 1)
for f in facets(mesh):
#If on boundary between two domains
if 4 in markers.array()[f.entities(dim)] and 7 in markers.array()[f.entities(dim)]:
#Compute facet area of that face (couldn't figure out what to do with FacetArea(mesh))
x_vert = np.array([])
y_vert = np.array([])
z_vert = np.array([])
norm_vector = f.normal()
for vertices in f.entities(dim - 3):
temp_point = Vertex(mesh, vertices).point()
x_vert = np.append(x_vert, temp_point.x())
y_vert = np.append(y_vert, temp_point.y())
z_vert = np.append(z_vert, temp_point.z())
a1 = np.array([x_vert, y_vert, ons])
a2 = np.array([y_vert, z_vert, ons])
a3 = np.array([z_vert, x_vert, ons])
a1 = np.linalg.det(a1)
a2 = np.linalg.det(a2)
a3 = np.linalg.det(a3)
face_area = 0.5 * math.sqrt(a1*a1 + a2*a2 + a3*a3)
#Measure coordinates and solution in each of the two neighboring cells
for cell in f.entities(dim):
if markers[cell] == 4:
volt_inside = u(Cell(mesh, cell).midpoint().x(), Cell(mesh, cell).midpoint().y(), Cell(mesh, cell).midpoint().z())
point_inside = Cell(mesh, cell).midpoint()
else:
volt_outside = u(Cell(mesh, cell).midpoint().x(), Cell(mesh, cell).midpoint().y(), Cell(mesh, cell).midpoint().z())
point_outside = Cell(mesh, cell).midpoint()
#Calculate derivate * surface area
diff_point = point_outside - point_inside
volt_diff1 = np.append(volt_diff1, ((volt_outside - volt_inside)) / abs(norm_vector.dot(diff_point)) * face_area)
elif 5 in markers.array()[f.entities(dim)] and 7 in markers.array()[f.entities(dim)]:
x_vert = np.array([])
y_vert = np.array([])
z_vert = np.array([])
norm_vector = f.normal()
for vertices in f.entities(dim - 3):
temp_point = Vertex(mesh, vertices).point()
x_vert = np.append(x_vert, temp_point.x())
y_vert = np.append(y_vert, temp_point.y())
z_vert = np.append(z_vert, temp_point.z())
a1 = np.array([x_vert, y_vert, ons])
a2 = np.array([y_vert, z_vert, ons])
a3 = np.array([z_vert, x_vert, ons])
a1 = np.linalg.det(a1)
a2 = np.linalg.det(a2)
a3 = np.linalg.det(a3)
face_area = 0.5 * math.sqrt(a1*a1 + a2*a2 + a3*a3)
for cell in f.entities(dim):
if markers[cell] == 5:
volt_inside = u(Cell(mesh, cell).midpoint().x(), Cell(mesh, cell).midpoint().y(), Cell(mesh, cell).midpoint().z())
point_inside = Cell(mesh, cell).midpoint()
else:
volt_outside = u(Cell(mesh, cell).midpoint().x(), Cell(mesh, cell).midpoint().y(), Cell(mesh, cell).midpoint().z())
point_outside = Cell(mesh, cell).midpoint()
diff_point = point_outside - point_inside
volt_diff2 = np.append(volt_diff2, ((volt_outside - volt_inside)) / abs(norm_vector.dot(diff_point)) * face_area)
print("Flux 1: %f" % np.sum(volt_diff1))
print("Flux 2: %f" % np.sum(volt_diff2))
I’m very suspicious of both the (’+’) that is needed for the surface integral, and the FacetNormal values.
In a few iterations, my manual calculation was asymmetric until I took the absolute value of the face normal vector, which makes me think that the direction the FacetNormal points is not always consistent. I would expect all of the normal vectors to point either inside the electrode or outside, but if they’re randomly mixed the solution is wrong.
I also do not understand how the (’+’) and (’-’) side of each face are determined and if they have a similar issue where the positive side is not always in the same subdomain.
If anyone has any ideas, it is very appreciated.