TPZSBFemVolume.h

Go to the documentation of this file.

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

#ifndef TPZSBFemVolume_hpp
#define TPZSBFemVolume_hpp

#include <stdio.h>
#include "pzcompel.h"
#include "pzelmat.h"
#include "pzinterpolationspace.h"

class TPZSBFemVolume : public TPZInterpolationSpace
{
    
    int64_t fElementGroupIndex = -1;
    
    TPZCompEl *fElementGroup = 0;
    
    int64_t fSkeleton = -1;
    
    TPZIntPoints *fIntRule = 0;
    
    TPZFNMatrix<30,std::complex<double> > fPhi;
    
    TPZManVector<std::complex<double> > fEigenvalues;
    
    TPZFNMatrix<30,std::complex<double> > fCoeficients;
    
    TPZManVector<int64_t> fLocalIndices;
    
    void ExtendShapeFunctions(TPZMaterialData &data1d, TPZMaterialData &data2d);
    
    REAL fDensity = 1.;

    void AdjustAxes3D(const TPZFMatrix<REAL> &axes2D, TPZFMatrix<REAL> &axes3D, TPZFMatrix<REAL> &jac3D, TPZFMatrix<REAL> &jacinv3D, REAL detjac);
public:
    
    TPZSBFemVolume(TPZCompMesh &mesh, TPZGeoEl *gel, int64_t &index);
    
    virtual ~TPZSBFemVolume()
    {
        Reference()->ResetReference();
    }
    
    void ComputeKMatrices(TPZElementMatrix &E0, TPZElementMatrix &E1, TPZElementMatrix &E2, TPZElementMatrix &M0);
    
    void SetSkeleton(int64_t skeleton);
    
    int64_t SkeletonIndex()
    {
        return fSkeleton;
    }
    
    virtual void InitMaterialData(TPZMaterialData &data);

    void SetElementGroupIndex(int64_t index);
    
    virtual TPZCompEl *Clone(TPZCompMesh &mesh) const
    {
        // till I remember how this works
        DebugStop();
        return 0;
    }
    
    int64_t ElementGroupIndex() const
    {
        return fElementGroupIndex;
    }
    virtual TPZCompEl *ClonePatchEl(TPZCompMesh &mesh,
                                    std::map<int64_t,int64_t> & gl2lcConMap,
                                    std::map<int64_t,int64_t> & gl2lcElMap) const
    {
        // till I remember how this works
        DebugStop();
        return 0;
    }

    virtual int NConnects() const
    {
        if (fElementGroup == 0) {
            return 0;
        }
        return fElementGroup->NConnects();
    }
    
    virtual int64_t ConnectIndex(int i) const
    {
        if (fElementGroup == 0) {
            DebugStop();
        }
        return fElementGroup->ConnectIndex(i);
    }
    virtual int Dimension() const
    {
        TPZGeoEl *reference = Reference();
        return reference->Dimension();
    }
    
    REAL Density()
    {
        return fDensity;
    }
    
    void SetDensity(REAL density)
    {
#ifdef PZDEBUG
        if(density <= 0.) DebugStop();
#endif
        fDensity = density;
    }
    
    virtual void BuildCornerConnectList(std::set<int64_t> &connectindexes) const;
    
    virtual void SetConnectIndex(int inode, int64_t index)
    {
        DebugStop();
    }
    virtual int NSideConnects(int iside) const
    {
        DebugStop();
                                return 0;
    }
    
    virtual int SideConnectLocId(int icon,int is) const
    {
        DebugStop();
                                return 0;
    }

    virtual int NShapeF() const
    {
        int nc = NConnects();
        int nshape = 0;
        for (int ic=0; ic<nc; ic++) {
            TPZConnect &c = Connect(ic);
            nshape += c.NShape();
        }
        return nshape;
    }
    
    virtual int NConnectShapeF(int icon, int order) const
    {
        TPZConnect &c = Connect(icon);
        return c.NShape();
    }
    
    void InitializeIntegrationRule()
    {
        if (fIntRule) {
            DebugStop();
        }
        int nsides = Reference()->NSides();
        fIntRule = Reference()->CreateSideIntegrationRule(nsides-1, 1);
    }
    
    virtual void SetIntegrationRule(int order);
    
    virtual const TPZIntPoints &GetIntegrationRule() const
    {
        if(!fIntRule) DebugStop();

                                return *fIntRule;
    }
    
    virtual TPZIntPoints &GetIntegrationRule()
    {
        if(!fIntRule) InitializeIntegrationRule();
        return *fIntRule;
    }
    
    virtual void PRefine ( int order );

    virtual void SetPreferredOrder ( int order );
    


    void SetPhiEigVal(TPZFMatrix<std::complex<double> > &phi, TPZManVector<std::complex<double> > &eigval);
    
    TPZFMatrix<std::complex<double> > Phi()
    {
        return fPhi;
    }
    
    TPZManVector<std::complex<double> > Eigenvalues()
    {
        return fEigenvalues;
    }
    
    TPZFMatrix<std::complex<double> > Coeficients()
    {
        return fCoeficients;
    }
    
    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 = fCoeficients.Rows(),cols = fCoeficients.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) = fCoeficients(i,j).real();
            }
        }
        return coefreal;
    }
    
    virtual void LoadCoef(TPZFMatrix<std::complex<double> > &coef);
    
    virtual void ComputeSolution(TPZVec<REAL> &qsi,
                                 TPZSolVec &sol, TPZGradSolVec &dsol,TPZFMatrix<REAL> &axes);
    
    void EvaluateError(std::function<void(const TPZVec<REAL> &loc,TPZVec<STATE> &val,TPZFMatrix<STATE> &deriv)> fp,
                       TPZVec<REAL> &errors,bool store_error);
    
    virtual void Shape(TPZVec<REAL> &qsi,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphidxi);
    
    virtual void ComputeShape(TPZVec<REAL> &intpoint, TPZVec<REAL> &X, TPZFMatrix<REAL> &jacobian, TPZFMatrix<REAL> &axes,
                              REAL &detjac, TPZFMatrix<REAL> &jacinv, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi, TPZFMatrix<REAL> &dphidx);
    

    virtual void ComputeSolution(TPZVec<REAL> &qsi, TPZMaterialData &data)    {
        ComputeSolution(qsi, data.sol, data.dsol, data.axes);
    }

    virtual void Solution(TPZVec<REAL> &qsi,int var,TPZVec<STATE> &sol);
    

    virtual void Print(std::ostream &out = std::cout) const
    {
        out << "Printing " << __PRETTY_FUNCTION__ << std::endl;
        TPZCompEl::Print(out);
        out << "Group Element Index " << fElementGroupIndex << std::endl;
        out << "Skeleton Element Index " << fSkeleton << std::endl;
        out << "Local Indices " << fLocalIndices << std::endl;
        fCoeficients.Print("Coef =",out,EMathematicaInput);
        fPhi.Print("Phi = ",out,EMathematicaInput);
        if (fCoeficients.Rows())
        {
            TPZManVector<std::complex<double>,5> prod(fPhi.Rows(),0.);
            for (int i=0; i<fPhi.Rows(); i++) {
                for (int j=0; j<fPhi.Cols(); j++) {
                    prod[i] += fPhi.GetVal(i,j)*fCoeficients.GetVal(j,0);
                }
            }
            out << "Values at border " << prod << std::endl;
        }
        for (int i=0; i<NConnects(); i++) {
            Connect(i).Print(*Mesh(),out);
        }
    }
    
    void CreateGraphicalElement(TPZGraphMesh &, int);

};


TPZCompEl * CreateSBFemCompEl(TPZGeoEl *gel,TPZCompMesh &mesh,int64_t &index);



#endif /* TPZSBFemVolume_hpp */