# Numerically integrating a 2-form over a triangulated MeshRegion surface

I'd like to integrate a 2-form of the form $$adx\wedge dy+bdy\wedge dz+cdz\wedge dx$$ over a MeshRegion, which is topologically a compact surface with boundary.

Can Mathematica compute it automatically? The only way in my mind is to compute a normal vector of unit length for each triangle and replace the 2-form with a multiple of the volume form of the triangle on each point in the triangle, and then sum up.

Besides, how should I tell Mathematica the orientation of the surface? The surface is obtained from TriangulateMesh and I have a specified orientation of the boundary, given by ordering points on the boundary.

I know that there is a function

Integrate[f, {x, y, z} \[Element] mesh]


which computes the integration of $$f(x,y,z)$$ times the volume form over the surface, but the $$f$$ here is too hard to compute.

• I guess that's not built in, one will have to implement themselves. To do that in differential form, I found Keenan Crane's CS 15-458/858: Discrete Differential Geometry is a very good starting point, where algorithms and data structures are discussed in detail, side-by-side along with corresponding exterior calculus concepts. Commented Apr 26, 2023 at 16:08
• @Silvia Thanks. I've already written such a function mentioned above and it works well…… Commented Apr 28, 2023 at 14:48
• Nice! You do know our site encourages people answering their own questions right? :) Commented May 1, 2023 at 4:33
• @Silvia I've posted the codes together with a small example :) Commented May 8, 2023 at 14:40

Here is my solution. The 2-form is $$\mu=\mu_1dy \wedge dz + \mu_2dz \wedge dx+\mu_3 dx \wedge dy$$, where $$\{\mu_1, \mu_2, \mu_3\}$$ is

$$\left\{\frac{1}{2} \left(\frac{x}{2 \left(y^2-1\right) \left(x^2+y^2-1\right)}+\frac{\tan ^{-1}\left(\frac{x}{\sqrt{y^2-1}}\right)}{2 \left(y^2-1\right)^{3/2}}\right),\frac{1}{2} \left(\frac{y}{2 \left(x^2-1\right) \left(x^2+y^2-1\right)}+\frac{\tan ^{-1}\left(\frac{y}{\sqrt{x^2-1}}\right)}{2 \left(x^2-1\right)^{3/2}}\right),0\right\}.$$

For an oriented triangle with vertices $$a,b,c$$, construct a parametrization

paratriangle[a_, b_, c_, co1_, co2_] :=
a + (b - a)*co1 + (c - a)*co2


$$\{\partial x/\partial u,\partial y/\partial u,\partial z/\partial u\}$$ is

D[paratriangle[vert1, vert2, vert3, u, v], u]


Now we pull back the forms and get

dy wedge dz = (D[paratriangle[vert1, vert2, vert3,u,v],u][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[3]] -D[paratriangle[vert1, vert2, vert3,u,v],v][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],u][[3]]) du wedge dv

dz wedge dx= (D[paratriangle[vert1, vert2, vert3,u,v],u][[3]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[1]] - D[paratriangle[vert1, vert2, vert3,u,v],u][[1]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[3]]) du wedge dv

dx wedge dy=(D[paratriangle[vert1, vert2, vert3,u,v],u][[1]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[2]] -D[paratriangle[vert1, vert2, vert3,u,v],u][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[1]]) du wedge dv

Assume the pullback of $$\mu$$ is $$coeff \ du \wedge dv$$. Then

coeff=mu1* (D[paratriangle[vert1, vert2, vert3,u,v],u][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[3]]
-D[paratriangle[vert1, vert2, vert3,u,v],v][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],u][[3]])
+ mu2*(D[paratriangle[vert1, vert2, vert3,u,v],u][[3]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[1]]
- D[paratriangle[vert1, vert2, vert3,u,v],u][[1]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[3]])
+mu3 *(D[paratriangle[vert1, vert2, vert3,u,v],u][[1]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[2]]
-D[paratriangle[vert1, vert2, vert3,u,v],u][[2]]*D[paratriangle[vert1, vert2, vert3,u,v],v][[1]])


Let

standardtri = Triangle[{{0, 0}, {1, 0}, {0, 1}}]


We need to Integrate[coeff,{u,v}[Element]standardtri].

coeff[vertex1_, vertex2_, vertex3_] := Module[{vert1 = vertex1
, vert2 = vertex2
, vert3 = vertex3
, mu1val, mu2val, mu3val},
mu1val = mu1 /. {x -> paratriangle[vert1, vert2, vert3, u, v][[1]],
y -> paratriangle[vert1, vert2, vert3, u, v][[2]],
z -> paratriangle[vert1, vert2, vert3, u, v][[3]]};
mu2val = mu2 /. {x -> paratriangle[vert1, vert2, vert3, u, v][[1]],
y -> paratriangle[vert1, vert2, vert3, u, v][[2]],
z -> paratriangle[vert1, vert2, vert3, u, v][[3]]};
mu3val = mu3 /. {x -> paratriangle[vert1, vert2, vert3, u, v][[1]],
y -> paratriangle[vert1, vert2, vert3, u, v][[2]],
z -> paratriangle[vert1, vert2, vert3, u, v][[3]]};
mu1val*(D[paratriangle[vert1, vert2, vert3, u, v], u][[2]]*
D[paratriangle[vert1, vert2, vert3, u, v], v][[3]] -
D[paratriangle[vert1, vert2, vert3, u, v], v][[2]]*
D[paratriangle[vert1, vert2, vert3, u, v], u][[3]])
+ mu2val*(D[paratriangle[vert1, vert2, vert3, u, v], u][[3]]*
D[paratriangle[vert1, vert2, vert3, u, v], v][[1]] -
D[paratriangle[vert1, vert2, vert3, u, v], u][[1]]*
D[paratriangle[vert1, vert2, vert3, u, v], v][[3]])
+ mu3val*(D[paratriangle[vert1, vert2, vert3, u, v], u][[1]]*
D[paratriangle[vert1, vert2, vert3, u, v], v][[2]] -
D[paratriangle[vert1, vert2, vert3, u, v], u][[2]]*
D[paratriangle[vert1, vert2, vert3, u, v], v][[1]])
]


Now let's work out an example checking the Stokes' theorem.

Take the tetrahedra with vertices {0,0,0},{0.5,0,0},{0,0.5,0},{0,0,0.5}.

Integrate[
coeff[{0.5, 0, 0}, {0, 0.5, 0}, {0, 0, 0.5}], {u, v} \[Element]
standardtri] -
Integrate[
coeff[{0, 0, 0}, {0, 0.5, 0}, {0, 0, 0.5}], {u, v} \[Element]
standardtri] +
Integrate[
coeff[{0, 0, 0}, {0.5, 0, 0}, {0, 0, 0.5}], {u, v} \[Element]
standardtri] -
Integrate[
coeff[{0, 0, 0}, {0.5, 0, 0}, {0, 0.5, 0}], {u, v} \[Element]
standardtri]


we get

0.0232213


Define volfunc[x, y, z] to be $$\frac{1}{\left(x^2+y^2-1\right)^2}$$.

Simplify[D[mu1, x] + D[mu2, y] + D[mu3, z]] == volfunc[x, y, z]


returns True.

NIntegrate[
volfunc[x, y, z], {x, y, z} \[Element]
Tetrahedron[{{0, 0, 0}, {0.5, 0, 0}, {0, 0.5, 0}, {0, 0, 0.5}}]]


we get

0.0232213


as expected.