How are 2D manifolds immersed/embedded in 3D space handled in dolfinx?

dolfinx
,
6

I’ve revived the code of Marie Rognes from 2013 (see supplementary material of: GMD - Automating the solution of PDEs on the sphere and other manifolds in FEniCS 1.2)

from mpi4py import MPI
import meshio
import ufl
import basix.ufl
import dolfinx
import scifem
import numpy as np

order = 2
tmp_data = meshio.dolfin.read("../hdiv-l2/meshes/sphere_ico3.xml")

mesh = dolfinx.mesh.create_mesh(
    comm=MPI.COMM_WORLD,
    cells=tmp_data.cells_dict["triangle"],
    x=tmp_data.points,
    e=ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(3,))),
)

V_exact = dolfinx.fem.functionspace(mesh, ("Lagrange", order + 1, (mesh.geometry.dim,)))

# Define global normal
x = ufl.SpatialCoordinate(mesh)


S = dolfinx.fem.functionspace(mesh, ("DG", order, (mesh.geometry.dim,)))
L = dolfinx.fem.functionspace(mesh, ("DG", order))
V = dolfinx.fem.functionspace(mesh, ("Lagrange", order + 1))
R = scifem.create_real_functionspace(mesh)
W = ufl.MixedFunctionSpace(*(S, L, V, R))

(sigma, l, u, r) = ufl.TrialFunctions(W)
(w, gamma, v, t) = ufl.TestFunctions(W)

g = x[0] * x[1] * x[2]
up = ufl.as_vector((x[0], x[1], x[2]))


a = (
    ufl.dot(sigma, w)
    - ufl.dot(ufl.grad(u), w)
    - ufl.dot(w, up) * l
    + ufl.dot(sigma, ufl.grad(v))
    + gamma * ufl.dot(sigma, up)
    + u * t
    + r * v
) * ufl.dx
L_blocked = [
    ufl.ZeroBaseForm((w,)),
    ufl.ZeroBaseForm((gamma,)),
    -v * g * ufl.dx,
    ufl.ZeroBaseForm((t,)),
]

sigma_h = dolfinx.fem.Function(S)
l_h = dolfinx.fem.Function(L)
u_h = dolfinx.fem.Function(V)
r_h = dolfinx.fem.Function(R)
wh = (sigma_h, l_h, u_h, r_h)

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
    "ksp_monitor": None,
}


problem = dolfinx.fem.petsc.LinearProblem(
    ufl.extract_blocks(a),
    L_blocked,
    u=wh,
    bcs=[],
    petsc_options=petsc_options,
    petsc_options_prefix="mixed_poisson_",
    kind="mpi",
)
problem.solve()

with dolfinx.io.XDMFFile(mesh.comm, "new_mesh.xdmf", "w") as xdmf_file:
    xdmf_file.write_mesh(mesh)


with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u_h]) as vtx_writer:
    vtx_writer.write(0.0)

u_exact = -x[0] * x[1] * x[2] / 12
print(dolfinx.fem.assemble_scalar(dolfinx.fem.form(u_exact**2 * ufl.dx)))
L2_error_u = dolfinx.fem.form(ufl.inner(u_h - u_exact, u_h - u_exact) * ufl.dx)
print(
    "L2 error in u:",
    np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error_u), op=MPI.SUM)),
)
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u_h]) as vtx_writer:
    vtx_writer.write(0.0)

yields:
L2 error in u: 0.0004837339846541759