$1/r^2$ gravitational force of triangles in 2D Suppose I have two triangles relatively close together (so they probably shouldn't really be treated as point masses).  I want to calculate the gravitational force (and potentially torque?) generated between the two bodies in the 2D plane.
For spheres/circles you can just treat them as point masses and go from there, but can you do that for arbitrary triangles (or tetrahedrons in 3D)?
I know the answer is probably to do a spatial integral across both triangles, but it's been a long time since I knew how to do that :)
The end goal is to be able to compute the gravitational force between arbitrary polygons/polyhedra.  I figured decomposing it in to triangles/tetrahedrons would be a good start.
...
UPDATE:
Okay, my multidimensional calculus is a bit rusty, but I think this is a promising direction:
Let:
$$\vec{f} = (a - c) \mu_1 + (b - c) \upsilon_1 - (x - z) \mu_2 - (y-z) \upsilon_2 - (z - c)$$
be the separation vector between two points on either triangle,
where $a, b, c$ are the vertices of triangle 1, and $\mu_1, \upsilon_1$ are the barycentric coordinates (corresponding to $a$ and $b$) for the point on triangle 1.  Likewise for $x, y, z$ and $\mu_2, \upsilon_2$ for triangle 2.  Using the barycentric coordinates let's us form the spatial integral to arrive at an answer.
So the (linear, non-torque) force between them is proportional to:
$$\vec{F_G} = \displaystyle\int_0^1 \int_0^{1-v_2} \int_0^{1} \int_0^{1-v_1} \! \frac{f}{||{f}||^3} \, \mathrm{d} \mu_1 \mathrm{d} \upsilon_1 \mathrm{d} \mu_2 \mathrm{d} \upsilon_2 $$
I think this admits a closed form solution, though I'm still wrestling with Mathematica.  I'd really be surprised if this integral hasn't been done somewhere before, though.
 A: First problem. How to parametrize the area enclosed by a triangle. 
Suppose the vertices are given by ${\bf r}_A$ , ${\bf r}_B$, and ${\bf r}_C$. Any position inside the triangle is described by two parameters as $$ {\bf r}(t,s) = (1-t) {\bf r}_A + t (1-s) {\bf r}_B + t s {\bf r}_C$$
A differential area at ${\bf r}(t,s)$ is
$$ {\rm d}A = \| \frac{\partial {\bf r}(t,s)}{\partial t} \times \frac{\partial {\bf r}(t,s)}{\partial s} \| = t \left( {\bf r}_A \times {\bf r}_B + {\bf r}_B \times {\bf r}_C + {\bf r}_C \times {\bf r}_A \right) {\rm d}s {\rm d}t $$
You can confirm that the area of the triangle is
$$ A = \int \limits_0^1  \int \limits_0^1 t\,\left( {\bf r}_A \times {\bf r}_B + {\bf r}_B \times {\bf r}_C + {\bf r}_C \times {\bf r}_A \right) {\rm d}s {\rm d}t = \frac{{\bf r}_A \times {\bf r}_B + {\bf r}_B \times {\bf r}_C + {\bf r}_C \times {\bf r}_A}{2}$$
The area differential is simplified by ${\rm d}A = 2 A t \,{\rm d}t{\rm d}s$. 
The mass of the triangle is easily found by integrating ${\rm d}m = \int \lambda\, {\rm d}A$ where $\lambda$ is the unit mass per area. This can be constant or not a function of position $\lambda(t,s)$. In general $$ m = \int \lambda 2 A t \,{\rm d}t {\rm d}s $$
Second Problem. Now we have two triangles ${\bf r}_1(t,s)$ and ${\bf r}_2(u,v)$ and we need to describe the differential gravity between them
$$ \begin{align} 
  {\rm d}{\bf F}_1 & = -G \frac{{\bf r}_1-{\bf r}_2}{\|{\bf r}_1-{\bf r}_2\|^3} {\rm d}m_1 {\rm d}m_2 \\
  {\rm d}{\bf F}_2 & = +G \frac{{\bf r}_1-{\bf r}_2}{\|{\bf r}_1-{\bf r}_2\|^3} {\rm d}m_1 {\rm d}m_2 
\end{align}$$
or in expanded form
$$ \begin{align} 
  {\bf F}_1 & = -G \int \int \int \int \frac{{\bf r}_1(t,s)-{\bf r}_2(u,v)}{\|{\bf r}_1(t,s)-{\bf r}_2(u,v)\|^3} (\lambda_1 2 A_1 t) (\lambda_2 2 A_2 u)\, {\rm d}t\,{\rm d}s\, {\rm d}u\, {\rm d}v \\
  {\bf F}_2 & = +G \int \int \int \int \frac{{\bf r}_1(t,s)-{\bf r}_2(u,v)}{\|{\bf r}_1(t,s)-{\bf r}_2(u,v)\|^3} (\lambda_1 2 A_1 t) (\lambda_2 2 A_2 u)\, {\rm d}t\,{\rm d}s\, {\rm d}u\, {\rm d}v 
\end{align}$$
Third problem, the integral requires the complete integrals of the first and second kind to be solved analytically.
Fourth problem, find any torque values. To get the moments about the origin (needs to be transformed to the center of mass for dynamics) is
$$ \begin{align} 
  {\rm d}{\bf M}_1 & = -G {\bf r}_1 \times \frac{{\bf r}_1-{\bf r}_2}{\|{\bf r}_1-{\bf r}_2\|^3} {\rm d}m_1 {\rm d}m_2 \\
  {\rm d}{\bf M}_2 & = +G {\bf r}_2 \times\frac{{\bf r}_1-{\bf r}_2}{\|{\bf r}_1-{\bf r}_2\|^3} {\rm d}m_1 {\rm d}m_2 
\end{align}$$
A: A straightforward way to test your answer is to use a Monte Carlo method to estimate it. The Matlab code below evenly picks random points within two triangles that are both symmetrical about the x-axis (simplifies code because the net gravitational force must act along the axis, but it could be easily adapted to allow other configurations). It calculates the gravitational force as 1/r^2 for each pair of point and the x-component of that force. It then simply gets the mean value of 1/r^2, and uses that to calculate the distance between the centres of action of the gravitational forces.
The code plots the points so you can see the two shapes. Change the value of x_offset to alter the distance between them, and change the conditions to alter the shapes of the triangles.
clear
N=100000;
%Pick N random points with x between 0 and 1 and y between -0.5 and 0.5
x1=rand(1,N);
y1=rand(1,N)-0.5;
%find which points are within a triangle and only retain those points in x1
%and y1
isintri1=find(abs(y1)<=0.5*x1);
x1=x1(isintri1);
y1=y1(isintri1);
plot(x1,y1)
hold on
x_offset = 2;
%Pick N random points with x between x_offset and x_offset + 1 and y between -0.5 and 0.5
x2=rand(1,N)+x_offset;
y2=rand(1,N)-0.5;
%find which points are within a triangle and only retain those points in x2
%and y2
isintri2=find(abs(y2)<=0.5-0.5*(x2-x_offset));
x2=x2(isintri2);
y2=y2(isintri2);
plot(x2,y2)
%Make sure that the two sets of points are the same length
min_len = min(length(x1),length(x2));
x2=x2(1:min_len);
x1=x1(1:min_len);
y2=y2(1:min_len);
y1=y1(1:min_len);
%Calculate distances between corresponding random points in the two triangles
distance_between_points = sqrt((x2-x1).^2+(y2-y1).^2);
%Force proportional to 1/r^2
force = 1./distance_between_points.^2;
%Calculate x-component of force
x_comp=force.*(x2-x1)/distance_between_points;
%Calculate mean x-component of force
mean_x_comp=mean(x_comp)
%Calculate equivalent ditance between point masses
equiv_dist = sqrt(1/mean_x_comp)
A: Google for triangle mass center (WP).
Concentrate the mass of each triangle in the locus of their 'center of mass' then apply $1/r^2$ or $1/r$ law.
Then substitute a pair of centers of mass by one center of mass equivalent, that lies in the line that conect them, you figure out the formula ;). Apply until only one centre of mass remains. This is the centre of mass of the system. 
After comments:  Yes, the objections are valid.
In particular I didn't realize that the triangles were close, Sorry.
The answer above serves only if they are far from each other.
If the objects are near we enter a tidal regime and we have to make an integration of 'all against all' .  
In the this general case I would divide the surfaces into small enough rectangles the same size/mass and will use the centers as atoms where the forces are applied.
I will choose to use the universal gravitational law (Newton).   
For each object: for each sub-object calculate the vector force due to all the sub-objects of the other object. For each object I would replace every two vectors for an equivalent resting on the right point on the line joining the surface elements. 
If the final resultant vector is is not centered in the center of mass then it will produce a torque that will make the object rotate.
The same reasoning would apply to other object.
In the NASA site we can find a fortran program to calculate the tidal force between the Moon and the Earth.  It uses discretization and 'all against all' .  
