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:
- Air and copper wires (with current sources $J_z$) - relative permeability of 1
- 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...")
a = (1/mu)*inner( grad(A_z), grad(v) )*dx
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()
ADDITIONAL REMARKS:
@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.
.mshfile (Ansys extension) and solution as binary file. When looking into.mshfile it seems that the progr $\endgroup$