pzmaterialdata.cpp

Go to the documentation of this file.

#include "pzmaterialdata.h"
#include "TPZMaterial.h"
#include "pzcompel.h"
#include "pzelmat.h"
#include <sstream>
#include "pzerror.h"
#include "TPZInterfaceEl.h"
#include "pzdiscgal.h"
#include "pzaxestools.h"
#include "pzlog.h"

#ifdef LOG4CXX
static LoggerPtr logger(Logger::getLogger("pz.matrix.tpzfmatrix"));
static LoggerPtr loggerCheck(Logger::getLogger("pz.checkconsistency"));
#endif

TPZMaterialData::TPZMaterialData() : TPZRegisterClassId(&TPZMaterialData::ClassId), fShapeType(EEmpty), numberdualfunctions(0){
    this->SetAllRequirements(false);
    this->fNeedsNormalVecFad = false;
    this->intLocPtIndex = -1;
    this->intGlobPtIndex = -1;
    this->NintPts = -1;
    this->sol.Resize(1);
    this->dsol.Resize(1);
    this->gelElId = -1;
    this->HSize = 0.;
    this->detjac = 0.;
    this->numberdualfunctions = 0;
    this->gelElId = -1;
    this->fDirectionsOnMaster = 0;
#ifdef _AUTODIFF
    this->fNormalVecFad = 0;
#endif
}

TPZMaterialData::TPZMaterialData( const TPZMaterialData &cp ) : 
TPZRegisterClassId(&TPZMaterialData::ClassId),
fShapeType(cp.fShapeType) {
    this->operator =(cp);
}

TPZMaterialData & TPZMaterialData::operator= (const TPZMaterialData &cp ){
    this->fShapeType = cp.fShapeType;
    this->fNeedsSol = cp.fNeedsSol;
    this->fNeedsNeighborSol = cp.fNeedsNeighborSol;
    this->fNeedsHSize = cp.fNeedsHSize;
    this->fNeedsNeighborCenter = cp.fNeedsNeighborCenter;
    this->fNeedsNormal = cp.fNeedsNormal;
    this->fNeedsNormalVecFad = cp.fNeedsNormalVecFad;
    this->phi = cp.phi;
    this->dphi = cp.dphi;
    this->dphix = cp.dphix;
    this->divphi = cp.divphi;
    this->axes = cp.axes;
    this->jacobian = cp.jacobian;
    this->jacinv = cp.jacinv;
    this->normal = cp.normal;
    this->x = cp.x;
    this->p = cp.p;
    this->sol = cp.sol;
    this->dsol = cp.dsol;
    this->divsol = cp.divsol;
    this->HSize = cp.HSize;
    this->detjac = cp.detjac;
    this->intLocPtIndex = cp.intLocPtIndex;
    this->intGlobPtIndex = cp.intGlobPtIndex;
    this->NintPts = cp.NintPts;
    this->XCenter = cp.XCenter;
    this->fDirectionsOnMaster = cp.fDirectionsOnMaster;
    this->fVecShapeIndex = cp.fVecShapeIndex;
    this->fNormalVec = cp.fNormalVec;
#ifdef _AUTODIFF
    this->fNormalVecFad = cp.fNormalVecFad;
#endif
    this->numberdualfunctions = cp.numberdualfunctions;
    this->gelElId = cp.gelElId;
    
    return *this;
}

TPZMaterialData::~TPZMaterialData(){
    //NOTHING TO BE DONE!
}

void TPZMaterialData::SetAllRequirements(bool set){
    this->fNeedsSol = set;
    this->fNeedsNeighborSol = set;
    this->fNeedsHSize = set;
    this->fNeedsNeighborCenter = set;
    this->fNeedsNormal = set;
}

bool TPZMaterialData::Compare(TPZSavable *copy, bool override)
{
    TPZMaterialData *comp = dynamic_cast<TPZMaterialData *>(copy);
    if(!comp) return false;
    bool result = true;
    bool locres;
    locres = phi.Compare(&comp->phi,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"phi different")
    }
    result = result && locres;
    locres = dphi.Compare(&comp->dphi,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"dphi different")
    }
    result = result && locres;
    locres = dphix.Compare(&comp->dphix,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"dphix different")
    }
    result = result && locres;
    locres = divphi.Compare(&comp->divphi,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"dphix different")
    }
    result = result && locres;
    locres = axes.Compare(&comp->axes,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"axes different")
    }
    result = result && locres;
    locres = jacobian.Compare(&comp->jacobian,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"jacobian different")
    }
    result = result && locres;
    locres = jacinv.Compare(&comp->jacinv,override);
    if(!locres)
    {
        LOGPZ_DEBUG(loggerCheck,"jacinv different")
    }
    result = result && locres;
    return result;
}

// Compare the object for identity with the object pointed to, eventually copy the object
/*
 * compare both objects bitwise for identity. Put an entry in the log file if different
 * overwrite the calling object if the override flag is true
 */
bool TPZMaterialData::Compare(TPZSavable *copy, bool override) const
{
    DebugStop();
    return true;
}

void TPZMaterialData::Print(std::ostream &out) const
{
    phi.Print("phi",out);
    dphi.Print("dphi",out);
    dphix.Print("dphix",out);
    divphi.Print("div phi",out);
    axes.Print("axes",out);
    jacobian.Print("jacobian",out);
    jacinv.Print("jacinv",out);
    out << "normal " << normal << std::endl;
    out << "x " << x << std::endl;
    out << "p " << p << std::endl;
    out << "sol " << sol << std::endl;
    int nsol = dsol.size();
    for (int is=0; is<nsol; is++) {
        dsol[is].Print("dsol",out);
        
    }
    out << "divsol " << divsol << std::endl;
    out << "HSize " << HSize << std::endl;
    out << "detjac " << detjac << std::endl;
    out << "XCenter " << XCenter << std::endl;
    out << "fDirectionsOnMaster" << fDirectionsOnMaster << std::endl;
    out << "intLocPtIndex " << intLocPtIndex << std::endl;
    out << "intGlobPtIndex " << intGlobPtIndex << std::endl;
    out << "NintPts " << NintPts << std::endl;
    out << "gelElId " << gelElId << std::endl;
}

void TPZMaterialData::PrintMathematica(std::ostream &out) const
{
    phi.Print("phi = ",out,EMathematicaInput);
    dphi.Print("dphi = ",out,EMathematicaInput);
    dphix.Print("dphix = ",out,EMathematicaInput);
    axes.Print("axes = ",out,EMathematicaInput);
    jacobian.Print("jacobian = ",out,EMathematicaInput);
    jacinv.Print("jacinv = ",out,EMathematicaInput);
    out << "normal = {" << normal << "};" << std::endl;
    out << "x = {" << x << "};" << std::endl;
    out << "p = " << p << ";" << std::endl;
    out << "sol = { " << sol << "};" << std::endl;
    int nsol=dsol.size();
    for (int is=0; is<nsol; is++) {
        std::stringstream sout;
        sout << "dsol" << is << " = ";
        dsol[is].Print(sout.str().c_str(),out,EMathematicaInput);
    }
    out << "divsol = { " << divsol << "};" << std::endl;

    out << "HSize = " << HSize << ";" << std::endl;
    out << "detjac = " << detjac << ";" << std::endl;
    out << "XCenter = {" << XCenter << "};" << std::endl;
    out << "fDirectionsOnMaster" << fDirectionsOnMaster << std::endl;
    out << "intLocPtIndex = " << intLocPtIndex << ";" <<std::endl;
    out << "intGlobPtIndex = " << intGlobPtIndex << ";" <<std::endl;
    out << "NintPts = " << NintPts << ";" <<std::endl;
    out << "gelElId = " << gelElId << ";" <<std::endl;
}

void TPZMaterialData::Write(TPZStream &buf, int withclassid) const
{
    int shapetype = fShapeType;
    buf.Write(&shapetype);
    phi.Write(buf,0);
    dphi.Write(buf,0);
    dphix.Write(buf,0);
    divphi.Write(buf, 0);
    axes.Write(buf,0);
    jacobian.Write(buf,0);
    jacinv.Write(buf,0);
    buf.Write(normal.begin(),normal.size());
    buf.Write(x.begin(),x.size());
    buf.Write(&p,1);
    int nsol = sol.size();
    buf.Write(&nsol);
    for (int is=0; is<nsol; is++) {
        buf.Write(sol[is].begin(),sol[is].size());
    }
    
    nsol = dsol.size();
    buf.Write(&nsol);
    for (int is=0; is<nsol; is++) {
        dsol[is].Write(buf,0);
    }
    nsol = divsol.size();
    buf.Write(&nsol);
    for (int is=0; is<nsol; is++) {
        buf.Write(divsol[is].begin(),divsol[is].size());
    }

    buf.Write(&HSize,1);
    buf.Write(&detjac,1);
    buf.Write(XCenter.begin(),XCenter.size());
    buf.Write(&intLocPtIndex,1);
    buf.Write(&intGlobPtIndex,1);
    buf.Write(&NintPts,1);
    buf.Write(&gelElId,1);
}

void TPZMaterialData::Read(TPZStream &buf, void *context)
{
    int shapetype;
    buf.Read(&shapetype);
    fShapeType = (MShapeFunctionType) shapetype;
    phi.Read(buf,0);
    dphi.Read(buf,0);
    dphix.Read(buf,0);
    divphi.Read(buf, 0);
    axes.Read(buf,0);
    jacobian.Read(buf,0);
    jacinv.Read(buf,0);
    buf.Read(normal);
    buf.Read(x);
    buf.Read(&p,1);
    int nsol;
    buf.Read(&nsol,1);
    for (int is=0; is<nsol; is++) {
        buf.Read(sol[is]);
    }
    buf.Read(&nsol,1);
    for (int is = 0; is<nsol; is++) {
        dsol[is].Read(buf,0);
    }
    buf.Read(&nsol,1);
    for (int is=0; is<nsol; is++) {
        buf.Read(divsol[is]);
    }

    buf.Read(&HSize,1);
    buf.Read(&detjac,1);
    buf.Read(XCenter);
    buf.Read(&intLocPtIndex,1);
    buf.Read(&intGlobPtIndex,1);
    buf.Read(&NintPts,1);
    buf.Read(&gelElId,1);
}

int TPZMaterialData::ClassId() const{
    return Hash("TPZMaterialData");
}

void TPZMaterialData::ComputeFluxValues(TPZFMatrix<REAL> & fluxes){
    
    if (fShapeType != EVecandShape) {
        std::cout << __PRETTY_FUNCTION__ << "works only for Vec and Shape type." << std::endl;
        return;
    }
    
    int n_shape = fVecShapeIndex.size();
    fluxes.Redim(3, n_shape);
    
    for (int i = 0; i < n_shape; i++) {
        int i_v = fVecShapeIndex[i].first;
        int i_s = fVecShapeIndex[i].second;
        
        for (int j = 0; j < 3; j++) {
            fluxes(j,i) = phi(i_s,0) * fNormalVec(j,i_v);
        }
    }
    
}

void TPZMaterialData::ComputeFunctionDivergence()
{
    
    // Getting test and basis functions
    TPZFMatrix<REAL> dphi_s       = dphi; // Derivative For H1  test functions
    
    int n_phi_v = fVecShapeIndex.NElements();
    divphi.Redim(n_phi_v,1);
    REAL det_jac = detjac;

    int i_vec = 0;
    int i_phi_s = 0;
    
    for (int iq = 0; iq < n_phi_v; iq++)
    {
        i_vec = fVecShapeIndex[iq].first;
        i_phi_s = fVecShapeIndex[iq].second;
        
        int n_dir = dphi_s.Rows();
        for (int k = 0; k < n_dir; k++) {
            divphi(iq,0) +=  dphi(k,i_phi_s)*fDirectionsOnMaster(k,i_vec)/detjac;
        }
    }
        

}



template class TPZRestoreClass<TPZMaterialData>;