Dear FEniCSx community,
I’m struggling with the computation of the curvature of a surface of a geometry. A great amount of work has been done so far, and I’ve been able to create a functioning code to compute the curvature of a sphere. However, I would need my code to be more general. In order to do so, the only thing I’m changing is the way of computing the normal of my surface.
Unfortunately, I haven’t figured out how to convert the dolfin code to a dolfinx one, that would be robust.
Here is my code :
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
import dolfinx.mesh as msh
from dolfinx import fem, plot
import ufl
import numpy as np
from dolfinx.fem import petsc
import pyvista
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.CharacteristicLengthMax", lc)
gmsh.model.mesh.generate(3)
model.mesh.setOrder(2)
# 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}')
print(f'facet_tags.values : {facet_tags.values}')
# Close gmsh
gmsh.finalize()
tdim = final_mesh.topology.dim
fdim = tdim - 1
xref = [0.0, 0.05, 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]
submesh_info = [submesh, entity_maps_mesh]
return mesh_info, submesh_info
def compute_normal_vector(submesh):
"""
Compute the normal vector field on a surface mesh.
For a 2D surface embedded in 3D, we compute the normal using the cross product
of tangent vectors derived from the parametrization.
Returns:
n_func: Function containing the normal vector at each point
"""
gdim = submesh.geometry.dim
tdim = submesh.topology.dim
V_vector = fem.functionspace(submesh, ("CG", 1, (gdim,)))
n_func = fem.Function(V_vector)
# Get mesh coordinates
x = ufl.SpatialCoordinate(submesh)
# For a 2D surface in 3D, we compute the normal as the cross product
# of the tangent vectors obtained from the Jacobian
# The cell normal can be computed using the cross product of coordinate gradients
# Create DG0 space to compute cell-wise normals
V_DG = fem.functionspace(submesh, ("DG", 0, (gdim,)))
n_dg = fem.Function(V_DG)
# Compute the Jacobian of the mapping
J = ufl.Jacobian(submesh)
# For a 2D manifold in 3D, J is a 3x2 matrix
# The two columns are tangent vectors to the surface
# Normal = tangent1 × tangent2
if tdim == 2 and gdim == 3:
t1 = ufl.as_vector([J[0, 0], J[1, 0], J[2, 0]])
t2 = ufl.as_vector([J[0, 1], J[1, 1], J[2, 1]])
# Cross product
n_cross = ufl.as_vector([
t1[1] * t2[2] - t1[2] * t2[1],
t1[2] * t2[0] - t1[0] * t2[2],
t1[0] * t2[1] - t1[1] * t2[0]
])
# Normalize
norm_n = ufl.sqrt(ufl.dot(n_cross, n_cross))
n_normalized = n_cross / norm_n
else:
raise ValueError(f"Expected 2D surface in 3D space, got tdim={tdim}, gdim={gdim}")
# Project to DG0 space first
u_dg = ufl.TrialFunction(V_DG)
v_dg = ufl.TestFunction(V_DG)
a_dg = ufl.inner(u_dg, v_dg) * ufl.dx
L_dg = ufl.inner(n_normalized, v_dg) * ufl.dx
problem_dg = petsc.LinearProblem(a_dg, L_dg, u=n_dg, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type":
"mumps"}, petsc_options_prefix="dg")
problem_dg.solve()
# Now interpolate from DG0 to CG1 (averaging at nodes)
u_cg = ufl.TrialFunction(V_vector)
v_cg = ufl.TestFunction(V_vector)
a_cg = ufl.inner(u_cg, v_cg) * ufl.dx
L_cg = ufl.inner(n_dg, v_cg) * ufl.dx
problem_cg = petsc.LinearProblem(a_cg, L_cg, u=n_func, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type":
"mumps"}, petsc_options_prefix="cg")
problem_cg.solve()
# Normalize the normal vector at each node
n_values = n_func.x.array.reshape(-1, gdim)
norms = np.sqrt(np.sum(n_values**2, axis=1))
norms[norms < 1e-10] = 1.0
n_values /= norms[:, np.newaxis]
n_func.x.array[:] = n_values.flatten()
return n_func
def harry_plotter(space, sol, str_value, show_edges = True):
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(space)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data[str_value] = sol.x.array.real
u_grid.set_active_scalars(str_value)
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=show_edges)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
u_plotter.show()
def compute_mean_curvature(submesh, n_func=None):
"""
Calculate mean curvature H = div_s(n) on an arbitrary surface.
The surface divergence is: div_s(n) = div(n) - n^T · (∇n) · n
Parameters:
submesh: The surface mesh (2D manifold in 3D space)
n_func: Optional pre-computed normal vector. If None, it will be computed.
Returns:
H_func: Function containing the mean curvature
n_func: Function containing the normal vector
"""
gdim = submesh.geometry.dim
V_scalar = fem.functionspace(submesh, ("CG", 1))
# Compute normal vector if not provided
if n_func is None:
print("Computing normal vector...")
n_func = compute_normal_vector(submesh)
# Calculate mean curvature H = div_s(n)
# Surface divergence: div_s(n) = trace((I - n⊗n) · ∇n)
# This equals: div(n) - (n · ∇n) · n
H_func = fem.Function(V_scalar)
# Test and trial functions
u = ufl.TrialFunction(V_scalar)
v = ufl.TestFunction(V_scalar)
# Compute gradient of normal vector
grad_n = ufl.grad(n_func)
# Surface divergence using the formula
# div_s(n) = div(n) - n^T · grad(n) · n
# where div(n) = trace(grad(n))
div_n = ufl.tr(grad_n)
normal_derivative = ufl.dot(n_func, ufl.dot(grad_n, n_func))
div_s_n = div_n - normal_derivative
# Variational form: project div_s(n) onto V_scalar
a = ufl.inner(u, v) * ufl.dx
L = ufl.inner(div_s_n, v) * ufl.dx
# Solve the linear system
problem = petsc.LinearProblem(a, L, u=H_func, petsc_options={"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_mat_solver_type":
"mumps"}, petsc_options_prefix="kappa")
problem.solve()
# Calculate statistics
H_values = H_func.x.array.real
H_mean = np.mean(H_values)
H_std = np.std(H_values)
H_min = np.min(H_values)
H_max = np.max(H_values)
harry_plotter(V_scalar, H_func, 'kappa', show_edges = True)
print("\n" + "="*60)
print("MEAN CURVATURE RESULTS")
print("="*60)
print(f"Average mean curvature: {H_mean:.6f}")
print(f"Standard deviation: {H_std:.6f}")
print(f"Min value: {H_min:.6f}")
print(f"Max value: {H_max:.6f}")
print("="*60 + "\n")
return H_func, n_func
if __name__ == "__main__":
# Create the curved cubic domain
print("Creating curved cubic mesh...")
#mesh_info, submesh_info = curvedcubic_domain(side_box=0.11, radius=0.1, lc=8e-3)
mesh_info, submesh_info = create_sphere_mesh(R=0.3, lc=0.02)
submesh = submesh_info[0]
print(f"Submesh (surface) has {submesh.topology.index_map(2).size_local} triangles")
print(f"Submesh (surface) has {submesh.topology.index_map(0).size_local} vertices")
# Compute mean curvature
print("\nComputing mean curvature on the surface...")
H_func, n_func = compute_mean_curvature(submesh)
If the n_func is defined thanks to
x = ufl.SpatialCoordinate(submesh)
r = ufl.sqrt(x[0]**2 + x[1]**2 + x[2]**2)
n_analytical = ufl.as_vector([x[0]/r, x[1]/r, x[2]/r])
n_func = fem.Function(V_vector)
n_expr = fem.Expression(n_analytical, V_vector.element.interpolation_points)
n_func.interpolate(n_expr)
then there is no problem. But this is specific to a sphere…
The first long code can be executed without error, and it helps at visualizing the value of the curvature on the structure too.
I give you the outputs:
============================================================
MEAN CURVATURE RESULTS
============================================================
Average mean curvature: -0.058697
Standard deviation: 16.334957
Min value: -63.747906
Max value: 71.704344
============================================================
Many thanks for the help.
All the best,
Pierre.