EMM_DisplacementFEMOperator.h

Go to the documentation of this file.

#ifndef EMM_DisplacementFEMOperator_H
#define EMM_DisplacementFEMOperator_H

#include "EMM_DisplacementOperator.h"

#include "EMM_MassMatrix.h"

#include "MATH_Solver.h"
#include "MATH_Vector.h"

DEFINE_SPTR(EMM_DisplacementFEMOperator);
class EMM_DisplacementFEMOperator : public EMM_DisplacementOperator   {
    
    SP_OBJECT(EMM_DisplacementFEMOperator);
    // ATTRIBUTES
public:

    static const tFlag CANONICAL_MASS_MATRIX=0;
    static const tFlag BLOCK_MASS_MATRIX=1;
    static const tFlag CONDENSED_MASS_MATRIX=2;
    
private:
   
    
    /* \brief compute \f$ \int_{\hat \omega} \frac{\partial \hat \Phi_i}{\partial r_j} d{\hat \omega}=|\hat \omega| . INT_GRAD_PHI[i*3+j]/ (H_j . INT_GRAD_PHI_FACTOR)\f$
     */
    static const tCInt INT_GRAD_PHI[];
    static const tUCInt INT_GRAD_PHI_FACTOR;

    // ASSOCIATION
 
private:

    tFlag mMassMatrixType;
    
    //accelerator matrix
    SP::EMM_MassMatrix mMassMatrix;

    
    
protected: 

    // METHODS
  
    // CONSTRUCTORS 
  
    EMM_DisplacementFEMOperator(void);
  
    // DESTRUCTORS 
    
    virtual ~EMM_DisplacementFEMOperator(void);   
  
  
public:
    
    //limit conditions input data
    //===========================
 
  
protected:
   
    virtual void buildDataOnBoundaryFaces(const EMM_Grid3D& mesh,
                                          const EMM_LimitConditionArray& limitConditionOnPoints,
                                          const EMM_RealField& U0,
                                          EMM_RealField& DnU0) {

        //temporaray vector to compute \f$ \int_{\partial \omega_1} C. d\sigma \f$
        SP::EMM_RealField C=DnU0.NewInstance();
        if (buildDataOnNeumannBoundaryFaces(mesh, limitConditionOnPoints,DnU0,*C.get())) {
            //not null constraint
            DnU0.copy(*C.get());
        } else {
            //null constraint
            DnU0.setSize(0);
        }
             
    }
    
private:
    tBoolean buildDataOnNeumannBoundaryFaces(const EMM_Grid3D& mesh,
                                             const EMM_LimitConditionArray& limitConditionOnPoints,
                                             const EMM_RealField& DnU0,
                                             EMM_RealField& C) const;
    
    //discretization mehods
    //======================
public:
    virtual tBoolean discretize(const EMM_LandauLifschitzSystem& system);

    virtual tBoolean resetToInitialState(const EMM_LandauLifschitzSystem& system);

    
    inline void setMassMatrixType(const tFlag& t) {
        mMassMatrixType=t;
    }
    inline const tFlag&  getMassMatrixType() const {
        return mMassMatrixType;
    }
    
    virtual SP::EMM_MassMatrix NewMassMatrix() const;
 
    
    inline const tUIndex* getPointsNumber() const {
        return mMassMatrix->getPointsNumber();
    }

    inline const tCellFlag& getPeriodicIndicator() const {
        return mMassMatrix->getPeriodicIndicator();
    }

    inline const EMM_CellFlagArray& getNeighborsIndicator() const {
        return mMassMatrix->getNeighborsIndicator();
    }
    
    virtual tULLInt getMemorySize() const {
        return EMM_DisplacementOperator::getMemorySize()+
            mMassMatrix->getMemorySize();
    }
    
    //output data
    //=============

    virtual tBoolean getDataFieldSpace(const tUSInt& index,tString& dataName,tFlag& supportType,tString& dFieldType,tUIndex& n,tDimension& dim) const {
        supportType=EMM_Grid3D::POINT;
        dFieldType=CORE_Object::getTypeName<tReal>();
        switch(index) {
        case 0:
            dim=getDisplacement().getDimension();
            dataName="U";
            n=getDisplacement().getSize();
            return true;
            break;
        case 1:
            dataName="dU_dt";
            dim=getVelocity().getDimension();
            n=getVelocity().getSize();
            return true;
            break;
        case 2:
            dataName="acc";
            dim=getAccelerator().getDimension();
            n=getAccelerator().getSize();
            return true;
            break;
        }
        return false;
    }
   
   

    
   

    //time evolutions of fields for unsteady state
    //=============================================

protected:

    
    virtual tBoolean solveAcceleratorSystem(EMM_RealField& V);
    
    virtual void spaceProjection(EMM_RealField& V) const;
    
    virtual void spaceProjection(const EMM_RealField& V0, EMM_RealField& V) const;
     
protected:
    virtual void spaceRelevant(EMM_RealField& V) const;
    
    //steady state method
    //===================

    
public:
    
    virtual  void computeEquilibriumMatrixDiagonalConditioner(MATH_Vector& D) const;
  
 

    
    
    //elastic tensor methods
    //======================
    
protected:
     
    virtual void computeElasticTensor(const tUIndex& nData,
                                      const tDimension& dim,
                                      const tReal* U,
                                      EMM_2PackedSymmetricTensors& eTensor) const;
    
    
    //stress methods
    //==============
protected:
    virtual void computeElasticStress(const tBoolean& withConstraints,
                                      const tReal& beta,
                                      const tUIndex& nData,const tDimension& dim, const tReal* U,tReal* D) const;
    
    
private:
    void addBoundaryElasticStress(const tReal& beta,const tUIndex& nPoints,const tDimension& dim,tReal* D) const;

protected:
    virtual void computeMagneticStress(const tReal& alpha,const tReal& beta,
                                       const tUIndex& nCells,const tDimension& dim, const EMM_RealArray& sigma,
                                       const tReal* M,
                                       const tUIndex& nData,tReal* D) const;

  
    

    //energy methods
    //===============
public:
    virtual tReal computeCineticEnergy(const EMM_RealField& V) const;
    
    /* ! \brief compute the  potential energy to the space variation of the displacement 
     * @param U: displacement field
     * @return potential energy value \f$ \int_\omega (\lambda^e:\varepsilon(u)):\varepsilon(u) d\omega \f$ 
     * 
     */
    virtual tReal computePotentialEnergy(const EMM_RealField& U) const;
    
    tReal  computeStressConstraintEnergy(const EMM_RealField& Uf) const;
    
    
    //STRING representation
    //======================
public:
    virtual tString toString() const {
        return EMM_DisplacementFEMOperator::toString();
    }
    
    
   
};

#endif