# Electromagnetism FEM (FEniCS) interpolation - leakage effect

As for the background of what is going on:

• I'm using FEniCS that is dedicated FEM solver
• The problem I'm solving is magnetostatic problem where the governing PDE is $$\bf{\nabla} \times \frac{1}{\mu} \bf{\nabla} \times \textbf{A} = \textbf{J}$$ since the domain is 2D it simplifies to Poisson's equation: $$\frac{1}{\mu} \bf{\nabla}^{2} \it{A_z} = J_z$$. The domain contains subdomains of different magnetic probabilities:

1. Air and copper wires (with current sources $$J_z$$) - relative permeability of 1
2. Ferromagnetic core - relative permeability of 10e3

Also, there is boundary condition: $$A_z = 0$$ at the outer edges of air subdomain encapsulating everything.

The computation was a success and I've obtained the distribution of magnetic vector potential $$A_z$$. Next, I'm evaluating magnetic flux distribution as: $$\bf{B} = \bf{\nabla} \times \bf{A}$$.

Now, if we try to visualize the solution as surface, which requires obviously interpolation inside the cells we obtain the following:

I’ve encountered this issue when visualizing data. I know that solutions are associated with nodes. There is some non-zero solution on nodes that are placed on the edge, the boundary between low and high permeability domains. For example, the very next node has a value 0 (zero), but the values are being interpolated inside the cell. We know, that magnetic flux doesn't behave this way and the transition between these domains should be steep.

The solution at nodes seems to be right, but the problem occurs when you try to interpolate data:

• for visualization purpose
• when reading data at any point and projecting solution to another mesh (when simulating multi physics)

How to address this problem? Do I need some additional boundary condition at the edge of ferromagnetic material (or any material for that matter)? I can't assume that the permeability will always be so high, sometimes it might be $$\mu_r = 10$$.

Edit: Here's the Python code

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import mshr
import copy
from scipy import constants
#############################################################################
import time
#############################################################################
from ProcessSubDomains import *
#############################################################################
mesh = Mesh('Mesh.xml')
# ---------------------------------------------------------------------------
markSubdomains('subdomains.xml', subdom_filename) # custom function that marks each subdomain with a number
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
File(subdom_filename+'.xml') >> materials

#############################################################################
# MATERIAL PROPERTIES
class MaterialProperty(UserExpression):
def __init__(self, materials, property, material_subdomains, property_index, **kwargs):
super(MaterialProperty, self).__init__(**kwargs)
self.materials = materials
self.property = property
self.material_subdomains = material_subdomains
self.property_index = property_index

def eval_cell(self, values, x, cell):
label = self.materials[cell.index]
for key in self.material_subdomains:
if label in self.material_subdomains[key]:
values[0] = self.property[ self.property_index[key] ]

print("Assigning material properties...")
material_subdomains = {'air': [0, 2], 'iron': [1], 'wire': [3, 4, 5, 6]} # numbers associeted with subdomains
property_index = {'air': 0, 'iron': 1, 'wire': 2} # property index in permeability array
permeability = constants.mu_0*np.array([1, 35e3, 0.9991])
mu = MaterialProperty(materials, permeability, material_subdomains, property_index, degree=0)

#############################################################################
# DEFINE FUNCTION SPACE
V = FunctionSpace(mesh, 'P', 1)

#############################################################################
# BOUNDRY CONDITIONS
tol = 1e-14
def boundry(x, on_boundry):
return on_boundry and ( near(abs(x[0]), 0.1, tol) or near(abs(x[1]), 0.1, tol) ) # checking when on Dirichlet BC
u_D = Constant(0) # value on Dirichlet BC
bc = DirichletBC(V, u_D, boundry)

#############################################################################
# REDEFINE INTEGRATION MEASURES
dx = Measure('dx', domain = mesh, subdomain_data = materials)

#############################################################################
# SOURCES
I = 400000.0
J_A = Constant(I)
J_B = Constant(I)

#############################################################################
# TRIAL AND TEST FUNCTIONS
A_z = TrialFunction(V)
v = TestFunction(V)

#############################################################################
# SOLVE VARIATIONAL PROBLEM
print("Solving variational form...")
L = J_A*v*dx(3) - J_A*v*dx(4) + J_B*v*dx(5) - J_B*v*dx(6)

A_z = Function(V)
solve( a == L, A_z, bc)

#############################################################################
# POSTPROCESSING
W = VectorFunctionSpace(mesh, 'P', 1)

Bx = A_z.dx(1)
By = -A_z.dx(0)
B = project( as_vector(( Bx, By )), W )

B_abs = np.power( Bx**2 + By**2, 0.5 )

plot(A_z)
plt.show()

plot(B, mode='glyphs')
plt.show()

plot(B_abs)
plt.show()



After finding solution for scalar Poisson's equation - values $$A_z$$ it would be the best to create corresponding vector field $$\textbf{A} = [0, 0, A_z]$$, but the following operation fails:

A_vec = project( as_vector([0, 0, A_z]), W )


and then perform:

B_vec = project(curl(A_vec), W)


But such operation doesn't seem possible due to 2D nature of geometry (probably).

EDIT (SOLUTION):

So I've found the solution! It's here: FEniCS course This requires somewhat solving the PDE once again (in reality it's special type of mapping). Let's say you've solved the equation and have the distribution of $$A_z$$.

The new formulation that has to be solved in Discontinouos Galerkin function space is:

$$\int_{\Omega} w_h v \mathrm{d}x = \int_{\Omega} f v \mathrm{d}x$$

where $$f$$ is your previous solution that you want to map, i.e. $$f = |\textbf{B}|$$

# POSTPROCESSING
# calculate derivatives
Bx = A_z.dx(1)
By = -A_z.dx(0)

B_abs = np.power( Bx**2 + By**2, 0.5 ) # compute length of vector

# define new function space as Discontinuous Galerkin
abs_B = FunctionSpace(mesh, 'DG', 0)
f = B_abs # obtained solution is "source" for solving another PDE

# make new weak formulation
w_h = TrialFunction(abs_B)
v = TestFunction(abs_B)

a = w_h*v*dx
L = f*v*dx

w_h = Function(abs_B)
solve(a == L, w_h)

# plot the solution
plot(w_h)
plt.show()


@AntonMenshov: The interpolation is done in ParaView on its own as well as interpolation by reading values at any point in script.

Yes, the visualization is inside the commercial software (Ansys Maxwell/Electronics). Unfortunately the documentation doesn't state anything about how it's done except for which PDE is solved. The second issue with it, is that it's basically impossible to export the data (mesh and solution associated with nodes) outside the software.

The mesh is saved as .msh file (Ansys extension) and solution as binary file. When looking into .msh file it seems that the program is using 2nd order elements instead of 1st order which is some sort of clue why this behaves in such a way. Maybe it could be done by reading position of each node, reconstructing mesh and reading values on those nodes. But this will requires some additional labor, a lot of it.

Another hint from the software is that it sets default Neumann boundary condition between every subdomain (if not stated otherwise), so that normal component of $$\bf{B}$$ and tangent component of $$\bf{H}$$ are continuous at the boundaries.

@Paul: Yes, I'm using here the Continuous Galerkin elements (Lagrange). I thought about adding additional layer of elements close to the boundary where the values will reach small values quickly which will reduce this glowing in post processing. I haven't tried it yet with Discontinuous Galerkin elements, which might be an option but requires tweaking the variational formulation.

When trying to interpolate values along some line for solution $$A_z$$ there is no problem in python/FEniCS. But when I'm trying to interpolate even the component of $$\bf{B}$$ which is calculated as derivative of course i.e.

B_x = A_z.dx(1)


And when tring to plot the component i get:

Traceback (most recent call last):
...
some lines
...
/anaconda3/lib/python3.6/site-packages/dolfin/function/function.py", line 257, in ufl_evaluate
assert derivatives == ()   # TODO: Handle derivatives


So, as for the FEniCS version 2018.1.0 such operation isn't possible. This also affects the results when visualizing it in ParaView or build-in FEniCS plot() function.

But there's a way out, at least for the visualization part. The solution that is correct is distribution of magnetic vector potential $$\bf{A}$$. We can then calculate $$\nabla \times \bf{A}$$ to get $$\bf{B}$$ which is possible thanks to calculator in ParaView.

The interpolation inside the python script still lacks this functionality. Maybe it'll be possible by adding custom code to perform something similar to what is shown above.

• first question: it seems that you are using ParaView to visualize the data from FEniCS. Does the interpolation happen inside ParaView or outside? Second question: it seems like the visualization from the commercial tool happens inside the commercial tool. Did you try exporting the "commercial" solution and visualizing it in ParaView using exactly the same steps? – Anton Menshov Jun 8 '19 at 19:56
• This behavior seems normal to me if you’re using continuous galerkin elements such as lagrange. If you really want to model sharp transitions close to that boundary, you’ll either need a lot of mesh refinement, higher order elements, or discontinuous galerkin elements. – Paul Jun 8 '19 at 19:58
• The solution $A_z$ should be continuous, but the quantity $\mathbf{B}$ is discontinuous across edges. I don't know anything about FEniCs, but if e.g. $A_z$ uses linear rooftop basis functions then $\mathbf{B}$ is constant within each triangle. That means that $\mathbf{B}$ is multivalued at each node, but you're plotting a single value at each node. How are you getting that single value for $\mathbf{B}$? – LedHead Jun 9 '19 at 11:57
• @AntonMenshov: The interpolation is done in ParaView on its own as well as interpolation by reading values at any point in script. Yes, the visualization is inside the commercial software (Ansys Maxwell/Electronics). Unfortunately the documentation doesn't state anything about how it's done except for which PDE is solved. The second issue with it, is that it's basically impossible to export the data (mesh and solution associated with nodes) outside the software. The mesh is saved as .msh file (Ansys extension) and solution as binary file. When looking into .msh file it seems that the progr – antagim Jun 9 '19 at 16:43
• By "multivalued" I did not mean that it is a vector quantity, of course it is, i meant that at each node $\mathbf{B}$ has a different value for each triangle that touches the node. You were plotting $|\mathbf{B}|$ at each node, which doesn't really make sense. – LedHead Jun 9 '19 at 18:10

Your problem is using project() in FEniCS. Here's the FEniCS documentation that discusses why you might want to use the project operator. Note that in that example, the exact flux is continuous, while the numerically computed one is discontinuous. That's not your situation. In your case, your $$\mathbf{B}$$ is actually discontinuous, but when you project it onto a continuous space, you incorrectly make it continuous.

What you want to do is to simply avoid all use of the project function. Your postprocessing step should be simplified to:

EDIT: OP says the code didn't work. Maybe for some reason FEniCS calculates Bx and By from A_z and averages the quantities at nodes, which would create "glowing" effect. I don't have FEniCS, so I'm really just guessing at how it works based on OP's discussions and code. Anyway, here's a second attempt that first projects the solution A_z onto a "DG" space that hopefully disassociates the solutions from nodes of the mesh, which should allow Bx and By to exhibit the proper discontinuities.

#############################################################################
# POSTPROCESSING

# project A_z onto a linear DG space
W = FunctionSpace(mesh, 'DG', 1)
A_z2 = project(A_z,W)

# now Bx and By should be associated with triangles instead of nodes
Bx = A_z2.dx(1)
By = -A_z2.dx(0)

B = as_vector(( Bx, By ))

B_abs = np.power( Bx**2 + By**2, 0.5 )

print("Plotting fields...")

File('results/A_z.xml') << A_z
File('results/A_z.pvd') << A_z
plot(A_z)
plt.show()

File('results/B_vec.pvd') << B
plot(B, mode='glyphs')
plt.show()

plot(B_abs)
plt.show()

plot(mesh)
plot(B_abs)
plt.show()

• This still creates "glowing" effect due to interpolation. The correct way is to project the solution to Discontinuous Galerkin elements as i've shown above. This way, magnetic field stays encapsulated inside high permeability region. – antagim Jun 10 '19 at 16:16
• @antagim There's no reason you should need to go through all those complications. If that works for you, then great, but it's far from the "correct" solution. (1) It's needlessly circular, because you're going from discontinuous to continous back to discontinuous, and (2) can you guarantee that this new discontinuous solution is actually the same as the original since you've gone through 2 projections??? Seems fishy. – LedHead Jun 10 '19 at 16:24
• @antagim I've modified the code as a complete guess at how to fix your problem the right-ish way. If you have time, see if it works for you. If not, I'll humbly accept defeat. Regardless, I hope you see my point about how all this business with projections really should be unnecessary. – LedHead Jun 10 '19 at 16:42

@LedHead: Yes, I've used project to calculate values of $$\bf{B}$$ which are shown above as solutions at nodes.

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
import mshr
import copy
from scipy import constants
#############################################################################
import time
#############################################################################
from ProcessSubDomains import *
#############################################################################
mesh = Mesh('Mesh.xml')
# ---------------------------------------------------------------------------
markSubdomains('subdomains.xml', subdom_filename) # custom function that marks each subdomain with a number
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
File(subdom_filename+'.xml') >> materials

#############################################################################
# MATERIAL PROPERTIES
class MaterialProperty(UserExpression):
def __init__(self, materials, property, material_subdomains, property_index, **kwargs):
super(MaterialProperty, self).__init__(**kwargs)
self.materials = materials
self.property = property
self.material_subdomains = material_subdomains
self.property_index = property_index

def eval_cell(self, values, x, cell):
label = self.materials[cell.index]
for key in self.material_subdomains:
if label in self.material_subdomains[key]:
values[0] = self.property[ self.property_index[key] ]

print("Assigning material properties...")
material_subdomains = {'air': [0, 2], 'iron': [1], 'wire': [3, 4, 5, 6]} # numbers associeted with subdomains
property_index = {'air': 0, 'iron': 1, 'wire': 2} # property index in permeability array
permeability = constants.mu_0*np.array([1, 35e3, 0.9991])
mu = MaterialProperty(materials, permeability, material_subdomains, property_index, degree=0)

#############################################################################
# DEFINE FUNCTION SPACE
V = FunctionSpace(mesh, 'P', 1)

#############################################################################
# BOUNDRY CONDITIONS
tol = 1e-14
def boundry(x, on_boundry):
return on_boundry and ( near(abs(x[0]), 0.1, tol) or near(abs(x[1]), 0.1, tol) ) # checking when on Dirichlet BC
u_D = Constant(0) # value on Dirichlet BC
bc = DirichletBC(V, u_D, boundry)

#############################################################################
# REDEFINE INTEGRATION MEASURES
dx = Measure('dx', domain = mesh, subdomain_data = materials)

#############################################################################
# SOURCES
I = 400000.0
J_A = Constant(I)
J_B = Constant(I)

#############################################################################
# TRIAL AND TEST FUNCTIONS
A_z = TrialFunction(V)
v = TestFunction(V)

#############################################################################
# SOLVE VARIATIONAL PROBLEM
print("Solving variational form...")
L = J_A*v*dx(3) - J_A*v*dx(4) + J_B*v*dx(5) - J_B*v*dx(6)

A_z = Function(V)
solve( a == L, A_z, bc)

#############################################################################
# POSTPROCESSING
W = VectorFunctionSpace(mesh, 'P', 1)

Bx = A_z.dx(1)
By = -A_z.dx(0)
B = project( as_vector(( Bx, By )), W )

B_abs = np.power( Bx**2 + By**2, 0.5 )

plot(A_z)
plt.show()

plot(B, mode='glyphs')
plt.show()

plot(B_abs)
plt.show()



After finding solution for scalar Poisson's equation - values $$A_z$$ it would be the best to create corresponding vector field $$\textbf{A} = [0, 0, A_z]$$, but the following operation fails:

A_vec = project( as_vector([0, 0, A_z]), W )


and then perform:

B_vec = project(curl(A_vec), W)


But such operation doesn't seem possible due to 2D nature of geometry (probably).

At this point, it's not a problem to extract solution at nodes, cell centers (using interpolation), normal vectors or whatever - it's more about the algorithm itself to calculate derivative of one field and ensure that it's constant across cell, as shown in derivative obtained in ParaView.

I've done something yo suggested from the tutorial:

V = A_z.function_space()
mesh = V.mesh()
degree = V.ufl_element().degree()

W = VectorFunctionSpace(mesh, 'P', degree)
B_vec = project( as_vector([Bx, By]), W)

plot(B_vec, title='flux field')
plt.show()


and as expected the vectors are associated with nodes.

But still plotting $$|\textbf{B}|$$ is the problem. :)

EDIT:

So I've found the solution! It's here: FEniCS course This requires somewhat solving the PDE once again (in reality it's special type of mapping). Let's say you've solved the equation and have the distribution of $$A_z$$.

The new formulation that has to be solved in Discontinouos Galerkin function space is:

$$\int_{\Omega} w_h v \mathrm{d}x = \int_{\Omega} f v \mathrm{d}x$$

where $$f$$ is your previous solution that you want to map, i.e. $$f = |\textbf{B}|$$

# POSTPROCESSING
W = VectorFunctionSpace(mesh, 'P', 1) # new function space for mapping B as vector

# calculate derivatives
Bx = A_z.dx(1)
By = -A_z.dx(0)

B = project( as_vector(( Bx, By )), W ) # project B as vector to new function space
B_abs = np.power( Bx**2 + By**2, 0.5 ) # compute length of vector

# plot B vectors
plot(B)
plt.show()

# define new function space as Discontinuous Galerkin
abs_B = FunctionSpace(mesh, 'DG', 0)
f = B_abs # obtained solution is "source" for solving another PDE

# make new weak formulation
w_h = TrialFunction(abs_B)
v = TestFunction(abs_B)

a = w_h*v*dx
L = f*v*dx

w_h = Function(abs_B)
solve(a == L, w_h)

# plot the solution
plot(w_h)
plt.show()


• You're doing a lot of work to end up back at the same place. You shouldn't be using project at all. Can't you just do B= as_vector((Bx,By)) and call it a day? – LedHead Jun 10 '19 at 15:45
• @LedHead: Sure, but now I can use the solution of $|\textbf{B}|$, or anything else, interpolate over the domain, read values at points and pass it further to compute something else in code interactively - whatever. Finally, i can simply visualize the solution and say - this is magnitude of $\textbf{B}$ vector in this domain! Without misunderstanding of audience. – antagim Jun 10 '19 at 15:52