Hi there ! Currently, I am dealing with an interface condition where the following relationship holds at the interface
-
Displacement continuity condition:
u_L = u_R \quad \text{on } \Gamma -
Normal derivative jump condition:
where n_L = (1,0) is the outward normal from \Omega_L at the interface, and n_R = (-1,0) is the outward normal from \Omega_R at the interface.
I used the Lagrange multiplier method to derive the variational formulation, but the expected results were not obtained. I would like to ask how to handle this type of boundary condition. Your help would be greatly appreciated.
1. Problem Description
This code solves an interface problem: in a rectangular domain \Omega = [0,2] \times [0,1], divided into left and right subdomains by x=1, namely \Omega_L = [0,1] \times [0,1] and \Omega_R = [1,2] \times [0,1], we solve a Poisson equation in each subdomain while enforcing specific continuity conditions at the interface.
1.1 Governing Equations
In the left subdomain \Omega_L:
In the right subdomain \Omega_R:
where \Delta is the Laplacian operator, and f_L=1 and f_R=2 are source terms.
1.2 Boundary Conditions
-
Dirichlet boundary conditions: at x=2 and y=0, y=1
u = 0 -
Neumann boundary condition: at x=0
\frac{\partial u_L}{\partial n} = g_N = 0where n is the outward normal vector to the boundary.
1.3 Interface Conditions
At the interface \Gamma = \{(x,y) | x=1, 0 \leq y \leq 1\}:
-
Displacement continuity condition:
u_L = u_R \quad \text{on } \Gamma -
Normal derivative jump condition:
where n_L = (1,0) is the outward normal from \Omega_L at the interface, and n_R = (-1,0) is the outward normal from \Omega_R at the interface.
2. Weak Form Derivation
2.1 Classical Weak Form
Applying the weak form to each subdomain:
For \Omega_L, multiply by test function v_L and integrate over the domain:
For \Omega_R, multiply by test function v_R and integrate over the domain:
2.2 Introduction of Lagrange Multipliers
To handle the interface conditions, we introduce two Lagrange multipliers:
- \lambda_1: for the displacement continuity condition
- \lambda_2: for the normal derivative jump condition
2.3 Application of Lagrange Multiplier Method
The displacement continuity constraint u_L = u_R is enforced via Lagrange multiplier \lambda_1:
The normal derivative jump constraint \frac{\partial u_L}{\partial n_L} - \frac{\partial u_R}{\partial n_R} = G is enforced via Lagrange multiplier \lambda_2:
3. Complete Variational Form
Combining all the terms above, we get the complete variational problem:
Find (u_L, u_R, \lambda_1, \lambda_2) \in V_L \times V_R \times V_{I1} \times V_{I2} such that for all (v_L, v_R, \mu_1, \mu_2) \in V_L \times V_R \times V_{I1} \times V_{I2}:
where a is the bilinear form:
and L is the linear form:
4. Theoretical Explanation of the Lagrange Multiplier Method
The code is as blew:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
# Define subdomains for boundary conditions
class Dirichlet(SubDomain):
"""Define Dirichlet boundary conditions on the right edge (x=2) and top/bottom edges (y=0,1)"""
def inside(self, x, on_boundary):
return on_boundary and (near(x[0], 2) or near(x[1], 0) or near(x[1], 1))
class Neumann(SubDomain):
"""Define Neumann boundary conditions on the left edge (x=0)"""
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0)
class Interface(SubDomain):
"""Define the interface between left and right subdomains (x=1)"""
def inside(self, x, on_boundary):
return near(x[0], 1)
# Create the rectangular mesh [0,2]×[0,1]
n = 32 # mesh refinement parameter
mesh = RectangleMesh(Point(0, 0), Point(2, 1), 2 * n, n)
# Mark subdomains: 0 for left subdomain (x<1), 1 for right subdomain (x>1)
domain = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
for cell in cells(mesh):
domain[cell] = 1 if cell.midpoint().x() > 1.0 else 0
# Mark the interface facets (x=1)
facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
Interface().mark(facets, 1)
# Create mesh views for each subdomain and the interface
mesh_L = MeshView.create(domain, 0) # Left subdomain mesh
mesh_R = MeshView.create(domain, 1) # Right subdomain mesh
mesh_I = MeshView.create(facets, 1) # Interface mesh
# Define function spaces using continuous Galerkin elements (CG) of order 1
V_L = FunctionSpace(mesh_L, "CG", 1) # Function space for left subdomain
V_R = FunctionSpace(mesh_R, "CG", 1) # Function space for right subdomain
V_I1 = FunctionSpace(mesh_I, "CG", 1) # Function space for displacement continuity Lagrange multiplier
V_I2 = FunctionSpace(mesh_I, "CG", 1) # Function space for flux jump Lagrange multiplier
# Create a mixed function space for all variables
W = MixedFunctionSpace(V_L, V_R, V_I1, V_I2)
# Mark boundary conditions on the left subdomain boundary
ff_L = MeshFunction("size_t", mesh_L, mesh_L.topology().dim() - 1, 0)
Dirichlet().mark(ff_L, 1) # Mark Dirichlet boundaries with tag 1
Neumann().mark(ff_L, 2) # Mark Neumann boundaries with tag 2
Interface().mark(ff_L, 3) # Mark interface with tag 3
# Mark boundary conditions on the right subdomain boundary
ff_R = MeshFunction("size_t", mesh_R, mesh_R.topology().dim() - 1, 0)
Dirichlet().mark(ff_R, 1) # Mark Dirichlet boundaries with tag 1
Interface().mark(ff_R, 3) # Mark interface with tag 3
# Define Dirichlet boundary conditions (u=0 on marked boundaries)
g_D = Constant(0.0) # Value for Dirichlet boundary condition
bc_L = DirichletBC(W.sub_space(0), g_D, ff_L, 1) # Apply to left subdomain
bc_R = DirichletBC(W.sub_space(1), g_D, ff_R, 1) # Apply to right subdomain
bcs = [bc_L, bc_R] # Collect boundary conditions
# Define problem parameters
f_L = Constant(1.0) # Source term in left subdomain (-Δu_L = f_L)
f_R = Constant(2.0) # Source term in right subdomain (-Δu_R = f_R)
g_N = Constant(0.0) # Neumann boundary condition (∂u/∂n = g_N on x=0)
G = Constant(10.0) # Interface flux jump condition (∂u_L/∂n_L - ∂u_R/∂n_R = G)
# Define measures for integration over different domains
dx_L = Measure("dx", domain=mesh_L) # Integration over left subdomain
ds_L = Measure("ds", domain=mesh_L, subdomain_data=ff_L) # Integration over left boundary
dx_R = Measure("dx", domain=mesh_R) # Integration over right subdomain
ds_R = Measure("ds", domain=mesh_R, subdomain_data=ff_R) # Integration over right boundary
ds_I = Measure("dx", domain=mesh_I) # Integration over interface
# Define trial and test functions
(u_L, u_R, lambda_1, lambda_2) = TrialFunctions(W)
(v_L, v_R, mu_1, mu_2) = TestFunctions(W)
# Define the normal vectors at the interface
n_L = Constant((1.0, 0.0)) # Outward normal from left subdomain (pointing right)
n_R = Constant((-1.0, 0.0)) # Outward normal from right subdomain (pointing left)
# Define the bilinear form 'a' for the variational problem
a = (
# Standard Poisson terms for both subdomains
inner(grad(u_L), grad(v_L)) * dx_L # ∫_Ω_L ∇u_L·∇v_L dx
+ inner(grad(u_R), grad(v_R)) * dx_R # ∫_Ω_R ∇u_R·∇v_R dx
# Interface terms for displacement continuity (u_L = u_R)
+ inner(lambda_1, v_R - v_L) * ds_I # ∫_Γ λ₁(v_R - v_L) ds
+ inner(u_R - u_L, mu_1) * ds_I # ∫_Γ (u_R - u_L)μ₁ ds
# Interface terms for flux jump condition (∂u_L/∂n_L - ∂u_R/∂n_R = G)
+ inner(lambda_2, dot(grad(v_L), n_L) - dot(grad(v_R), n_R)) * ds_I # ∫_Γ λ₂(∇v_L·n_L - ∇v_R·n_R) ds
+ inner(dot(grad(u_L), n_L) - dot(grad(u_R), n_R), mu_2) * ds_I # ∫_Γ (∇u_L·n_L - ∇u_R·n_R)μ₂ ds
)
# Define the linear form 'L' for the variational problem
L = (
# Source terms
inner(f_L, v_L) * dx_L # ∫_Ω_L f_L·v_L dx
+ inner(f_R, v_R) * dx_R # ∫_Ω_R f_R·v_R dx
# Neumann boundary condition term
+ inner(g_N, v_L) * ds_L(2) # ∫_∂Ω_L∩{x=0} g_N·v_L ds
# Interface flux jump condition term
+ inner(G, mu_2) * ds_I # ∫_Γ G·μ₂ ds
)
# Solve the variational problem
u = Function(W)
solve(a == L, u, bcs, solver_parameters={"linear_solver": "direct"})
# Extract solution components
u_L, u_R, lambda_1, lambda_2 = u.split()
# Create individual functions for evaluation
u_L_func = Function(V_L)
u_R_func = Function(V_R)
lambda_1_func = Function(V_I1)
lambda_2_func = Function(V_I2)
# Assign solution components to individual functions
u_L_func.assign(u_L)
u_R_func.assign(u_R)
lambda_1_func.assign(lambda_1)
lambda_2_func.assign(lambda_2)
# Project gradients for flux evaluation
grad_u_L = project(grad(u_L), VectorFunctionSpace(mesh_L, "CG", 1))
grad_u_R = project(grad(u_R), VectorFunctionSpace(mesh_R, "CG", 1))
# Define interface sampling points
num_points = 20
y_points = np.linspace(0.01, 0.99, num_points)
# Initialize arrays to store sampled values
u_L_values = [] # Left solution values at interface
u_R_values = [] # Right solution values at interface
flux_L_values = [] # Left flux values at interface
flux_R_values = [] # Right flux values at interface
# Sample solution and flux values at the interface
for y in y_points:
# Define points slightly to the left and right of the interface
point_L = Point(0.9999, y) # Point just left of interface
point_R = Point(1.0001, y) # Point just right of interface
# Evaluate solution values
u_L_val = u_L_func(point_L)
u_R_val = u_R_func(point_R)
# Evaluate gradient values
grad_L_val = grad_u_L(point_L)
grad_R_val = grad_u_R(point_R)
# Calculate normal fluxes (derivatives in x-direction)
flux_L = grad_L_val[0] # Left flux (∂u_L/∂x at x=1⁻)
flux_R = -grad_R_val[0] # Right flux (negative of ∂u_R/∂x at x=1⁺ due to opposite normal)
# Store values for analysis
u_L_values.append(u_L_val)
u_R_values.append(u_R_val)
flux_L_values.append(flux_L)
flux_R_values.append(flux_R)
# Calculate errors in interface conditions
continuity_error = np.array(u_L_values) - np.array(u_R_values) # u_L - u_R should be 0
flux_jump_error = np.array(flux_L_values) - np.array(flux_R_values) - 10.0 # flux_L - flux_R should be G=10
# Print statistical information about errors
print("Interface continuity check:")
print(f" Maximum error: {np.max(np.abs(continuity_error)):.6e}")
print(f" Average error: {np.mean(np.abs(continuity_error)):.6e}")
print("Flux jump condition check:")
print(f" Maximum error: {np.max(np.abs(flux_jump_error)):.6e}")
print(f" Average error: {np.mean(np.abs(flux_jump_error)):.6e}")
Got wrong output at the interface
Interface continuity check:
Maximum error: 1.612573e-02
Average error: 3.787325e-03
Flux jump condition check:
Maximum error: 2.121465e+02
Average error: 5.768207e+01
Any hint will be greatly appreciated. I’m not fixated on using the Lagrange multiplier method. I welcome any feasible methods that can help me solve this problem of handling the interface condition with the given relationship at the interface to obtain the expected results from the variational formulation.
