Hi again,
Here is a MWE that could work if it was possible to integrate over the edge of a surface of a volume domain. Actually the geometry is the one I use, and might not be relevant for everyone.. And I hope this is not too much case specific. The error provided by this code comes after.
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
import ufl
from ufl import (
TestFunction,
TrialFunction,
Measure,
inner,
grad,
)
from basix.ufl import element
from dolfinx.fem import functionspace, form, petsc
import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio
import gmsh
from mpi4py import MPI
import dolfinx.io.gmshio as gmshio
import dolfinx.mesh as dmesh
import numpy as np
def cubic_domain(side_box = 0.11, radius = 0.1, lc = 8e-3, model_name = "Cubic"):
gmsh.initialize()
comm = MPI.COMM_WORLD
model_rank = 0
model = gmsh.model()
gmsh.model.add(model_name)
# Definition of the points
p1 = gmsh.model.geo.addPoint(0, 0, 0, lc)
p2 = gmsh.model.geo.addPoint(side_box, 0, 0, lc)
p3 = gmsh.model.geo.addPoint(side_box, side_box, 0, lc)
p4 = gmsh.model.geo.addPoint(0, side_box, 0, lc)
p5 = gmsh.model.geo.addPoint(0, 0, side_box, lc)
p6 = gmsh.model.geo.addPoint(side_box, 0, side_box, lc)
p7 = gmsh.model.geo.addPoint(side_box, side_box, side_box, lc)
p8 = gmsh.model.geo.addPoint(0, side_box, side_box, lc)
p9 = gmsh.model.geo.addPoint(side_box, radius, 0, lc)
p10 = gmsh.model.geo.addPoint(side_box, 0, radius, lc)
# Definition of the lines
l1 = gmsh.model.geo.addLine(p1, p4)
l2 = gmsh.model.geo.addLine(p4, p3)
l3 = gmsh.model.geo.addLine(p3, p9)
l4 = gmsh.model.geo.addLine(p9, p2)
l5 = gmsh.model.geo.addLine(p2, p1)
l6 = gmsh.model.geo.addLine(p5, p6)
l7 = gmsh.model.geo.addLine(p6, p7)
l8 = gmsh.model.geo.addLine(p7, p8)
l9 = gmsh.model.geo.addLine(p8, p5)
l10 = gmsh.model.geo.addLine(p2, p10)
l11 = gmsh.model.geo.addLine(p10, p6)
l12 = gmsh.model.geo.addLine(p3, p7)
l13 = gmsh.model.geo.addLine(p4, p8)
l14 = gmsh.model.geo.addLine(p1, p5)
# definition of the quarter circle
c15 = gmsh.model.geo.addCircleArc(p9, p2, p10)
# Curve loops
cl1 = [l1, l2, l3, l4, l5]
cl2 = [l6, l7, l8, l9]
cl3 = [-l3, l12, -l7, -l11, -c15]
cl4 = [c15, -l10, -l4]
cl5 = [-l2, l13, -l8, -l12]
cl6 = [-l5, l10, l11, -l6, -l14]
cl7 = [-l1, l14, -l9, -l13]
# Surfaces
s1 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl1)])
s2 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl2)])
s3 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl3)])
s4 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl4)])
s5 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl5)])
s6 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl6)])
s7 = gmsh.model.geo.addPlaneSurface([gmsh.model.geo.addCurveLoop(cl7)])
# Surface loop
sl1 = gmsh.model.geo.addSurfaceLoop([s1, s2, s5, s3, s4, s6, s7])
# Definition of the volume
v1 = gmsh.model.geo.addVolume([sl1])
gmsh.model.geo.synchronize()
# Definition of the physical identities
gmsh.model.addPhysicalGroup(1, [l12, -l7, -l6, -l14, l1, l2], tag=111)
gmsh.model.addPhysicalGroup(2, [s4], tag=1) # Neumann
gmsh.model.addPhysicalGroup(2, [s7, s5, s2], tag=3)
gmsh.model.addPhysicalGroup(3, [v1], tag=1) # Volume physique
# Generation of the 3D mesh
gmsh.model.mesh.generate(3)
# Save the mesh
gmsh.write(model_name + ".msh")
# Affect the mesh, volumic tags and surfacic tag to variables
mesh_data = gmshio.model_to_mesh(model, comm, model_rank)
#print(f'mesh_data : {mesh_data}')
final_mesh = mesh_data.mesh
cell_tags = mesh_data.cell_tags
facet_tags = mesh_data.facet_tags
ridge_tags = mesh_data.ridge_tags
# Close gmsh
gmsh.finalize()
tdim = final_mesh.topology.dim
fdim = tdim - 1
xref = [side_box, 0, 0]
submesh, entity_map = dmesh.create_submesh(final_mesh, fdim, facet_tags.find(3))[0:2]
entity_maps_mesh = [entity_map]
#mesh_info = [final_mesh, cell_tags, facet_tags, xref]
mesh_info = [final_mesh, cell_tags, facet_tags, xref, ridge_tags] # this line is new
submesh_info = [submesh, entity_maps_mesh]
return mesh_info, submesh_info
# Example usage:
if __name__ == "__main__":
# Generate a cube mesh and retrieve mesh data
mesh_info, submesh_info = cubic_domain()
final_mesh = mesh_info[0]
cell_tags = mesh_info[1]
facet_tags = mesh_info[2]
xref = mesh_info[3]
ridge_tags = mesh_info[4]
submesh = submesh_info[0]
entity_maps_mesh = submesh_info[1]
# Define volume (dx) and surface (ds) measures
dx = Measure("dx", domain=final_mesh)
ds = Measure("ds", domain=final_mesh, subdomain_data=facet_tags)
dr = Measure("dr", domain=final_mesh, subdomain_data=ridge_tags)
dx1 = Measure("dx", domain=submesh)
ds1 = Measure("ds", domain=submesh)
# Define a Lagrange element of degree 3
P1 = element("Lagrange", final_mesh.basix_cell(), 2)
P = functionspace(final_mesh, P1)
Q1 = element("Lagrange", submesh.basix_cell(), 2)
Q = functionspace(submesh, Q1)
# Define trial and test functions
p, v = TrialFunction(P), TestFunction(P)
u, w = TrialFunction(Q), TestFunction(Q)
# Define the weak form
a_00_0 = inner(p, v) * dx
a_01_0 = inner(u, v) * ds(3)
a_10_0 = inner(grad(p)[0] +grad(p)[1] + grad(p)[2], w) * dr(111)
a_10_1 = inner(grad(p), grad(w)) * ds(3)
a_11_0 = inner(grad(u)[0] +grad(u)[1] + grad(u)[2], w) * ds1
a_11_1 = inner(grad(u), grad(w)) * dx1
# Assemble the total form
A_00 = form(a_00_0)
A_01 = form(a_01_0, entity_maps=entity_maps_mesh)
A_10 = form(a_10_0 + a_10_1, entity_maps=entity_maps_mesh)
A_11 = form(a_11_0 + a_11_0)
print("end")
The error :
File "/Users/pierremariotti/Documents/PhD/FEniCSx/main-branch/root/ABC_FEniCSx_classical/mwe_FEniCSx.py", line 179, in <module>
A_10 = form(a_10_0 + a_10_1, entity_maps=entity_maps_mesh)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 438, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 430, in _create_form
return _form(form)
^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/fem/forms.py", line 339, in _form
ufcx_form, module, code = jit.ffcx_jit(
^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/jit.py", line 61, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-dolfinx-main-iqvj2zu4tzntqaicsg2qyaehn6q5jb4r/lib/python3.12/site-packages/dolfinx/jit.py", line 216, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 226, in compile_forms
raise e
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 206, in compile_forms
impl = _compile_objects(
^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/codegeneration/jit.py", line 331, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/ir/representation.py", line 189, in compute_ir
_compute_integral_ir(
File "/Users/pierremariotti/spack/opt/spack/darwin-m1/py-fenics-ffcx-main-5i36ozqb5qxae2mt5gy3l4ovafzniqlc/lib/python3.12/site-packages/ffcx/ir/representation.py", line 262, in _compute_integral_ir
entity_type = _entity_types[itg_data.integral_type]
~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^
KeyError: 'ridge'