How to get the area of a hemi sphere

mesh
1

Hello, I am trying to build a hemi sphere and calculate the area by following the tutorial Integration measures, but the code does not give the correct result. Here is my code.

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

R = 1.0
grid_size = 0.2
order = 2

gmsh.initialize()
gmsh.model.add("hemi_sphere")

sphere = gmsh.model.occ.add_sphere(0.0, 0.0, 0.0, R)
box = gmsh.model.occ.add_box(-R-0.1, -R-0.1, 0, 2*R+0.1, 2*R+0.1, R+0.1)

result = gmsh.model.occ.intersect([(3, sphere)], [(3, box)])
gmsh.model.occ.synchronize()

surfaces = gmsh.model.getEntities(dim=2)

for surface in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
    if np.allclose(com, [0, 0, 0]): 
        plane = surface
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], 10, "circle")
    elif np.allclose(com, [0.0, 0.0, R*0.5]):
        hemi_sphere = surface
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], 20, "hemi_sphere")

# plane should be used here. There is an extra curve in hemi_sphere
boundary = gmsh.model.getBoundary([(plane[0], plane[1])])
gmsh.model.addPhysicalGroup(1, [b[1] for b in boundary], 1, "boundary")

gmsh.option.setNumber("Mesh.CharacteristicLengthMax", grid_size)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(order)
mesh, cell_marker, facet_marker = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()

dx = ufl.Measure("dx", domain=mesh, subdomain_data=cell_marker)

area = 1 * dx(20)
compiled_area = dolfinx.fem.form(area)
local_area = dolfinx.fem.assemble_scalar(compiled_area)
global_area = mesh.comm.allreduce(local_area, op=MPI.SUM)

print(global_area)

If I integrate 1*dx(20), the result should be close to 2\pi, but the result is 3.14, it is uncorrect. If I integrate 1*dx, it is 6.28, but it should be 3\pi.
Does anyone have some idea ? Thanks.

2

Here you are setting the geometrical dimension to two, but it should be 3, as you are making a 2D manifold embedded in 3D space, i.e.

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

which gives the correct area.
Please note that this is super easy to debug, if you had written the mesh to file with:

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "hemi_sphere.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(cell_marker, mesh.geometry)

and looked at the output in paraview, as it would be:

image
image1763×1791 277 KB

2 Likes
3

Thanks for your help, may I ask another question ? I’m trying to calculate the integral of vector field curl on the hemi sphere. But my code produces the error:

ValueError: Cannot determine geometric dimension from expression.

The code to generate the hemi-sphere is same as above, here is the code that calculates the integral of curl.

dx = ufl.Measure("dx", domain=mesh, subdomain_data=cell_marker)
x = ufl.SpatialCoordinate(mesh)
F = ufl.as_vector([-x[1], x[0], 0.0])
cell_normal_vector = ufl.CellNormal(mesh)
# this line causes error.
form = dolfinx.fem.form(ufl.dot(ufl.curl(F), cell_normal_vector) * dx(20))
4

If you change your code

to:

zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
F = ufl.as_vector([-x[1], x[0], zero])
cell_normal_vector = ufl.CellNormal(mesh)
form = dolfinx.fem.form(ufl.dot(ufl.curl(F), cell_normal_vector) * dx(20))

it will run.
The problem is doing curl of 0., as it is not associated with a the mesh.

1 Like
5

Sorry to bother you again, but how should I calculate the integral of vector \vec{F} along the boundary ? (I am trying to verify Stokes theorem) Is there any function like FacetNormal or CellNormal to get the unit vector along the boundary? Thanks.

6

As you are working on a manifold, do you mean the boundary of each cell (a line), and thus you would like to know the tangent t\in \mathbb{R}^3 which is perpendicular to n=FacetNormal(mesh)?

You could try using : CellEdgeVectors: ufl package — Unified Form Language (UFL) 2021.1.0 documentation

7

Maybe not? I want the unit vectors along the boundary as the figures shows

屏幕截图 2025-05-30 080805
屏幕截图 2025-05-30 080805349×492 19.7 KB

I tried CellEdgeVectors, it does not work

# The line integral of F on the boundary ?
ds = ufl.Measure("ds", subdomain_data=facet_marker)
cell_edge_vector = ufl.classes.CellEdgeVectors(mesh)
form = dolfinx.fem.assemble_scalar(ufl.dot(cell_edge_vector, F) * ds )

The code produces the errors:

Traceback (most recent call last):
  File "/mnt/e/workdir/FEM/fenics/bb2.py", line 48, in <module>
    form = dolfinx.fem.assemble_scalar(ufl.dot(cell_edge_vector, F) * ufl.ds )
                                       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~
  File "/mnt/d/software_install/fenics-0.9/lib/python3.13/site-packages/ufl/measure.py", line 382, in __rmul__
    raise ValueError(
    ...<3 lines>...
    )
ValueError: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3,) and free indices with labels ().

Here is the complete code:

import gmsh 
import dolfinx, ufl
from mpi4py import MPI
import numpy as np
import ufl.classes 
R = 1.0
grid_size = 0.2
order = 2

gmsh.initialize()
gmsh.model.add("hemi_sphere")
sphere = gmsh.model.occ.add_sphere(0.0, 0.0, 0.0, R)
box = gmsh.model.occ.add_box(-R-0.1, -R-0.1, 0, 2*R+0.1, 2*R+0.1, R+0.1)
result = gmsh.model.occ.intersect([(3, sphere)], [(3, box)])
gmsh.model.occ.synchronize()

surfaces = gmsh.model.getEntities(dim=2)
for surface in surfaces:
    com = gmsh.model.occ.getCenterOfMass(surface[0], surface[1])
    if np.allclose(com, [0, 0, 0]): # 平面
        plane = surface
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], 10, "circle")
    elif np.allclose(com, [0.0, 0.0, R*0.5]):
        hemi_sphere = surface
        gmsh.model.addPhysicalGroup(surface[0], [surface[1]], 20, "hemi_sphere")

boundary = gmsh.model.getBoundary([(plane[0], plane[1])])
gmsh.model.addPhysicalGroup(1, [b[1] for b in boundary], 1, "boundary")

gmsh.option.setNumber("Mesh.CharacteristicLengthMax", grid_size)
gmsh.model.mesh.generate(2)
gmsh.model.mesh.setOrder(order)
mesh, cell_marker, facet_marker = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, 0, gdim=3)
gmsh.finalize()

dx = ufl.Measure("dx", domain=mesh, subdomain_data=cell_marker)
x = ufl.SpatialCoordinate(mesh)
zero = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0))
F = ufl.as_vector([-x[1], x[0], zero])
cell_normal_vector = ufl.CellNormal(mesh)

form2 = dolfinx.fem.form(ufl.dot(ufl.curl(F), cell_normal_vector)* dx(10)) 
print(dolfinx.fem.assemble_scalar(form2))

# The line integral of F on the boundary ?
ds = ufl.Measure("ds", subdomain_data=facet_marker)
cell_edge_vector = ufl.classes.CellEdgeVectors(mesh)
form = dolfinx.fem.assemble_scalar(ufl.dot(cell_edge_vector, F) * ds )