Flux recovery in H(div) conforming finite element space

8 days ago by
Hi guys,

I'm currently having some trouble when trying to recover the flux in a H(div) conforming finite element space. As a benchmark problem I have adapted the magnetostatics tutorial. But when I'm visualizing the recovered magnetic flux in paraview it seems like that the recovered flux lies in continous Lagrange finite element space.
But I only want the normal component of the flux to be continous.

Here is a short working example to reproduce the error:
from __future__ import print_function
from fenics import *
from mshr import *
from math import sin, cos, pi

a = 1.0   # inner radius of iron cylinder
b = 1.2   # outer radius of iron cylinder
c_1 = 0.8 # radius for inner circle of copper wires
c_2 = 1.4 # radius for outer circle of copper wires
r = 0.1   # radius of copper wires
R = 5.0   # radius of domain
n = 10    # number of windings

# Define geometry for background
domain = Circle(Point(0, 0), R)

# Define geometry for iron cylinder
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)

# Define geometry for wires (N = North (up), S = South (down))
angles_N = [i*2*pi/n for i in range(n)]
angles_S = [(i + 0.5)*2*pi/n for i in range(n)]
wires_N = [Circle(Point(c_1*cos(v), c_1*sin(v)), r) for v in angles_N]
wires_S = [Circle(Point(c_2*cos(v), c_2*sin(v)), r) for v in angles_S]

# Set subdomain for iron cylinder
domain.set_subdomain(1, cylinder)

# Set subdomains for wires
for (i, wire) in enumerate(wires_N):
    domain.set_subdomain(2 + i, wire)
for (i, wire) in enumerate(wires_S):
    domain.set_subdomain(2 + n + i, wire)

# Create mesh
mesh = generate_mesh(domain, 100)

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

# Define boundary condition
bc = DirichletBC(V, Constant(0), 'on_boundary')

# Define subdomain markers and integration measure
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

# Define current densities
J_N = Constant(1.0)
J_S = Constant(-1.0)

# Define magnetic permeability
class Permeability(Expression):
    def __init__(self, markers, **kwargs):
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 0:
            values[0] = 4*pi*1e-7 # vacuum
        elif self.markers[cell.index] == 1:
            values[0] = 1e-5      # iron (should really be 6.3e-3)
            values[0] = 1.26e-6   # copper

mu = Permeability(markers, degree=1)

# Define variational problem
A_z = TrialFunction(V)
v = TestFunction(V)
a = (1 / mu)*dot(curl(A_z), curl(v))*dx
L_N = sum(J_N*v*dx(i) for i in range(2, 2 + n))
L_S = sum(J_S*v*dx(i) for i in range(2 + n, 2 + 2*n))
L = L_N + L_S

# Solve variational problem
A_z = Function(V)
solve(a == L, A_z, bc)

# Compute magnetic field (B = curl A)
W = FunctionSpace(mesh, 'BDM', 1)
B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)

# Save solution to file
vtkfile_B = File('magnetostatics/field.pvd')
vtkfile_B << B​

Does anybody have experience with recovering the flux in aH(div) conforming finite element space?

I will be grateful for any help.

Best regards,

Community: FEniCS Project

1 Answer

22 hours ago by
There is some support for visualising H(div) fields using: Paraview 5.5, FEniCS 2018.1.0 and the ``XDMFFile.write_checkpoint`` method.

There is nothing wrong with your field in FEniCS, it's just Paraview had very little support for visualising this type of field until very recently.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »