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

Go to the documentation of this file.

#ifndef OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_H
#define OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_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>

namespace OpenTissue
{
  namespace fem
  {
    namespace detail
    {
      template < typename fem_mesh, typename real_type >
      inline void dynamics_assembly(
        fem_mesh & mesh,
        real_type const & mass_damping,
        real_type const & dt
        )
      {
        typedef typename fem_mesh::vector3_type                  vector3_type;
        typedef typename fem_mesh::matrix3x3_type                matrix3x3_type;
        typedef typename fem_mesh::node_iterator                 node_iterator;
        typedef typename fem_mesh::node_type::matrix_iterator    matrix_iterator;

        node_iterator begin = mesh.node_begin();
        node_iterator end   = mesh.node_end();
        for (node_iterator n_i = begin ; n_i != end; ++n_i)
        {
          unsigned int i     =  n_i->idx();
          vector3_type & b_i =  n_i->m_b;
          real_type & m_i    =  n_i->m_mass;

          b_i.clear();

          matrix_iterator Kbegin = n_i->Kbegin();
          matrix_iterator Kend   = n_i->Kend();
          for (matrix_iterator K = Kbegin; K != Kend;++K)
          {
            unsigned int     j    = K->first;
            matrix3x3_type & K_ij = K->second;
            node_iterator    n_j  = mesh.node(j);
            vector3_type &   x_j  = n_j->m_coord;
            matrix3x3_type & A_ij = n_i->A(j);

            A_ij = K_ij * (dt*dt);
            b_i -= K_ij * x_j;
            if (i == j)
            {
              real_type c_i = mass_damping*m_i;
              real_type tmp = m_i + dt*c_i;
              A_ij(0,0) += tmp; A_ij(1,1) += tmp;  A_ij(2,2) += tmp;
            }
          }
          b_i -= n_i->m_f0;
          b_i += n_i->m_f_external;
          b_i *= dt;
          b_i += n_i->m_velocity * m_i;
        }
      }

    } // namespace detail
  } // namespace fem
} // namespace OpenTissue

//OPENTISSUE_DYNAMICS_FEM_FEM_DYNAMICS_ASSEMBLY_H
#endif