Mesh refinement gives solution with zero values

mesh
1

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

2

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

1 Like
3

Yes , will do that. Thanks.

4

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.

5

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.

6

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.
However, if you use the adapt function to redefine your markers:

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):

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
1 Like
7

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