HexahedronFEMForceField.h

Go to the documentation of this file.

/******************************************************************************
*       SOFA, Simulation Open-Framework Architecture, version 1.0 beta 3      *
*                (c) 2006-2008 MGH, INRIA, USTL, UJF, CNRS                    *
*                                                                             *
* This library is free software; you can redistribute it and/or modify it     *
* under the terms of the GNU Lesser General Public License as published by    *
* the Free Software Foundation; either version 2.1 of the License, or (at     *
* your option) any later version.                                             *
*                                                                             *
* This library is distributed in the hope that it will be useful, but WITHOUT *
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or       *
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License *
* for more details.                                                           *
*                                                                             *
* You should have received a copy of the GNU Lesser General Public License    *
* along with this library; if not, write to the Free Software Foundation,     *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.          *
*******************************************************************************
*                               SOFA :: Modules                               *
*                                                                             *
* Authors: The SOFA Team and external contributors (see Authors.txt)          *
*                                                                             *
* Contact information: contact@sofa-framework.org                             *
******************************************************************************/
#ifndef SOFA_COMPONENT_FORCEFIELD_HEXAHEDRONFEMFORCEFIELD_H
#define SOFA_COMPONENT_FORCEFIELD_HEXAHEDRONFEMFORCEFIELD_H

#include <sofa/core/componentmodel/behavior/ForceField.h>
#include <sofa/component/MechanicalObject.h>
#include <sofa/core/componentmodel/topology/BaseMeshTopology.h>
#include <sofa/component/topology/SparseGridTopology.h>
#include <sofa/helper/vector.h>
#include <sofa/defaulttype/Vec.h>
#include <sofa/defaulttype/Mat.h>

namespace sofa
{

  namespace component
  {

    namespace forcefield
    {

      using namespace sofa::defaulttype;
      using sofa::helper::vector;

      template<class DataTypes>
      class HexahedronFEMForceField : virtual public core::componentmodel::behavior::ForceField<DataTypes>, public virtual core::objectmodel::BaseObject
      {
      public:
        typedef typename DataTypes::VecCoord VecCoord;
        typedef typename DataTypes::VecDeriv VecDeriv;
        typedef VecCoord Vector;
        typedef typename DataTypes::Coord Coord;
        typedef typename DataTypes::Deriv Deriv;
    typedef typename Coord::value_type Real;

        typedef core::componentmodel::topology::BaseMeshTopology::index_type Index;
#ifdef SOFA_NEW_HEXA
        typedef core::componentmodel::topology::BaseMeshTopology::Hexa Element;
        typedef core::componentmodel::topology::BaseMeshTopology::SeqHexas VecElement;
#else
        typedef core::componentmodel::topology::BaseMeshTopology::Cube Element;
        typedef core::componentmodel::topology::BaseMeshTopology::SeqCubes VecElement;
#endif

        static const int LARGE = 0;   
        static const int POLAR = 1;   
        
        
        

      protected:
        //component::MechanicalObject<DataTypes>* object;
          

        typedef Vec<24, Real> Displacement;     

        typedef Mat<6, 6, Real> MaterialStiffness;  
        typedef vector<MaterialStiffness> VecMaterialStiffness;         
        VecMaterialStiffness _materialsStiffnesses;                 

        typedef Mat<24, 24, Real> ElementStiffness;
        typedef helper::vector<ElementStiffness> VecElementStiffness;
        Data<VecElementStiffness> _elementStiffnesses;
                        
        typedef Mat<3, 3, Real> Mat33;
        typedef Mat33 Transformation; 


        typedef std::pair<int,Real> Col_Value;
        typedef vector< Col_Value > CompressedValue;
        typedef vector< CompressedValue > CompressedMatrix;
        CompressedMatrix _stiffnesses;
        double m_potentialEnergy;


        sofa::core::componentmodel::topology::BaseMeshTopology* _mesh;
        topology::SparseGridTopology* _sparseGrid;
        const VecElement *_indexedElements;
        Data< VecCoord > _initialPoints; 
        

              Mat<8,3,int> _coef; 
#ifndef SOFA_NEW_HEXA
              static const int _indices[8]; 
#endif

public:
              
    int method;
        Data<std::string> f_method; 
        Data<Real> f_poissonRatio;
        Data<Real> f_youngModulus;
        Data<bool> f_updateStiffnessMatrix;
        Data<bool> f_assembling;


        HexahedronFEMForceField()
            : _elementStiffnesses(initData(&_elementStiffnesses,"stiffnessMatrices", "Stiffness matrices per element (K_i)"))
            , _mesh(NULL)
            , _sparseGrid(NULL)
            , _indexedElements(NULL)
            , _initialPoints(initData(&_initialPoints,"initialPoints", "Initial Position"))
            , f_method(initData(&f_method,std::string("large"),"method","\"large\" or \"polar\" displacements"))
            , f_poissonRatio(initData(&f_poissonRatio,(Real)0.45f,"poissonRatio",""))
            , f_youngModulus(initData(&f_youngModulus,(Real)5000,"youngModulus",""))
            , f_updateStiffnessMatrix(initData(&f_updateStiffnessMatrix,false,"updateStiffnessMatrix",""))
            , f_assembling(initData(&f_assembling,false,"assembling",""))
        {
            _coef[0][0]=-1;
            _coef[1][0]=1;
            _coef[2][0]=1;
            _coef[3][0]=-1;
            _coef[4][0]=-1;
            _coef[5][0]=1;
            _coef[6][0]=1;
            _coef[7][0]=-1;
            _coef[0][1]=-1;
            _coef[1][1]=-1;
            _coef[2][1]=1;
            _coef[3][1]=1;
            _coef[4][1]=-1;
            _coef[5][1]=-1;
            _coef[6][1]=1;
            _coef[7][1]=1;
            _coef[0][2]=-1;
            _coef[1][2]=-1;
            _coef[2][2]=-1;
            _coef[3][2]=-1;
            _coef[4][2]=1;
            _coef[5][2]=1;
            _coef[6][2]=1;
            _coef[7][2]=1;
            
            _alreadyInit=false;

        }

        void parse(core::objectmodel::BaseObjectDescription* arg);

        void setPoissonRatio(Real val) { this->f_poissonRatio.setValue(val); }

        void setYoungModulus(Real val) { this->f_youngModulus.setValue(val); }

                void setMethod(int val) 
                { 
                  method = val; 
                  switch(val)
                  {
                    case POLAR: f_method.setValue("polar"); break;
                    default   : f_method.setValue("large"); 
                  };
                }

        void setUpdateStiffnessMatrix(bool val) { this->f_updateStiffnessMatrix.setValue(val); }

        void setComputeGlobalMatrix(bool val) { this->f_assembling.setValue(val); }

        //  component::MechanicalObject<DataTypes>* getObject() { return object; }

        virtual void init();
        virtual void reinit();

        virtual void addForce (VecDeriv& f, const VecCoord& x, const VecDeriv& v);

        virtual void addDForce (VecDeriv& df, const VecDeriv& dx);

    virtual double getPotentialEnergy(const VecCoord& x);

        void addKToMatrix(sofa::defaulttype::BaseMatrix *mat, SReal k, unsigned int &offset);


        void draw();

      protected:

         
    virtual void computeElementStiffness( ElementStiffness &K, const MaterialStiffness &M, const helper::fixed_array<Coord,8> &nodes, const int elementIndice, double stiffnessFactor=1.0);
        Mat33 integrateStiffness( int signx0, int signy0, int signz0, int signx1, int signy1, int signz1, const Real u, const Real v, const Real w, const Mat33& J_1  );

        void computeMaterialStiffness(int i);
        
        void computeForce( Displacement &F, const Displacement &Depl, const ElementStiffness &K );


        vector<helper::fixed_array<Coord,8> > _rotatedInitialElements;   
        vector<Transformation> _rotations;
        void initLarge(int i, const Element&elem);
        void computeRotationLarge( Transformation &r, Coord &edgex, Coord &edgey);
        virtual void accumulateForceLarge( Vector& f, const Vector & p, int i, const Element&elem  );

        void initPolar(int i, const Element&elem);
        void computeRotationPolar( Transformation &r, Vec<8,Coord> &nodes);
        virtual void accumulateForcePolar( Vector& f, const Vector & p, int i, const Element&elem  );
        

        
        bool _alreadyInit;
      };

    } // namespace forcefield

  } // namespace component

} // namespace sofa

#endif