Skip to content

Instantly share code, notes, and snippets.

@finsberg

finsberg/create_mesh.py

Last active August 15, 2025 13:34
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save finsberg/d219d8c8d1d057f984bada3ced710222 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save finsberg/d219d8c8d1d057f984bada3ced710222 to your computer and use it in GitHub Desktop.
Download ZIP
Gradient precision - dolfinx_adjoint vs dolfin_adjoint
Raw
create_mesh.py
# Run in dolfinx-nightly container
from pathlib import Path
from mpi4py import MPI
import dolfinx
import gmsh
def generate_mesh(
path: Path = Path("annulus.msh"),
c=(0.0, 0.0, 0.0), # origin
r_outer: float = 50,
r_inner: float = 10,
char_length=10.0,
verbose=True,
):
gmsh.initialize()
if not verbose:
gmsh.option.setNumber("General.Verbosity", 0)
gmsh.model.add("annulus")
outer_id = gmsh.model.occ.addSphere(*c, r_outer)
inner_id = gmsh.model.occ.addSphere(*c, r_inner)
outer = [(3, outer_id)]
inner = [(3, inner_id)]
annulus, _ = gmsh.model.occ.cut(outer, inner, removeTool=False)
surfaces = gmsh.model.occ.getEntities(dim=2)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(
dim=surfaces[0][0],
tags=[surfaces[0][1]],
tag=1,
name="INNER",
)
gmsh.model.addPhysicalGroup(
dim=surfaces[1][0],
tags=[surfaces[1][1]],
tag=2,
name="OUTER",
)
gmsh.model.add_physical_group(
dim=3, tags=[t[1] for t in annulus], tag=3, name="WALL"
)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", char_length)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", char_length)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Netgen")
gmsh.write(path.as_posix())
gmsh.finalize()
def main():
c = (0.0, 0.0, 0.0)
r_outer: float = 50
r_inner: float = 10
char_length = 15.0
msh_file = Path("annulus.msh")
comm = MPI.COMM_WORLD
if not msh_file.exists() and comm.rank == 0:
generate_mesh(
msh_file,
r_inner=r_inner,
r_outer=r_outer,
c=c,
verbose=False,
char_length=char_length,
)
comm.barrier()
msh = dolfinx.io.gmshio.read_from_msh(msh_file, comm=comm)
msh.facet_tags.name = "facet_tags"
with dolfinx.io.XDMFFile(comm, "mesh.xdmf", "w") as xdmf:
xdmf.write_mesh(msh.mesh)
xdmf.write_meshtags(msh.facet_tags, msh.mesh.geometry)
if __name__ == "__main__":
main()
Raw
main_dolfin.py
# Use container ghcr.io/scientificcomputing/fenics-gmsh:2024-05-30
# Run pip install dolfin-adjoint
from pathlib import Path
from dolfin import *
from dolfin_adjoint import *
import numpy as np
import ufl_legacy as ufl
comm = MPI.comm_world
OUTER = 2
INNER = 1
# ----- Read mesh -----
mesh = Mesh(comm)
with XDMFFile("mesh.xdmf") as xdmf:
xdmf.read(mesh)
ffun_val = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mesh.xdmf") as xdmf:
xdmf.read(ffun_val, "facet_tags")
ffun = MeshFunction("size_t", mesh, ffun_val)
# ----- Define measure -----
quad_degree = 6
dx = dx(mesh, metadata={"quadrature_degree": quad_degree})
ds = ds(
domain=mesh,
subdomain_data=ffun,
metadata={"quadrature_degree": quad_degree},
)
# ------ Define control --------
# V_control = FunctionSpace(mesh, "R", 0)
V_control = FunctionSpace(mesh, "DG", 0)
v_control = Function(V_control, name="Control")
v_control.vector()[:] = 0.5
p_control = v_control
mu = Constant(500.0, name="lambda")
lmbda = Constant(10_000.0, name="lambda")
# ------ Define variational problem --------
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
I = Identity(3)
eps = lambda u: 0.5 * (ufl.grad(u) + grad(u).T)
n = FacetNormal(mesh)
sigma = 2 * mu * eps(u) + lmbda * div(u) * I
a = inner(sigma, eps(v)) * dx
L = inner(-p_control * 1000 * n, v) * ds(INNER)
# ------ Define Dirichlet BC --------
bcs = DirichletBC(V, (0.0, 0.0, 0.0), ffun, OUTER)
# ------ Solve problem --------
uh = Function(V, name="Displacement Field")
solve(a == L, uh, bcs=bcs, solver_parameters={"linear_solver": "mumps"})
# ---- Set up cost function ----
J_symbolic = 0.5 * inner(uh, uh) * dx
J = assemble(J_symbolic)
control = Control(v_control)
Jhat = ReducedFunctional(J, control)
dJ = Jhat.derivative()
print(dJ.vector().get_local().max())
# 96.76507197055057
print(np.linalg.norm(dJ.vector().get_local()))
# 262.0681691322662
Raw
main_dolfinx.py
# In dolfinx nightly container, run pip install git+https://github.com/jorgensd/dxa.git
from pathlib import Path
from mpi4py import MPI
import dolfinx
import numpy as np
import pyadjoint
# import scifem
import ufl
import dolfinx_adjoint
comm = MPI.COMM_WORLD
OUTER = 2
INNER = 1
# ----- Read mesh -----
with dolfinx.io.XDMFFile(comm, "mesh.xdmf", "r") as xdmf:
mesh = xdmf.read_mesh()
mesh.topology.create_connectivity(2, 3)
facet_tags = xdmf.read_meshtags(mesh, "facet_tags")
# ----- Define measure -----
quad_degree = 6
dx = ufl.dx(mesh, metadata={"quadrature_degree": quad_degree})
ds = ufl.ds(
domain=mesh,
subdomain_data=facet_tags,
metadata={"quadrature_degree": quad_degree},
)
# ------ Define control --------
# V_control = scifem.create_real_functionspace(mesh, value_shape=())
V_control = dolfinx.fem.functionspace(mesh, ("DG", 0))
v_control = dolfinx_adjoint.Function(V_control, name="Control")
v_control.x.array[:] = 0.5
p_control = v_control
mu = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(500))
lmbda = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(10_000.0))
# ------ Define variational problem --------
V = dolfinx.fem.functionspace(mesh, ("CG", 1, (3,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
I = ufl.Identity(3)
eps = lambda u: 0.5 * (ufl.grad(u) + ufl.grad(u).T)
n = ufl.FacetNormal(mesh)
sigma = 2 * mu * eps(u) + lmbda * ufl.div(u) * I
a = ufl.inner(sigma, eps(v)) * dx
L = ufl.inner(-p_control * 1000 * n, v) * ds(INNER)
# ------ Define Dirichlet BC --------
outer_dofs = dolfinx.fem.locate_dofs_topological(
V, facet_tags.dim, facet_tags.find(OUTER)
)
bcs = [dolfinx.fem.dirichletbc(np.zeros(3), outer_dofs, V)]
# ------ Solve problem --------
uh = dolfinx_adjoint.Function(V, name="Displacement Field")
problem = dolfinx_adjoint.LinearProblem(
a,
L,
u=uh,
bcs=bcs,
petsc_options={
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type": "mumps",
},
)
problem.solve()
# ---- Set up cost function ----
J_symbolic = 0.5 * ufl.inner(uh, uh) * dx
J = dolfinx_adjoint.assemble_scalar(J_symbolic)
control = pyadjoint.Control(v_control)
Jhat = pyadjoint.ReducedFunctional(J, control)
dJ = Jhat.derivative()
print(dJ.x.array.max())
# 96.76482374990631
print(np.linalg.norm(dJ.x.array))
# 262.0680227282903
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment