Compute length/area of equipotential lines/surfaces

variational formulation
1

Hello,

I would like to know if there is a way to extract/compute the lengths of equipotential lines (or areas of equipotential surfaces in 3D)?

I am solving Poisson’s equation with Dirichlet BC. The solution is then restricted to a conductor of arbitrary shape, and I would like to extract the length of the equipotential lines inside this conductor. A finite number of lines could be explicitly computed and the rest interpolated between the previously computed values.

(I know that in practice the conductivity should be accounted for in the PDE for conduction problems but let’s assume here that it is =1 everywhere).

Here is a minimal working example of the problem:

from dolfin import *
from dolfin_adjoint import *

from matplotlib import pyplot as plt

mesh = UnitSquareMesh(50,50)
U = FunctionSpace(mesh,'CG', 2)

class top_left(SubDomain):
    def inside(self,x,on_boundary):
        return (near(x[1], 1.0, DOLFIN_EPS) and (x[0]<=0.4+DOLFIN_EPS) and on_boundary)
                
class top_right(SubDomain):
    def inside(self,x,on_boundary):
        return (near(x[1], 1.0, DOLFIN_EPS) and (x[0]>=0.6+DOLFIN_EPS) and on_boundary)

bcs = []
bc1 = DirichletBC(U,2.0,top_left())
bc2 = DirichletBC(U,1.0,top_right())

bcs.append(bc1)
bcs.append(bc2)

# Poisson's equation
f = interpolate(Constant(1.0e-5), U)  # the source term for the PDE for elec
u = Function(U, name="Potential")
v = TestFunction(U)

F = inner(grad(v),grad(u))*dx - f*v*dx
solve(F == 0, u, bcs, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7, "maximum_iterations": 20}})

# Conductor distribution (expression not known beforehand)
expr = Expression('(pow(x[0]-0.5,2) + pow(x[1]-1,2) > pow(0.25,2)) && (pow(x[0]-0.5,2) + pow(x[1]-1,2) <= pow(0.5,2)) ? 1.0 : 0.0',degree=3)
conductor = Function(U, name="Conductor")
conductor.assign(project(expr,U))

# Function of interest (compute the length of equipotential lines)
v_c = u*conductor

# Plots for visualization
pl, ax = plt.subplots();
fig = plt.gcf()
fig.set_size_inches(16, 4)
plt.subplot(1, 2, 1); p = plot(u, title='potential everywhere'); p.set_cmap("coolwarm"); cbar = plt.colorbar(p);
plt.subplot(1, 2, 2); p = plot(v_c, title='potential in conductor',mode='color',vmin=0.0); p.set_cmap("magma"); cbar = plt.colorbar(p);

image

Is there a way to compute these lengths/areas with a PDE or another operation?

(If this operation can easily be accounted for with dolfin_adjoint to compute its gradient afterwards, this would be perfect!)

Thank you for your time

2

Here’s a possible direction (with some potential problems that I’ll point out): If you note that the normal to the surface u = C is \mathbf{n} = \nabla u / \vert\nabla u\vert, then you can formally apply the divergence theorem to say that

\vert\partial\Omega\vert = \int_{\partial\Omega}1\,d\mathbf{s}= \int_{\partial\Omega} \left(\frac{\nabla u}{\vert\nabla u\vert}\right)\cdot\mathbf{n}\,d\mathbf{s} = \int_\Omega \nabla\cdot\left(\frac{\nabla u}{\vert\nabla u \vert}\right)\,d\mathbf{x}\text{ ,}

where \Omega is the set of points on which u < C (integration over which can be implemented by multiplying with a characteristic function defined with conditional), so that \nabla u points outward.

The problems:

  • \vert\nabla u\vert may be small, which could cause precision issues. One workaround might be to replace it with \sqrt{\nabla u\cdot\nabla u + \epsilon} for some small \epsilon.
  • The integrand involves second derivatives of u, and directly calling assemble with it ignores distributional contributions from interior facets in a C^0 finite element solution. Expanding \nabla\cdot(\nabla u/\vert\nabla u\vert), you get a Laplacian term that can be replaced with -f, but there are still other second derivatives in \nabla(1/\vert\nabla u\vert). It may be okay to just naively assemble the integrand over element interiors if you’re using higher-order elements (e.g., CG2, as in the example), but it seems kind of iffy to me.
2 Likes
4

Thank you very much for your answer @kamensky.
However, I am encountering some issues.

  • First, the proposed approach would not lead to desired results if I want to compute the length of the equipotential lines when the field is restricted to a conductor (v_c as shown in the right picture of my original post). When building the \Omega set of points, the resulting contour does not correspond to the equipotential line.
    –
    image
    –

  • Then, if the field corresponding to the potential everywhere is studied (u), the actual value doesn’t seem to be correct. I choose to compute the length of the equipotential line for potential = 1.5. As the domain and boundary conditions are placed with a certain symmetry, this equipotential line corresponds to a straight vertical line in the middle of the domain of length l=1.0.
    –
    image
    –
    However, when computing the length with the proposed approach it is actually returning l=0.59 … Do you think this could come from the second problem you mentioned?
    Here is the completed code:

from dolfin import *
from dolfin_adjoint import *

from matplotlib import pyplot as plt
from ufl import sqrt

mesh = UnitSquareMesh(50,50)
U = FunctionSpace(mesh,'CG', 2)

class top_left(SubDomain):
    def inside(self,x,on_boundary):
        return (near(x[1], 1.0, DOLFIN_EPS) and (x[0]<=0.4+DOLFIN_EPS) and on_boundary)
                
class top_right(SubDomain):
    def inside(self,x,on_boundary):
        return (near(x[1], 1.0, DOLFIN_EPS) and (x[0]>=0.6+DOLFIN_EPS) and on_boundary)

bcs = []
bc1 = DirichletBC(U,2.0,top_left())
bc2 = DirichletBC(U,1.0,top_right())

bcs.append(bc1)
bcs.append(bc2)

# Poisson's equation
f = interpolate(Constant(1.0e-5), U)  # the source term for the PDE for elec
u = Function(U, name="Potential")
v = TestFunction(U)

F = inner(grad(v),grad(u))*dx - f*v*dx
solve(F == 0, u, bcs, solver_parameters={"newton_solver": {"absolute_tolerance": 1.0e-7, "maximum_iterations": 20}})

# Conductor distribution (expression not known beforehand)
expr = Expression('(pow(x[0]-0.5,2) + pow(x[1]-1,2) > pow(0.25,2)) && (pow(x[0]-0.5,2) + pow(x[1]-1,2) <= pow(0.5,2)) ? 1.0 : 0.0',degree=3)
conductor = Function(U, name="Conductor")
conductor.assign(project(expr,U))

# Function of interest (compute the length of equipotential lines)
v_c = u*conductor

# Plots for visualization
pl, ax = plt.subplots();
fig = plt.gcf()
fig.set_size_inches(16, 4)
plt.subplot(1, 2, 1); p = plot(u, title='potential everywhere'); p.set_cmap("coolwarm"); cbar = plt.colorbar(p);
plt.subplot(1, 2, 2); p = plot(v_c, title='potential in conductor',mode='color',vmin=0.0); p.set_cmap("magma"); cbar = plt.colorbar(p);

# -------------------- Compute the length of th equipotential/iso - line
isovalue = 1.5
epsilon = 0.001

u_c = conditional(lt(u,isovalue),u,0.0)

length_isoline = assemble(div( grad(u_c) /sqrt(dot(grad(u_c),grad(u_c))+epsilon) ) *dx)
print("The length of the isoline is l=",length_isoline)
5

So, the method I described really only makes sense as-is for closed isocontours. In the case of an isocontour that terminates on the boundaries of the mesh, like your second example, you can simply subtract off the contribution from the mesh boundary, which seems to be pretty accurate:

#length_isoline = assemble(div( grad(u_c) /sqrt(dot(grad(u_c),grad(u_c))+epsilon) ) *dx)
n = FacetNormal(mesh)
n_u = grad(u_c)/sqrt(dot(grad(u_c),grad(u_c))+DOLFIN_EPS)
length_isoline = assemble(div(n_u)*dx - dot(n,n_u)*ds)

I’m not sure what would be the best way to eliminate contributions from subsets of an isocontour, though.

1 Like
6

Thank you very much @kamensky!
You’re right I forgot to subtract the contribution coming from the mesh boundary. It now works for the case of closed isocontours.

I’ll try to figure something out for the other cases and post it here if conclusive.