TPZSBFemElementGroup.h

Go to the documentation of this file.

//
//  TPZSBFemElementGroup.hpp
//  PZ
//
//  Created by Philippe Devloo on 4/4/16.
//
//

#ifndef TPZSBFemElementGroup_hpp
#define TPZSBFemElementGroup_hpp

#include <stdio.h>

#include "pzelementgroup.h"
#include "TPZSBFemVolume.h"


class TPZSBFemElementGroup : public TPZElementGroup
{
    
public:
    enum EComputationMode {EStiff, EOnlyMass, EMass};
    
private:
    
    TPZFMatrix<std::complex<double> > fPhi;
    
    TPZFMatrix<std::complex<double> > fPhiInverse;
    
    TPZManVector<std::complex<double> > fEigenvalues;
    
    TPZFMatrix<std::complex<double> > fCoef;
    
    TPZFMatrix<STATE> fMassMatrix;
    
    EComputationMode fComputationMode = EStiff;
    
    REAL fMassDensity = 1.;
    
    REAL fDelt = 1.;
    
    void ComputeMassMatrix(TPZElementMatrix &M0);
    
public:
    
    TPZSBFemElementGroup() : TPZElementGroup()
    {
        
    }
    
    TPZSBFemElementGroup(TPZCompMesh &mesh, int64_t &index) : TPZElementGroup(mesh,index)
    {
        
    }
    
    virtual void AddElement(TPZCompEl *cel);
    

    void ComputeMatrices(TPZElementMatrix &E0, TPZElementMatrix &E1, TPZElementMatrix &E2, TPZElementMatrix &M0);
    
    virtual void CalcStiff(TPZElementMatrix &ek,TPZElementMatrix &ef);
    
    void SetDensity(REAL density)
    {
        fMassDensity = density;
    }
    void SetComputeOnlyMassMatrix()
    {
        if(fMassMatrix.Rows() == 0)
        {
            DebugStop();
        }
        fComputationMode = EOnlyMass;
    }
    
    void SetComputeTimeDependent(REAL delt)
    {
        fDelt = delt;
        fComputationMode = EMass;
    }
    
    void SetComputeStiff()
    {
        fComputationMode = EStiff;
    }
    virtual void Print(std::ostream &out = std::cout) const
    {
        out << __PRETTY_FUNCTION__ << std::endl;
        TPZElementGroup::Print(out);
        int nel = fElGroup.size();
        out << "Element indexes of the volume elements ";
        for (int el=0; el<nel; el++) {
            out << fElGroup[el]->Index() << " ";
        }
        out << std::endl;
        out << "Indices of the associated computational skeleton elements\n";
        for (int el=0; el<nel; el++) {
            TPZCompEl *cel = fElGroup[el];
            TPZSBFemVolume *vol = dynamic_cast<TPZSBFemVolume *>(cel);
            if(!vol) DebugStop();
            out << vol->SkeletonIndex() << " ";
        }
        out << std::endl;
        out << "Connect indexes of the contained elements\n";
        for (int el=0; el<nel; el++) {
            TPZCompEl *cel = fElGroup[el];
            int nc = cel->NConnects();
            for (int ic=0; ic<nc; ic++) {
                out << cel->ConnectIndex(ic) << " ";
            }
            out << std::endl;
        }

/*
        out << "EigenVectors for displacement\n";
        fPhi.Print("Phi = ",out,EMathematicaInput);
        out << "Inverse EigenVectors\n";
        fPhiInverse.Print("PhiInv = ",out,EMathematicaInput);
        out << "EigenValues " << fEigenvalues << std::endl;
        out << "Mass Matrix\n";
        fMassMatrix.Print("Mass = ",out);
        out << "Solution Coeficients\n";
        fCoef.Print("Coef ",out);
        for (int el=0; el<nel; el++) {
            fElGroup[el]->Print(out);
        }
 */
        out << "End of " << __PRETTY_FUNCTION__ << std::endl;
    }
    

    
    virtual void CalcResidual(TPZElementMatrix &ef)
    {
        TPZElementMatrix ek(Mesh(),TPZElementMatrix::EK);
        CalcStiff(ek,ef);
    }
    

    virtual void LoadSolution();

    virtual void LoadElementReference()
    {
        for (int64_t i = 0; i < fElGroup.size(); i++) {
            fElGroup[i]->LoadElementReference();
        }
    }
    


    int64_t NumEigenValues()
    {
        return fEigenvalues.size();
    }
    
    void LoadEigenVector(int64_t eig);
    
    TPZVec<STATE> MultiplyingCoeficients()
    {
        int nel = fCoef.Rows();
        TPZVec<STATE> result(nel,0.);
        for (int ir=0; ir<nel; ir++) {
            result[ir] = fCoef(ir,0).real();
        }
        return result;
    }
    
    TPZVec<std::complex<double> > &EigenValues()
    {
        return fEigenvalues;
    }
    
    TPZFMatrix<std::complex<double> > &Phi()
    {
        return fPhi;
    }
    
    TPZFMatrix<std::complex<double> > &PhiInverse()
    {
        return fPhiInverse;
    }
    
    TPZFMatrix<STATE> &MassMatrix()
    {
        return fMassMatrix;
    }
    
    TPZFMatrix<std::complex<double> > Coeficients()
    {
        return fCoef;
    }
    
    TPZFMatrix<double> PhiReal()
    {
        int64_t rows = fPhi.Rows(),cols = fPhi.Cols();
        TPZFMatrix<double> phireal(rows,cols);
        for(int64_t i=0; i<rows; i++)
        {
            for(int64_t j=0; j<cols; j++)
            {
                phireal(i,j) = fPhi(i,j).real();
            }
        }
        return phireal;
    }
    
    TPZManVector<double> EigenvaluesReal()
    {
        int64_t nel = fEigenvalues.NElements();
        TPZManVector<double> eig(nel);
        for(int64_t el=0; el<nel; el++)
        {
            eig[el] = fEigenvalues[el].real();
        }
        return eig;
    }
    
    TPZFMatrix<double> CoeficientsReal()
    {
        int64_t rows = fCoef.Rows(),cols = fCoef.Cols();
        TPZFMatrix<double> coefreal(rows,cols);
        for(int64_t i=0; i<rows; i++)
        {
            for(int64_t j=0; j<cols; j++)
            {
                coefreal(i,j) = fCoef(i,j).real();
            }
        }
        return coefreal;
    }

};

#endif /* TPZSBFemElementGroup_hpp */