In this essay, I implemented a Linear Finite Element Method (FEM) to solve the Laplace-Beltrami operator eigenproblem on a given triangular mesh surface. The method computes eigenvalues and eigenvectors corresponding to frequencies and modes of the Laplace-Beltrami operator, commonly used in shape analysis and geometry processing. The implementation is based on piecewise linear basis functions, meaning higher-order terms are not included. It converts an input surface mesh into a triangular mesh and processes elements. Local mass and stiffness matrices are constructed for each triangle, using a transformation to 2D for computation. These local matrices are then assembled into global stiffness and mass matrices using sparse matrix operations. The code supports different boundary conditions, including Dirichlet conditions that zero out boundary nodes and Neumann conditions which are left as default. Since this is only a linear FEM approach, it does not include higher-order elements such as quadratic or cubic, which limits its accuracy for larger eigenstates, esp for highly curved surfaces.
An intriguing problem is determining the eigensystem of the Laplacian operator on 2D manifolds embedded in 3D such as a gyroid surface. As of now, you cannot feed those surfaces into built-in functionalities of Mathematica, because the embedding dimension is not the same as the geometric dimension.
In[]:=
NDEigensystem[-Laplacian[u[x,y,z],{x,y,z}],u,{x,y,z}∈DiscretizeRegion[ImplicitRegion[Cos[x]+Cos[y]+Cos[z]==0&&0<=x<=2Pi&&0<=y<=2Pi&&0<=z<=2Pi,​​{x,y,z}]],1]
NDEigensystem
:The current version of NDSolve cannot solve equations over boundaries or surfaces. Please specify a region where the embedding dimension is the same as the dimension.
NDEigensystem-
(0,0,2)
u
[x,y,z]-
(0,2,0)
u
[x,y,z]-
(2,0,0)
u
[x,y,z],u,{x,y,z}∈
,1
There is a function in the Wolfram Function Repository that partially addresses this problem: LaplaceBeltramiOperator. However, as noted in the “Details” section, it is not efficient and does not support Dirichlet boundary conditions. In this essay, I explore a different approach which is more efficient and closer to the original FEM ideas.

Schwarz P surface

The Schwarz P surface is a triply periodic minimal surface defined by the implicit equation cos(x) + cos(y) + cos(z) = 0 within the domain [0, 2π] for x, y, and z. It is similar to the gyroid in that it is triply periodic and has zero mean curvature at every point, but it is simpler due to its more intuitive cubic symmetry and non-chiral nature. This surface has important applications in materials science, where it appears in self-assembled block copolymers and porous materials, as well as in physics, particularly in triply periodic electron density distributions and photonic crystals. It is also used in architecture and structural engineering to design lightweight yet strong structures, often incorporated into 3D-printed materials for mechanical optimization. Its simplicity and geometric properties make it a fundamental example of minimal surfaces with widespread practical significance.
Define an implicit region in 3D where the gyroid equation holds:
In[]:=
ir=ImplicitRegion[Cos[x]+Cos[y]+Cos[z]==0&&0<=x<=2Pi&&0<=y<=2Pi&&0<=z<=2Pi,​​{x,y,z}](*ImplicitRegion[Sin[x]Cos[y]+Sin[y]Cos[z]+Sin[z]Cos[x]==0&&0<=x<=2Pi&&0<=y<=2Pi&&0<=z<=2Pi,​​{x,y,z}]*);
Convert the implicit region into a discrete mesh region:
In[]:=
mr=DiscretizeRegion[ir,MaxCellMeasure->.005]
Out[]=
In[]:=
numEigenToCompute=10;
In[]:=
{eigenValues,eigenStates}=LaplaceBeltramiOperatorEigensystem[mr,10,"BoundaryCondition"->"Dirichlet"];
In[]:=
ListPlot[eigenValues,PlotRange->All,AspectRatio->1/2,GridLines->Automatic]
Out[]=
Plot eigenstates:
In[]:=
schwarzViz=With[{pnts=MeshCoordinates[mr],faces=MeshCells[mr,2]},Table[Graphics3D[{EdgeForm[],GraphicsComplex[pnts,faces,VertexColors->(ColorData["Rainbow"]/@Rescale[eigenStates[[e]]])]},Lighting->"Neutral",Boxed->False],{e,1,numEigenToCompute}]]
Out[]=

Gyroid surface

A gyroid surface is a type of triply periodic minimal surface, meaning it extends infinitely in three directions while maintaining a repeating structure. It has zero mean curvature, making it a minimal surface, and lacks both straight lines and mirror symmetries. The gyroid is found in natural systems like butterfly wings, biological membranes, and self-assembling polymers. Its equation, defines a complex, labyrinth-like structure often used in material science and photonics.
Define an implicit region in 3D where the gyroid equation holds:
In[]:=
ir=ImplicitRegion[Sin[x]Cos[y]+Sin[y]Cos[z]+Sin[z]Cos[x]==0&&0<=x<=2Pi&&0<=y<=2Pi&&0<=z<=2Pi,​​{x,y,z}];
Convert the implicit region into a discrete mesh region:
In[]:=
mr=DiscretizeRegion[ir,MaxCellMeasure->.01]
Out[]=
Number of eigenstates to compute:
In[]:=
numEigenToCompute=13;
Find eigenstates:
In[]:=
{eigenValues,eigenStates}=LaplaceBeltramiOperatorEigensystem[mr,numEigenToCompute,"BoundaryCondition"->"Dirichlet"];
Plot eigenvalues:
In[]:=
eigenValues[[;;3]]
Out[]=
{1.,1.,1.}
In[]:=
ListPlot[{eigenValues[[4;;]],eigenValues[[;;3]]},PlotMarkers->{Automatic,"X"},PlotRange->All,AspectRatio->1/2,GridLines->Automatic,PlotLabel->"Eigenvalues (three of them are trivial ones)",PlotLegends->{"Nontrivial","Trivial"}]
Out[]=
Nontrivial
Trivial
Since the first 3 eigenpairs are trivial (from the boundary nodes, artifact of calculation) and we will drop them.
Plot eigenstates:

Shortcoming of LaplaceBeltramiOperatorEigensystem

The primary limitation of LaplaceBeltramiOperatorEigensystem is that it uses a linear FEM, excluding higher-order terms. As a result, the accuracy is more dependent on mesh size, requiring a very fine mesh for improved results, which significantly increases computational cost.
User should note that we have considered the case of
Find eigenvalues with disk as the region (analytical solutions are available):
Find the eigenvalues from linear FEM:
Calculate %error:

Visualization

LaplaceBeltramiOperatorEigensystem

CITE THIS NOTEBOOK

A linear FEM to find eigensystem over a 2D manifold embedded in 3D​
by Mohammad Bahrami​
Wolfram Community, STAFF PICKS, March 6, 2025
​https://community.wolfram.com/groups/-/m/t/3410080