 # Anisotropic material definition and results issue

Hi guys!
I’m trying to solve a very simple PDE with anisotropic material property: The PDE is fairly normal magnetostatic equation that solves for magnetic vector potential

\nabla \times \upsilon \nabla \times \bf{A} = \bf{J}

where the reluctivity matrix can be defined as inverse of permeability matrix:

\upsilon = \mu ^{-1} = \begin{bmatrix} \upsilon_{x} & 0 & 0 \\ 0 & \upsilon_{y} & 0 \\ 0 & 0 & \upsilon_{z} \end{bmatrix} = \begin{bmatrix} {1}/{\mu_{x}} & 0 & 0 \\ 0 & {1}/{\mu_{y}} & 0 \\ 0 & 0 & {1}/{\mu_{z}} \end{bmatrix}

As for electromagnetism, the higher the magnetic permeability (lower reluctivity) the higher the magnetic should be. In this case, where there is such high anisotropy, the flux should align itself along X axis.

But the thing is, I’m getting totally opposite effect where flux aligns itself along Y axis! I’m not sure where does this effect come from… I’ve simplified example to 2D problem.

Am I doing something wrong here or I do not understand how FEniCS operates on gradients, tensors, coordinates?

from dolfin import *
from mshr import *
from scipy import constants
import matplotlib.pyplot as plt

domain = Rectangle(Point(-10, -10), Point(10, 10))
mesh = generate_mesh(domain, 64)

plot(mesh)
plt.show()

# function space
V = FunctionSpace(mesh, 'P', 1)

# boundary conditions
walls = 'on_boundary && (near(abs(x), 10.0) || near(abs(x), 10.0))'
bc = DirichletBC(V, Constant(0.0), walls)

tol = 1e-14

# Wire
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x**2 + x**2 <= 4 + tol

# Space
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x**2 + x**2 > 4 + tol

materials = MeshFunction('size_t', mesh, mesh.topology().dim())

subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

plot(materials)
plt.show()

dx = Measure('dx', domain = mesh, subdomain_data = materials)

A_z = Function(V) # magnetic vector potential
v = TestFunction(V)

J = 5e6; # current source

# anisotropic material parameters
reluctivity = as_matrix( ( \
( 1/(constants.mu_0*10000), 0 ), \
( 0, 1/(constants.mu_0*1) ) \
) )
"""
reluctivity = 1/constants.mu_0
"""

solve(F == 0, A_z, bc)

plot(A_z)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

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

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

plot(B)
plt.xlabel('x')
plt.ylabel('y')
plt.show()

Might the problem be that the weak form in the code (definition of F) is using the gradient, rather than the curl? Note that, although the UFL curl operator is only defined for 3D vectors, you can conveniently define your own 2D scalar-to-vector curl as a Python function

def curl2D(v):
return as_vector((v.dx(1),-v.dx(0)))


Good thing you’ve pointed it out! You’re absolutely correct!
The curl doesn’t equal to gradient, components are being swapped (and even mixed) in curl.

So either way, i can swap my permeability matrix (which is messy and unreadable) or redefine this 2D curl as you’ve shown.

Thank you!