### How to accurately apply solution of one mesh to an adjacent mesh as a Dirichlet boundary condition?

264
views
4
6 months ago by

Hello all,
I have two adjacent meshes that that have coincident nodes along their shared boundary.
I've solved for the velocity solution on Vector_space_1 defined on mesh_1 with "CG" elements of degree 1.

u1, v1 = U.split(deepcopy = True)

The trouble arises in applying the velocity solution, v1, as a Dirichlet boundary condition to the shared boundary of the neighboring mesh, mesh_2.

bcu = DirichletBC(Vector_space_2, v1, facets_2, 20)

where Vector_space_2 is the vector function space defined on mesh_2 and 'facets_2, 20' identifies the appropriate facets. Vector_space_2 is also "CG" elements of degree 1.

I expect the solution on mesh_2 on the boundary nodes to be very close to v1. However the relative norm comparing v1 and the solved solution on mesh_2 along the shared boundary is on the order of 10^{-4} as opposed to the desired and expected 10^{-16}. Evidently the DirichletBC reads the right values and applies them to the correct nodes, but incurs some error in the process.

The slight disparity persists if I first project v1 onto the vector space of mesh_2

v2 = project(v1, V_space_2 , solver_type = "mumps",\
form_compiler_parameters = {"cpp_optimize" : True, "representation" : "uflacs"} )

Is there a robust way to take the solution of one mesh and ensure the values along a given boundary are prescribed exactly as values on the same boundary in a neighboring mesh?

Here's a minimum working example. The crux is at the end where the values on the boundary do not match between the two meshes. I think it's likely I'm misusing the project function. Ultimately, I need values in the 2nd function space just to be accurate on that boundary for a Dirichlet BC.

from fenics import *
from mshr import *

# Two adjacent meshes that have coincident nodes on shared bouandary:
box_1 = Rectangle(Point(0, 0), Point(0.1, 0.1))
box_2 = Rectangle(Point(0, 0.1), Point(0.1, 0.2))

# Generate mesh
mesh_1 = generate_mesh(box_1, 4)
mesh_2 = generate_mesh(box_2, 4)

# Set up a vector function space on each mesh
Vector_space_1 = VectorFunctionSpace(mesh_1, "CG", 1)
Vector_space_2 = VectorFunctionSpace(mesh_2, "CG", 1)

# find the indices that correspond to matching nodes on the shared boundary.
# Independent meshes are designed so that these nodes do match.

# degree of freedom coordinates
dofs_1 = Vector_space_1.tabulate_dof_coordinates().reshape((Vector_space_1.dim(),-1))
dofs_2 = Vector_space_2.tabulate_dof_coordinates().reshape((Vector_space_2.dim(),-1))

# Locate coordinates on shared boundary.
indices_1 = np.where((dofs_1[:,1] <= 0.1 + DOLFIN_EPS) & (dofs_1[:,1] >= 0.1 - DOLFIN_EPS))[0]
indices_2 = np.where((dofs_2[:,1] <= 0.1 + DOLFIN_EPS) & (dofs_2[:,1] >= 0.1 - DOLFIN_EPS))[0]

# Treat x and y components seperately
indices_1_x = indices_1[::2]
indices_1_y = indices_1[1::2]
indices_2_x = indices_2[::2]
indices_2_y = indices_2[1::2]

# reorder indices
indices_1_x = indices_1_x[dofs_1[indices_1_x][:,0].argsort()]
indices_1_y = indices_1_y[dofs_1[indices_1_y][:,0].argsort()]
indices_2_x = indices_2_x[dofs_2[indices_2_x][:,0].argsort()]
indices_2_y = indices_2_y[dofs_2[indices_2_y][:,0].argsort()]

# concatenate
indices_1 = np.concatenate((indices_1_x,indices_1_y),axis = 0)
indices_2 = np.concatenate((indices_2_x,indices_2_y),axis = 0)

# Given some solution on mesh1/ vector_function_space_1

v_expression = Expression(("sin(pi*x[0]) + sin(pi*x[1])","sin(pi*x[0]) + sin(pi*x[1])"), degree = 1)
v_1 = interpolate(v_expression, Vector_space_1)

# Project solution from mesh 1 onto mesh 2
v_2 = project(v_1, Vector_space_2 , solver_type = "mumps",\
form_compiler_parameters = {"cpp_optimize" : True, "representation" : "uflacs"} )

# Retrieve values from each mesh that lie on boundary.
#Indices go from left to right in x and then y. #dofs_1[indices_1]

v_1_boundary = v_1.vector()[indices_1]
v_2_boundary = v_2.vector()[indices_2]

# Compare values on boundary:
print "relative error boundary values v_1 and v_2", np.linalg.norm(v_1_boundary - v_2_boundary)/np.linalg.norm(v_1_boundary)


I would like the relative error calclated at the end to be very small, so that these values can be enforced exactly on the boundary.
Thanks for the help!

Community: FEniCS Project
Suggest you add a minimum working example to help us understand exactly what you are doing. How are you defining variables on two meshes? You might check
https://fenicsproject.org/qa/1605/2-unknowns-solved-in-2-different-subdomains/
for some info. It might be an issue with how you are implementing the two variables (e.g.: If you are applying a DirichletBC over the domain in order to eliminate one from the solution then continuity at the boundary might be messing with your BC along the facets...

I'm also not sure why you are using 'deepcopy=True'...
written 6 months ago by Mike Welland

Hi Mike, appreciate the help,

I've added a minimum working example. In essence I have two independent meshes, and the issue I have is communicating values on the boundary between them.

Fair query on 'deepcopy = True'. I don't believe this changes the outcome so I'll leave that alone for now.

Hope the code clarifies things. My method to print values on the boundary is a little cumbersome, so if there's a more elegant solution there that'd be great too!

written 6 months ago by felixnewberry
Thanks, what you wrote is a good idea... although I think you mean to project it on Vector_space_2

# Project solution from mesh 1 onto mesh 2
v_2 = project(v_1, Vector_space_2 , solver_type = "mumps",\
form_compiler_parameters = {"cpp_optimize" : True, "representation" : "uflacs"} )​

I don't have a good answer for you, but my understanding is that the project operator tries to solve (v_1-f_1) (where f_1 is a function  on Vector_space_2) over the full domain... which I don't think is what you want. Therefore I think your DirichletBC approach would be the way to go, but I'm not sure why you still see an error...
written 6 months ago by Mike Welland

You're right, that should read Vector_space_2. I've updated it accordingly.

Yes, I think that DirichletBC must project or do something similar under the hood. To take the values from one function space to another.

written 6 months ago by felixnewberry

3
6 months ago by
If I'm understanding the problem correctly, the following example seems to get the job done:

from fenics import *

# mshr generates non-regular meshes, so I've just used RectangleMesh()
# to make sure the nodes line up as expected.
mesh_1 = RectangleMesh(Point(0.0,0.0),Point(0.1,0.1),3,3)
mesh_2 = RectangleMesh(Point(0.0,0.1),Point(0.1,0.2),3,3)

# Set up a vector function space on each mesh
V1 = VectorFunctionSpace(mesh_1, "CG", 1)
V2 = VectorFunctionSpace(mesh_2, "CG", 1)

# Create some non-trivial function on mesh 1. Define expr with higher
# degree, so that, when checking u2 vs expr, we have some non-zero error
# (that converges as the mesh is refined). (This mainly serves to check that
# we're not just getting zero by integrating over a wrong/empty domain.)
expr = Expression(("cos(x[0])","exp(x[0])"),degree=3)
u1 = interpolate(expr,V1)

# allow extrapolation on that function
u1.set_allow_extrapolation(True)

# use u1 as a BC on mesh 2
def boundary2(x, on_boundary):
return near(x[1],0.1)
bc2 = DirichletBC(V2, u1, boundary2)

# apply BC on mesh 2
u2 = Function(V2)
u2.vector().zero()
bc2.apply(u2.vector())

# now check the accuracy:

x1 = Expression("x[1]",degree=1)

# Look at L2 error on the shared boundary

# this should be zero, since u1 should equal u2 exactly on the shared
# boundary
print("Error between u2 and u1 (should be machine zero): "
# (conditional isolates bottom boundary of mesh_2)
+str(assemble(conditional
(gt(0.1+DOLFIN_EPS,x1),
inner(u2-u1, u2-u1),
Constant(0.0))*ds(domain=mesh_2))))

# this should be nonzero, but converging with mesh refinement, since u1
# (and hence u2) is not exactly equal to the higher-degree expression expr
print("Error between u2 and expr (should be nonzero but convergent): "
+str(assemble(conditional
(gt(0.1+DOLFIN_EPS,x1),
inner(u2-expr, u2-expr),
Constant(0.0))*ds(domain=mesh_2))))
​

Thank you! I think this resolves the issue.

I'm hoping you could parse the conditional statement a little for me. My take is that the first conditional applies inner(u2-u1, u2-u1) when 0.1+DOLFIN_EPS > x1, and is set to zero elsewhere. I thought that gt(0.1+DOLFIN_EPS,x1) would exclude the boundary on 0.1.

Thanks again.

written 6 months ago by felixnewberry
Your interpretation of the conditional() sounds right, but, since DOLFIN_EPS is positive, the case x1==0.1 satisfies 0.1+DOLFIN_EPS > x1, and the integration should follow the first branch of the conditional on that side of the domain.  (I tried using a SubDomain and flagging facets first, actually, but somehow this didn't seem to reliably select only the facets on the desired boundary; maybe I just mixed up my definition of the SubDomain in some way.)
written 6 months ago by David Kamensky
Ah yes, that checks out. Thanks!
written 6 months ago by felixnewberry