MultiMesh: Mixed space over different meshes

3 months ago by
I am trying to implement a mixed system of the following form:

Find \((r,u) \in V \times U\) such that
\((J(r), v) +(Bu, v) = (f,v)\) for all \(v \in V\)
\(\phantom{(J(r), v) +} (Bw, r) = 0   \) for all \(w \in U\),
where J is a non-linear, B a linear operator, and U and V are finite element spaces.

The finite element spaces I would like to use are the following:
mesh_1 = UnitSquareMesh(N,N)
mesh_2 = UnitSquareMesh(2*N,2*N)

U      = FunctionSpace(mesh_1, 'Lagrange',n)
V      = FunctionSpace(mesh_2, 'Lagrange',m)​

Here m>n. If I had the same mesh for both, this would easily be possible using MixedElement.
Is it possible to do this with MultiMesh?

Let's simplify things a little and consider the following coupled linear system:
\[ - \Delta u_1 + u_2 = f_1 \text{ in } \Omega\\
-\Delta u_2 +u_1 = f_2 \text{ in } \Omega\]
Note that both equations are defined on the same domain.

Here is my attempt at implementing the system:

# coupled linear system
# -(dxx u1 + dyy u1) + u2 = f1
# -(dxx u2 + dyy u2) + u1 = f2
# u1 = u2 = sin(2 pi x)sin(2 pi y)

from fenics import *

#Create multimesh consisting of two overlapping meshes
N = 10
mesh_1 = UnitSquareMesh(N,N)
shift = 1e-15
mesh_2 =RectangleMesh(Point(shift, shift), Point(1-shift, 1-shift), 2*N, 2*N)
multimesh = MultiMesh()

#Create function space over multimesh
n = 1
W = MultiMeshFunctionSpace(multimesh, "Lagrange", n)

#Create test and trial functions corresponding to each of the meshes
u = TrialFunction(W)
v =  TestFunction(W)

nu = FacetNormal(multimesh)

sig = 1000
#Variational form integrating over overlap
a = inner(grad(u('+')),grad(v('+')))*dO + u('+')*v('-')*dO\
    +inner(grad(u('-')), grad(v('-')))*dO + u('-')*v('+')*dO\
    +inner(grad(u),grad(v))*dC + u*v*dC\
    +inner(grad(u),grad(v))*dX + u*v*dX\
    -dot(grad(u('+')),nu('+'))*v('+')*dI - u('+')*dot(grad(v('+')),nu('+'))*dI\
    -dot(grad(u('-')),nu('-'))*v('-')*dI - u('-')*dot(grad(v('-')),nu('-'))*dI\

f1 = Expression("(8*pow(pi,2)+1)*sin(2*x[1]*pi)*sin(2*x[0]*pi)", degree=2)
f2 = Expression("(8*pow(pi,2)+1)*sin(2*x[1]*pi)*sin(2*x[0]*pi)", degree=2)
l = f1*v('+')*dO+f2*v('-')*dO+f2*v*dC+f2*v*dX

A = assemble_multimesh(a)
b = assemble_multimesh(l)

#bc = MultiMeshDirichletBC(W, Constant(0),  DomainBoundary())
# Compute solution
w = MultiMeshFunction(W)
solve(A, w.vector(), b, "cg")

File('u.pvd') << w.part(0)
File('r.pvd') << w.part(1)

I still have the following issues:

1) The integrals over the overlap that are defined using the measure dO would have to be over the whole domain instead of only a certain neighbourhood of the interface to express the mixed terms as well as the terms on the first mesh properly.

2) I seem to have a problem with MultiMeshDirichletBC in general. I am using dolfin 2017.2.0. I don’t know if this issue has been resolved in a development version of dolfin. In the code above I am applying boundary conditions using Nitsche’s method. Even this was only possible after making the domain for the top mesh slightly smaller (see definition of mesh_2) so that the domains aren’t exactly aligned.

3) Is it possible to have different polynomial degrees for different parts?

Is there a way to address the above issues?
Is there a different way of approaching my problem?

Community: FEniCS Project

1 Answer

3 months ago by
I dont have much idea of the multimesh that you are referring but if you want to solve PDE on one mesh and then project the solution to the other then you  may try

You can project from one mesh to another using  interpolate_nonmatching_mesh_any   after installing fenicstools

This works in fenics latest version.
Interpolation is not projection.
written 3 months ago by Nate  
Thank you for the suggestion, but this is unfortunately not what I am trying to do. I am trying to solve on two meshes simultaneously. I have a coupled system of equations that I cannot decouple and each of the variables has to live in a different finite element space.
written 3 months ago by Sarah R  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »