Mesh refinement gives solution with zero values

mesh
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