Add manifold example with real spaces

[jorgensd]

Code based on @meg-simula's paper, now using GMSH to generate higher order grids:

from mpi4py import MPI
import gmsh
import numpy as np
import ufl
import dolfinx
import scifem

order = 3
exact_case = 2
res = 0.2

gmsh.initialize()
outer_sphere = gmsh.model.occ.addSphere(0, 0, 0, 1)
gmsh.model.occ.synchronize()
boundary = gmsh.model.getBoundary([(3, outer_sphere)], recursive=False)
gmsh.model.addPhysicalGroup(boundary[0][0], [boundary[0][1]], tag=1)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", res)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(order)

mesh_data = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)


mesh = mesh_data.mesh


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


V = dolfinx.fem.functionspace(mesh, ("RT", order))
Q = dolfinx.fem.functionspace(mesh, ("DG", order - 1))
R = scifem.create_real_functionspace(mesh)
W = ufl.MixedFunctionSpace(*(V, Q, R))

(sigma, u, r) = ufl.TrialFunctions(W)
(tau, v, t) = ufl.TestFunctions(W)

global_orientation = ufl.sign(ufl.dot(x, ufl.CellNormal(mesh)))
sigma_ = global_orientation * sigma
tau_ = global_orientation * tau
# Choose exact solution based on user preference
if exact_case == 1:
    g = x[2]
    u_exact = -0.5 * x[2]
elif exact_case == 2:
    g = x[0] * x[1] * x[2]
    u_exact = -x[0] * x[1] * x[2] / 12
    sigma_exact = ufl.as_vector(
        (
            -(
                x[1] * x[2]
                - 3
                * x[0]
                * x[0]
                * x[1]
                * x[2]
                / (x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
            )
            / 12,
            -(
                x[0] * x[2]
                - 3
                * x[0]
                * x[1]
                * x[1]
                * x[2]
                / (x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
            )
            / 12,
            -(
                x[0] * x[1]
                - 3
                * x[0]
                * x[1]
                * x[2]
                * x[2]
                / (x[0] * x[0] + x[1] * x[1] + x[2] * x[2])
            )
            / 12,
        )
    )
elif exact_case == 3:
    g = 3 * x[2] * x[2] - 1
    u_exact = -(3 * x[2] * x[2] - 1) / 6
else:
    raise Exception("Unrecognized exact case (%d)" % exact_case)

a = (
    ufl.inner(sigma_, tau_) + ufl.div(sigma_) * v + ufl.div(tau_) * u + r * v + t * u
) * ufl.dx
L = [ufl.ZeroBaseForm((tau,)), g * v * ufl.dx, ufl.ZeroBaseForm((t,))]

petsc_options = {
    "ksp_type": "preonly",
    "pc_type": "lu",
    "pc_factor_mat_solver_type": "mumps",
    "ksp_error_if_not_converged": True,
}
problem = dolfinx.fem.petsc.LinearProblem(
    ufl.extract_blocks(a),
    L,
    bcs=[],
    petsc_options=petsc_options,
    petsc_options_prefix="mixed_poisson_",
    kind="mpi",
)
sigma_h, u_h, r_h = problem.solve()

V_out = dolfinx.fem.functionspace(mesh, ("DG", order))
v_out = dolfinx.fem.Function(V_out)
v_out.interpolate(u_h)
v_out.name = "u_h"
with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [v_out]) as vtx_writer:
    vtx_writer.write(0.0)

L2_error_u = dolfinx.fem.form(ufl.inner(u_h - u_exact, u_h - u_exact) * ufl.dx)
sigma_approx = global_orientation * sigma_h
L2_error_sigma = dolfinx.fem.form(
    ufl.inner(sigma_approx - sigma_exact, sigma_approx - sigma_exact) * ufl.dx
)
Hdiv_error_sigma = dolfinx.fem.form(
    ufl.inner(ufl.div(sigma_approx) - g, ufl.div(sigma_approx) - g) * ufl.dx
    + ufl.inner(sigma_approx - sigma_exact, sigma_approx - sigma_exact) * ufl.dx
)

E_u = np.sqrt(mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error_u), op=MPI.SUM))
E_sigma = np.sqrt(
    mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_error_sigma), op=MPI.SUM)
)
E_sigma_hdiv = np.sqrt(
    mesh.comm.allreduce(dolfinx.fem.assemble_scalar(Hdiv_error_sigma), op=MPI.SUM)
)
print("L2 error in u:", E_u)
print("L2 error in sigma:", E_sigma)
print("H(div) error in sigma:", E_sigma_hdiv)

Output:

L2 error in u: 6.449031074513135e-07
L2 error in sigma: 1.8679131149781764e-06
H(div) error in sigma: 1.408483169816556e-05

Solution visualized in Paraview:

[jorgensd]

Maybe add the other example:

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)

based of
https://fenicsproject.discourse.group/t/how-are-2d-manifolds-immersed-embedded-in-3d-space-handled-in-dolfinx/18066/6?u=dokken