pzshapeextend.h

Go to the documentation of this file.

// $Id: pzshapeextend.h,v 1.1 2008-10-08 02:44:18 phil Exp $
#ifndef SHAPEEXTENDHPP
#define SHAPEEXTENDHPP

#include "pzfmatrix.h"
#include "pzstack.h"
#include "pztrnsform.h"
//#include "pzgeoelside.h"


namespace pzshape {
                
                template<class TFather>
                class SPr : public TFather {
                                
                public:
                                
                                typedef typename TFather::TMem FatMem;
                                typedef typename TFather::Top FatTop;
                                
                                typedef  pztopology::Pr<FatTop> Top;
                                
                                class TMem
                                {
                                                FatMem fUpper;
                                                FatMem fLower;
                                                FatMem fExtension;
                                                int fExtSideOrders[TFather::Top::NSides];
                                                TPZManVector<int,TFather::Top::NSides> fSign;
                                                TPZManVector<int,TFather::Top::NSides> fNShapeBefore;
                                                TPZManVector<int,TFather::Top::NSides> fNShapeAfter;
                                                TPZManVector<REAL> fExtShapes[TFather::Top::NSides];
                                                TPZManVector<int,TFather::Top::NSides> fNShape;
                                                int fMaxOrder;
                                                TPZManVector<REAL> fShape;
                                                
                                };
                                
                                static void Shape(TPZVec<REAL> &pt, TPZVec<int> &id, TPZVec<int> &order, TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi, TMem &memory);
                                static void SideShape(int side, TPZVec<REAL> &pt, TPZVec<int> &id, TPZVec<int> &order, TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi);
                                
                                
                                static void ShapeCorner(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi);
                                
                                static int NSideShapeF(int side, int order);
                                
                                //static int NShapeF(TPZVec<int> &order);
                                
                                static int NSideShapeF(int side, const TMem &mem);
                                
                                static int NShapeF(TPZVec<int> &order);
                                static void ShapeInternal(TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi,
                                                                                                                                  TPZFMatrix<REAL> &dphi);//,int quad_transformation_index
                                
                                
                };
                
                template<class TFather>
                inline void SPr<TFather>::Shape(TPZVec<REAL> &pt, TPZVec<int> &id, TPZVec<int> &order, TPZFMatrix<REAL> &shape,TPZFMatrix<REAL> &dshape, TMem &mem)
                {
                                if(!mem.IsInitialized())
                                {
                                                mem.Initialize(id,order);
                                }
                                else if(pt == mem.fpt)
                                {
                                                shape = mem.fphi;
                                                dshape = mem.fdphi;
                                                return;
                                }
                                TPZFNMatrix<100> shapeupper,dshapeupper;
                                TPZFNMatrix<100> shapelower,dshapelower;
                                TPZFNMatrix<100> shapeextend,dshapeextend;
                                // compute (or copy) the shape functions of the upper and lower sides
                                TFather::Shape(pt,shapeupper,dshapeupper,TMem::fUpper);
                                TFather::Shape(pt,shapelower,dshapelower,TMem::fLower);
                                // compute (or copy) the shape function values of the shapefunctions which will be extended
                                TFather::Shape(pt,shapeextend,dshapeextend,TMem::fExtension);
                                
                                REAL ExtendPoint = pt[Top::Dimension-1];
                                TPZFNMatrix<100> shape1d(mem.fMaxOrder,1),dshape1d(1,mem.fMaxOrder);
                                if(ExtendPoint == mem.fExtendPoint)
                                {
                                                shape1d = mem.fshape1d;
                                                dshape1d = mem.fdshape1d;
                                }
                                else
                                {
                                                // compute the onedimensional shapefunctions
                                }
                                // compute the shapefunctions of the lower and upper sides
                                int i,j;
                                int imax0 = shapelower.Rows();
                                REAL fac = (1-ExtendPoint)/2.;
                                for(i=0; i<imax0; i++)
                                {
                                                shape(i,0) = shapelower(i,0)*fac;
                                                for(j=0; j< Top::Dimension-1; j++)
                                                {
                                                                dshape(j,i) = dshapelower(j,i)*fac;
                                                }
                                                dshape(j,i) = -shapelower(i,0)/2;
                                }
                                int imax1 = imax0+shapeupper.Rows();
                                fac = (1+ExtendPoint)/2.;
                                for(i=imax0; i<imax1; i++)
                                {
                                                shape(i,0) = shapeupper(i-imax0,0)*fac;
                                                for(j=0; j< Top::Dimension-1; j++)
                                                {
                                                                dshape(j,i) = dshapelower(j,i-imax0)*fac;
                                                }
                                                dshape(j,i) = shapelower(i-imax0,0)/2;
                                }
                                // compute the contributions of the extended side shape functions
                                
//                              REAL bubble = (1.+ExtendPoint)*(1.-ExtendPoint);
                                // for each side
        int is;
                                for(is = 0; is < TFather::Top::NSides; is++)
                                {
                                                int ishape, ishapebef, ishapeaft, ishapeext;
                                                int next = mem.fExtShapes[is].NElements();
                                                for(ishape=0; ishape<mem.fNShape[is]; ishape++)
                                                {
                                                                ishapebef = ishape%mem.fNShapeBefore[is];
                                                                ishapeext = ishape/mem.fNShapeBefore[is];
                                                                ishapeaft = ishape/(mem.fNShapeBefore[is]*next);
                                                                int orgind = ishapebef+ishapeaft*mem.fNShapeBefore[is];
                                                                // this could be lower too
                                                                shape(ishape,0) = mem.fExtShapes[is][orgind]*shapelower[ishapeext];
                                                }
                                }
                                // project the point to the side
                                // compute the internal shape functions
                                // compute the extended derivatives of the internal shape functions
                                // compute the extension function for the side
                                // multiply both
                }
                
                template<class TFather>
    inline void SPr<TFather>::ShapeCorner(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
                {
                                TPZManVector<REAL,FatTop::NCornerNodes> phifather;
                                TPZFNMatrix<TFather::TOP::Dimension*TFather::TOP::NCornerNodes> dphifather(FatTop::Dimension,FatTop::NCornerNodes);
                                TFather::ShapeCorner(pt,phifather,dphifather);
                                phi.Resize(Top::NCornerNodes,1);
                                dphi.Resize(Top::Dimension,Top::NCornerNodes);
                                REAL faclower = -(pt[Top::Dimension-1]-1.)/2.;
                                REAL facupper = (pt[Top::Dimension-1]+1.)/2.;
                                int i;
                                for(i=0; i< FatTop::NCornerNodes; i++)
                                {
                                                phi(i,0) = faclower*phifather(i,0);
                                                phi(i+FatTop::NCornerNodes,0) = facupper*phifather(i,0);
                                                dphi(Top::Dimension-1,i) = -0.5*phifather(i,0);
                                                dphi(Top::Dimension-1,i+FatTop::NCornerNodes) = 0.5*phifather(i,0);
                                                int j;
                                                for(j=0; j<FatTop::Dimension; j++)
                                                {
                                                                dphi(j,i) = faclower*dphifather(j,i);
                                                                dphi(j,i+FatTop::NCornerNodes) = facupper*dphifather(j,i);
                                                }
                                }
                }
                
};

#endif