Thanks, Stein. That’s exactly what I want. I am new to fenics and by using two fields defined over two submeshes is cited from topic_example
Also I have found topic is similar to my problem ,and I have rerite the code
import matplotlib.pyplot as plt
import dolfinx.fem.petsc
import ufl
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, div, FacetNormal, inner, avg
import numpy as np
from dolfinx import default_scalar_type
from petsc4py import PETSc
from mpi4py import MPI
from dolfinx import fem, mesh, io, plot
import dolfinx
def V_exact(mode):
"""Define exact solution function"""
return lambda x: -mode.cos(mode.pi * x[1])
V_numpy = V_exact(np) # Used for interpolation
V_ufl = V_exact(ufl) # Used for defining source term
def solve_poisson(N, degree=1):
"""
Solve Poisson equation with interface conditions
"""
# Create unit square mesh
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])],
[N, N], mesh.CellType.triangle)
# Why connectivity is needed:
# - Boundary conditions are defined on edges
# - Finite element degrees of freedom are associated with cells
# - To know which DOFs are affected by boundary conditions, the code needs to know which cells connect to boundary edges
domain.topology.create_connectivity(domain.topology.dim - 1, domain.topology.dim)
# Create discontinuous Galerkin space for representing discontinuous material properties
Q = fem.functionspace(domain, ("DG", 0))
# Define subdomains
def Omega_0(x):
return x[1] <= 0.50
def Omega_1(x):
return x[1] >= 0.50
# Set spatial coordinates
x = SpatialCoordinate(domain)
# Create conductivity function and assign values
sigma = fem.Function(Q)
cells_0 = mesh.locate_entities(domain, domain.topology.dim, Omega_0)
cells_1 = mesh.locate_entities(domain, domain.topology.dim, Omega_1)
def anode_conductivity(T):
return 1. / (5.929e-5 - T * 1.235e-8)
# Set different conductivity values for different regions
sigma.x.array[cells_0] = np.full_like(cells_0, 210, dtype=default_scalar_type)
sigma.x.array[cells_1] = np.full_like(cells_1, anode_conductivity(800), dtype=default_scalar_type)
# Calculate source term based on exact solution
f = div(-sigma * grad(V_ufl(x)))
# Create first-order Lagrange function space for solution
W = fem.functionspace(domain, ("Lagrange", 1))
V = TrialFunction(W)
phi = TestFunction(W)
# Variational form left side: bilinear form
a = dot(sigma * grad(V), grad(phi)) * dx
# Mark boundaries and interface
boundaries = [(1, lambda x: np.isclose(x[0], 0)), # Left boundary
(2, lambda x: np.isclose(x[0], 1)), # Right boundary
(3, lambda x: np.isclose(x[1], 0)), # Bottom boundary
(4, lambda x: np.isclose(x[1], 1)), # Top boundary
(5, lambda x: np.isclose(x[1], 0.5))] # Interface
# Create boundary markers
facet_indices, facet_markers = [], []
fdim = domain.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(domain, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
# Define interior interface integration measure
dS = ufl.Measure("dS", domain=domain, subdomain_data=facet_tag)
# Calculate flux jump at the interface
n = FacetNormal(domain)
# Calculate flux jump at the interface based on exact solution
g = ufl.pi * (sigma('+') - sigma('-')) # Because cos(pi*y) has zero derivative at y=0.5
# Variational form right side: linear form
L = f * phi * dx + g * avg(phi) * dS(5) # Using average of test function
# Set Dirichlet boundary conditions
facets_bottom = facet_tag.find(3)
dofs_bottom = fem.locate_dofs_topological(W, fdim, facets_bottom)
facets_top = facet_tag.find(4)
dofs_top = fem.locate_dofs_topological(W, fdim, facets_top)
BCs = [
fem.dirichletbc(PETSc.ScalarType(-1), dofs_bottom, W), # Bottom boundary condition
fem.dirichletbc(PETSc.ScalarType(1), dofs_top, W) # Top boundary condition
]
# Solve linear problem
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=BCs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
return problem.solve(), V_ufl(x)
def error_L2_func(Vh, V_ex, degree_raise=3):
"""Calculate L2 error norm"""
# Create higher-order function space
degree = 1
family = Vh.function_space.ufl_element().family_name
mesh = Vh.function_space.mesh
Q = fem.functionspace(mesh, (family, degree + degree_raise))
# Interpolate numerical solution
V_W = fem.Function(Q)
V_W.interpolate(Vh)
# Interpolate exact solution
V_ex_W = fem.Function(Q)
if isinstance(V_ex, ufl.core.expr.Expr):
u_expr = fem.Expression(V_ex, Q.element.interpolation_points)
V_ex_W.interpolate(u_expr)
else:
V_ex_W.interpolate(V_ex)
# Calculate error
e_W = fem.Function(Q)
e_W.x.array[:] = V_W.x.array - V_ex_W.x.array
# Integrate error
error = fem.form(ufl.inner(e_W, e_W) * ufl.dx)
error_local = fem.assemble_scalar(error)
error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
return np.sqrt(error_global)
# Convergence test
N = [40, 80, 160, 320, 640]
error_L2 = []
error_H1 = []
h = []
mpi_rank = 5
for i in range(len(N)):
Vh, Vex = solve_poisson(N[i])
comm = Vh.function_space.mesh.comm
# Calculate L2 error
error_L2 += [error_L2_func(Vh, V_numpy)]
# Calculate H1 seminorm error
eh = Vh - Vex
error_H10 = fem.form(dot(grad(eh), grad(eh)) * dx)
E_H10 = np.sqrt(comm.allreduce(fem.assemble_scalar(error_H10), op=MPI.SUM))
error_H1 += [E_H10]
h += [1. / N[i]]
if comm.rank == 0:
mpi_rank = comm.rank
print(f"h: {h[i]:.2e} Error L2: {error_L2[i]:.2e}")
print(f"h: {h[i]:.2e} Error H1: {error_H1[i]:.2e}")
# Plot convergence curves
if mpi_rank == 0:
plt.figure(figsize=(10, 6))
plt.loglog(N, error_L2, 'o-', label='$L^2$ error')
plt.loglog(N, error_H1, 's-', label='$H^1$ error')
# Reference lines
plt.loglog(N, h, '--', label='$O(h)$')
h_square = [x ** 2 for x in h]
plt.loglog(N, h_square, '-.', label='$O(h^2)$')
plt.xlabel('Mesh resolution N')
plt.ylabel('Error')
plt.title('Error Convergence Analysis')
plt.legend()
plt.grid(True)
plt.show()
But please forgive my ignorance as a beginner. I still want to know the implementation process of using weak forms such as the Lagrange multiplier method or the Nitsche method. I’ve found two similar examples, but I don’t know how to modify them into the format of interface flux jump to meet the requirements of my problem.
Nitsche method example:
# In this demo, we implement a domain decomposition scheme for
# the Poisson equation based on Nitche's method. The scheme can
# be found in "Mortaring by a method of J. A. Nitsche" by Rolf
# Stenberg. See also "A finite element method for domain
# decomposition with non-matching grids" by Becker et al.
from dolfinx import mesh, fem, io
from mpi4py import MPI
import ufl
from ufl import inner, grad, avg, div
import numpy as np
from petsc4py import PETSc
from utils import norm_L2, convert_facet_tags, interface_int_entities
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from meshing import create_square_with_circle
def u_e(x, module=np):
"A function to represent the exact solution"
return module.exp(-((x[0] - 0.5) ** 2 + (x[1] - 0.5) ** 2) / (2 * 0.05**2)) + x[0]
# Set some parameters
comm = MPI.COMM_WORLD
h = 0.05 # Maximum cell diameter
k_0 = 1 # Polynomial degree in omega_0
k_1 = 3 # Polynomial degree in omega_1
# Create mesh and sub-meshes
msh, ct, ft, vol_ids, surf_ids = create_square_with_circle(comm, h)
tdim = msh.topology.dim
submesh_0, sm_0_to_msh = mesh.create_submesh(msh, tdim, ct.find(vol_ids["omega_0"]))[:2]
submesh_1, sm_1_to_msh = mesh.create_submesh(msh, tdim, ct.find(vol_ids["omega_1"]))[:2]
# Define function spaces on each submesh
V_0 = fem.functionspace(submesh_0, ("Lagrange", k_0))
V_1 = fem.functionspace(submesh_1, ("Lagrange", k_1))
W = ufl.MixedFunctionSpace(V_0, V_1)
# Test and trial functions
u = ufl.TrialFunctions(W)
v = ufl.TestFunctions(W)
# We use msh as the integration domain, so we require maps from cells
# in msh to cells in submesh_0 and submesh_1. These can be created
# as follows:
cell_imap = msh.topology.index_map(tdim)
num_cells = cell_imap.size_local + cell_imap.num_ghosts
msh_to_sm_0 = np.full(num_cells, -1)
msh_to_sm_0[sm_0_to_msh] = np.arange(len(sm_0_to_msh))
msh_to_sm_1 = np.full(num_cells, -1)
msh_to_sm_1[sm_1_to_msh] = np.arange(len(sm_1_to_msh))
# Compute integration entities for the interface integral
fdim = tdim - 1
interface_facets = ft.find(surf_ids["interface"])
domain_0_cells = ct.find(vol_ids["omega_0"])
domain_1_cells = ct.find(vol_ids["omega_1"])
# Create interface integration entities and modify msh_to_sm maps
interface_entities, msh_to_sm_0, msh_to_sm_1 = interface_int_entities(
msh, interface_facets, msh_to_sm_0, msh_to_sm_1
)
# Create entity maps using the modified msh_to_sm maps
entity_maps = {submesh_0: msh_to_sm_0, submesh_1: msh_to_sm_1}
# Create integration measures
dx = ufl.Measure("dx", domain=msh, subdomain_data=ct)
dS = ufl.Measure("dS", domain=msh, subdomain_data=[(surf_ids["interface"], interface_entities)])
x = ufl.SpatialCoordinate(msh)
kappa = [1.0 + 0.1 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1]) for _ in range(2)]
# Penalty parameter (including harmonic mean on kappa on interface)
# TODO Add k dependency
gamma = 10 * 2 * kappa[0] * kappa[1] / (kappa[0] + kappa[1])
h = ufl.CellDiameter(msh)
n = ufl.FacetNormal(msh)
def jump_i(v, n):
return v[0]("+") * n("+") + v[1]("-") * n("-")
def grad_avg_i(v, kappa):
return kappa[1] / (kappa[0] + kappa[1]) * kappa[0] * grad(v[0]("+")) + kappa[0] / (
kappa[0] + kappa[1]
) * kappa[1] * grad(v[1]("-"))
a = (
inner(kappa[0] * grad(u[0]), grad(v[0])) * dx(vol_ids["omega_0"])
+ inner(kappa[1] * grad(u[1]), grad(v[1])) * dx(vol_ids["omega_1"])
- inner(grad_avg_i(u, kappa), jump_i(v, n)) * dS(surf_ids["interface"])
- inner(jump_i(u, n), grad_avg_i(v, kappa)) * dS(surf_ids["interface"])
+ gamma / avg(h) * inner(jump_i(u, n), jump_i(v, n)) * dS(surf_ids["interface"])
)
# Compile LHS forms
a = fem.form(ufl.extract_blocks(a), entity_maps=entity_maps)
# Define right-hand side forms
f = [-div(kap * grad(u_e(ufl.SpatialCoordinate(msh), module=ufl))) for kap in kappa]
# L = inner(f[0], v[0]) * dx(vol_ids["omega_0"]) + inner(f[1], v[1]) * dx(vol_ids["omega_1"])
#
G = 1.0 #
#
L = (inner(f[0], v[0]) * dx(vol_ids["omega_0"]) +
inner(f[1], v[1]) * dx(vol_ids["omega_1"]) +
0.5 * G * (v[0]("+") + v[1]("-")) * dS(surf_ids["interface"]))
# Compile RHS forms and set block structure
L = fem.form(ufl.extract_blocks(L), entity_maps=entity_maps)
# Apply boundary conditions. We require the DOFs of V_0 on the domain
# boundary. These can be identified via that facets of submesh_0 that
# lie on the domain boundary.
ft_sm_0 = convert_facet_tags(msh, submesh_0, sm_0_to_msh, ft)
bound_facets_sm_0 = ft_sm_0.find(surf_ids["boundary"])
submesh_0.topology.create_connectivity(fdim, tdim)
bound_dofs = fem.locate_dofs_topological(V_0, fdim, bound_facets_sm_0)
u_bc_0 = fem.Function(V_0)
u_bc_0.interpolate(u_e)
bc_0 = fem.dirichletbc(u_bc_0, bound_dofs)
bcs = [bc_0]
# Assemble linear system of equations
A = assemble_matrix_block(a, bcs=bcs)
A.assemble()
b = assemble_vector_block(L, a, bcs=bcs)
# Set-up solver
ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("superlu_dist")
# Compute solution
x = A.createVecRight()
ksp.solve(b, x)
# Recover solution
u_0, u_1 = fem.Function(V_0), fem.Function(V_1)
offset = V_0.dofmap.index_map.size_local * V_0.dofmap.index_map_bs
u_0.x.array[:offset] = x.array_r[:offset]
u_1.x.array[: (len(x.array_r) - offset)] = x.array_r[offset:]
u_0.x.scatter_forward()
u_1.x.scatter_forward()
# Write solution to file
with io.VTXWriter(msh.comm, "u_0.bp", u_0, "BP4") as f:
f.write(0.0)
with io.VTXWriter(msh.comm, "u_1.bp", u_1, "BP4") as f:
f.write(0.0)
# Compute error in solution
e_L2_0 = norm_L2(msh.comm, u_0 - u_e(ufl.SpatialCoordinate(submesh_0), module=ufl))
e_L2_1 = norm_L2(msh.comm, u_1 - u_e(ufl.SpatialCoordinate(submesh_1), module=ufl))
e_L2 = np.sqrt(e_L2_0**2 + e_L2_1**2)
if msh.comm.rank == 0:
print(e_L2)
Lagrange multiplier example from github_link (I want to implement it solely using FEniCS)
# Primal, multiscale
from dolfin import *
from xii import *
from emi.utils import H1_norm, L2_norm, subdomain_interpolate
import emi.fem.common as common
from xii.assembler.trace_matrix import trace_mat_no_restrict
def setup_problem(n, mms, params):
'''Multi-dimensional primal formulation'''
base_mesh = UnitSquareMesh(mpi_comm_self(), *(n,) * 2)
# Marking
inside = ['(0.25-tol<x[0])', '(x[0] < 0.75+tol)', '(0.25-tol<x[1])', '(x[1] < 0.75+tol)']
inside = CompiledSubDomain(' && '.join(inside), tol=1E-10)
mesh_f = MeshFunction('size_t', base_mesh, base_mesh.topology().dim(), 0)
inside.mark(mesh_f, 1)
inner_mesh = SubMesh(base_mesh, mesh_f, 1) # Inside
outer_mesh = SubMesh(base_mesh, mesh_f, 0) # Ouside
interface_mesh = BoundaryMesh(inner_mesh, 'exterior')
# Spaces
V1 = FunctionSpace(outer_mesh, 'CG', 1)
V = FunctionSpace(inner_mesh, 'CG', 1)
Q = FunctionSpace(interface_mesh, 'CG', 1)
W = [V1, V, Q]
u1, u, p = map(TrialFunction, W)
v1, v, q = map(TestFunction, W)
Tu1, Tv1 = (Trace(f, interface_mesh) for f in (u1, v1))
Tu, Tv = (Trace(f, interface_mesh) for f in (u, v))
# Mark subdomains of the interface mesh (to get source terms therein)
subdomains = mms.subdomains[1] #
marking_f = MeshFunction('size_t', interface_mesh, interface_mesh.topology().dim(), 0)
[subd.mark(marking_f, i) for i, subd in enumerate(map(CompiledSubDomain, subdomains), 1)]
# The line integral
dx_ = Measure('dx', domain=interface_mesh, subdomain_data=marking_f)
kappa, epsilon = map(Constant, (params.kappa, params.eps))
a = block_form(W, 2)
a[0][0] = kappa * inner(grad(u1), grad(v1)) * dx
a[0][2] = inner(Tv1, p) * dx_
a[1][1] = inner(grad(u), grad(v)) * dx
a[1][2] = -inner(Tv, p) * dx_
a[2][0] = inner(Tu1, q) * dx_
a[2][1] = -inner(Tu, q) * dx_
a[2][2] = -epsilon * inner(p, q) * dx_
# Data
# Source for domains, outer boundary data, source for interface
f1, f, gBdry, gGamma, hGamma = mms.rhs
L = block_form(W, 1)
L[0] = inner(f1, v1) * dx
L[0] += sum(inner(hi, Tv1) * dx_(i) for i, hi in enumerate(hGamma, 1))
# Iface contribution
L[1] = inner(f, v) * dx
L[2] = sum(inner(gi, q) * dx_(i) for i, gi in enumerate(gGamma, 1))
A, b = map(ii_assemble, (a, L))
subdomains = mms.subdomains[0] #
facet_f = MeshFunction('size_t', outer_mesh, outer_mesh.topology().dim() - 1, 0)
[subd.mark(facet_f, i) for i, subd in enumerate(map(CompiledSubDomain, subdomains), 1)]
V1_bcs = [DirichletBC(V1, gi, facet_f, i) for i, gi in enumerate(gBdry, 1)]
bcs = [V1_bcs, [], []]
A, b = apply_bc(A, b, bcs)
return A, b, W
setup_mms = common.setup_mms
def setup_error_monitor(mms_data, params):
'''Compute function mapping numerical solution to errors'''
# Error of the solution ...
exact = mms_data.solution
subdomains = mms_data.subdomains[1]
def get_error(wh, subdomains=subdomains, exact=exact, mms=mms_data, params=params):
u1h, uh, ph = wh
sigma_exact, u_exact, p_exact, I_exact = exact
# Mutliplier error
mesh = ph.function_space().mesh()
Q = FunctionSpace(mesh, 'CG', 3)
cell_f = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
# Spaces on pieces
for tag, subd in enumerate(subdomains, 1):
CompiledSubDomain(subd, tol=1E-10).mark(cell_f, tag)
dx = Measure('dx', domain=mesh, subdomain_data=cell_f)
p, q = TrialFunction(Q), TestFunction(Q)
a = inner(p, q) * dx
L = sum(inner(I_exact_i, q) * dx(i) for i, I_exact_i in enumerate(I_exact, 1))
I_exact = Function(Q)
A, b = map(assemble, (a, L))
solve(A, I_exact.vector(), b)
V = uh.function_space()
mesh = V.mesh()
#
# Recover the transmembrane potential by postprocessing
#
V = uh.function_space()
mesh = V.mesh()
interface_mesh = BoundaryMesh(mesh, 'exterior')
V = FunctionSpace(interface_mesh, 'CG', 1)
Tu1 = PETScMatrix(trace_mat_no_restrict(u1h.function_space(), V)) * u1h.vector()
Tu = PETScMatrix(trace_mat_no_restrict(uh.function_space(), V)) * uh.vector()
vh = Function(V, Tu1 - Tu)
# Now using flux we have
Q = ph.function_space()
mesh = Q.mesh()
cell_f = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
# Spaces on pieces
for tag, subd in enumerate(subdomains, 1):
CompiledSubDomain(subd, tol=1E-10).mark(cell_f, tag)
dx = Measure('dx', domain=mesh, subdomain_data=cell_f)
p, q = TrialFunction(Q), TestFunction(Q)
gGamma = mms.rhs[3]
a = inner(p, q) * dx
L = inner(params.eps * ph, q) * dx + sum(inner(gi, q) * dx(i) for i, gi in enumerate(gGamma, 1))
A, b = map(assemble, (a, L))
vh_P = Function(Q) # Since we project
solve(A, vh_P.vector(), b)
vh_I = subdomain_interpolate(zip(gGamma, subdomains), Q)
vh_I.vector().axpy(params.eps, ph.vector())
# Simply by interpolation
return (sqrt(H1_norm(u_exact[0], u1h) ** 2 + H1_norm(u_exact[1], uh) ** 2),
L2_norm(p_exact, vh),
L2_norm(p_exact, vh_P),
L2_norm(p_exact, vh_I))
error_types = ('|u|_1', '|v|_0', '|v|_{0P}', '|v|_{0I}')
return get_error, error_types
# --------------------------------------------------------------------
# How is the problem parametrized
PARAMETERS = ('kappa', 'eps')
I’d be extremely grateful for any hints and assistance. Currently, I’m working on a semiconductor device simulation project based on FEniCS. Those who are interested are also welcome to communicate and cooperate with me.
Thank you again.