Для начала нужно обзавестись этим самым "mesh"-ем поверхности, или триангуляцией поверхности, полигональной сеткой, разбиением двумерного многообразия. В данном случае работа будет вестись именно с треугольной сеткой, но все ниже представленные формулы и код (если немного модифицировать), будет работать с сеткой состоящей из любых полигонов. Главное, чтобы они были малые, от этого зависит точность, чем меньше - тем лучше.

Код здесь написан языке Python, mesh импортируется из .stl файла. Откуда он будет импортироваться совершенно неважно, Python дальше будет работать с данными, которые имеют тип массивов numpy.

# Треугольники (список треугольников, в каждый из которых входят 3 точки)
[[[-5.0500002e+00  0.0000000e+00  0.0000000e+00]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-4.9604182e+00  2.5664667e-02 -4.8957970e-03]]

 [[-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-5.0500002e+00  0.0000000e+00  0.0000000e+00]]

 [[-4.9604182e+00  2.6127456e-02  0.0000000e+00]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-5.0500002e+00  0.0000000e+00  0.0000000e+00]]

 ...

 [[-4.8888855e+00  4.6034887e-02 -2.8962698e-03]
  [-4.9604182e+00  2.5921434e-02 -3.2746387e-03]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]]

 [[-4.8888855e+00  4.6034887e-02 -2.8962698e-03]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.8888855e+00  4.6125907e-02  0.0000000e+00]]

 [[-4.8888855e+00  4.6125907e-02  0.0000000e+00]
  [-4.9604182e+00  2.6075900e-02 -1.6405566e-03]
  [-4.9604182e+00  2.6127456e-02  0.0000000e+00]]]
# Нормали к треугольникам (вектора)
[[ 0.27986652 -0.94821906  0.15018432]
 [ 0.27986658 -0.95577824  0.09034719]
 [ 0.27986655 -0.95956516  0.03015531]
 ...
 [ 0.26912326 -0.95883155  0.09063575]
 [ 0.26912323 -0.96263045  0.03025224]
 [ 0.26912326 -0.96263045  0.03025234]]

# pip install numpy-stl  # чтобы установить библиотеку
my_mesh = mesh.Mesh.from_file('название_файла.STL')
Normals = my_mesh.normals 
Triangles = my_mesh.vectors

# и сразу нормализуем нормали, на всякий пожарный
Nnorm = np.linalg.norm(Normals, axis=1) 
Nnorm = np.tile(Nnorm.reshape(len(Nnorm), 1), (1, 3))
Normals = Normals/Nnorm

Теперь у нас есть треугольники и единичные нормали к ним. Но для дальнейшей работы еще понадобится список площадей треугольников и их центров масс.

A, B, C = Triangles[:, 0], Triangles[:, 1], Triangles[:, 2] 
R1, R2 = C-A, B-A

# список центров масс
R = 1/3*(A+B+C)

# список площадей 
S = np.cross(R1, R2)
S = np.linalg.norm(S, axis=1)/2

Теперь у нас есть:

• S - список площадей треугольников

• R - список центров масс треугольников

• Normals - список нормалей к треугольникам

• Triangles - список треугольников (треугольник содержит 3 точки)

Теперь следует напомнить как вычисляются (чисто математически) массы, объемы, площади, центры масс, моменты инерции и что вообще это такое и зачем вообще нужно.

А те кто в теме и сомневаются, что все эти вещи можно вычислять для объемных тел используя только 2d сетку, скажу: это делать можно. И нужно, если лень разбивать 3d область пространства на маленькие элементы - тетраэдры там, призмы... Да, это делать возможно, и мы это сделаем. Но всё по порядку.

Масса - вычисляется как интеграл по объему (площади, линии):

здесь под интегралом (тройным) - плотность вещества. Объемная плотность вещества, имеющая размерность кг/м^3. Но, если в наличии имеется поверхностная плотность (кг/м^2), или линейная (?? кг/м), то интегралы соответственно будут двойные и простые.

Объяснять зачем нужна масса мало кому нужно, вернее об это знают не только лишь все. Масса на ускорение - это сила. Второй закон Ньютона, известный всем еще с седьмого класса.

Для большинства приложений в инженерии достаточно работать с плотностями, которые постоянны в пространстве и во времени =) Потому мы упростим себе жизнь и вынесем за все интегралы с которыми будем работать все плотности. И даже работать мы будем не с массами, а с площадями и объемами. А когда будем вычислить центры масс, плотности вообще сами по себе сократятся, и неважны.

Кстати, центр масс:

Ну или момент первого порядка, или как там его называют в статистике.

В инженерном деле (конструкторском) центр масс важен, например, для описания движения твердых (недеформируемых) тел. Конструктора помещают туда связанную с телом координатную систему (как-то располагая её относительно тела, под каким-то углом в смысле) и с этих пор тело и связанная кс неразлучны друг с другом и неподвижны друг относительно друга.

Зачем помещают связанную систему в центр масс? Ну так принято, пойдем дальше.

Вот мы пришли к моментам инерции. Момент инерции - это то, что влияет на выбор конструктора, когда перед ним встает задача: как именно закрепить на теле связанную координатную систему.

Удачный выбор закрепления будет влиять вполне конкретно на уравнения движения этого самого тела. Момент инерции - это аналог массы во втором законе Ньютона для точки. Еще говорят "тензор инерции" - он участвует в описании вращательного движения тела, но это далеко не скалярная величина, это не просто число. Тензор инерции имеет валентность матрицы, то есть таблицы порядка 2 на 2 или 3 на 3, в зависимости от размерности пространства.

Ну а закрепляют кс исходя из симметрии тела (корабли, самолеты, дирижабли.. ) имеют зачастую симметрию относительно вертикальной плоскости. Вот поэтому одна из координатных плоскостей кс совпадает с ней. А дальше у нас остается только одна степень свободы как расположить - угол. И обычно одну из осей кс, находящихся в вертикальной плоскости, направляют горизонтально, параллельно земле. Тогда тензор инерции будет иметь почти идеальный диагональный вид, и, кстати, его можно было бы сделать диагональным в любом случае, для любого тела, но всё же кс располагают так как я описал. Ну чтобы диагональные элементы имели более ясный, практичный физический смысл. Что-то далеко я отошел от темы, вычисляется момент инерции так:

это момент инерции для точки, для тела бы записалось как сумма по всем его точкам, или интеграл. А точнее, куча интегралов, в точности 6 различных штук, ведь тензор инерции - симметричная матрица (и еще её всегда можно диагонализировать - повернуть в пространстве координатную систему, относительно которой и вычисляется момент, так, чтобы матрица приняла диагональный вид - все элементы 0, которые не на диагонали).

А еще тензор инерции можно пересчитывать в другие координатные системы смещенные относительно центра масс по теореме Гюйгенса. То есть, есть формула, связывающая моменты инерции в разных координатных системах, и зная в одном можно найти в другом имея только лишь вектор смещения. Но важно! теорема Гюйгенса работает только когда одна из кс расположена в центре масс. Если есть две кс вне центра масс, в формуле будут появляться дополнительные члены (моменты первого порядка).

Итак, достаточно. Всё что нужно напомнилось (было напомнено), короче, вспомнено. Восполнено, ладно, клоунада как клоунада, она иногда надо.

А давайте посчитаем все эти вещи просто для поверхности, чтобы набить немного руку?

Масса поверхности получается как сумма площадей всех треугольников умноженная на поверхностную плотность:

def get_m_obolochka(S, rho):
    m = np.sum(S) * rho
    return m

Центр масс поверхности (оболочки):

def get_cm_obolochka(S, R):
    S_sum = np.sum(S)
    cm = np.einsum('ij,i', R, S) / S_sum
    return cm

Момент инерции оболочки:

где матрица вычисляется по приведенной где-то выше формуле, а еще здесь она не умноженная на i-тую массу.

def get_J_for_ob(rho, R, S):
	x, y, z = R[:, 0], R[:, 1], R[:, 2]
	x2, y2, z2 = x**2, y**2, z**2

	Jx, Jy, Jz = y2+z2, x2+z2, x2+y2
	Jxy, Jyz, Jxz = -x*y, -y*z, -x*z

	J = np.array([	[Jx, Jxy, Jxz],
									[Jxy, Jy, Jyz], 
									[Jxz, Jyz, Jz]
							])
	J = J.T # это нужно, потому что если подумать, то нам нужен список J[i] для всех треугольников
	# короче делаем выше правильный shape =)
  J = np.einsum('ijk,i->jk', J, S)*rho
	return J

Теперь, когда размялись, можно приступать к основному, тому, что написано в названии статьи. Вычислим все эти характеристики для объемного тела, которое ограничено mesh-ем, используя только mesh.

Для этого понадобится, внимание... теорема Остроградского-Гаусса:

На мой взгляд это самая красивая теорема математики... на самом деле она вроде как есть частный случай операций над дифференциальными формами, но я не слишком математик, поверхностно знаком ;) Да... там и теорема Гаусса-Бонне и все в этом духе, это топология, детка. Это всё очень красиво, и может быть я этим когда-нибудь займусь, когда меня освободят из рабства.

Формула позволяет переходить от вычислений по объему, к вычислениям по поверхности, ограничивающей заданный объем. Это можно интерпретировать как проекцию. Мы первое что видим, когда рождаемся - это проекцию трехмерной мамы на нашу двухмерную сетчатку глаза....

Да и куда ни глянь, всё есть проекцией, всё является лишь тенью более сложной "настоящей" вещи, вещи из "реального" мира.. а она в свою очередь - еще одна проекция из следующей итерации, следующего под-уровня реальности, под-Матрицы Вачовской.

Пренебрегая философией, и подбирая векторное поле следующим образом:

мы видим, что дивергенция становится равна:

Это означает, что теперь зная лишь mesh поверхности, мы можем вычислить объем:

Можно сказать, что здесь:

(Центры масс треугольников, напоминаю)

Код:

def Volume(S, R, Normals):
	V = np.einsum('ij,ij,i', R, Normals, S)*1/3 # einsum как работает можно найти в документации
	return V

Но можно обнаглеть, и пойти дальше объема:

Используя ту же теорему и задав вектор F так (для вычисления координаты x):

(чего-то индекс x не туда залез, я его хотел чутка ниже поместить)

Его дивергенция:

Координата центра масс тела:

Аналогичным образом вычисляются остальные две координаты. Компактно это можно записать если ввести операцию покоординатного умножения (использующегося в numpy):

где

Воплощение в коде:

def get_cm_V(S, R, Normals, V_t):
	F = R**2/2
	cm = np.einsum('ij,ij,i->j', F, Normals, S)
	return cm/V_t

Для объемного тела подсчитать тензор инерции можно используя только триангуляцию оболочки, аналогично как делалось для подсчета объема и центра масс:

Где вектор F можно, например, задать так:

Для компоненты должно выполняться:

Соответственно можно задать F:

Аналогично и с остальными компонентами.

def get_J_for_V(S, R, Normals, rho):
	x, y, z = R[:, 0], R[:, 1], R[:, 2]
	x2, y2, z2 = x**2, y**2, z**2

	Fx, Fy, Fz = x*(y2+z2), y*(x2+z2), z*(x2+y2)
	Fxyz = -x*y*z

	Fx, Fy, Fz = Fx[:, None], Fy[:, None], Fz[:, None]
	Fxyz = Fxyz[:, None]

	N = len(x)
	zero = np.zeros((N, 1))

	Fxx, Fyy, Fzz = np.hstack((Fx, zero, zero)), np.hstack((zero, Fy, zero)), np.hstack((zero, zero, Fz))
	Fxy, Fxz, Fyz = np.hstack((zero, zero, Fxyz)), np.hstack((zero, Fxyz, zero)), np.hstack((Fxyz, zero, zero))

	Jx, Jy, Jz = np.einsum('ij,ij,i', Fxx, Normals, S), np.einsum('ij,ij,i', Fyy, Normals, S), np.einsum('ij,ij,i', Fzz, Normals, S)
	Jxy, Jyz, Jxz = np.einsum('ij,ij,i', Fxy, Normals, S), np.einsum('ij,ij,i', Fyz, Normals, S), np.einsum('ij,ij,i', Fxz, Normals, S)

	J = np.array([	[Jx, Jxy, Jxz],
					[Jxy, Jy, Jyz], 
					[Jxz, Jyz, Jz]
				])*rho

	return J
# не так элегантно, но работает вроде.. 

На этом всё, спасибо за внимание.