Вычислительный тензор инерции в 2D

Для любой фигуры вам нужно разбить ее на треугольники и обрабатывать каждый треугольник отдельно. Затем в итоге объединили результаты, используя следующие правила

В общем и целом

% Combined total area of all triangles
total_area = SUM( area(i), i=1:n )
total_mass = SUM( mass(i), i=1:n )
% Combined centroid (center of mass) coordinates
combined_centroid_x = SUM( mass(i)*centroid_x(i), i=1:n)/total_mass
combined_centroid_y = SUM( mass(i)*centroid_y(i), i=1:n)/total_mass
% Each distance to triangle (squared)
centroid_distance_sq(i) = centroid_x(i)*centroid_x(i)+centroid_y(i)*centroid_y(i)
% Combined mass moment of inertia
combined_mmoi = SUM(mmoi(i)+mass(i)*centroid_distance_sq(i), i=1:n)

Теперь для каждого треугольника.

Рассмотрим три угловые вершины с векторными координатами, точки , В а также С

a=[ax,ay]
b=[bx,by]
c=[cx,cy]

и следующие комбинации точек и перекрестных произведений (скалярных)

a·a = ax*ax+ay*ay
b·b = bx*bx+by*by
c·c = cx*cx+cy*cy
a·b = ax*bx+ay*by
b·c = bx*cx+by*cy
c·a = cx*ax+cy*ay
a×b = ax*by-ay*bx
b×c = bx*cy-by*cx
c×a = cx*ay-cy*ax

Свойства треугольника (с t(i) толщина и rho массовая плотность)

area(i) = 1/2*ABS( a×b + b×c + c×a )
mass(i) = rho*t(i)*area(i)
centroid_x(i) = 1/3*(ax + bx + cx)
centroid_y(i) = 1/3*(ay + by + cy)
mmoi(i) = 1/6*mass(i)*( a·a + b·b + c·c + a·b + b·c + c·a )

По компоненту выше

area(i) = 1/2*ABS( ax*(by-cy)+ay*(cx-bx)+bx*cy-by*cx)
mmoi(i) = mass(i)/6*(ax^2+ax*(bx+cx)+bx^2+bx*cx+cx^2+ay^2+ay*(by+cy)+by^2+by*cy+cy^2)

аппендикс

Немного теории здесь. Площадь каждого треугольника определяется с помощью

Area = 1/2 * || (b-a) × (c-b) ||

где × является векторным перекрестным произведением, и || .. || является векторной нормой (функцией длины).

Треугольник параметризован двумя переменными t а также s такой, что двойной интеграл A = INT(INT(1,dx),dy) дает общую площадь

% position r(s,t) = [x,y]
[x,y] = [ax,ay] + t*[bx-ax, by-zy] + t*s*[cx-bx,cy-by]

% gradient directions along s and t
(dr/dt) = [bx-ax,by-ay] + s*[cx-bx,cy-by]
(dr/ds) = t*[cx-bx,cy-by]

% Integration area element
dA = || (dr/ds)×(dr/dt) || = (2*A*t)*ds*dt
%
%   where A = 1/2*||(b-a)×(c-b)||

% Check that the integral returns the area
Area = INT( INT( 2*A*t,s=0..1), t=0..1) = 2*A*(1/2) = A

% Mass moment of inertia components

/  /  /  | y^2+z^2  -x*y    -x*z   |
I = 2*m*|  |  | t*|  -x*y   x^2+z^2  -y*z   | dz ds dt
/  /  /   |  -x*z    -y*z   x^2+y^2 |

% where [x,y] are defined from the parametrization