Detailed Description

Implementation of an implicit backward euler solver compatible with non-linear materials.

We are trying to solve to following

$\begin{eqnarray*} \mat{M} \ddot{\vect{x}} + \mat{C} \dot{\vect{x}} + \vect{R}(\vect{x}) = \vect{P} \end{eqnarray*}$

Where $\mat{M}$ is the mass matrix, $\mat{C}$ is the damping matrix, $\vect{R}$ is the (possibly non-linear) internal elastic force residual and $\vect{P}$ is the external force vector (for example, gravitation force or surface traction).

We first transform this second-order differential equation to a first one by introducing two independant variables:

$\begin{eqnarray*} \vect{v} &= \dot{\vect{x}} \\ \vect{a} &= \ddot{\vect{x}} \end{eqnarray*}$

Using the backward Euler scheme, we pose the following approximations:

$\begin{align} \vect{x}_{n+1} &= \vect{x}_{n} + h \vect{v}_{n+1} \label{eq:backward_euler_x_step} \\ \vect{v}_{n+1} &= \vect{v}_{n} + h \vect{a}_{n+1} \label{eq:backward_euler_v_step} \end{align}$

where is the delta time between the steps and .

Substituting inside gives

$\begin{eqnarray*} \vect{x}_{n+1} &= \vect{x}_{n} + h \left[ \vect{v}_{n} + h \vect{a}_{n+1} \right] \\ &= \vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1} \end{eqnarray*}$

And the problem becomes:

$\begin{eqnarray*} \mat{M} \vect{a}_{n+1} + \mat{C} \left[ \vect{v}_{n} + h \vect{a}_{n+1} \right] + \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}) = \vect{P}_n \end{eqnarray*}$

where $\vect{a}_{n+1}$ is the vector of unknown accelerations.

Finally, assuming $\vect{R}$ is non-linear in $\vect{x}_{n+1}$ , we iteratively solve for $\vect{a}_{n+1}$ using the Newton-Raphson method and using the previous approximations to back propagate it inside the velocity and position vectors.

Let be the Newton iteration number for a given time step , we pose

$\begin{align*} \vect{F}(\vect{a}_{n+1}^i) &= \mat{M} \vect{a}_{n+1}^i + \mat{C} \left[ \vect{v}_{n} + h \vect{a}_{n+1}^i \right] + \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}^i) - \vect{P}_n \\ \mat{J} = \frac{\partial \vect{F}}{\partial \vect{a}_{n+1}} \bigg\rvert_{\vect{a}_{n+1}^i} &= \mat{M} + h \mat{C} + h^2 \mat{K}(\vect{a}_{n+1}^i) \end{align*}$

where $\vect{x}_{n}$ and $\vect{x}_{n}$ are the position and velocity vectors at the beginning of the time step, respectively, and remains constant throughout the Newton iterations.

We then solve for $\vect{a}_{n+1}^{i+1}$ with

$\begin{align*} \mat{J} \left [ \Delta \vect{a}_{n+1}^{i+1} \right ] &= - \vect{F}(\vect{a}_{n+1}^i) \\ \vect{a}_{n+1}^{i+1} &= \vect{a}_{n+1}^{i} + \Delta \vect{a}_{n+1}^{i+1} \end{align*}$

And propagate back the new acceleration using and .

In addition, this component implicitly adds a Rayleigh's damping matrix $\mat{C}_r = r_m \mat{M} + r_k \mat{K}(\vect{x}_{n+1})$ . We therefore have

$\begin{align*} \vect{F}(\vect{a}_{n+1}^i) &= \mat{M} \vect{a}_{n+1}^i + \mat{C} \left[ \vect{v}_{n} + h \vect{a}_{n+1}^i \right] + \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}^i) - \vect{P}_n \\ &= \mat{M} \vect{a}_{n+1}^i + (\mat{C}_r+\mat{C}) \left[ \vect{v}_{n} + h \vect{a}_{n+1}^i \right] + \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}^i) - \vect{P}_n \\ &= \mat{M} \vect{a}_{n+1}^i + (r_m\mat{M}+r_k\mat{K}) \left[ \vect{v}_{n} + h \vect{a}_{n+1}^i \right] + \mat{C} \left[ \vect{v}_{n} + h \vect{a}_{n+1}^i \right] + \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}^i) - \vect{P}_n \\ &= \left[ (1 + hr_m)\mat{M} + h\mat{C} + hr_k\mat{K} \right] \vect{a}_{n+1}^i + \left[ r_m\mat{M} + \mat{C} + r_k \mat{K} \right] \vect{v}_n + \left[ \vect{R}(\vect{x}_{n} + h \vect{v}_{n} + h^2 \vect{a}_{n+1}^i) - \vect{P}_n \right] \\ \mat{J} = \frac{\partial \vect{F}}{\partial \vect{a}_{n+1}} \bigg\rvert_{\vect{a}_{n+1}^i} &= (1 + hr_m)\mat{M} + h \mat{C} + h(h+r_k) \mat{K}(\vect{a}_{n+1}^i) \end{align*}$

#include <BackwardEulerODESolver.h>

Inheritance diagram for SofaCaribou::ode::BackwardEulerODESolver:

Public Types
template<typename T >
using	Data = sofa::core::objectmodel::Data< T >

Public Types inherited from SofaCaribou::ode::NewtonRaphsonSolver
template<typename T >
using	Data = sofa::core::objectmodel::Data< T >

template<typename T >
using	Link = sofa::core::objectmodel::SingleLink< NewtonRaphsonSolver, T, sofa::core::objectmodel::BaseLink::FLAG_STRONGLINK >

Public Member Functions
	SOFA_CLASS (BackwardEulerODESolver, NewtonRaphsonSolver)

void	solve (const sofa::core::ExecParams *params, SReal dt, sofa::core::MultiVecCoordId x_id, sofa::core::MultiVecDerivId v_id) override

Public Member Functions inherited from SofaCaribou::ode::NewtonRaphsonSolver
	SOFA_CLASS (NewtonRaphsonSolver, sofa::core::behavior::OdeSolver)

void	init () override

void	solve (const sofa::core::ExecParams *params, SReal dt, sofa::core::MultiVecCoordId x_id, sofa::core::MultiVecDerivId v_id) override

auto	iteration_times () const -> const std::vector< UNSIGNED_INTEGER_TYPE > &
	List of times (in nanoseconds) that each Newton-Raphson iteration took to compute in the last call to Solve().

auto	squared_residuals () const -> const std::vector< FLOATING_POINT_TYPE > &
	The list of squared residual norms (\|\|r\|\|^2) of every newton iterations of the last solve call.

auto	squared_initial_residual () const -> const FLOATING_POINT_TYPE &
	The initial squared residual (\|\|r0\|\|^2) of the last solve call.

The documentation for this class was generated from the following files:

/home/jnbrunet/sources/caribou/src/Plugin/Ode/BackwardEulerODESolver.h
/home/jnbrunet/sources/caribou/src/Plugin/Ode/BackwardEulerODESolver.cpp

Detailed Description

Public Types

Public Member Functions