How to correctly use MeshView/Submeshes


I am trying to use use fenics for multiphysics problem but i do not understand how to properly use submeshes. How do I define variational form that solves two problem simultaniously?
And perhaps how do I create NonlinearVariationalProblem and Jacobian having mixed function spaces (potentialy with different number of variables) in each submesh? I always get errors and I cannot find a consistent way of doing such stuff.

Are there any good tutorials that use only fenics? Should i use multiphemics instead?

Consider this MWE:

from mshr import *
from fenics import *
import matplotlib.pyplot as plt
tol  = 1E-14

# define two parts of the complete-domain
omega1 = Rectangle(Point(0, 0), Point(0.5, 1))
omega2 = Rectangle(Point(0.5, 0), Point(1, 1))
domain = omega1 + omega2

# set left and right domains
domain.set_subdomain(1, omega1)
domain.set_subdomain(2, omega2)

# generate mesh fuction that marks domains
mesh = generate_mesh(domain, 32)
marker = MeshFunction("size_t", mesh, 2,

# create submeshes
submesh1 = MeshView.create(marker, 1)
submesh2 = MeshView.create(marker, 2)

CG = FiniteElement("CG", mesh.ufl_cell(), 1) # element type
element_1 = MixedElement([CG, CG])         # mixed element for left domain
element_2 = MixedElement([CG, CG])         # mixed element for right
W1 = FunctionSpace(submesh1, element_1)    # function space for left domain
W2 = FunctionSpace(submesh2, element_2)    # function space for right domain
W = MixedFunctionSpace(W1, W2)

dx_1 = Measure("dx", domain=submesh1)   # so we can integrate only over left domain
dx_2 = Measure("dx", domain=submesh2)   # so we can integrate only over right domain

u = Function(W)
u1, u2 = u.split()  # mixed solutions of each domain
T1, C1 = u1.split() # individual fields of the let domain
T2, C2 = u2.split() # individual fields of the roght

v1, v2 = TestFunctions(W)
T1_v, C1_v = v1
T2_v, C2_v = v2

lmbd_1 = 1
lmbd_2 = 2
# poisson for left side
PDE = inner(lmbd_1*nabla_grad(T1), nabla_grad(T1_v))*dx_1

# poisson for right side 
PDE += inner(lmbd_2*nabla_grad(T2), nabla_grad(T2_v))*dx_2

# hello C1 be constant on the left
PDE += C1*C1_v*dx_1 - Constant(1.0)*C1_v*dx_1

# hello C2 be constant on the right but higher
PDE += C2*C2_v*dx_2 - Constant(2.0)*C2_v*dx_2

# poisson needs a boundary conditions
def boundary_left(x, on_boundary):
    return on_boundary and near(x[0], 0., tol)

def boundary_mid(x, on_boundary):
    # How does on_boundary behave with meshview?
    return on_boundary and near(x[0], 0.5, tol)

def boundary_right(x, on_boundary):
    return on_boundary and near(x[0], 1., tol)

# not sure about this
bc1 = DirichletBC(W.sub_space(0).sub(0), Constant(1.0), boundary_left)
bc2 = DirichletBC(W.sub_space(0).sub(0), Constant(2.0), boundary_mid)
bc3 = DirichletBC(W.sub_space(1).sub(0), Constant(1.0), boundary_mid)
bc4 = DirichletBC(W.sub_space(1).sub(0), Constant(2.0), boundary_right)
bcs = [bc1, bc2, bc3, bc4]

solve(PDE == 0, u, bcs) #KSP error internal 76

sol_u1, sol_u2 = u.split() # domain solutions
sol_T1, sol_C1 = sol_u1.split() # left side solution
c = plot(sol_T1)

Which gives this error:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

How do I make the MWE above work?

Thank you for any reply.

Have you considered the supplementary material: for the paper on MeshView (View article)

1 Like

Yes I have seen those. I did not read the paper though, so I am going to read it. Thank you for the link!

I was unable to find how to treat the case where I need mixed elements in both submeshes (there is example that uses VectorFunctionSpace, but I would need to index into the individual components in the form - I am not sure if that is the best way to do it). It also does not show any example where derivative such as Jacobian is needed. Or where we want to define NonlinearVaratiationalProblem (for problems in form of F == 0 defined via Function instead of TrialFunction), the examples show only how to use solve directly. It focuses on examples in form a == L ( I understood that there is fundamental difference between the use of Function and TrialFunction in definition of the form).

Or maybe it shows that and I just do not understand how to derive that from those examples. Maybe a do not understand the limitations of fenics in this regard.

I guess my question is:
Is it even meant to be used as I am trying to do in the MWE (having different mixed function spaces in each domain with different physics, being perhaps coupled through some internal boundary)?

Honestly I am just confused right now.

I have figured out that SubMesh is not ment to be used in parallel and although MeshView is supposed to be used that way, it is currently not working. (See BoundaryMesh of MeshView in Parallel - #3 by LiborKudela).

The MWE provided at the top does not work in serial either… I still do not know why. It might not be meant for nonlinear variational problem.

Anyway I am going to try to find a workaround.

MeshView is meant to work in parallel, but there seems to be some bugs.
As I’m not the author of MeshView, I cannot comment to much on what is working and what is not.

There is a fix for some issues in the branch Bitbucket
which there is a docker image for at fenics versions · scientificcomputing · GitHub
I would suggest @LiborKudela and you to test it and see if it resolves any issues.

1 Like

Also note this is legacy DOLFIN and no longer actively developed (for ~3-4 years now). I’d recommend DOLFINx.