 # Asymmetric flux calculation with a symmetric solution

#1

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*x + x*x) > 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 = k_fascicle
elif self.markers[cell.index] == 2:
values = k_perineurium
elif self.markers[cell.index] == 3:
values = k_epineurium
elif self.markers[cell.index] == 6:
values = k_cuff
elif self.markers[cell.index] == 7:
values = k_saline
else:
values = 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.

1 Like
#2

Hi, this thread could shed some light on ‘+’, ‘-’ orientation and how to compute the surface integrals.

1 Like
#3

Thank you!
That was exactly what I needed.

Changing from

flux_1 = assemble(dot(grad_u, n)(’+’)*ds(1))
flux_2 = assemble(dot(grad_u, n)(’+’)*ds(2))

To

flux_1 = assemble(dot(grad_u, n)(’+’)*ds(1) + Constant(0)*dx(domain=mesh, subdomain_data=markers))
flux_2 = assemble(dot(grad_u, n)(’+’)*ds(2) + Constant(0)*dx(domain=mesh, subdomain_data=markers))

made everything symmetric again.