How to solve a coupled multiphysics problem using one mesh with subdomains?

4 months ago by
Dear FEniCS community,

I am a new user of FEniCS trying to solve a multiphysics problem with its help. I have generated a simple mesh from two subdomains with a conforming interface between them, and created a FunctionSpace on the whole mesh using a MixedElement (composed of a Vector- and FiniteElement). Now I want to solve for two dependent variables (scalar and vector) in the first subdomain and a single scalar dependent variable in the second subdomain (the scalar variable is common for the two subdomains=the same variable), with a coupling condition on the interface. How can I do this in FEniCS? If I add any (non-physical, unwanted) component with the vector variable (even multiplied by a very small number) to the bilinear form in the second subdomain, I am able to get some (wrong) results. Otherwise, the solution is nan everywhere except where the boundary conditions were applied. Eliminating superfluous DOFs by setting
DirichletBC(W.sub(0), Constant((0.,)*ndim), SecondDomainIdentifier, method='pointwise')​
does not work correctly either. What is the proper way for solving such problems? What is worth mentioning here is that I have found a quite recent presentation (Link, page 5, 7th slide) tackling similar issue to some extent, but the key function FunctionSpaceProduct that in my opinion will serve the discussed purpose is not available in my distribution of dolfin (2017.2.0, ubuntu). Am I doing something wrong?

A minimal working example is attached below.

Thank you very much for your time and help! Best regards.

# run using dolfin v2017.2.0
from dolfin import *
from mshr import *

def build_mesh(l1, l2, h, nl1, nl2, nh):
	domain_1 = Rectangle(Point(.0, .0), Point(l1, h)) #
	domain_2 = Rectangle(Point(l1, h), Point(l1+l2, .0)) 
	domain_tot = domain_1 + domain_2

	domain_tot.set_subdomain(1, domain_1) #
	domain_tot.set_subdomain(2, domain_2)
	m = generate_mesh(domain_tot, 25) # FIXME: should use nl1, nl2, nh to control the mesh resolution

	return m

#===== BORDERS =====
# Geometric parameters
H = 4.E-3       # height
L2 = 24.E-3     # second thickness
L1 = L2         # first thickness

#===== MESH =====
# Create mesh
mesh = build_mesh(L1, L2, H, 10, 10, 8)
#n = FacetNormal(mesh) # FIXME: this normal vector is not consistent on interior facets as it is probably defined cell-wise
V = VectorFunctionSpace(mesh, "CG", mesh.topology().dim(), dim=2)
n = project(Constant((-1., 0.)), V) # normal vectors pointing in the negative x-direction (outside the second domain at the interface)
ndim = mesh.geometry().dim()

#===== SUBDOMAINS =====
# Define subdomain markers
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim(),

# Define integration measure for domain integrals
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

#===== FUNCTION SPACES =====
Uh = VectorElement('P', mesh.ufl_cell(), 2)
Ph = FiniteElement('P', mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Uh, Ph]))
(u, p) = TrialFunctions(W)
(du, dp) = TestFunctions(W)

u_bc = Constant((0.,)*ndim)
p_bc = Constant(1.)

# Define boundaries
tol = 1.E-14
def Left(x, on_boundary):
	return on_boundary and x[0] < tol
def Right(x, on_boundary):
	return on_boundary and near(x[0], L1 + L2, tol)
def Symmetry(x, on_boundary):
	return on_boundary and near(x[1], 0., tol) or near(x[1], H, tol)
def Interface(x):
	return near(x[0], L1, tol)

class SecondDomain(SubDomain): #
	def inside(self, x, on_boundary):
		return x[0] >= L1 - tol
SecondDomainIdentifier = SecondDomain()

bcs = [DirichletBC(W.sub(0), u_bc, Left), DirichletBC(W.sub(0).sub(1), Constant(0.), Symmetry), DirichletBC(W.sub(1), p_bc, Right)\
#	, DirichletBC(W.sub(0), Constant((0.,)*ndim), SecondDomainIdentifier, method='pointwise')\

# Create mesh function over cell facets
boundary_subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundary_subdomains.set_all(0) # 0 for interior cell boundaries
interface_boundary = AutoSubDomain(Interface)
interface_boundary.mark(boundary_subdomains, 1) # 1 for interface boundary

# Define integration measure for boundary integrals
dS = Measure('dS', domain=mesh, subdomain_data=boundary_subdomains) # measure for internal boundaries

a = p*dp*dx(2)
#a += inner(u, du)*dx(2) # TODO: non-physical component
a += ( p*dp + inner(u, du) )*dx(1)
a_interface = ( avg(dp)*inner(u('+'),n) + avg(p)*inner(du('+'),n) )*dS(1)
L = Constant(0.)*dp*dx(1)

A = assemble(a) # component matrix
A_interface = assemble(a_interface) # interface coupling matrix
A += A_interface # governing matrix
b = assemble(L) # right-hand-side vector

for bc in bcs: bc.apply(A, b)

#===== SOLUTION =====
up = Function(W)
solve(A, up.vector(), b, 'lu')​
Community: FEniCS Project
Please login to add an answer/comment or follow this question.

Similar posts:
Search »