Flux recovery in H(div) conforming finite element space
8 days ago by
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 = 4*pi*1e-7 # vacuum elif self.markers[cell.index] == 1: values = 1e-5 # iron (should really be 6.3e-3) else: values = 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.
Community: FEniCS Project
22 hours ago by
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.