Hello everyone!
I am trying to solve the following biharmonic problem with FEniCSx 0.9 and dolfinx_mpc with a \nabla u\cdot n=0 condition on interior facets \Gamma_{I}:
\begin{align}
\Delta^2\,u &= f &\text{on}\quad\Omega\tag{1.1}\\
u &= u_0 &\text{on}\quad\partial\Omega\tag{1.2}\\
\nabla u\cdot n &= g_0 &\text{on}\quad\partial\Omega\tag{1.3}\\
u &= u_I &\text{on}\quad\Gamma_I\tag{1.4}\\
\nabla u\cdot n &= 0 &\text{on}\quad\Gamma_I\tag{1.5}
\end{align}
My geometry is given in the figure below. The bottom and top boundaries \Gamma_{YL} and \Gamma_{YH} are periodic such that \partial\Omega consists only of the left and right parts:
\begin{align}
\partial\Omega = \Gamma_{XL}\cup\Gamma_{XH}\tag{2}
\end{align}
The line \Gamma_I is located at x=0.75. The whole domain is a unit square with \Omega = [0,1]^2.
I chose a mixed finite element approach with the following definitions:
\begin{align}
\sigma &:= \nabla u\tag{3.1}\\
\alpha &:= \nabla\cdot\sigma\tag{3.2}
\end{align}
With this the problem in equation (1) transforms into:
\begin{align}
\Delta\alpha &= f &\text{on}\quad\Omega\tag{4.1}\\
\sigma - \nabla u&=0 &\text{on}\quad\Omega\tag{4.2}\\
\alpha - \nabla\cdot\sigma &= 0&\text{on}\quad\Omega\tag{4.3}\\
u &= u_0 &\text{on}\quad\partial\Omega\tag{4.4}\\
\sigma\cdot n &= g_0 &\text{on}\quad\partial\Omega\tag{4.5}\\
u &= u_I &\text{on}\quad\Gamma_I\tag{4.6}\\
\sigma\cdot n &= 0 &\text{on}\quad\Gamma_I\tag{4.7}
\end{align}
Multiplying equations (4.1), (4.2) and (4.3) with test functions \nu, \tau and \eta and integrating gives:
\begin{align}
\int\limits_\Omega\Delta\alpha\cdot\nu\;dx\ &= \int\limits_\Omega f\cdot\nu\;dx \tag{5.1}\\
\int\limits_\Omega\sigma\cdot\tau\;dx\ - \int\limits_\Omega\nabla u\cdot\tau\;dx\ &=\ 0 \tag{5.2}\\
\int\limits_\Omega\alpha\cdot\eta\;dx\ - \int\limits_\Omega(\nabla\cdot\sigma)\cdot\eta\;dx\ &=\ 0\tag{5.3}
\end{align}
Integrating equation (5.1) by parts gives:
\begin{align}
\int\limits_\Omega\nabla\alpha\cdot\nabla\nu\;dx\ - \int\limits_{\partial\Omega}(\nabla\alpha\cdot n)\cdot\nu\;ds\ &=\ -\int\limits_\Omega f\cdot\nu\;dx \tag{6}
\end{align}
where the second term on the lhs vanishes because of \nu|_{\partial\Omega}=0, such that we get:
\begin{align}
\int\limits_\Omega\nabla\alpha\cdot\nabla\nu\;dx\ &=\ -\int\limits_\Omega f\cdot\nu\;dx \tag{7}
\end{align}
Applying the divergence theorem on the second term of equation (5.3) transforms this equation into:
\begin{align}
\int\limits_\Omega\alpha\cdot\eta\;dx\ + \int\limits_\Omega\nabla\eta\cdot\sigma\;dx\ &= \int\limits_{\partial\Omega}(\sigma\cdot n)\cdot\eta\;ds\tag{8}
\end{align}
where the boundary condition from equation (4.5) can be inserted:
\begin{align}
\int\limits_\Omega\alpha\cdot\eta\;dx\ + \int\limits_\Omega\nabla\eta\cdot\sigma\;dx\ &= \int\limits_{\partial\Omega}g_0\cdot\eta\;ds\tag{9}
\end{align}
I did not touch equation (5.2) further. In summary (and using eq. (2)) the weak form of my mixed biharmonic problem looks as follows:
\begin{align}
\int\limits_\Omega\nabla\alpha\cdot\nabla\nu\;dx\ &=\ -\int\limits_\Omega f\cdot\nu\;dx \tag{10.1}\\
\int\limits_\Omega\sigma\cdot\tau\;dx\ - \int\limits_\Omega\nabla u\cdot\tau\;dx\ &=\ 0 \tag{10.2}\\
\int\limits_{\Omega}\alpha\cdot\eta\;dx\ + \int\limits_{\Omega}\nabla\eta\cdot\sigma\;dx\ &= \int\limits_{\Gamma_{XL}}g_0\cdot\eta\;ds + \int\limits_{\Gamma_{XH}}g_0\cdot\eta\;ds\tag{10.3}
\end{align}
where the dirichlet boundary conditions on u must be enforced strongly.
The problem now is that I want to enforce equations (4.6) and (4.7):
\begin{align}
u &= u_I &\text{on}\quad\Gamma_I\tag{4.6}\\
\sigma\cdot n &= 0 &\text{on}\quad\Gamma_I\tag{4.7}
\end{align}
For eq. (4.6) it works to just apply a dirichlet BC on the internal facets just like a normal dirichlet BC on the boundary. For eq. (4.7) I tried to use dolfinx_mpc and apply a slip constraint. However if I do that the result explodes and gives me values \sim10^{38}.
Can anybody help me to solve that problem? I do not care if it is using dolfinx_mpc or any other method. My goal is just to solve the biharmonic equation (1) with u=u_I and \nabla u\cdot n=0 on \Gamma_{I}.
My code is given in the following. It is not quite a MWE because the mesh generation is not very simple, so I splitted it into a separate file. To run everything you need to create the two files create_mesh.py and mixed_biharmonic_mpc.py in one folder and run them via the commands >>> python create_mesh.py and >>> python mixed_biharmonic_mpc.py.
In the mixed_biharmonic_mpc.py file you can activate/deactivate the slip constraint in line 41 by setting add_slip_on_sigma = True/False. When you inspect the generated xdmf result files in the results folder you will see the following:
Without slip condition (add_slip_on_sigma = False):
With slip condition (add_slip_on_sigma = True):
Here is the code for create_mesh.py:
import gmsh
import os
import pickle
def create_mesh(periodic: str = "",
name: str = "mesh",
folder: str = "mesh") -> None:
dim: int = gmsh.model.getDimension()
allowed_directions: set = {"X", "Y", "Z"} if dim == 3 else {"X", "Y"}
set_mesh_properties()
directions: set = set(periodic) & allowed_directions
for direction in directions:
set_periodicity(direction)
gmsh.model.mesh.generate(dim)
# Save mesh
if not os.path.exists(folder):
os.makedirs(folder)
gmsh.write(folder + "//" + name + ".msh")
# Save physical names
physical_names = get_physical_names_dict()
with open(folder + '//' + name + '.pkl', 'wb') as physical_file:
pickle.dump(physical_names, physical_file)
def create_geometry(H: float = 0.5,
h: float = 0.5,
l: float = 0.25) -> None:
gmsh.model.occ.addPoint(-H, -H, 0, tag = 1)
gmsh.model.occ.addPoint(-H, +H, 0, tag = 2)
gmsh.model.occ.addPoint(+l, +H, 0, tag = 3)
gmsh.model.occ.addPoint(+H, +H, 0, tag = 4)
gmsh.model.occ.addPoint(+H, -H, 0, tag = 5)
gmsh.model.occ.addPoint(+l, -H, 0, tag = 6)
gmsh.model.occ.addPoint(l, +h/2, 0, tag = 7)
gmsh.model.occ.addPoint(l, -h/2, 0, tag = 8)
gmsh.model.occ.synchronize()
all_points = gmsh.model.getEntities()
gmsh.model.occ.translate(all_points, H, H, 0)
gmsh.model.occ.synchronize()
gmsh.model.occ.addLine(1, 2, tag = 1)
gmsh.model.occ.addLine(2, 3, tag = 2)
gmsh.model.occ.addLine(3, 4, tag = 3)
gmsh.model.occ.addLine(4, 5, tag = 4)
gmsh.model.occ.addLine(5, 6, tag = 5)
gmsh.model.occ.addLine(6, 1, tag = 6)
gmsh.model.occ.addLine(3, 7, tag = 7)
gmsh.model.occ.addLine(7, 8, tag = 8)
gmsh.model.occ.addLine(8, 6, tag = 9)
tag_loop_left = gmsh.model.occ.addCurveLoop([1,2,7,8,9,6])
tag_loop_right = gmsh.model.occ.addCurveLoop([4,5,9,8,7,3])
gmsh.model.occ.addPlaneSurface([tag_loop_left], tag = 1)
gmsh.model.occ.addPlaneSurface([tag_loop_right], tag = 2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag = 1, name = "S_L")
gmsh.model.addPhysicalGroup(2, [2], tag = 2, name = "S_R")
gmsh.model.addPhysicalGroup(1, [1], tag = 1, name = "GAMMA_XL")
gmsh.model.addPhysicalGroup(1, [2,3], tag = 2, name = "GAMMA_YH")
gmsh.model.addPhysicalGroup(1, [4], tag = 3, name = "GAMMA_XH")
gmsh.model.addPhysicalGroup(1, [5,6], tag = 4, name = "GAMMA_YL")
gmsh.model.addPhysicalGroup(1, [8], tag = 5, name = "GAMMA_I")
gmsh.model.occ.synchronize()
def get_entity_tags_from_physical_groups(dim: int = 1) -> list[int]:
physical = gmsh.model.getPhysicalGroups(dim = dim)
entity_tags = []
for entity_dim_tag in physical:
entity_tags += list(gmsh.model.getEntitiesForPhysicalGroup(dim = dim, tag = entity_dim_tag[1]))
return entity_tags
def get_physical_names_dict(dim=-1) -> dict:
dim_tags = gmsh.model.getPhysicalGroups(dim=dim)
names = dict()
for dim_tag in dim_tags:
name = gmsh.model.getPhysicalName(*dim_tag)
names[name] = dim_tag[1]
return names
def set_mesh_properties(min_lenght: float = 0.002,
max_lenght: float = 0.05,
min_dist: float = 0.01,
max_dist: float = 0.15) -> None:
dim: int = gmsh.model.getDimension()
curves_or_surfaces = "SurfacesList" if dim == 3 else "CurvesList"
entities = get_entity_tags_from_physical_groups(dim = dim - 1)
distance = gmsh.model.mesh.field.add("Distance")
gmsh.model.mesh.field.setNumbers(distance, curves_or_surfaces, entities)
gmsh.model.mesh.field.setNumber(distance, "Sampling", 100)
threshold = gmsh.model.mesh.field.add("Threshold")
gmsh.model.mesh.field.setNumber(threshold, "InField", distance)
gmsh.model.mesh.field.setNumber(threshold, "SizeMin", min_lenght)
gmsh.model.mesh.field.setNumber(threshold, "SizeMax", max_lenght/2)
gmsh.model.mesh.field.setNumber(threshold, "DistMin", min_dist)
gmsh.model.mesh.field.setNumber(threshold, "DistMax", max_dist)
minimum = gmsh.model.mesh.field.add("Min")
gmsh.model.mesh.field.setNumbers(minimum, "FieldsList", [threshold])
gmsh.model.mesh.field.setAsBackgroundMesh(minimum)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
def get_sorted_entities(dim_tags: list[tuple], direction: str) -> list[tuple]:
dim: int = gmsh.model.getDimension()
tags: list[int] = [tag for (dim, tag) in dim_tags]
if direction == "X":
index_1 = 1
index_2 = 2
if direction == "Y":
index_1 = 0
index_2 = 2
if direction == "Z":
index_1 = 0
index_2 = 1
sorted_entities = []
for entity in tags:
com = gmsh.model.occ.getCenterOfMass(dim - 1, entity)
com = [com[index_1], com[index_2]] if dim == 3 else [com[index_1]]
sorted_entities.append((entity, com))
sorted_entities = sorted(sorted_entities, key = lambda t: t[1])
return sorted_entities
def get_translation(direction: str,
tags_min: list[tuple],
tags_max: list[tuple]) -> list[float]:
dim: int = gmsh.model.getDimension()
source = [tag[0] for tag in tags_min]
destin = [tag[0] for tag in tags_max]
dl = gmsh.model.occ.getDistance(dim - 1, tags_min[0][0], dim - 1, tags_max[0][0])[0]
if direction == "X":
translation = [1, 0, 0, dl,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1]
if direction == "Y":
translation = [1, 0, 0, 0,
0, 1, 0, dl,
0, 0, 1, 0,
0, 0, 0, 1]
if direction == "Z":
translation = [1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, dl,
0, 0, 0, 1]
return destin, source, translation
def set_periodicity(direction: str) -> None:
"""
Parameters
----------
direction : str
The direction in which periodicity should be applied. Must be X, Y or Z.
Raises
------
ValueError
direction must be X, Y or Z
Returns
-------
None.
"""
dim: int = gmsh.model.getDimension()
boundary_type: str = "S" if dim == 3 else "GAMMA"
try:
dim_tags_low: list[tuple] = gmsh.model.getEntitiesForPhysicalName(boundary_type + "_" + direction + "L")
dim_tags_high: list[tuple] = gmsh.model.getEntitiesForPhysicalName(boundary_type + "_" + direction + "H")
except:
raise ValueError("direction must be X, Y or Z")
tags_low_sorted = get_sorted_entities(dim_tags_low, direction = direction)
tags_high_sorted = get_sorted_entities(dim_tags_high, direction = direction)
destin, source, translation = get_translation(direction, tags_low_sorted, tags_high_sorted)
gmsh.model.mesh.setPeriodic(dim - 1, destin, source, translation)
gmsh.initialize()
create_geometry()
create_mesh(periodic = "XY")
gmsh.fltk.run()
gmsh.finalize()
Here is the code for mixed_biharmonic_mpc.py:
import pickle
from petsc4py import PETSc
import dolfinx
from mpi4py import MPI
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh, default_scalar_type, default_real_type
from ufl import Measure, SpatialCoordinate, TestFunctions, TrialFunctions, grad, div, exp, inner, FacetNormal
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps"
}
def get_physical_names(file_name: str = "mesh.pkl", folder: str = "mesh"):
with open(folder + "//" + file_name, 'rb') as f:
physical_names = pickle.load(f)
return physical_names
def periodic_boundary(x):
return np.isclose(x[1], 1, atol=250 * np.finfo(default_scalar_type).resolution)
def periodic_relation(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]
out_x[1] = 1 - x[1]
return out_x
# Change this to True if you want to compute the solution with slip condition.
# Unfortunately the = True case does not work currently.
# Solutions for True/False will be written to separate .xdmf files, so you can
# compare the cases and see that the slip condition does not work.
add_slip_on_sigma = False
# Read mesh and physical markers from gmsh
msh, cell_markers, facet_markers = io.gmshio.read_from_msh("mesh//mesh.msh", MPI.COMM_WORLD, gdim=2)
physical_names = get_physical_names(file_name = "mesh.pkl", folder = "mesh")
# Define mixed element space
VE = element("CG", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
SE = element("CG", msh.basix_cell(), 1, dtype=default_real_type)
ME = mixed_element([VE, SE, SE])
V = fem.functionspace(msh, ME)
# Subspaces
V0_sub = V.sub(0)
V1_sub = V.sub(1)
V2_sub = V.sub(2)
V0_sub_collapse, _ = V.sub(0).collapse()
V2_sub_collapse, _ = V.sub(2).collapse()
# Trial and Test functions
(sigma, alpha, u) = TrialFunctions(V)
(tau, eta, v) = TestFunctions(V)
# Definition of source and boundary values
x = SpatialCoordinate(msh)
dx_ = x[0] - 0.5
dy_ = x[1] - 0.9
f = 160 * exp(-(dx_ * dx_ + dy_ * dy_) / 0.02) # Source term
g_0 = 0.0 # Boundary value for inner(sigma, n)
u_0 = 0.0 # Boundary value for u
u_I = 0.0 # Value for u on internal line
# Define variational form
dx = Measure("dx", msh)
a = inner(grad(alpha), grad(v))*dx # eq. 10.1 (only lhs)
a += inner(sigma, tau)*dx - inner(grad(u), tau)*dx # eq. 10.2
a += inner(alpha, eta)*dx + inner(grad(eta), sigma)*dx # eq. 10.3
L = - inner(f, v) * dx # eq. 10.1 (only rhs)
# Neumann Boundary Conditions
ds = Measure("ds", msh, subdomain_data=facet_markers)
L += inner(g_0, eta)*ds(physical_names["GAMMA_XL"]) # eq. 10.3 (rhs term for XL)
L += inner(g_0, eta)*ds(physical_names["GAMMA_XH"]) # eq. 10.3 (rhs term for XH)
# Dirichlet Boundary Conditions
dofs_XL = fem.locate_dofs_topological((V2_sub, V2_sub_collapse), 1, facet_markers.find(physical_names["GAMMA_XL"]))
dofs_XH = fem.locate_dofs_topological((V2_sub, V2_sub_collapse), 1, facet_markers.find(physical_names["GAMMA_XH"]))
bc_XL = fem.dirichletbc(u_0, dofs_XL[0], V2_sub)
bc_XH = fem.dirichletbc(u_0, dofs_XH[0], V2_sub)
bcs = [bc_XL, bc_XH]
# Condition on internal line
dofs_I = fem.locate_dofs_topological((V2_sub, V2_sub_collapse), 1, facet_markers.find(physical_names["GAMMA_I"]))
bc_I = fem.dirichletbc(u_I, dofs_I[0], V2_sub)
bcs += [bc_I]
# Periodicity in y-direction:
mpc = MultiPointConstraint(V)
mpc.create_periodic_constraint_geometrical(V0_sub, periodic_boundary, periodic_relation, bcs) # sigma
mpc.create_periodic_constraint_geometrical(V1_sub, periodic_boundary, periodic_relation, bcs) # alpha
mpc.create_periodic_constraint_geometrical(V2_sub, periodic_boundary, periodic_relation, bcs) # u
# Internal slip condition: inner(sigma, n) = 0 <==> inner(grad(u), n) = 0
if add_slip_on_sigma:
n_space = V0_sub_collapse
normals = dolfinx.fem.Function(n_space)
# This does not work:
mpc.create_slip_constraint(V0_sub, (facet_markers, physical_names["GAMMA_I"]), normals, bcs=bcs)
mpc.finalize()
# Solve problem
problem = LinearProblem(a, L, mpc, bcs, petsc_options=petsc_options)
w = problem.solve()
sigma_s, alpha_s, u_s = w.split()
# Write solution files
if add_slip_on_sigma:
file_name = "u_slip"
else:
file_name = "u_no_slip"
with io.XDMFFile(msh.comm, "results//"+file_name+".xdmf", "w") as file:
file.write_mesh(msh)
u_s.name = "u"
file.write_function(u_s)
