SemiImplicit.h

Go to the documentation of this file.

/*  Copyright (C) 2003-2021 Free Electron Organization
    Any use of this software requires a license.  If a valid license
    was not distributed with this file, visit freeelectron.org. */

/** @file */

#ifndef __solve_SemiImplicit_h__
#define __solve_SemiImplicit_h__

#include "signal/signal.h"
#include "datatool/datatool.h"
#include "shape/shape.h"

#define PARTICLE_DEBUG
//#define USE_FILTERS

namespace fe
{
namespace ext
{

typedef FE_UWORD            FE_UWORD;


/** Semi Implicit time integration

    @copydoc SemiImplicit_info
    */
class FE_DL_EXPORT SemiImplicit : public Handled<SemiImplicit>
{
    public:
                SemiImplicit(void);
virtual         ~SemiImplicit(void);

virtual void    compile(sp<RecordGroup> rg_input);
virtual void    initialize(sp<Scope> a_spScope);
virtual void    extract(sp<RecordGroup> rg_output);
virtual void    prestep(void);
virtual void    step(t_solve_real a_timestep, t_solve_v3 &a_totalConstraintForce);

virtual void    setGravity(t_solve_v3 a_g) { m_gravity = a_g; }
virtual void    setDirectSolve(bool a_flag) { m_direct_solve = a_flag; }
virtual void    set2D(bool a_flag) { m_2D = a_flag; }
virtual void    setRayleighDamping(bool a_flag) { m_rayleigh_damping = a_flag; }
virtual void    setRayleighDamping(t_solve_real a_stiffness, t_solve_real a_mass)
                {
                    m_rayleigh_stiffness = a_stiffness;
                    m_rayleigh_mass = -a_mass;
                }

    class Particle
    {
        public:
            t_solve_v3          m_location;
            t_solve_v3          m_estimated;
            t_solve_v3          m_velocity;
            t_solve_v3          m_force;
            t_solve_v3          m_force_weak;
            t_solve_v3          m_force_external;
            t_solve_real        m_mass;
            t_solve_v3          m_prev_location;
            t_solve_v3          m_prev_velocity;
            t_solve_v3          m_constraint_force;
#ifdef PARTICLE_DEBUG
            Color                   m_color;
#endif
    };

        std::vector<Particle>   &particles(void) { return m_particles; }
        bool                    lookupIndex(unsigned int &a_particle, Record &r_particle);

    class FE_DL_EXPORT CompileMatrix
    {
        public:
                                            CompileMatrix(void);
                                            ~CompileMatrix(void);

            typedef t_solve_matrix                          **t_ppBlock;
            typedef std::vector<t_ppBlock>                  t_dfdx_array;
            typedef std::vector<t_ppBlock>                  t_dfdv_array;
            typedef std::pair<t_dfdx_array, t_dfdv_array>   t_entry;
            typedef std::map<unsigned int, t_entry>         t_row;

            std::vector<t_row>              m_rows;

            void                            clear(void);
            void                            setRows(unsigned int a_count);
            unsigned int                    rows(void);
            t_row                           &row(unsigned int a_index);
            t_entry                         &entry( unsigned int a_i,
                                                    unsigned int a_j);

            typedef std::set<unsigned int>                  t_nonzero_set;
            typedef std::vector< t_nonzero_set >            t_nonzero_pattern;

            void                            symbolicFill(void);
    };

    //typedef std::size_t t_size;
    typedef unsigned int t_size;
    typedef std::pair<t_size, t_size> t_pair;
    typedef std::vector<t_pair> t_pairs;
    class Force : public Counted, public CastableAs<Force>
    {
        public:
                    Force(void){}
        virtual     ~Force(void){}

        virtual     void clear(void){}
        virtual     void accumulate(void){}
        virtual     bool validate(void){ return true; }

        virtual     void compile(   sp<RecordGroup> rg_input,
                                    std::vector<Particle> &a_particles,
                                    CompileMatrix &a_compileMatrix){}
        virtual     void precompile(sp<RecordGroup> rg_input) {}
        virtual     void pairs(     sp<RecordGroup> rg_input,
                                    t_pairs &a_pairs) {}
        virtual     void constrain( std::vector< FilterConstraint<t_solve_real> > &a_filters,
                                    std::vector<t_solve_v3> &a_dv,
                                    std::vector<t_solve_v3> &a_dx,
                                    t_solve_real a_dt){}
    };

virtual void    addForce(sp<Force> a_force, bool a_add_damping=false);



    private:

        void reorder(std::vector<unsigned int> &a_order, t_pairs &a_pairs);

        std::vector<Particle>       m_particles;

        AsParticle              m_asParticle;
        AsSolverParticle        m_asSolverParticle;
        AsLineConstrained       m_asLineConstrained;
        AsPlaneConstrained      m_asPlaneConstrained;
        AsComponent             m_asComponent;
        AsForcePoint            m_asForcePoint;
        AsTemporal              m_asTemporal;
        AsValidate              m_asValidate;
        AsAccumulate            m_asAccumulate;
        AsClear                 m_asClear;
        AsUpdate                m_asUpdate;
#ifdef PARTICLE_DEBUG
        AsColored               m_asColored;
#endif
        AsForceFilter           m_asForceFilter;
        hp<SignalerI>           m_hpSignaler;

        unsigned int                    m_n;
        unsigned int                    m_n_sim;
        std::vector< sp<Force> >        m_forces_add_damping;
        std::vector< sp<Force> >        m_forces_as_is;

        t_solve_matrix                  m_dummy_block;

        sp< UpperTriangularBCRS<t_solve_real> >         m_dfdx;
        sp< UpperTriangularBCRS<t_solve_real> >         m_dfdv;
        sp< UpperTriangularBCRS<t_solve_real> >         m_lhs;
        sp< UpperTriangularBCRS<t_solve_real> >         m_lhs_snapshot;

        std::vector<t_solve_v3>     m_rhs;
        std::vector<t_solve_v3>     m_dv;
        std::vector<t_solve_v3>     m_tmp;

        t_solve_v3                  m_gravity;

        std::map<FE_UWORD, unsigned int>    m_recordToParticle;
        //std::map<unsigned int, unsigned int>  m_recordToParticle;

        BlockDiagonalPreconditioner<t_solve_real>       m_preconditioner;
        //OverlappedInversionPreconditioner m_preconditioner2;
        BlockPCG<t_solve_real>                      m_solver;

        t_solve_real                    m_dv2dxRatio;
        t_solve_real                    m_dxImplicitness;
        t_solve_real                    m_dvImplicitness;

        t_solve_real                    m_ratio;
        unsigned int                    m_subdivcnt;
        unsigned int                    m_subdivsz;
        t_solve_real                    m_subdivmult;

        bool                            m_direct_solve;
        bool                            m_2D;

        bool                            m_rayleigh_damping;
        t_solve_real                    m_rayleigh_stiffness;
        t_solve_real                    m_rayleigh_mass;
        std::vector<t_solve_v3>         m_perturb;
#ifndef USE_FILTERS
        std::vector<FilterConstraint< t_solve_real> > m_dummy_filters;
#endif
};

} /* namespace */
} /* namespace */


#endif /* __solve_SemiImplicit_h__ */