This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

0 votes

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 ?
asked Jun 1, 2017 by ALimkilde FEniCS Novice (120 points)

1 Answer

0 votes

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

answered Jul 13, 2017 by lhdamiani FEniCS User (2,580 points)

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

...