# General Solution for the Gravity Between Two 3D Triangles

I would like to find the general solution for the gravity between two (flat) triangles in 3D, including the location $(x,y,z)$ where this force should be applied (in order to later account for rotational forces.)

Given two triangles $t_1$ and $t_2$, which are defined by the points $t_{1v_1}$, $t_{1v_2}$, $t_{1v_3}$, $t_{2v_1}$, $t_{2v_2}$, and $t_{2v_3}$. Each of these points has $x$, $y$, and $z$ coordinates. For example, $t_{1v_2}y$ would be triangle 1 $\to$ vertex 2 $\to$ y coordinate.

Note: This will ultimately be used in the physics engine for a 3D space game (used with JBullet).

My attempt at a solution:

The formula for gravity for point masses is:

$$F_{\text{gravity}}=G\frac{\text{mass}_1 \times \text{mass}_2}{\text{radius}^2}$$

Which is easily expanded out to:

$$F_{\text{gravity}}=G\frac{\text{mass}_1 \times \text{mass}_2}{(x_1-x_2)^2+(y_1-y_2)^2+(z_1-z_2)^2}$$

With unit direction: *which is not considered in the following equations, but ultimately needs to be

$$\text{dir} = \text{normal}(\text{point}_1-\text{point}_2)$$

We must consider the entirety of each triangle:

$$\int_0^1\int_0^1(t_{1v_3}x \times q+t_{1v_2}x \times (1-q)) \times r+t_{1v_1}x \times (1-r)\ d_q\ d_r$$

By combining the equation for gravity with the integral of a single variable over a triangle, we get:

$$\int_0^1\int_0^1\int_0^1\int_0^1 \left[ \frac{(\text{mass}_1 \times \text{mass}_2 \times G)}{(((((t_{1v_3}x \times q+t_{1v_2}x \times (1-q)) \times r+t_{1v_1}x \times (1-r))-((t_{2v_3}x \times s+t_{2v_2}x \times (1-s)) \times t+t_{2v_1}x \times (1-t)) ))^2+((((t_{1v_3}y \times q+t_{1v_2}y \times (1-q)) \times r+t_{1v_1}y \times (1-r))-((t_{2v_3}y \times s+t_{2v_2}y \times (1-s)) \times t+t_{2v_1}y \times (1-t)) ))^2+((((t_{1v_3}z \times q+t_{1v_2}z \times (1-q)) \times r+t_{1v_1}z \times (1-r))-((t_{2v_3}z \times s+t_{2v_2}z \times (1-s)) \times t+t_{2v_1}z \times (1-t)) ))^2 )} \right] \ d_q\ d_r\ d_s\ d_t$$

Is this correct? What is the final solution? (including the position and/or angle of the force) Alternatively... Is this problem unsolvable?

• $1.$ Forces are vectors; they have direction, so you still need to take that into account. $2.$ The force due to one triangle acts at every point of the other triangle, which will lead to torques as well as translations. $3.$ It is almost certainly easier to use the potential formulation if you want an exact answer. $4.$ The easiest thing to do would be to just consider both of the triangles as point masses located at their centers of mass for the purposes of gravity (this is an approximation, however.) – Mr. G Jan 2 '15 at 4:50
• Referring to the answer by John Hughes, in the Form Factor theory (Radiative Heat Transfer) it is possible to replace the double-double integral by a line-line integral, by applying Green's theorem for each of the two surfaces (if I remember well). If such is possible, then it would reduce the computational complexity of the problem quadratically. – Han de Bruijn Jan 8 '15 at 12:57

This computation is essentially the same as computing the "Form Factor" in rendering (for computer graphics). It turns out that the integrals (once you include the direction of the force as you should) are not nice (i.e., involve dilogarithms). Here's a paper that describes the computation in some detail.

• By the way, following up on user7530's answer, there's been a good deal of effort expended in graphics in rapidly approximating this value. The simplest/cheapest is probably the "hemicube" method. Also: you might want to do those numerical computations and see how the results, when the distance between triangles is more than say, 5 to 10 times the size (e.g., longest edge) of the larger. I'll bet it looks pretty close to $Gm_1 m_2/r^2$ where the distance $r$ is taken to be the distance between the centers of mass. This may save you a great deal of work. – John Hughes Dec 29 '14 at 15:27

I'm also interested in whether this integral has a nice closed form (I'm a bit skeptical). It will almost certainly involve choosing a clever set of coordinates (such as slicing one triangle intro strips parallel to the other) rather than brute-force barycentric interpolation. When searching for a solution to this integral, you may want to also search for the net Coulomb interaction between two triangular plates, and/or the integral of the gradient of the Green's function of the Laplacian over two triangles (which might have been addressed somewhere in the Boundary Element Method (BEM) literature).

For your purposes I'm going to suggest an alternative approach, though: instead of trying for an exact formula, approximate the integral using numerical quadrature. This PDF lists, towards the end, good quadrature weights and quadrature point locations for triangles. Numerical quadrature will work well provided that the two triangles are well-separated: if they touch or overlap you are going to need more sophisticated methods.

The coordinate transformation employed in this answer has a couple of references in:

Suppose that the two triangles, with vertices numbered as $(0,1,2)$, are successively described by: $$\begin{array}{lll} (r) & : &\vec{r}(\xi,\eta) = \vec{r_0}+\xi\,(\vec{r_1}-\vec{r_0})+\eta\,(\vec{r_2}-\vec{r_0})\\ (\rho) & : &\vec{\rho}(\alpha,\beta) = \vec{\rho_0}+\alpha\,(\vec{\rho_1}-\vec{\rho_0})+\beta\,(\vec{\rho_2}-\vec{\rho_0}) \end{array}$$ with $\,0 \le \xi \le 1$ , $0 \le \eta \le 1-\xi\,$ and $\,0 \le \alpha \le 1$ , $0 \le \beta \le 1-\alpha$ .

Then the force upon an infinitesimal part $(d\alpha\,d\beta)$ of triangle $(\rho)$ being attracted by an infinitesimal part $(d\xi\,d\eta)$ of triangle $(r)$ is: $$d\vec{F} = G\cdot m_r\cdot m_\rho \frac{\vec{r}(\xi,\eta)-\vec{\rho}(\alpha,\beta)} {\left|\vec{r}(\xi,\eta)-\vec{\rho}(\alpha,\beta)\right|^{3/2}} 2\,d\xi\,d\eta\cdot 2\,d\alpha\,d\beta$$ Where the factors $2$ comes from the fact that the parent triangles have area $1/2$ instead of $1$ : $$2\int_{\xi=0}^{\xi=1} d\xi \int_{\eta=0}^{\eta=1-\xi} d\eta = 2\int_{\alpha=0}^{\alpha=1} d\alpha \int_{\beta=0}^{\beta=1-\alpha} d\beta = 1$$ The force as a whole with $\vec{r} = (x,y,z)$ , $\vec{\rho} = (a,b,c)$ is: $$\left[ \begin{array}{c} F_x\\F_y\\F_z \end{array} \right] = \left[\begin{array}{c} \int_{\xi=0}^{\xi=1}\int_{\eta=0}^{\eta=1-\xi}\int_{\alpha=0}^{\alpha=1}\int_{\beta=0}^{\beta=1-\alpha} \frac{4\,G\cdot m_r\cdot m_\rho\cdot (x-a)\cdot d\xi\,d\eta\cdot d\alpha\,d\beta} {\left|(x-a)^2+(y-b)^2+(z-c)^2\right|^{3/2}} \\ \int_{\xi=0}^{\xi=1}\int_{\eta=0}^{\eta=1-\xi}\int_{\alpha=0}^{\alpha=1}\int_{\beta=0}^{\beta=1-\alpha} \frac{4\,G\cdot m_r\cdot m_\rho\cdot (y-b)\cdot d\xi\,d\eta\cdot d\alpha\,d\beta} {\left|(x-a)^2+(y-b)^2+(z-c)^2\right|^{3/2}} \\ \int_{\xi=0}^{\xi=1}\int_{\eta=0}^{\eta=1-\xi}\int_{\alpha=0}^{\alpha=1}\int_{\beta=0}^{\beta=1-\alpha} \frac{4\,G\cdot m_r\cdot m_\rho\cdot (z-c)\cdot d\xi\,d\eta\cdot d\alpha\,d\beta} {\left|(x-a)^2+(y-b)^2+(z-c)^2\right|^{3/2}} \end{array} \right]$$ Where: $$\left\{\begin{array}{l} x = x_0+\xi\,(x_1-x_0)+\eta\,(x_2-x_0)\\ y = y_0+\xi\,(y_1-y_0)+\eta\,(y_2-y_0)\\ z = z_0+\xi\,(z_1-z_0)+\eta\,(z_2-z_0)\end{array}\right. \qquad \qquad \left\{\begin{array}{l} a = a_0+\alpha\,(a_1-a_0)+\beta\,(a_2-a_0)\\ b = b_0+\alpha\,(b_1-b_0)+\beta\,(b_2-b_0)\\ c = c_0+\alpha\,(c_1-c_0)+\beta\,(c_2-c_0)\end{array}\right.$$ A short excercise with MAPLE learns that even calculating the innermost integrals is a disaster.
And it doesn't finish at all with calculating the innermost double integrals. Therefore a "closed form" seems to be out of reach anyway, unless someone can devise some clever trick.
It has been suggested that it would be easier to calculate the potential energy of a triangle $(r)$ with respect to a triangle $(\rho)$ : $$V = \int_{\xi=0}^{\xi=1}\int_{\eta=0}^{\eta=1-\xi}\int_{\alpha=0}^{\alpha=1}\int_{\beta=0}^{\beta=1-\alpha} \frac{4\,G\cdot m_r\cdot m_\rho\cdot d\xi\,d\eta\cdot d\alpha\,d\beta} {\sqrt{(x-a)^2+(y-b)^2+(z-c)^2}}$$ But, apart from calculating the integrals, I have no idea how to proceed from there to obtain the force.

### Numerically

As suggested in a comment by Mr. G , The easiest thing to do would be to just consider both of the triangles as point masses located at their centers of mass for the purposes of gravity (this is an approximation, however) : $$\vec{F} = G\cdot m_r\cdot m_\rho \frac{\vec{r}(\xi,\eta)-\vec{\rho}(\alpha,\beta)} {\left|\vec{r}(\xi,\eta)-\vec{\rho}(\alpha,\beta)\right|^{3/2}} \qquad \mbox{with} \qquad \xi = \eta = \alpha = \beta = \frac{1}{3}$$ With a subdivision of both triangles, like in the above picture on the right, the approximation can be made more accurate. Define $N^2$ midpoints in one triangle and $M^2$ in the other: $$\begin{array}{l} \vec{r}(i,j) = \vec{r_0}+\frac{i+(1/3\, \vee\, 2/3)}{N}\,(\vec{r_1}-\vec{r_0})+\frac{j+(1/3\,\vee\, 2/3)}{N}\,(\vec{r_2}-\vec{r_0})\\ \vec{\rho}(k,l) = \vec{\rho_0}+\frac{k+(1/3\,\vee\, 2/3)}{M}\,(\vec{\rho_1}-\vec{\rho_0})+\frac{l+(1/3\,\vee\, 2/3)}{M}\,(\vec{\rho_2}-\vec{\rho_0}) \end{array}$$ For suitable integers $\,0 \le i+j < N(-1) \; ; \; 0 \le k+l < M(-1)$ , namely such that all midpoints are inside the main triangles. But I don't think that implementing such a subdivision is intended by the OP, because the gravitational force between the two main triangles is meant as simple case in itself. But anyway, then the net force would become: $$\vec{F} = \frac{G\cdot m_r\cdot m_\rho}{N^2 M^2} \sum_i \sum_j \sum_k \sum_l \frac{\vec{r}(i,j)-\vec{\rho}(k,l)} {\left|\vec{r}(i,j)-\vec{\rho}(k,l)\right|^{3/2}}$$

• That's a beautiful answer, thank you. I will probably use a similar triangle subdivision algorithm in my final code (with recursive subdivision though). – Kent Jan 18 '15 at 20:18