pzgradientreconstruction.h

Go to the documentation of this file.

//
//  pzgradientreconstruction.h
//  PZ
//
//  Created by Agnaldo Farias on 4/10/13.
//
//

#ifndef __PZ__pzgradientreconstruction__
#define __PZ__pzgradientreconstruction__

#include <iostream>

#include "pzcompel.h"
#include "pzvec.h"
#include "pzcmesh.h"
#include "tpzautopointer.h"
#include "pzfunction.h"

class TPZGradientReconstruction
{
    
    /*
     *@brief Class to calculate the gradient and slope limiter
     */
    class TPZGradientData
    {
    public:
        /*
         *@param cel: computational element to reconstruct the gradient
         *@param useweight [in]: if true indicates the choice of using weights in the distorted mesh
         *@param paramk[in]: value equal 1 or 2
         */
        TPZGradientData();
        
         ~TPZGradientData();
        
        TPZGradientData(const TPZGradientData &cp);
        
        TPZGradientData &operator=(const TPZGradientData &copy);
        
        void SetCel(TPZCompEl *cel, bool useweight, REAL paramK);
        
        
//        struct TNeighcell
//        {
//            int fIsBoundary;
//            STATE fNeighSol;
//            TPZManVector<REAL,3> fCenterNeighbour;
//            TPZManVector<REAL,3> fCenterInterface;
//        };
//        
//        TPZStack<TNeighcell> fData;
        
        void Print(std::ostream &out) const;
        
        void GetCenterPointAndCellAveraged(TPZCompEl *cel, TPZManVector<REAL,3> &xcenter, STATE &solcel);
        
        void InitializeGradData(TPZCompEl *cel);
        
        /*
         *@brief Method to reconstruction gradient by Least Squares and apply the slope limiter
         */
        void ComputeGradient();
        
        void QRFactorization(TPZFMatrix<REAL> &matA,TPZFMatrix<REAL> &vecb);
        
        //para problema linear
        void ComputeSlopeLimiter();
        
        //para problema linear
        void ComputeSlopeLimiter2();
        
        //Para problema nao-linear
        void ComputeSlopeLimiter3();
        
        void ComputeWeights(REAL paramk);
        
        void NodeCloserCenterX(TPZManVector<REAL,3> &nodecelX);
        
        void InsertWeights(TPZFMatrix<REAL> &DeltaH, TPZFMatrix<REAL> &DifSol);
        
        
        /*
         *@brief Returns data from the computational element K (K cell)
         */
        void GetData(TPZManVector<REAL,3> &centerPoint, TPZManVector<STATE,3> &grad, STATE &cellAverage, STATE &slopeLimiter)
        {
            centerPoint = fCenterPointCellAndNeighbors[0];
            grad = fGradient;
            cellAverage = fSolCellAndNeighbors[0];
            slopeLimiter = fSlopeLimiter;
        }
        
        TPZAutoPointer<TPZFunction<STATE> > fForcingFunctionExact;
        
        void SetForcingFunctionExact(TPZAutoPointer<TPZFunction<STATE> > fp)
        {
            fForcingFunctionExact = fp;
        }
        
        bool fUseForcinfFuncion;
        
        void EnableForcinFucnction(){
            fUseForcinfFuncion=true;
        }
        void DisableForcinFucnction(){
            fUseForcinfFuncion=false;
        }
        int HasForcingFunctionExact() {
            if(fUseForcinfFuncion){
                return (fForcingFunctionExact != 0);
            }
            else{
                fForcingFunctionExact = 0;
                return 0;
            }
        }
        
        void UseGhostsNeighbors(TPZVec<REAL> LxLyLz, TPZVec<int> MatIdBC, TPZManVector<TPZVec<REAL> > coordmin, TPZManVector<TPZVec<REAL> > coordmax){
            
            fGhostNeighbor = true;
            
            fLxLyLz = LxLyLz;
            fMatIdBC = MatIdBC;
            fcoordminBC = coordmin;
            fcoordmaxBC = coordmax;
        }
        
        void CreateGhostsNeighbors(TPZCompEl *cel);
        
        
    protected:
        
        TPZStack<STATE> fSolCellAndNeighbors;
        
        TPZStack<TPZManVector<REAL,3> > fCenterPointCellAndNeighbors;
        
        TPZStack<TPZManVector<REAL,3> > fCenterPointInterface;
        
        TPZManVector<STATE,3> fGradient;
        
        STATE fSlopeLimiter;
        
        int fdim;
        
        TPZStack<TPZCompEl *> fCelAndNeighbors;
        
        TPZManVector<REAL,3> fWeightsGrad;
        
        bool fUseWeight;
        
        REAL fparamK;
        
        
        //parameter the ghosts neighbors
        bool fGhostNeighbor;
        
        TPZVec<REAL> fLxLyLz;
        TPZVec<int> fMatIdBC;
        TPZManVector< TPZVec<REAL> > fcoordminBC;
        TPZManVector< TPZVec<REAL> > fcoordmaxBC;
    };
    
public:
    
    /*
     *@param distmesh: Parameter indicates if the mesh is distorted
     *@param param: parameter used in the harmonic, equal 1 or 2, to calculate the weights
     */
    TPZGradientReconstruction(bool distmesh, REAL paramK);
    
    ~TPZGradientReconstruction();
    
    TPZGradientReconstruction(const TPZGradientReconstruction &cp);
    
    TPZGradientReconstruction &operator=(const TPZGradientReconstruction &copy);

    
    /*
     *@brief Method to replace the solution by finite element method from the reconstruction of the gradient.
     *Using the L2 projection
     *@param cmesh [in]: Computational mesh
     *@param datagradients [in]: Matrix of data about gradient in each element.
     *@param matidl2proj [in]: Id of the l2 projection material
     */
    void ProjectionL2GradientReconstructed(TPZCompMesh *cmesh, int matidl2proj);
    
    //change material id current to material id of the L2Projection
    void ChangeMaterialIdIntoCompElement(TPZCompEl *cel, int oldmatid, int newmatid);
    
    //assemble pos l2 projection
    void AssembleGlobalMatrix(TPZCompEl *el, TPZElementMatrix &ek, TPZElementMatrix &ef,TPZMatrix<STATE> & stiffmatrix, TPZFMatrix<STATE> &rhs);
    
    void GetDataDistortedMesh(bool &useweight, REAL &paramK){
        useweight = fDistortedMesh;
        paramK = fparam;
    }
    
    /*
     *@brief Method to insert the data of the boundaries (No flux) that there are ghosts neighbors.
     *@param LxLyLz [in]: Dimension of the domain: Lx X Ly X Lz
     *@param MatIdBC [in]: Indices of the Neumann type boundary conditions (No flux).
     *@param coordmin [in]: The minimum coordinate for each boundary (xmin,ymin,zmin).
     *@param coordmax [in]: The maximum coordinate for each boundary (xmax,ymax,zmax).
     */
    void SetDataGhostsNeighbors(TPZVec<REAL> LxLyLz, TPZVec<int> MatIdBC, TPZManVector< TPZVec<REAL> > coordmin, TPZManVector< TPZVec<REAL> > coordmax);
    
    TPZGradientData * fGradData;
    
protected:
    
    bool fDistortedMesh;
    
    REAL fparam;
};

#endif /* defined(__PZ__pzgradientreconstruction__) */