# Mesh refinement gives solution with zero values

Hi everyone!
I’m trying to refine one specific area of a mesh in the magnetostatics example from the tutorial. And when plotting I got zero values for the results and encountered with this error:

home/manuel/anaconda3/envs/fenics2018/lib/python3.7/site-packages/matplotlib/quiver.py:686: RuntimeWarning: divide by zero encountered in double_scalars
length = a * (widthu_per_lenu / (self.scale * self.width))
/home/manuel/anaconda3/envs/fenics2018/lib/python3.7/site-packages/matplotlib/quiver.py:686: RuntimeWarning: invalid value encountered in multiply
length = a * (widthu_per_lenu / (self.scale * self.width))
/home/manuel/anaconda3/envs/fenics2018/lib/python3.7/site-packages/matplotlib/quiver.py:738: RuntimeWarning: invalid value encountered in less
short = np.repeat(length < minsh, 8, axis=1)
/home/manuel/anaconda3/envs/fenics2018/lib/python3.7/site-packages/matplotlib/quiver.py:751: RuntimeWarning: invalid value encountered in less
tooshort = length < self.minlength

the code is basically the magnetostatics example with a mesh refinement condition:

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

a = 1.0
b = 1.2
c_1 = 0.8
c_2 = 1.4
r = 0.1
R = 5.0
n = 10

domain = Circle(Point(0, 0), R)
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)
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]

domain.set_subdomain(1, cylinder)

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)

mesh = generate_mesh(domain, 32)

cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if  2 < p.x() < 3 and  2.5 < p.y() < 3.5 :
cell_markers[cell] = True

mesh = refine(mesh, cell_markers)

V = FunctionSpace(mesh, 'P', 1)

bc = DirichletBC(V, Constant(0), 'on_boundary')

markers = MeshFunction('size_t', mesh, 2, mesh.domains())

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

J_N = Constant(1.0)
J_S = Constant(-1.0)

class Permeability(UserExpression): #cambio
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*1
elif self.markers[cell.index] == 1:
values[0] = 1
else:
values[0] = 1

mu = Permeability(markers, degree=1)

A_z = TrialFunction(V)
v = TestFunction(V)
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

A_z = Function(V)
solve(a == L, A_z, bc)

W = VectorFunctionSpace(mesh, 'P', 1)
B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)

thanks for reading and any help will be appreciated

Please format your code using 3 times `. Also please remove all lines of code not needed to reproduce the error.

1 Like

Yes , will do that. Thanks.

To me the problem seems to be the permability, as all values are set to 1.
You can confirm this by removing the refine command and still getting a 0 solution.

I removed the refine command and got a vectors and values for B. But when activate the refine command, I got no vectors of B and when plotting in paraview the results of B are 0. Even if the iron permeability is set to 1500, still got 0 values of B when the refine command is active.

I strongly doubt that, as you can clearly see from just printing the markers before and after refinement:

from dolfin import *
from mshr import *
from math import sin, cos, pi

a = 1.0
b = 1.2
c_1 = 0.8
c_2 = 1.4
r = 0.1
R = 5.0
n = 10

domain = Circle(Point(0, 0), R)
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)
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]

domain.set_subdomain(1, cylinder)

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)

mesh = generate_mesh(domain, 32)

cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if  2 < p.x() < 3 and  2.5 < p.y() < 3.5 :
cell_markers[cell] = True

markers = MeshFunction('size_t', mesh, 2, mesh.domains())
print(markers.array(), min(markers.array()), max(markers.array()))
mesh = refine(mesh)
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
print(markers.array(), min(markers.array()), max(markers.array()))

which returns

[20 21 21 ...,  0  0  0] 0 21
[18446744073709551615 18446744073709551615 18446744073709551615 ...,
18446744073709551615 18446744073709551615 18446744073709551615] 18446744073709551615 18446744073709551615

This means that your mesh.domains() are not carried over when you refine the mesh.

everything should work.
Full working code (where I have also changed the space you are projecting B into to something suitable):

from dolfin import *
from mshr import *
from math import sin, cos, pi

a = 1.0
b = 1.2
c_1 = 0.8
c_2 = 1.4
r = 0.1
R = 5.0
n = 10

domain = Circle(Point(0, 0), R)
cylinder = Circle(Point(0, 0), b) - Circle(Point(0, 0), a)
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]

domain.set_subdomain(1, cylinder)

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)

mesh = generate_mesh(domain, 32)

cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
cell_markers.set_all(False)
for cell in cells(mesh):
p = cell.midpoint()
if  2 < p.x() < 3 and  2.5 < p.y() < 3.5 :
cell_markers[cell] = True

markers = MeshFunction('size_t', mesh, 2, mesh.domains())

# Refine mesh in given area
mesh = refine(mesh, cell_markers)

V = FunctionSpace(mesh, 'P', 1)
bc = DirichletBC(V, Constant(0), 'on_boundary')
dx = Measure('dx', domain=mesh, subdomain_data=markers)

J_N = Constant(1.0)
J_S = Constant(-1.0)

class Permeability(UserExpression): #cambio
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*1
elif self.markers[cell.index] == 1:
values[0] = 1
else:
values[0] = 1

mu = Permeability(markers, degree=1)
perm = project(mu, FunctionSpace(mesh, "DG", 0))

A_z = TrialFunction(V)
v = TestFunction(V)
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

A_z = Function(V)
solve(a == L, A_z, bc)

W = VectorFunctionSpace(mesh, 'DG', 0)
B = project(as_vector((A_z.dx(1), -A_z.dx(0))), W)
File("B.pvd") << B
1 Like

Thank you very much for the code example, worked great.
I’ll continue working and studying.
Cheers