TPZMatLaplacian.h

Go to the documentation of this file.

#ifndef MATLAPLACDH
#define MATLAPLACDH

#include <iostream>
#include "pzdiscgal.h"
#include "pzfmatrix.h"

class TPZMatLaplacian : public TPZDiscontinuousGalerkin {
                
                protected :
                
                STATE fXf;
                
                int fDim;
                
                STATE fK;
    
    TPZFNMatrix<9,STATE> fTensorK, fInvK;
                                
                REAL fSymmetry;
                
                enum EPenaltyType {ENoPenalty = 0, EFluxPenalty = 1, ESolutionPenalty, EBoth};
                
                EPenaltyType fPenaltyType;
                
    TPZAutoPointer<TPZFunction<STATE> > fPermeabilityFunction;

public:
                
                REAL fPenaltyConstant;
                
                void SetNoPenalty(){ this->fPenaltyType = ENoPenalty;}

                void SetFluxPenalty(){ this->fPenaltyType = EFluxPenalty; }

                void SetSolutionPenalty(){ this->fPenaltyType = ESolutionPenalty; }
                
                void SetBothPenalty(){ this->fPenaltyType = EBoth; }

                TPZMatLaplacian(int matid, int dim);
    
  TPZMatLaplacian(int matid)
    : TPZRegisterClassId(&TPZMatLaplacian::ClassId), 
    TPZDiscontinuousGalerkin(matid), fXf(0.), fDim(1), fK(1.), fTensorK(1,1,1.),
    fInvK(1,1,1.),
     fSymmetry(0.), fPenaltyType(ENoPenalty), fPenaltyConstant(0.)
  {

  }

                TPZMatLaplacian();

                TPZMatLaplacian(const TPZMatLaplacian &copy);

                virtual ~TPZMatLaplacian();

                TPZMatLaplacian &operator=(const TPZMatLaplacian &copy);

                void SetSymmetric(){
                                this->fSymmetry = -1.0;
                }

                void SetNonSymmetric() {
                                this->fSymmetry = +1.0;
                }

                bool IsSymetric(){
                                if (fSymmetry == -1.0) return true;
                                if (fSymmetry == +1.0) return false;
                                PZError << __PRETTY_FUNCTION__ << "\n Comparacao de numeros reais da errado\n";
                                return false;
                }

                virtual TPZMaterial * NewMaterial() override {
                                return new TPZMatLaplacian(*this);
                }

    virtual void FillDataRequirements(TPZMaterialData &data) override;
    
    virtual void FillDataRequirementsInterface(TPZMaterialData &data) override;

    virtual void FillBoundaryConditionDataRequirement(int type,TPZMaterialData &data) override
    {
        data.SetAllRequirements(false);
        if (type == 50) {
            data.fNeedsSol = true;
        }
        if (type == 3) {
            data.fNeedsNormal = true;
        }
    }


                virtual int Dimension() const override { return fDim;}

                virtual int NStateVariables() const override;

                void SetParameters(STATE diff, STATE f);

    void GetParameters(STATE &diff, STATE &f) const
    {
        diff = fK;
        f = fXf;
    }
    
    void SetPermeability(REAL perm) {
        fK = perm;
        fTensorK.Zero();
        fInvK.Zero();
        for (int i=0; i<fDim; i++) {
            fTensorK(i,i) = perm;
            fInvK(i,i) = 1./perm;
        }
    }
    
    void SetValPenaltyConstant(REAL penalty)
    {
        fPenaltyConstant = penalty;
    }
    //Set the permeability tensor and inverser tensor    
    void SetPermeabilityFunction(TPZAutoPointer<TPZFunction<STATE> > fp)
    {
        fPermeabilityFunction = fp;
    }
    

    //void GetParameters(STATE &diff, STATE &f);

  void SetDimension(int dim)
  {
      if(dim<0 || dim >3)
      {
          DebugStop();
      }
      fDim = dim;
      fTensorK.Redim(dim,dim);
      fInvK.Redim(dim,dim);
      for (int i=0; i<dim; i++) {
          fTensorK(i,i) = fK;
          fInvK(i,i) = STATE(1.)/fK;
      }
  }


                virtual void Print(std::ostream & out) override;

                virtual std::string Name() override { return "TPZMatLaplacian"; }

    virtual void Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef) override;

    virtual void Contribute(TPZMaterialData &data,REAL weight, TPZFMatrix<STATE> &ef) override;
    
    virtual void Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) override {
        TPZDiscontinuousGalerkin::Contribute(datavec,weight,ek,ef);
    }


                virtual void ContributeBCHDiv(TPZMaterialData &data,REAL weight,
                                                                                                                                  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc);

                virtual void ContributeHDiv(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef);

                virtual void ContributeBC(TPZMaterialData &data,REAL weight,
                                                                                                                  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) override;

    virtual void ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) override;

                virtual void ContributeBCInterface(TPZMaterialData &data, TPZMaterialData &dataleft, REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) override;

                virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight,
                                                                                                                                                 TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef) override;

                virtual int VariableIndex(const std::string &name) override;

                virtual int NSolutionVariables(int var) override;

                virtual int NFluxes() override { return 3;}

protected:

                virtual void Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout) override;

  public:

                virtual void Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout) override;
    virtual void Solution(TPZVec<TPZMaterialData> &datavec, int var, TPZVec<STATE> &Solout) override {
        DebugStop();
    }

                virtual void Flux(TPZVec<REAL> &x, TPZVec<STATE> &Sol, TPZFMatrix<STATE> &DSol, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux) override;

                virtual void Errors(TPZVec<REAL> &x,TPZVec<STATE> &u,
                                                                TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux,
                                                                TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) override;
    
                void ErrorsHdiv(TPZMaterialData &data,TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) override;


                virtual int NEvalErrors()  override {return 3;}

                virtual REAL ComputeSquareResidual(TPZVec<REAL>& X, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol) override;


                void InterfaceErrors(TPZVec<REAL> &/*x*/,
                                                                                                 TPZVec<STATE> &leftu, TPZFMatrix<STATE> &leftdudx, /* TPZFMatrix<REAL> &leftaxes,*/
                                                                                                 TPZVec<STATE> &rightu, TPZFMatrix<STATE> &rightdudx, /* TPZFMatrix<REAL> &rightaxes,*/
                                                                                                 TPZVec<STATE> &/*flux*/,
                                                                                                 TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values,
                                                                                                 TPZVec<STATE> normal, STATE elsize);

                virtual void BCInterfaceJump(TPZVec<REAL> &x, TPZSolVec &leftu,TPZBndCond &bc,TPZSolVec & jump) override;

                virtual int IsInterfaceConservative() override { return 1;}

    public:
virtual int ClassId() const override;


                virtual void Write(TPZStream &buf, int withclassid) const override;

                virtual void Read(TPZStream &buf, void *context) override;

};

#endif