Объем 3d-формы с использованием численного интегрирования со Scipy

Я написал функцию для вычисления объема пересечения куба и полупространства, и сейчас я пишу тесты для него.

Я попытался вычислить объем численно следующим образом:

integral = scipy.integrate.tplquad(lambda z, y, x: int(Vector(x, y, z).dot(normal) < distance),
                                   -0.5, 0.5,
                                   lambda x: -0.5, lambda x: 0.5,
                                   lambda x, y: -0.5, lambda x, y: 0.5,
                                   epsabs=1e-5,
                                   epsrel=1e-5)

... в основном я интегрирую по всему кубу, и каждая точка получает значение 1 или 0 в зависимости от того, находится ли она внутри полупространства. Это становится очень медленным (более нескольких секунд на вызов) и продолжает давать мне такие предупреждения, как

scipy.integrate.quadpack.IntegrationWarning: The integral is probably divergent, or slowly convergent

Есть ли лучший способ рассчитать этот объем?

Источник

3 ответа

Решение

интеграция

Интеграция прерывистой функции проблематична, особенно в многомерном. Необходима некоторая предварительная работа, сводящая задачу к интегралу от непрерывной функции. Здесь я определяю высоту (сверху вниз) как функцию от x и y и использую dblquad для этого: он возвращается в 36.2 ms,

Я выражаю уравнения плоскости как a*x + b*y + c*z = distance, Нужна некоторая осторожность со знаком c, так как плоскость может быть частью вершины или низа.

from scipy.integrate import dblquad
distance = 0.1
a, b, c = 3, -4, 2  # normal
zmin, zmax = -0.5, 0.5  # cube bounds

# preprocessing: make sure that c > 0
# by rearranging coordinates, and flipping the signs of all if needed

height = lambda y, x: min(zmax, max(zmin, (distance-a*x-b*y)/c)) - zmin
integral = dblquad(height, -0.5, 0.5, 
                   lambda x: -0.5, lambda x: 0.5,
                   epsabs=1e-5, epsrel=1e-5)

Методы Монте-Карло

Выбор случайных точек выборки (метод Монте-Карло) позволяет избежать проблем с разрывом: точность для прерывистости примерно одинакова, как для непрерывных функций, ошибка уменьшается со скоростью 1/sqrt(N) где N - количество точек выборки.

Многогранный пакет использует его внутри. С этим вычисление может пойти как

import numpy as np
import polytope as pc
a, b, c = 3, 4, -5  # normal vector
distance = 0.1 
A = np.concatenate((np.eye(3), -np.eye(3), [[a, b, c]]), axis=0)
b = np.array(6*[0.5] + [distance])
p = pc.Polytope(A, b)
print(p.volume)

Здесь A и B кодируют полупространства как Ax<=b: первые строки исправлений для граней куба, последние для плоскости.

Чтобы лучше контролировать точность, либо внедрите метод Монте-Карло самостоятельно (легко), либо используйте mcint пакет (примерно так же просто).

Объем многогранника: задача для линейной алгебры, а не для интеграторов

Вы хотите вычислить объем многогранника, выпуклого тела, образованного пересекающимися полупространствами. Это должно иметь алгебраическое решение. SciPy имеет класс HalfspaceIntersection для них, но пока (1.0.0) не реализует поиск объема такого объекта. Если бы вы могли найти вершины многогранника, то класс ConvexHull мог бы использоваться для вычисления объема. Но, как есть, кажется, что пространственный модуль SciPy не поможет. Возможно в будущей версии SciPy...

Источник

Интегрирование характеристической функции математически правильно, но не практично. Гораздо более подходящий подход состоит в том, чтобы сначала создать дискретизированную версию домена, а затем просто суммировать объемы маленьких тетраэдров.

Дискретизация в 3D может быть сделана, например, с помощью pygalmesh (проект сопряжения шахты CGAL). Код ниже дискретизирует отрезанный куб до

Затем файл читается с помощью meshio (другой мой проект), суммирование объемов тетраэдра дает

0.8050894798758184

Вы можете увеличить точность, уменьшив cell_size и / или edge_size, но создание сетки займет больше времени. Кроме того, вы можете указать "ребра объектов" для точного разрешения ребер пересечения. Это дало бы вам точный результат даже с самым крупным размером ячейки.

import pygalmesh
import numpy
import meshio


c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1])
h = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 4.0, 10.0)
u = pygalmesh.Intersection([c, h])
pygalmesh.generate_mesh(u, 'out.mesh', cell_size=3.0e-2, edge_size=3.0e-2)

points, cells, *_ = meshio.read('out.mesh')

def compute_tet_volumes(vertices, tets):
    cell_coords = vertices[tets]
    a = cell_coords[:, 1, :] - cell_coords[:, 0, :]
    b = cell_coords[:, 2, :] - cell_coords[:, 0, :]
    c = cell_coords[:, 3, :] - cell_coords[:, 0, :]
    # omega = <a, b x c>
    omega = numpy.einsum('ij,ij->i', a, numpy.cross(b, c))
    # https://en.wikipedia.org/wiki/Tetrahedron#Volume
    return abs(omega) / 6.0

print(numpy.sum(compute_tet_volumes(points, cells['tetra'])))

Посмотрите на pygalmesh (проект мой). Это позволит вам с легкостью создавать сложные трехмерные сетки, и, как только вы получите сетку, нужно просто суммировать объемы тетраэдров.

Источник

Если мы предположим, что граница полупространства задается как $\{(x, y, z) \mid ax + + cz + d = 0 \}$ с $c \not= 0$, и что половина пространство интереса состоит в том, что ниже плоскости (в $ z $ -направлении) ваш интеграл определяется

scipy.integrate.tplquad(lambda z, y, x: 1,
                        -0.5, 0.5,
                        lambda x: -0.5, lambda x: 0.5,
                        lambda x, y: -0.5, lambda x, y: max(-0.5, min(0.5, -(b*y+a*x+d)/c)))

Поскольку хотя бы одно из $ a $, $ b $ и $ c $ должно быть ненулевым, случай $c = 0$ может быть обработан путем изменения координат.

Источник
Другие вопросы по тегам