Как рассчитать массу неоднородной сферы?

Question

Я хочу рассчитать массу сферы, основанную на трехмерном дискретном неоднородном распределении плотности. Допустим, что набор 3х3х3 кубов разной плотности вписан сферой. Каков самый быстрый способ суммировать секционированные массы с помощью Python?

Я пытался вычислить объем под математическим уравнением для сферы: х^2 + у^2 + г^2 = R^2 для диапазона одного из кубов с использованием scipy.integrate.dblquad. Однако результат действителен только в том случае, если границы меньше радиуса сферы и повторяющиеся вычисления позволяют предположить, что 50 000 сфер с 27 кубами будут довольно медленными.

С другой стороны, обычное уравнение для CoM-калорий не может использоваться, на мой взгляд, из-за довольно грубого и дискретного распределения масс.

Answer 1

Я не могу получить точный смысл вписанной сферой. Также я не пробовал scipy.integrate. Однако, вот некоторые из них:

Установите куб 3x3x3 на единицу плотности. Затем возьмите интеграцию для каждого куба соответственно, поэтому вы должны иметь объемный куб V_ijk здесь. Теперь для каждой сферы вы можете получить массу каждой сферы, суммируя V_ijk*D_ijk, где D_ijk - плотность сферы.

Это должно быть намного быстрее, потому что теперь вам не нужно выполнять интеграцию.

Answer 2

Вы можете получить аналитическую формулу для пересекающегося объема между кубом (или прямоугольной призмой) и сферой. Это будет нелегко, но это должно быть возможно. Я сделал это для произвольного треугольника и круга в 2D. Основная идея состоит в том, чтобы разложить пересечение на более простые части, такие как тетраэдры и объемные сферические сектора треугольника, для которых известны относительно простые формулы объема. Главная трудность заключается в рассмотрении всех возможных случаев пересечений. К счастью, оба объекта выпуклые, поэтому вам гарантирован один выпуклый объем пересечения.

Приблизительный метод может состоять в том, чтобы просто подразделить кубы до тех пор, пока не будет выполнен приблизительный алгоритм численного интегрирования; это должно быть относительно быстрым. Вы знаете о Pick's Theorem? Это работает только в 2D, но есть, я считаю, 3D generalizations.

Answer 3

Timing Эксперимент

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

Без оптимизации каждый интеграл в сферических координатах может быть оценен в среднем за 0,005 секунды на стандартном ноутбуке, если плотность кубов является прямолинейной функцией.

Подобно тому, как ссылки, это программа в Mathematica:

Clear@f; 
(* Define a cuboid as density function *) 
iP = IntegerPart; 
f[{x_, y_, z_}, {lx_, ly_, lz_}] := iP[x - lx] + iP[y - ly] + iP[z - lz] /; 
    lx <= x <= lx + 3 && ly <= y <= ly + 3 && lz <= z <= lz + 3; 

f[{x_, y_, z_}, {lx_, ly_, lz_}] := Break[] /; True; 

Timing[Table[s = RandomReal[{0, 3}, 3]; (*sphere center random*) 
    sphereRadius = Min[Union[s, 3 - s]]; (*max radius inside cuboid *) 
    NIntegrate[(f[{x, y, z} - s, -s] /. (*integrate in spherical coords *) 
     {x -> r Cos@th Sin@phi, 
     y -> r Sin@th Sin@phi, 
     z -> r Cos@phi}) r^2 Sin@phi, 
     {r, 0, sphereRadius}, {th, 0, 2 Pi}, {phi, 0, Pi}], 
     {10000}]][[1]] 

В результате 52 секунд для 10^4 итераций.

Так что, возможно, вам не нужно оптимизировать много ...