### Solution discontinuous over subdomains (continuous solution needed)

53
views
0
28 days ago by

Hi everyone,

I'm working on a magnetic field simulation with the code attached below and I'm getting very strange results.

The 3D-geometry is divided into various subdomains, some of which have currents inducing the field. The solution should be a continuous function over the whole geometry, first for the vector potential A, then for its curl B. Unfortunately, that's not what I get. The PDE-solution A is computed as constant where $L\ne0$L0, and as zero everywhere else (where L = 0). So it has discontinuities on each of the subdomain boundaries. It looks as if the solution was computed for each subdomain separately, instead of connecting all the conditions to one field.

My guess would be that something is wrong with the integration measures, but I have no idea what. What's especially confusing is that the (probably) exact same approach worked fine on a 2D model. I'd really appreciate any help finding out what's wrong here. Thanks in advance.

P. S. Don't be surprised because I split up all the vectors - that's to get a better overview over the components while working and looking for the error. I've had the same problem with the vectorial formulation.

Community: FEniCS Project

Sorry, I posted before adding the code. This is it:

from __future__ import print_function
from fenics import *
from mshr import *
from math import sin, cos, pi

import matplotlib.pyplot as plt
import numpy as np

# some constants
mu_ri = 1.0e6
n = 20000

# coil currents in Ampere
I_left = 10000.0
I_right = 10000.0

# following current densities
J_left = I_left*n/(109.4*27.88)
J_right = I_right*n/(109.4*27.88)

mesh = Mesh("Feldgenerator.xml")
subdomains = MeshFunction("size_t", mesh, "Feldgenerator_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "Feldgenerator_facet_region.xml")

V = FunctionSpace(mesh, "CG", 2)

class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
b = bool(near(x[0], -800) or near(x[0], 1200) or near(x[1], -1200) or near(x[1], 1200) or near(x[2], -800) or near(x[2], 800))
return b

b0 = Constant(0.0)
dbc = DirichletBoundary()
bc = DirichletBC(V, b0, dbc)

dV = Measure("dx")(domain = mesh, subdomain_data = subdomains)

# define relative permeability for subdomains
class Permeability(Expression):
def __init__(self, subdomains, **kwargs):
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 1:
values[0] = mu_ri
else:
values[0] = 1.0

# initialise permeability
Mu = 4.0*pi*1.0e-7*Permeability(subdomains, degree = 0)

# set mu for whole mesh
mu = project(Mu, V)

# define current densities in x-direction for subdomains
class CurrentDensity_x(Expression):
def __init__(self, subdomains, **kwargs):
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] > 1 and self.subdomains[cell.index] < 18:
if self.subdomains[cell.index] < 10:
J_value = J_left
else:
J_value = J_right
if (self.subdomains[cell.index] == 2 or self.subdomains[cell.index] == 14):
values[0] = J_value
elif (self.subdomains[cell.index] == 3 or self.subdomains[cell.index] == 13):
values[0] = 1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 4 or self.subdomains[cell.index] == 12):
values[0] = 0
elif (self.subdomains[cell.index] == 5 or self.subdomains[cell.index] == 11):
values[0] = -1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 6 or self.subdomains[cell.index] == 10):
values[0] = -J_value
elif (self.subdomains[cell.index] == 7 or self.subdomains[cell.index] == 17):
values[0] = -1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 8 or self.subdomains[cell.index] == 16):
values[0] = 0
else:                                                                           # subdomains 9 and 15
values[0] = 1/sqrt(2)*J_value
else:
values[0] = 0

# define current densities in x-direction for subdomains
class CurrentDensity_z(Expression):
def __init__(self, subdomains, **kwargs):
self.subdomains = subdomains
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] > 1 and self.subdomains[cell.index] < 18:            # check if subdomain is coil part
if self.subdomains[cell.index] < 10:
J_value = J_left
else:
J_value = J_right
if (self.subdomains[cell.index] == 2 or self.subdomains[cell.index] == 14):
values[0] = 0
elif (self.subdomains[cell.index] == 3 or self.subdomains[cell.index] == 13):
values[0] = -1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 4 or self.subdomains[cell.index] == 12):
values[0] = -J_value
elif (self.subdomains[cell.index] == 5 or self.subdomains[cell.index] == 11):
values[0] = -1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 6 or self.subdomains[cell.index] == 10):
values[0] = 0
elif (self.subdomains[cell.index] == 7 or self.subdomains[cell.index] == 17):
values[0] = 1/sqrt(2)*J_value
elif (self.subdomains[cell.index] == 8 or self.subdomains[cell.index] == 16):
values[0] = J_value
else:                                                                           # subdomains 9 and 15
values[0] = 1/sqrt(2)*J_value
else:
values[0] = 0

# initialise current densities
JX = CurrentDensity_x(subdomains, degree = 1)
JZ = CurrentDensity_z(subdomains, degree = 1)

Jx = project(JX, V)
Jz = project(JZ, V)

# define variational problems for x, y, z
Ax = TrialFunction(V)
vx = TestFunction(V)
Ay = Function(V)
vy = TestFunction(V)
Az = TrialFunction(V)
vz = TestFunction(V)

# left hand sides of variational problem

# right hand sides of variational problems
Lx = mu*Jx*vx*dV
#Ly = 0
Lz = mu*Jz*vz*dV

# solve variational problems
Ax = Function(V)
solve(ax == Lx, Ax, bc)
#Ay = Function(Vy)
solve(ay == 0, Ay, bc)
Az = Function(V)
solve(az == Lz, Az, bc)

plot(Ax)
plt.show()
plot(Ay)
plt.show()
plot(Az)
plt.show()

# compute magnetic field (B = curl(A))
W = VectorFunctionSpace(mesh, "CG", 2)
A = project(as_vector((Ax, Ay, Az)), W)

B = project(1.0e6*curl(A), W)

Bx, By, Bz = split(B)

plot(Bx)
plt.show()
plot(By)
plt.show()
plot(Bz)
plt.show()

# print solutions on screen
for v in vertices(mesh):
x = v.point().x()
y = v.point().y()
z = v.point().z()
print(x, y, z, A(Point(x, y, z))[0], A(Point(x, y, z))[1], A(Point(x, y, z))[2])

# save solutions A and B for plotting in paraview
vtkfile_A = File("solutions/scalar-03/potential.pvd")
vtkfile_A << A
vtkfile_B =  File("solutions/scalar-03/field.pvd")
vtkfile_B << B
written 28 days ago by Rebecca Pflanze