Surface integration over a portion of an ellipsoid

I would like to perform a surface integration over a portion $D$ of an ellipsoid. A plane arbitrarily intersects the ellipsoid forming two sections, of which one is $D$. I do not know how I can discretize $D$ in order to perform the numerical integration.

Attempt at a Solution

I parameterize the ellipsoid in spherical coordinates, then for each grid point I check whether that point lies in $D$ or not using the equation of the intersecting plane. If it does, I compute its contribution to a Riemann integral over a rectangular mesh.

Problems with Solution

1. Variable mesh size due to parameterization using spherical coordinates (big squares near center, tiny ones near poles).
2. Slow 1st order convergence (not to mention the if loop in python)
3. Inaccuracies near the plane-ellipsoid boundary due to use of square mesh

Desired Solution

I would like to use triangles to approximate the surface. I am aware of a well-developed literature on triangle meshes, but I cannot find how to deal with the irregular domain I have.

Please note that this is a revised version of a question I posted on Math Stack Exchange a week ago:

https://math.stackexchange.com/questions/851334/numerical-integration-over-two-regions-of-an-ellipsoid

• Have you thought about accurately solving for the boundary curve (in the spherical coordinate system)? I assume once the solution of this 3D intersection problem is found accurately, then the resultant surface integral can also be computed accurately since it now becomes the calculation of a planar integral bounded by a curve. In particular, you should be able to rotate your coordinate system so that it becomes two 1D integrals (a top branch and a lower branch?) – TSGM Jul 10 '14 at 10:48
• Why not using ellipsoidal coordinates? Maybe that parametrization is better for integration. – nicoguaro Aug 6 '14 at 18:43

Step 1: Map ellipsoid to unit sphere. For a point $r$ on the ellipsoid, you can compute the map $r'=Mr$ where $\left|r'\right|=1$. This should be a linear map (assuming the ellipsoid was centered at the origin). $M$ is a $3\times 3$ matrix that is a composition of scaling and rotation. Note that the cut plane also transforms into another plane. We can choose the mapping to make the plane parallel to the $xy$ plane.

Suppose your ellipsoid is defined by $r^T A r=1$, where $r\in\mathbb{R}^3$, and your cut plane is $n^T r = 1$, then the mapping $r' = A^{1/2}r$ is sufficient to get you to a sphere, and $r'=QA^{1/2}r$ will also map the cut plane appropriately. The pure rotation matrix $Q$ can be computed as follows. Compute $z'=\mathrm{normalize}\left[(A^T)^{-1/2}n\right]$. Apply a basis construction algorithm (Gram-Schmidt or QR factorization) to compute an orthonormal triad $x',y',z'$. The matrix $Q$ has $z'$ as the first row, and the other two rows are $x'$ and $y'$ in the order to make the determinant $+1$ instead of $-1$.

Note: if you've never seen this notation for an ellipsoid, then notice that for an axis aligned ellipsoide with semi-major axes $a,b,c$, the matrix $a$ is simply diagonal with diagonal elements $a^{-2},b^{-2},c^{-2}$. It should generally be symmetric positive definite. Also note that $n$ is not necessarily a unit vector.

Step 2: Okay, at this point you have your mapping $M = QA^{1/2}$, and $M^{-T} n$ should be parallel to the $+z$ axis, with a length that is less than 1 if your plane actually intersects the ellipsoid. Some basic analytical geometry well tell you the equation of the circle of intersection between the transformed plane and sphere, and you can apply whatever numerical integration technique you want to this spherical cap, which is a much nicer domain. The only problem that remains is the Jacobian of the transformation.

The Jacobian of the back transformation is $\det M^{-1}$, but we need it restricted to the tangent space. Fortunately, we know the tangent space; it is the orthogonal complement of the normal vector at any point on the ellipsoid, so for any point $r'$ on the unit sphere, the original point was $r = M^{-1}r'$ and the normal vector there is in the direction $Ar$. Using that vector and doing a full QR decomposition (or Gram-Schmidt) to make an orthonormal triad gives you a pair of orthonormal basis vectors that span the tangent space. Let $U\in\mathbb{R}^{3\times 2}$ be the matrix with those two vectors as columns. The Jacobian is then the projection of $M$ into this space: $J=\det(U^TM^{-1}U)$.

Step 3: Do the integration. Let your function on the ellipsoid be $f(r)$. Apply a quadrature routine to the spherical cap you obtained in step 2 to evaluate: $$\int_\textrm{your stupid domain} f(r) dr = \int_\textrm{spherical cap} f(M^{-1}r') J \,dr'$$

As for how to integrate over the spherical cap, you can do another projection into the plane to integrate in a circular disc, where you'll have to work out the Jacobian of the projection (it's pretty nasty if you're near the equator). Alternatively, there are good semi-analytic triangulations you can apply now.

• This is great. Prior to having read this, I used Delaunay triangulation to mesh the surface, and for each triangle checked whether its center was on side A or side B of the plane. This makes for a jagged boundary of the domain, but I find it converges pretty quickly. Your solution looks much better; if I experience any convergence problems as I proceed I'll certainly implement this. – Eric Kightley Aug 9 '14 at 20:56

My answer below assumes your integrand is nonsingular over the domain of integration. If you have singularities, you will need more complicated machinery.

A high-level "straightforward" solution I can think of for "simple" domains is to come up with a coordinate transformation (really, a diffeomorphism) to map your domain $D \subset \mathbb{R}^{3}$ (which, from what you've described, is a manifold) to another domain $D' \subset \mathbb{R}^{2}$). Assuming you can translate your description of $D'$ to something like a collection of cubic splines, you should be able to mesh it with triangles or quads using open source meshing software (Gmsh and Triangle come to mind). Good meshing software for finite elements will discretize $D'$ such that the mesh is relatively uniform, and no particular element is skewed, so as long as your integrand is nonsingular over your domain, this approach will probably work. Then you can integrate the surface of $D'$ using quadrature, and use your coordinate transformation to map the integration over $D'$ back to integration over $D$. Obviously, the implementation is tricky, and requires some pencil and paper analytical calculations, but the quadrature and meshing parts should be within the ability of anyone with a little experience in finite element methods

Another approach that would work would be to mesh the surface directly, again, using finite element method meshing software. Something I searched for that might work for you would be something like CGAL, although it looks a little involved. You'll need a fine enough mesh that it approximates the surface well. The integration should again be doable for anyone with a finite element background.

In both cases above, improving the order of convergence with respect to the fineness of the discretization (really, the diameter of the elements) means using higher-order mesh elements to approximate your domain, and higher-order quadrature rules.

Let's assume for a moment that your plane intersects the ellipsoid in such a way that less than one half of the surface of the ellipsoid forms your domain $D$. Then it should not be very difficult to find a transformation that maps a half-sphere (that is, one half of the surface of a ball) to your domain $D$. Now, we know how to triangulate a half-sphere: with triangles, use one half of the 8 triangles that form an octahedron; with quadrilaterals, take the six sides of a cube and refine them into 24 smaller quadrilaterals by pulling the mid-points out to the surface, then take one half of them for the half-sphere. In either of these cases, you will end up with triangles or quadrilaterals that are all of the same size. To mesh $D$, you then just have to map these triangles or quadrilaterals from the half-sphere to $D$.

In case your plane intersects the ellipsoid in a way that $D$ forms more than one half of the surface, one can employ a similar technique.

• This is clearly a valid approach, the trick is in the implementation. – Eric Kightley Jul 9 '14 at 17:55