No interface condition between SubMeshes in mixed-dimensional branch?

1

Hi @cdaversin,

thanks for your work on mixed-dimensional branch. Currently i’ m eploring it with some demos. What i found is in the solution with mixed function space defined with two submeshes, there seems no interface condition (continuity), but rather the natural boundary conditions for each part. What I tried is to modify the diffusion factor on one of those two part from the demo

from dolfin import *

#Create mesh and define function space
mesh = UnitSquareMesh(30, 30)

marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
    marker[c] = c.midpoint().y() < 0.5
    #marker[c] = c.midpoint().x() < 0.5

submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 0)

# Define Dirichlet boundary
def boundarySub1(x):
    return x[1] < DOLFIN_EPS
    #return x[0] < DOLFIN_EPS

def boundarySub2(x):
    return x[1] > 1.0 - DOLFIN_EPS
    #return x[0] > 1.0 - DOLFIN_EPS

#element2D = FiniteElement("Lagrange", triangle, 1)
W1 = FunctionSpace(submesh1, "Lagrange", 1)
W2 = FunctionSpace(submesh2, "Lagrange", 1)
# Define the mixed function space
V = MixedFunctionSpace( W1, W2 )

# Define boundary conditions
u0 = Constant(0.0)
# Subdomain 1
bc1 = DirichletBC(V.sub_space(0), u0, boundarySub1)
# Subdomain 2
bc2 = DirichletBC(V.sub_space(1), u0, boundarySub2)
bcs = [bc1,bc2]

# Define variational problem
# Use directly TrialFunction and TestFunction on the product space
# Subdomain 1
u1 = TrialFunction(W1)
v1 = TestFunction(W1)
# Subdomain 2
u2 = TrialFunction(W2)
v2 = TestFunction(W2)
# Mixed
(u1_m,u2_m) = TrialFunctions(V)
(v1_m,v2_m) = TestFunctions(V)

f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)

# Weak form - 1
a1 = 10*inner(grad(u1), grad(v1))*dx
L1 = f*v1*dx
u1 = Function(W1)
solve(a1 == L1, u1, bc1)
# Weak form - 1
a2 = inner(grad(u2), grad(v2))*dx
L2 = f*v2*dx
u2 = Function(W2)
solve(a2 == L2, u2, bc2)

# Mixed weak form
a = 10*inner(grad(u1_m), grad(v1_m))*dx + inner(grad(u2_m), grad(v2_m))*dx
L = f*v1_m*dx + f*v2_m*dx
# Solve the problem
sol = Function(V)
solve(a == L, sol, bcs, solver_parameters={"linear_solver":"direct"})

sol1 = sol.sub(0)
sol2 = sol.sub(1)

assert len(u1.vector()) == len(sol1.vector())
for i in range(len(u1.vector())):
    assert abs(sol1.vector()[i] - u1.vector()[i]) < 1e-10
    
assert len(u2.vector()) == len(sol2.vector())
for i in range(len(u2.vector())):
    assert abs(sol2.vector()[i] - u2.vector()[i]) < 1e-10

## Export result
encoding = XDMFFile.Encoding.HDF5 if has_hdf5() else XDMFFile.Encoding.ASCII

out_u1 = XDMFFile(MPI.comm_world, "block-assembly-2D2D-subdomain1-ref.xdmf")
out_u2 = XDMFFile(MPI.comm_world, "block-assembly-2D2D-subdomain2-ref.xdmf")
out_sub1 = XDMFFile(MPI.comm_world, "block-assembly-2D2D-subdomain1.xdmf")
out_sub2 = XDMFFile(MPI.comm_world, "block-assembly-2D2D-subdomain2.xdmf")

if MPI.size(MPI.comm_world) > 1 and encoding == XDMFFile.Encoding.ASCII:
    print("XDMF file output not supported in parallel without HDF5")
else:
    out_u1.write(u1, encoding)
    out_u2.write(u2, encoding)
    out_sub1.write(sol1, encoding)
    out_sub2.write(sol2, encoding)

The solution that I suppose is ‘u’ should still continuous at the interface? Please correct me if i made a mistake.

So my question is, if that continuous ‘u’ at the interface is what i need, did i implement it in an incorrect way or this is not the mixed-dimensional branch aimed at.

Thanks a lot for your help!

2

Hi @yang-sjtu,

Indeed there is no interface condition in this demo, it is only intended to illustrate the mixed function space and block assembly features. As you can see, the weak forms a1 = L1 and a2 = L2 are solved separately on submesh1 and submesh2, as reference solutions (no use of mixed-dimensional features here). Then we solve these uncoupled problems monolithically (a = L) using the mixed-dimensional features and check that we reproduce the reference solutions.

If you need some additional condition at the interface, you can possibly define a Lagrange multiplier living on the interface (defined as a lower submesh), you can have a look to the numerical examples here and the corresponding code

3

Hi @cdaversin,

Thank you so much for your explanation and providing the reference material. After reviewing the demo and Internal boundary condition with mixed_dimensional branch, I tried to implement the continuity condition as following, which is not working yet

  1. mesh, I assume the lower mesh for the interface is correctly defined at this step, as well as the measure dxGamma.
mesh = UnitSquareMesh(30, 30)

marker = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for c in cells(mesh):
    marker[c] = c.midpoint().y() < 0.5

submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 0)

gamma = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
for f in facets(mesh):
    if 0.5 - DOLFIN_EPS < f.midpoint().y() < 0.5 + DOLFIN_EPS:
        gamma[f] = 1
    
interface_mesh = MeshView.create(gamma, 1)

dxOmega = Measure("dx", domain=mesh)
dxGamma = Measure('dx', domain=interface_mesh)
  1. define the functional space, I got confused at this step, as in the demo, the mixed functional space including also the functional space defined on the interface lower mesh, while in Internal boundary condition with mixed_dimensional branch, the mixed functional space only mixed by the upper and lower subdomains. Please help me figure out which way is corret.
W1 = FunctionSpace(submesh1, "Lagrange", 1)
W2 = FunctionSpace(submesh2, "Lagrange", 1)
W3 = FunctionSpace(interface_mesh, "Lagrange", 1)
# Define the mixed function space
#V = MixedFunctionSpace( W1, W2, W3 )
V = MixedFunctionSpace( W1, W2 )
# Use directly TrialFunction and TestFunction on the product space
# Subdomain 1
u1 = TrialFunction(W1)
v1 = TestFunction(W1)
# Subdomain 2
u2 = TrialFunction(W2)
v2 = TestFunction(W2)
# Mixed
(u1_m,u2_m) = TrialFunctions(V)
(v1_m,v2_m) = TestFunctions(V)
#(u1_m,u2_m, lambda_) = TrialFunctions(V)
#(v1_m,v2_m, beta_) = TestFunctions(V)
  1. specifying the continuity at the interface. I tried F_coupling = va*ub*ds_int from Internal boundary condition with mixed_dimensional branch, but it ends up with error below:
Traceback (most recent call last):
  File "test_+Lagrange_multiplier.py", line 89, in <module>
    solve(a == L, sol, bcs, solver_parameters={"linear_solver":"direct"})
  File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/solving.py", line 259, in _solve_varproblem
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/problem.py", line 126, in __init__
    cpp.fem.MixedLinearVariationalProblem.__init__(self, a_list, L_list, u_comps, bcs)
RuntimeError: Cannot find common parent mesh

so is this way of implementing the continuity condition at the interface correct? If so, how should I solve this RuntimeError issue? Thanks a lot!

4

Hello @yang-sjtu,

From your snippets is not clear what ds_int is, but if you defined it as a “ds” type of measure, that might be the cause of that error. You should use the measure dxGamma that you defined at the end of your first snippet, which is created in based of a proper MeshView object.

1 Like