/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/dynamics/fem/fem_simulate.h

Go to the documentation of this file.

#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
#define OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
//
// OpenTissue Template Library
// - A generic toolbox for physics-based modeling and simulation.
// Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
//
// OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
//
#include <OpenTissue/configuration.h>

#include <OpenTissue/dynamics/fem/fem_clear_stiffness_assembly.h>
#include <OpenTissue/dynamics/fem/fem_update_orientation.h>
#include <OpenTissue/dynamics/fem/fem_reset_orientation.h>
#include <OpenTissue/dynamics/fem/fem_stiffness_assembly.h>
#include <OpenTissue/dynamics/fem/fem_add_plasticity_force.h>
#include <OpenTissue/dynamics/fem/fem_dynamics_assembly.h>
#include <OpenTissue/dynamics/fem/fem_conjugate_gradients.h>
#include <OpenTissue/dynamics/fem/fem_position_update.h>

namespace OpenTissue
{
  namespace fem
  {
    template < typename fem_mesh, typename real_type >
    inline void simulate(
      fem_mesh & mesh
      , real_type const & time_step
      , bool use_stiffness_warping
      )
    {
      // Some theory:
      //
      //
      // Implicit discretization of
      //
      //   M d^2x/dt^2 + C dx/dt + K (x-x0) = f_ext
      //
      // Evaluate at (i+1), and use d^2x/dt^2 = (v^{i+1} - v^i)/timestep and  dx/dt = v^{i+1}
      //
      //     M (v^{i+1} - v^i)/timestep + C v^{i+1} + K (x^{i+1}-x0) = f_ext
      //
      // and    x^{i+1} = x^i + v^{i+1}*timestep
      //
      //     M (v^{i+1} - v^i)/timestep + C v^{i+1} + K ( (x^i + v^{i+1}*timestep) -x0) = f_ext
      //
      //     M (v^{i+1} - v^i)/timestep + C v^{i+1} + K*x^i + timestep*K*v^{i+1} - K x0 = f_ext
      //
      //     M v^{i+1} - M v^i + timestep*C v^{i+1} + timestep*timestep*K*v^{i+1}  = timestep * (f_ext - K*x^i  + K x0)
      //
      //     (M  + timestep*C  + timestep*timestep*K) v^{i+1}  =  M v^i + timestep * (f_ext - K*x^i  + K x0)
      //
      // let f0 = -K x0
      //
      //     (M  + timestep*C  + timestep*timestep*K) v^{i+1}  =  M v^i + timestep * (f_ext - K*x^i  -f0)
      //
      //     (M  + timestep*C  + timestep*timestep*K) v^{i+1}  =  M v^i - timestep * (K*x^i  + f0 - f_ext)
      //
      // so we need to solve A v^{i+1} = b for v^{i+1}, where
      //
      //     A = (M  + timestep*C  + timestep*timestep*K)
      //     b =  M v^i - timestep * (K*x^i  + f0 - f_ext)
      //
      // afterwards we do a position update
      //
      //     x^{i+1} = x^i + v^{i+1}*timestep
      //
      // Notice that a fully implicit scheme requres K^{i+1}, however for linear elastic materials K is constant.

      detail::clear_stiffness_assembly(mesh.node_begin(),mesh.node_end());
      if(use_stiffness_warping)
        detail::update_orientation(mesh.tetrahedron_begin(),mesh.tetrahedron_end());
      else
        detail::reset_orientation(mesh.tetrahedron_begin(),mesh.tetrahedron_end());

      detail::stiffness_assembly(mesh.tetrahedron_begin(),mesh.tetrahedron_end());

      detail::add_plasticity_force(mesh.tetrahedron_begin(),mesh.tetrahedron_end(),time_step);

      real_type mass_damping = 2.0;  // TODO: Should be user controllable

      detail::dynamics_assembly(mesh,mass_damping,time_step);

      unsigned int min_iterations = 20;   // TODO: Should be user controllable
      unsigned int max_iterations = 20;   // TODO: Should be user controllable

      detail::conjugate_gradients(mesh, min_iterations, max_iterations);
      detail::position_update(mesh,time_step);
    }

  } // namespace fem
} // namespace OpenTissue

//OPENTISSUE_DYNAMICS_FEM_FEM_SIMULATE_H
#endif