Физическое моделирование разрушения с помощью opengl / c ++

Question

Я пытаюсь реализовать идеи из этой статьи для моделирования перелома: http://graphics.berkeley.edu/papers/Obrien-GMA-1999-08/index.html

Я застрял в какой-то момент (по сути, страница 4 ...) и был бы очень признателен за любую помощь Часть, на которой я застрял, включает в себя деформацию тетраэдра (с использованием FEM).

У меня есть один тетраэдр, определяемый четырьмя узлами (каждый узел имеет положение x, y, z), в котором я вычисляю следующие матрицы из:

• u: каждый столбец - это вектор, содержащий материальные координаты (x, y, z, 1) для каждого узла (итого 4 столбца) матрица 4х4

• B: обратный (и), он называет это базисной матрицей, матрицей 4х4

• P: каждый столбец - это вектор, содержащий координаты реального мира (x, y, z) для каждого узла я устанавливаю P изначально равным u, так как объект не деформируется в состоянии покоя, матрица 3х4

• V: дают некоторые начальные скорости для (x, y, z) в каждом узле, так что 3x4 матрица

• дельта: в основном единичная матрица, {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}, {0, 0, 0}}

Я получаю x(u) = P*B*u и v(u) = V*B*u, но не уверен, где их использовать ...

Кроме того, я получаю dx = P*B*delta и dv = V*B*delta

Я тогда получаю деформацию по тензору деформации Грина, epsilon = 1/2(dx+transpose(dx)) - Identity_3x3

А потом стресс, sigma = lambda*trace(epsilon)*Identity_3x3 + 2*mu*epsilon

Я получаю силу упругости по уравнению (24) на странице 4 статьи. Это просто большое суммирование.

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

Проблема, однако, в том, что сила невероятно мала ... что-то х 10 ^ -19 и т. Д. Итак, c ++ обычно округляется до 0. Я прошел вычисления и не могу понять, почему.

Я знаю, что я что-то здесь упускаю, просто не могу понять, что. Какое обновление я делаю неправильно?

Answer 1

Общая причина, по которой сила мала, состоит в том, что модуль Юнга (лямбда) слишком мал.Если вы используете шкалу метров, объект макромасштаба может иметь 10 ^ 5 модулей Юнга и коэффициент Пуассона от 0,3 до 4.

Похоже, что ваш тет все еще находится в остальной конфигурации.При отсутствии деформации деформация будет равна нулю, и, следовательно, напряжение и сила также будут примерно равны нулю.Вы можете нарушить вершины различными способами и убедиться, что ваше напряжение (эпсилон) вычисляется правильно.Один простой тест - масштабирование на 2 относительно центроида, которое должно дать вам положительное напряжение.Если вы масштабируете на 0,5 относительно центроида, вы получите отрицательное напряжение.Если вы переводите вершины равномерно, вы не получите никаких изменений в деформации (общий инвариант FEM).Если вы повернете их, вы, вероятно, получите изменение, но составная модель совместного вращения не будет.

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

Вам определенно не нужно использовать арифметику произвольной точности для примеров в статье.Фактически, для этих типов симуляции достаточно поплавков.

Answer 2

Я могу ошибаться, но двойные числа в c ++ идут только до 15 десятичных знаков (по крайней мере, так говорит мой std :: numeric_limits).Таким образом, вы далеко от точности.

Таким образом, вам может понадобиться библиотека для произвольной точности арифметики, например, http://gmplib.org/