C++ TPZFMatrix类代码示例

inline void TPZQuadraticQuad::GradX(const TPZFMatrix<REAL> &nodes,TPZVec<T> &loc, TPZFMatrix<T> &gradx){
    
    gradx.Resize(3,2);
    gradx.Zero();
    int nrow = nodes.Rows();
    int ncol = nodes.Cols();
#ifdef PZDEBUG
    if(nrow != 3 || ncol  != 8){
        std::cout << "Objects of incompatible lengths, gradient cannot be computed." << std::endl;
        std::cout << "nodes matrix must be 3x8." << std::endl;
        DebugStop();
    }
    
#endif
    TPZFNMatrix<3,T> phi(NNodes,1);
    TPZFNMatrix<6,T> dphi(2,NNodes);
    TShape(loc,phi,dphi);
    for(int i = 0; i < NNodes; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            gradx(j,0) += nodes.GetVal(j,i)*dphi(0,i);
            gradx(j,1) += nodes.GetVal(j,i)*dphi(1,i);
        }
    }
    
}

/// calcula a contribuicao para a matriz de rigidez
void TPBrSteamFlux::CalcStiff(TPZVec<REAL> &leftstate, TPZVec<REAL> &rightstate, TPZVec<REAL> &interfacestate, REAL delx, REAL area, REAL delt, 
			   TPZFMatrix<REAL> &ek, TPZFMatrix<REAL> &ef)
{
	const int totaleq = 2*TPBrCellConservation::NumCellEq+TPBrSteamFlux::NumFluxEq;
	TPZManVector<TFad<totaleq,REAL> , TPBrCellConservation::NumCellEq> leftcellfad(TPBrCellConservation::NumCellEq),rightcellfad(TPBrCellConservation::NumCellEq);
	TPZManVector<TFad<totaleq,REAL> , TPBrSteamFlux::NumFluxEq> interfacefad(TPBrSteamFlux::NumFluxEq);
	
	TPBrCellConservation::Initialize<totaleq>(leftstate,leftcellfad,0);
	TPBrSteamFlux::Initialize<totaleq>(interfacestate,interfacefad,TPBrCellConservation::NumCellEq);
	TPBrCellConservation::Initialize<totaleq>(rightstate,rightcellfad,TPBrCellConservation::NumCellEq+TPBrSteamFlux::NumFluxEq);
	
	TPZManVector<TFad<totaleq,REAL> , NumFluxEq> cellresidualfad(NumFluxEq); 
	
	FluxResidual(leftcellfad, rightcellfad, interfacefad, delx, area, delt, cellresidualfad );
	
	ek.Redim(NumFluxEq, totaleq);
	ef.Redim(NumFluxEq, 1);
	int i,j;
	for (i=0; i<NumFluxEq; i++) 
	{
		ef(i,0) = cellresidualfad[i].val();
		for (j=0; j<totaleq; j++) {
			ek(i,j) = cellresidualfad[i].d(j);
		}
	}
#ifdef LOG4CXX
	{
		std::stringstream sout;
		ek.Print("Flux stiffness",sout);
		LOGPZ_DEBUG(logger,sout.str())
	}
#endif
}

void TPZMaterialTest3D::Errors( TPZVec<REAL> &/*x*/,TPZVec<REAL> &u,TPZFMatrix<REAL> &dudx,
                                TPZFMatrix<REAL> &axes, TPZVec<REAL> &/*flux*/,TPZVec<REAL> & u_exact,
                                TPZFMatrix<REAL> &du_exact,TPZVec<REAL> &values)
{
    //TPZVec<REAL> sol(1),dsol(3);
    TPZManVector<REAL> sol(1),dsol(3);
    Solution(u,dudx,axes,1,sol);
    Solution(u,dudx,axes,2,dsol);
    if(dudx.Rows()<3)
    {
        REAL dx = du_exact(0,0)*axes(0,0)+du_exact(1,0)*axes(0,1);
        REAL dy = du_exact(0,0)*axes(1,0)+du_exact(1,0)*axes(1,1);
        REAL parc1 = fabs(dx-dudx(0,0));
        REAL parc2 = fabs(dy-dudx(1,0));
        //Norma L2
        values[1] = pow(fabs(u[0] - u_exact[0]),(REAL)2.0);
        //seminorma
        values[2] = pow(parc1,(REAL)2.)+pow(parc2,(REAL)2.);
        //Norma Energia
        values[0] = values[1]+values[2];
        return;
    }
    //values[1] : eror em norma L2
    values[1]  = pow(sol[0] - u_exact[0],(REAL)2.0);
    //values[2] : erro em semi norma H1
    values[2]  = pow(dsol[0] - du_exact(0,0),(REAL)2.0);
    if(dudx.Rows()>1) values[2] += pow(dsol[1] - du_exact(1,0),(REAL)2.0);
    if(dudx.Rows()>2) values[2] += pow(dsol[2] - du_exact(2,0),(REAL)2.0);
    //values[0] : erro em norma H1 <=> norma Energia
    values[0]  = values[1]+values[2];
}

void TPZMatRed<TVar,TSideMatrix>::SetF(const TPZFMatrix<TVar> & F)
{
	
	int64_t FCols=F.Cols(),c,r,r1;
	
	fF0.Redim(fDim0,FCols);
	fF1.Redim(fDim1,FCols);
	
	for(c=0; c<FCols; c++){
		r1=0;
		for(r=0; r<fDim0; r++){
			fF0.PutVal( r,c,F.GetVal(r,c) ) ;
		}
		//aqui r=fDim0
		for( ;r<fDim0+fDim1; r++){
			fF1.PutVal( r1++,c,F.GetVal(r,c) );
		}
	}
#ifdef LOG4CXX
    if (logger->isDebugEnabled()) {
        std::stringstream sout;
        F.Print("F Input",sout);
        fF0.Print("fF0 Initialized",sout);
        fF1.Print("fF1 Initialized",sout);
        LOGPZ_DEBUG(logger, sout.str())
    }

/**
 * @brief Computes a contribution to the stiffness matrix and load vector at one integration point to multiphysics simulation
 * @param data [in]
 * @param dataleft [in]
 * @param dataright [in]
 * @param weight [in]
 * @param ek [out] is the stiffness matrix
 * @param ef [out] is the load vector
 * @since June 5, 2012
 */
void TPZLagrangeMultiplier::ContributeInterface(TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, TPZVec<TPZMaterialData> &dataright, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
{
    TPZFMatrix<REAL> *phiLPtr = 0, *phiRPtr = 0;
    for (int i=0; i<dataleft.size(); i++) {
        if (dataleft[i].phi.Rows() != 0) {
            phiLPtr = &dataleft[i].phi;
            break;
        }
    }
    for (int i=0; i<dataright.size(); i++) {
        if (dataright[i].phi.Rows() != 0) {
            phiRPtr = &dataright[i].phi;
            break;
        }
    }
    
    if(!phiLPtr || !phiRPtr)
    {
        DebugStop();
    }
    TPZFMatrix<REAL> &phiL = *phiLPtr;
    TPZFMatrix<REAL> &phiR = *phiRPtr;
    
    
    int nrowl = phiL.Rows();
    int nrowr = phiR.Rows();
    static int count  = 0;

    if((nrowl+nrowr)*fNStateVariables != ek.Rows() && count < 20)
    {
        std::cout<<"ek.Rows() "<< ek.Rows()<<
        " nrowl " << nrowl <<
        " nrowr " << nrowr << " may give wrong result " << std::endl;
        count++;
    }

    int secondblock = ek.Rows()-phiR.Rows()*fNStateVariables;
    int il,jl,ir,jr;
    
    // 3) phi_I_left, phi_J_right
    for(il=0; il<nrowl; il++) {
        for(jr=0; jr<nrowr; jr++) {
            for (int ist=0; ist<fNStateVariables; ist++) {
                ek(fNStateVariables*il+ist,fNStateVariables*jr+ist+secondblock) += weight * fMultiplier * (phiL(il) * phiR(jr));
            }
        }
    }
    
    //	// 4) phi_I_right, phi_J_left
    for(ir=0; ir<nrowr; ir++) {
        for(jl=0; jl<nrowl; jl++) {
            for (int ist=0; ist<fNStateVariables; ist++) {
                ek(ir*fNStateVariables+ist+secondblock,jl*fNStateVariables+ist) += weight * fMultiplier * (phiR(ir) * phiL(jl));
            }
        }
    }

}

void TPZEulerConsLawDEP::ContributeBC(TPZMaterialData &data,REAL weight,
                                   TPZFMatrix &ek,TPZFMatrix &ef,TPZBndCond &bc) {
	TPZFMatrix dphi = data.dphix;
	TPZFMatrix dphiL = data.dphixl;
	TPZFMatrix dphiR = data.dphixr;
	TPZFMatrix phi = data.phi;
	TPZFMatrix phiL = data.phil;
	TPZFMatrix phiR = data.phir;
	TPZManVector<REAL,3> normal = data.normal;
	TPZManVector<REAL,3> x = data.x;
	// int POrder=data.p;
	// int LeftPOrder=data.leftp;
	// int RightPOrder=data.rightp;
	TPZVec<REAL> sol=data.sol;
	TPZVec<REAL> solL=data.soll;
	TPZVec<REAL> solR=data.solr;
	TPZFMatrix dsol=data.dsol;
	TPZFMatrix dsolL=data.dsoll;
	TPZFMatrix dsolR=data.dsolr;
	// REAL faceSize=data.HSize;
	
	int phr = phi.Rows();
	short in,jn,i,j;
	int nstate = NStateVariables();
	REAL v2[5];//m�imo nstate
	for(i=0;i<nstate;i++) v2[i] = bc.Val2()(i,0);
	
	switch (bc.Type()) {
		case 0 :// Dirichlet condition
			for(in = 0 ; in < phr; in++) {
				for(i = 0 ; i < nstate; i++)
					ef(in*nstate+i,0) += gBigNumber * weight * v2[i] * phi(in,0);
				for (jn = 0 ; jn < phr; jn++) {
					for(i = 0 ; i < nstate; i++)
						ek(in*nstate+i,jn*nstate+i) += gBigNumber * weight * phi(in,0) * phi(jn,0);
				}
			}
			break;
		case 1 :// Neumann condition
			for(in = 0 ; in < phi.Rows(); in++) {
				for(i = 0 ; i < nstate; i++)
					ef(in*nstate+i,0) += v2[i] * phi(in,0) * weight;
			}
			break;
		case 2 :// condi�o mista
			for(in = 0 ; in < phi.Rows(); in++) {
				for(i = 0 ; i < nstate; i++)
					ef(in*nstate+i, 0) += weight * v2[i] * phi(in, 0);
				for (jn = 0 ; jn < phi.Rows(); jn++) {
					for(i = 0 ; i < nstate; i++) for(j = 0 ; j < nstate; j++)
						ek(in*nstate+i,jn*nstate+j) += weight * bc.Val1()(i,j) * phi(in,0) * phi(jn,0);
				}
			}
	}
}

void TPZBiharmonic::Solution(TPZVec<REAL> &Sol,TPZFMatrix<REAL> &DSol,TPZFMatrix<REAL> &/*axes*/,
							 int var,TPZVec<REAL> &Solout){
	if(var == 0 || var == 1) Solout[0] = Sol[0];//function
	if(var == 2) {
		Solout.Resize(DSol.Rows());
		int id;
		for(id=0 ; id  < DSol.Rows(); id++) {
			Solout[id] = DSol(id,0);//derivate
		}
	}
}

void TElemento0d::CalcStiff(TMalha &malha, TPZFMatrix& stiff, TPZFMatrix& rhs)
{
  stiff.Redim(1,1);
  rhs.Redim(1,1);
  std::vector<double> phi(1,1.);
  TPZFMatrix dphi(0,1);
  std::vector<double> point(0);
  Shape(point,phi,dphi);
  double weight = 1.;
  TMaterial *mat = malha.getMaterial(this->fMaterialId);
  mat->Contribute(point,weight,phi,dphi,stiff,rhs);
}

REAL TPZIncNavierStokesKEps::Dot(TPZFMatrix<REAL> &A, TPZFMatrix<REAL> &B){
	REAL sum = 0.;
	int i, j, rows, cols;
	rows = A.Rows();
	cols = A.Cols();
	for(i = 0; i < rows; i++){
		for(j = 0; j < cols; j++){
			sum += A(i,j) * B(i,j);
		}
	}
	return sum;
}

/** This function creates a Visualization Tool Kit (VTK) file that allow to visualization of the value of a matrix passed as parameter */
void VisualMatrixVTK(TPZFMatrix<REAL> & matrix, const std::string &outfilename)
{
	const int nelx = matrix.Cols();
	const int nely = matrix.Rows();
	const int neltotal = nelx * nely;
	int i,j;
	ofstream out(outfilename.c_str());
	out << "# vtk DataFile Version 3.0\n";
	out << "Generated by PZ\n";
	out << "ASCII\n";
	out << "DATASET RECTILINEAR_GRID\n";
	out << "DIMENSIONS " << (nelx+1) << " " <<  (nely+1) << " 1\n";
	out << "X_COORDINATES " << nelx+1 << " float\n";
	for (i=0; i<=nelx; i++) {
		out << i << " ";
	}
	out << std::endl;
	out << "Y_COORDINATES " << nely+1 << " float\n";
	for (j=0; j<=nely; j++) {
		out << j << " ";
	}
	out << std::endl;
	out << "Z_COORDINATES " << 1 << " float\n0.\n";
	out << "CELL_DATA " << nelx*nely << std::endl;
	out << "SCALARS mat_value float 1\n";
	out << "LOOKUP_TABLE default\n";
	const REAL *elem = &matrix(0,0);
	for (i=0; i<neltotal; i++) {
		out << *(elem+i) << std::endl;
	}
	/*
	 # vtk DataFile Version 3.0
	 Cube example
	 ASCII
	 DATASET RECTILINEAR_GRID
	 DIMENSIONS 3 3 1
	 X_COORDINATES 3 float
	 0. 1. 2.
	 Y_COORDINATES 3 float
	 0. 1. 2.
	 Z_COORDINATES 1 float
	 0. 
	 
	 CELL_DATA 4
	 SCALARS toto float 1
	 LOOKUP_TABLE default
	 1
	 2
	 3
	 4
	 */
}

/// calcula a contribuicao para a matriz de rigidez das equacoes de entrada
void TPBrSteamFlux::InletCalcStiff(TPZVec<REAL> &rightstate, TPZVec<REAL> &interfacestate, REAL delx, REAL area, REAL delt, 
					TPZFMatrix<REAL> &ek, TPZFMatrix<REAL> &ef)
{
	const int totaleq = NumInletVars+NumFluxEq+TPBrCellConservation::NumCellEq;
//	TPZManVector<TFad<totaleq,REAL> , NumInletVars> inletfad(NumInletVars);
	TPZManVector<TFad<totaleq,REAL> , TPBrSteamFlux::NumFluxEq+TPBrSteamFlux::NumInletVars> interfacefad(TPBrSteamFlux::NumFluxEq+TPBrSteamFlux::NumInletVars);
	TPZManVector<TFad<totaleq,REAL> , TPBrCellConservation::NumCellEq> rightcellfad(TPBrCellConservation::NumCellEq);
	
	//TPBrSteamFlux::InitializeInlet<totaleq>(inletstate,inletfad,0);
	TPBrSteamFlux::InitializeInlet<totaleq>(interfacestate,interfacefad,0);
	TPBrCellConservation::Initialize<totaleq>(rightstate,rightcellfad,NumInletVars+TPBrSteamFlux::NumFluxEq);
#ifdef LOG4CXX
	{
		std::stringstream sout;
        sout << "before calling inlet flux residual\n";
//		sout << "inletfad " << inletfad << std::endl;
		sout << "interfacefad " << interfacefad << std::endl;
        sout << "rightcellfad " << rightcellfad << std::endl;
		LOGPZ_DEBUG(logger,sout.str())
	}
#endif
	
	TPZManVector<TFad<totaleq,REAL> , NumFluxEq+NumInletVars> cellresidualfad(NumFluxEq+NumInletVars); 
	
	InletFluxResidual(rightcellfad, interfacefad, delx, area, delt, cellresidualfad );
	
#ifdef LOG4CXX
    {
        std::stringstream sout;
        sout << "cellresidual " << cellresidualfad;
        LOGPZ_DEBUG(logger, sout.str())
    }
#endif
	ek.Redim(NumFluxEq+NumInletVars, totaleq);
	ef.Redim(NumFluxEq+NumInletVars, 1);
	int i,j;
	for (i=0; i<NumFluxEq+NumInletVars; i++) 
	{
		ef(i,0) = cellresidualfad[i].val();
		for (j=0; j<totaleq; j++) {
			ek(i,j) = cellresidualfad[i].d(j);
		}
	}
#ifdef LOG4CXX
	{
		std::stringstream sout;
		ek.Print("Inlet stiffness",sout);
		LOGPZ_DEBUG(logger,sout.str())
	}
#endif
}

void TPZQuadraticQuad::X(const TPZFMatrix<REAL> &nodes,TPZVec<T> &loc,TPZVec<T> &x){
    
    TPZFNMatrix<4,T> phi(NNodes,1);
    TPZFNMatrix<8,T> dphi(2,NNodes);
    TShape(loc,phi,dphi);
    int space = nodes.Rows();
    
    for(int i = 0; i < space; i++) {
        x[i] = 0.0;
        for(int j = 0; j < NNodes; j++) {
            x[i] += phi(j,0)*nodes.GetVal(i,j);
        }
    }
    
}

void TPZFrontNonSym::ExtractFrontMatrix(TPZFMatrix<REAL> &front)
{
	// Extend the front with the non initialized rigid body modes
	int ieq;
	int maxeq = fLocal.NElements();
	for (ieq = fNextRigidBodyMode; ieq< maxeq; ieq++) {
		int ilocal = Local(ieq);
		Element(ilocal, ilocal) = 1.;
	}
	
	int mineq = 0;
	for(mineq=0; mineq<maxeq; mineq++) if(fLocal[mineq] != -1) break;
	int numeq = maxeq-mineq;
	front.Redim(numeq,numeq);
	int jeq;
	for(ieq=mineq;ieq<maxeq;ieq++) {
		if(fLocal[ieq] == -1) continue;
		int il = ieq-mineq;
		for(jeq=0;jeq<maxeq;jeq++) {
			if(fLocal[jeq] == -1) continue;
			int jl = jeq-mineq;
			front(il,jl) = this->Element(fLocal[ieq],fLocal[jeq]);
		}
	}
	
}

void TPZDarcyAnalysis::PrintLS(TPZAnalysis *an)
{
    TPZAutoPointer< TPZMatrix<REAL> > KGlobal;
    TPZFMatrix<STATE> FGlobal;
    KGlobal =   an->Solver().Matrix();
    FGlobal =   an->Rhs();
    
#ifdef PZDEBUG
    #ifdef LOG4CXX
        if(logger->isDebugEnabled())
        {
            std::stringstream sout;
            KGlobal->Print("KGlobal = ", sout,EMathematicaInput);
            FGlobal.Print("FGlobal = ", sout,EMathematicaInput);
            LOGPZ_DEBUG(logger,sout.str())
        }

void TPZStructMatrixCS::Assemble(TPZFMatrix<STATE> & rhs,TPZAutoPointer<TPZGuiInterface> guiInterface){
    ass_rhs.start();
    if(fEquationFilter.IsActive())
    {
        int64_t neqcondense = fEquationFilter.NActiveEquations();
        int64_t neqexpand = fEquationFilter.NEqExpand();
        if(rhs.Rows() != neqexpand || Norm(rhs) != 0.)
        {
            DebugStop();
        }
        TPZFMatrix<STATE> rhsloc(neqcondense,1,0.);
        if(this->fNumThreads)
        {
            this->MultiThread_Assemble(rhsloc,guiInterface);
        }
        else
        {
            this->Serial_Assemble(rhsloc,guiInterface);
        }
        fEquationFilter.Scatter(rhsloc,rhs);
    }
    else
    {
        if(this->fNumThreads){
            this->MultiThread_Assemble(rhs,guiInterface);
        }
        else{
            this->Serial_Assemble(rhs,guiInterface);
        }
    }
    ass_rhs.stop();
}