### Flux recovery in H(div) conforming finite element space

29
views
0
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)
else:
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,

Dustin
Community: FEniCS Project

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