# Mesh refinement gives solution with zero values
**URL:** https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783
**Category:** mesh
**Created:** [May 11, 2021, 11:40pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783 "2021-05-11T23:40:47Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Manuel\_Pineda](https://fenicsproject.discourse.group/u/Manuel_Pineda)
#### Post date: [May 11, 2021, 11:40pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/1 "2021-05-11T23:40:47Z")
**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):
super(Permeability, self). __init__ (**kwargs) #agregado
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)
a = (1 / mu)*dot(grad(A_z), grad(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
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**
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [May 12, 2021, 5:50am UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/2 "2021-05-12T05:50:39Z")
Please format your code using 3 times `. Also please remove all lines of code not needed to reproduce the error.
---
### Author:  [@Manuel\_Pineda](https://fenicsproject.discourse.group/u/Manuel_Pineda)
#### Post date: [May 12, 2021, 3:50pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/3 "2021-05-12T15:50:11Z")
Yes , will do that. Thanks.
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [May 12, 2021, 4:09pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/4 "2021-05-12T16:09:35Z")
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.
---
### Author:  [@Manuel\_Pineda](https://fenicsproject.discourse.group/u/Manuel_Pineda)
#### Post date: [May 12, 2021, 4:18pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/5 "2021-05-12T16:18:41Z")
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.
---
### Author:  [@dokken](https://fenicsproject.discourse.group/u/dokken)
#### Post date: [May 12, 2021, 4:57pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/6 "2021-05-12T16:57:33Z")
I strongly doubt that, as you can clearly see from just printing the markers before and after refinement:
```auto
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
```auto
[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.
However, if you use the adapt function to redefine your markers:
```auto
markers = adapt(markers, mesh)
```
everything should work.
Full working code (where I have also changed the space you are projecting `B` into to something suitable):
```auto
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)
markers = adapt(markers, mesh)
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):
super(Permeability, self). __init__ (**kwargs) #agregado
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)
a = (1 / mu)*dot(grad(A_z), grad(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
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
```
---
### Author:  [@Manuel\_Pineda](https://fenicsproject.discourse.group/u/Manuel_Pineda)
#### Post date: [May 12, 2021, 5:19pm UTC](https://fenicsproject.discourse.group/t/mesh-refinement-gives-solution-with-zero-values/5783/7 "2021-05-12T17:19:07Z")
Thank you very much for the code example, worked great.
I’ll continue working and studying.
Cheers