Difference in solutions using asseble_system compared to bc.apply() after refinement??

Question

I have encountered a strange thing when solving a variational problem after refinement of the functionspace.

I am solving a problem given by

 < r,v >_{H^1(\Omega)}  = < g,v >_{L^2(\Omega)}- \int_\Gamma f v ds

where r is a vectorfunction. I have a dirichlet BC on the left part of the rectangular domain, but only in the x-direction (V.sub(0)). I refine the function space, boundary conditions etc. before i solve the problem.

But trying both

K_h1,F = assemble_system(a_h1,R,bcs = bcs2)
solve(K_h1,r.vector(),F)

and

K_h1 = assemble(a_h1)
F = assemble(R)
for bc in bcs2:
    bc.apply(K_h1)
    bc.apply(F)
solve(K_h1,r.vector(),F)

I get two different results, and it seems like the Dirichlet BC are not satisfied when using assemble_system(...). It only happens after refinement of the mesh. Any ideas why?

I am using FEniCS 2016.2.0 with Anaconda. Here is a minimal example:

from dolfin import*
import matplotlib.pyplot as plt

set_log_level(30)

# Define problem mesh
L = 3; H = 1;
mesh = RectangleMesh(Point(0,0), Point(L, H), 40, 40,"left/right")

# Function spaces
V = VectorFunctionSpace(mesh,"CG",1)

# Subdomains
TopBoundary = CompiledSubDomain("on_boundary && near(x[1],H) ",H = H)
Left = CompiledSubDomain("near(x[0],0) && on_boundary")

# Point load approximated with delta function
f_bdry = Expression(("0.0", "-1/(ep*sqrt(pi))*exp(-(x[0]*x[0])/(ep*ep))"),H=H,ep=0.01,degree=8)

# Dirichlet boundary conditions
BC_left     = DirichletBC(V.sub(0),Constant(0.0),Left)
bcs         = [BC_left]

# Mark the TopBoundary
markers = FacetFunctionSizet(mesh, 1)
markers.set_all(0)
TopBoundary.mark(markers, 1);
ds = Measure('ds', domain=mesh, subdomain_data=markers)

# Refine
cell_markers = CellFunction("bool", V.mesh())
cell_markers.set_all(False)
ref_mesh = refine(V.mesh(),cell_markers)
V2 = adapt(V,ref_mesh)

bcs2 = []
for bc in bcs:
   bcnew = adapt(bc,V2.mesh(),V2)
   bcs2.append(bcnew)

m2 = FacetFunctionSizet(V2.mesh(), 1)
m2.set_all(0)
TopBoundary.mark(m2,1)
ds2 = Measure('ds', domain=V2.mesh(), subdomain_data=m2)

# Setup a problem
r_ = TrialFunction(V2)
v_ = TestFunction(V2)
r = Function(V2)

a_h1 = inner(r_,v_)*dx + inner(grad(r_),grad(v_))*dx

g = Expression(('x[0]','x[1]'),degree = 1)
R =  inner(g,v_)*dx(V2.mesh()) - inner(f_bdry,v_)*ds2(1) # Residual form

K_h1,F = assemble_system(a_h1,R,bcs = bcs2)
solve(K_h1,r.vector(),F)

plt.figure(1)
fig = plot(r,backend = 'matplotlib')
plt.colorbar(fig)
plt.savefig('r_ass_sys.png') # The BC's are not satisfied?

r2 = Function(V2)
K_h1 = assemble(a_h1)
F = assemble(R)
for bc in bcs2:
    bc.apply(K_h1)
    bc.apply(F)
solve(K_h1,r2.vector(),F)

plt.figure(2)
fig = plot(r2,backend = 'matplotlib')
plt.colorbar(fig)
plt.savefig('r_ass_and_apply.png') # The BC's are satisfied here!!

plt.figure(3)
fig = plot(abs(r-r2),backend = 'matplotlib')
plt.colorbar(fig)
plt.savefig('r_difference.png') # There is a difference ?

Answer 1

Hi there,
I think the problem is not with assemble_system() and bc.apply()...

I say that because if you change the way you define your bcs2 to sth like:

bcs2 = []
bcd = DirichletBC(V2, (1.,1.), "x[1]<DOLFIN_EPS")
bcs2.append(bcd)

and then work with this version of bcs2 the both ways of applying the boundaries will produce the same result..

So, my guess is that something is not 100% correct when refining your mesh and adapting bcs to bcs2....

Regards,
Leonardo Hax Damiani

Comments

Does this also work if we take

bcd = DirichletBC(V2.sub(0), (1.,1.), "x[1]<DOLFIN_EPS")

I only observed the different results, when applying the Dirichlet BC in only one direction ( .sub(0) )?

---

For a subdomain, you need a scalar boundary value... if you use:

  bcd = DirichletBC(V2.sub(0), (1.), "x[1]<DOLFIN_EPS")

It works and the results are the same as when applying to whole V2 as I suggested before.

Regards,
Leonardo