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[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.

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.