Dear community,
I face a problem that I cannot explain scientifically, so either I’m wrong (highly possible) or something is wrong in my implementation. My problem is quite large to provide a relevant MWE, so I’m going to start small by a non-physical problem that includes my concern.
Let \Omega be a domain and \Gamma the boundary closed boundary of \Omega, e.g \Omega is a sphere. Let suppose the following weak formulation
(p,q)\in(\Omega\times\Gamma)
(v,u)\in(L_2(\Omega)\times L_2(\Gamma))
\begin{cases}
\int_{\Omega}\nabla p\cdot \nabla vd\Omega +\int_\Gamma q\cdot vd\Gamma = \int_\Omega g\cdot vd\Omega
\\
\int_{\Gamma}\nabla p\cdot \nabla ud\Gamma+\int_\Gamma q\cdot u d\Gamma=0
\end{cases}
My concern is about the implementation of the term \int_{\Gamma}\nabla p\cdot \nabla ud\Gamma.
Here is my full implementation leading to no error, followed by my configuration :
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
import dolfinx.mesh as msh
import numpy as np
from ufl import (TestFunction, TrialFunction, TestFunctions, TrialFunctions,
Measure,
FacetNormal, CellNormal,
inner, grad, div,
Dn,
MixedFunctionSpace,
extract_blocks)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc
def create_sphere_mesh(R=0.3, lc=0.03):
model_name = "sphere"
gmsh.initialize()
comm = MPI.COMM_WORLD
model_rank = 0
model = gmsh.model()
gmsh.model.add(model_name)
# Create a sphere of radius R centered at the origin
gmsh.model.occ.addSphere(0, 0, 0, R)
gmsh.model.occ.synchronize()
# Physical groups: volume (3D) and external surface (2D)
gmsh.model.addPhysicalGroup(3, [1], tag=1)
gmsh.model.addPhysicalGroup(2, [1], tag=2) # external surface
# Mesh generation settings
# gmsh.option.setNumber("Mesh.CharacteristicLengthMin", lc) # optional
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
# Save the mesh to file
gmsh.write(model_name + ".msh")
# Import mesh into DOLFINx
mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
final_mesh = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
print(f'facet_tags.indices : {facet_tags.indices}') # boundary entity indices
print(f'facet_tags.values : {facet_tags.values}') # corresponding tags
# Close gmsh
gmsh.finalize()
tdim = final_mesh.topology.dim
fdim = tdim - 1
#xref = [-0.07, 0.007, 0.007]
xref = [0.0, 0.2, 0.0]
submesh, entity_map = msh.create_submesh(final_mesh, fdim, facet_tags.find(2))[0:2]
entity_maps_mesh = [entity_map]
mesh_info = [final_mesh, cell_tags, facet_tags, xref] # this line is new
submesh_info = [submesh, entity_maps_mesh]
return mesh_info, submesh_info
mesh_info, submesh_info = create_sphere_mesh(R=0.3, lc=0.03)
final_mesh, cell_tags, facet_tags, xref = mesh_info
submesh, entity_maps_mesh = submesh_info
dx = Measure("dx", domain=final_mesh)
ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)
dx1 = Measure("dx", domain=submesh)
# Define a Lagrange element of degree 2
P1 = element("Lagrange", final_mesh.basix_cell(), 2)
P = functionspace(final_mesh, P1)
Q1 = element("Lagrange", submesh.basix_cell(), 2)
Q = functionspace(submesh, Q1)
p, v = TrialFunction(P), TestFunction(P)
q, u = TrialFunction(Q), TestFunction(Q)
# Define the weak form
a_00 = inner(p, v) * dx
a_01 = inner(q, v) * ds(2)
a_10 = inner(grad(p), grad(u)) * ds(2)
a_11 = inner(q, u) * dx1
# Assemble the total form
A_00 = form(a_00)
A_01 = form(a_01, entity_maps=entity_maps_mesh)
A_10 = form(a_10, entity_maps=entity_maps_mesh)
A_11 = form(a_11)
If it’s relevant this is my configuration (installed with spack) :
========== Versions FEniCSx ==========
dolfinx version : 0.10.0.0
ufl version : 2025.2.0.dev0
basix version : 0.10.0.0
mpi4py version : 4.0.1
======================================
I don’t really know how to present you the results I face without giving you thousands of lines of code, but it looks like the line
a_10 = inner(grad(p), grad(u)) * ds(2)
leads to a term equals to 0…
Would that trigger anything from you ? have you ever faced something like this ? would it be possible that the following line is wrongly expressed ?
a_10 = inner(grad(p), grad(u)) * ds(2)
Do not hesitate to ask for any details.
Thank you,
Pierre.