4
$\begingroup$

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.

enter image description here

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – Silvia
    Apr 26, 2023 at 16:08
  • $\begingroup$ @Silvia Thanks. I've already written such a function mentioned above and it works well…… $\endgroup$
    – Qing Lan
    Apr 28, 2023 at 14:48
  • $\begingroup$ Nice! You do know our site encourages people answering their own questions right? :) $\endgroup$
    – Silvia
    May 1, 2023 at 4:33
  • 1
    $\begingroup$ @Silvia I've posted the codes together with a small example :) $\endgroup$
    – Qing Lan
    May 8, 2023 at 14:40

1 Answer 1

4
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.