C++ TPZFMatrix::Redim方法代码示例

/// 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 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);
}

/// 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 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 StateAd(const TPZVec<REAL>&pt,TPZVec<STATE> &res, TPZFMatrix<STATE> & disp){
    disp.Redim(2,1);
    res[0]=0.;
    double x=pt[0];
    double y=pt[1];
    res[0]=10.*x*y*(1-x)*(1-y);
    
}

void TPZNonDarcyAnalysis::InitializeFirstSolution(TPZFMatrix<STATE> &Pressure, REAL &pini)
{

	Pressure.Redim(fSolution.Rows(),fSolution.Cols());
	TPZBlock<STATE> &block = this->Mesh()->Block();
	int nel = this->Mesh()->NElements();
	for (int iel = 0 ; iel < nel ; iel++){
		TPZCompEl *cel = this->Mesh()->ElementVec()[iel];
		if (!cel) continue;
		TPZInterpolatedElement *intel = dynamic_cast <TPZInterpolatedElement *> (cel);
		if (!intel) DebugStop();
		if (intel->Reference()->Dimension() != 2) continue; //soh elem 2d.
		int ncc = intel->NCornerConnects();
		//if (ncc != 4) DebugStop(); // I expect only quad element, althought it would work for triang
		for (int icc = 0 ; icc < ncc ; icc++){
			TPZConnect *con = &intel->Connect(icc);
			int conseq = con->SequenceNumber();
			int pos = block.Position(conseq);
			Pressure(pos,0) = pini;
		}
	}
}

void TPZSkylMatrix<TVar>::MultAdd(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,
							const TVar alpha,const TVar beta ,const int opt,const int stride ) const {
	// Computes z = beta * y + alpha * opt(this)*x
	//          z and x cannot overlap in memory
	if ((!opt && this->Cols()*stride != x.Rows()) || this->Rows()*stride != x.Rows())
		TPZMatrix<TVar>::Error(__PRETTY_FUNCTION__," <matrixs with incompatible dimensions>" );
	if(z.Rows() != x.Rows() || z.Cols() != x.Cols()) z.Redim(x.Rows(),x.Cols());
	if(x.Cols() != y.Cols() || x.Cols() != z.Cols() || x.Rows() != y.Rows() || x.Rows() != z.Rows()) {
		cout << "x.Cols = " << x.Cols() << " y.Cols()"<< y.Cols() << " z.Cols() " << z.Cols() << " x.Rows() " << x.Rows() << " y.Rows() "<< y.Rows() << " z.Rows() "<< z.Rows() << endl;
		TPZMatrix<TVar>::Error(__PRETTY_FUNCTION__," incompatible dimensions\n");
	}
	this->PrepareZ(y,z,beta,opt,stride);
	int rows = this->Rows();
	int xcols = x.Cols();
	int ic, r;
	for (ic = 0; ic < xcols; ic++) {
		for( r = 0 ; r < rows ; r++ ) {
			int offset = Size(r);
			TVar val = 0.;
			const TVar *p = &x.g((r-offset+1)*stride,ic);
			TVar *diag = fElem[r] + offset-1;
			TVar *diaglast = fElem[r];
			while( diag > diaglast ) {
				val += *diag-- * *p;
				p += stride;
			}
			if( diag == diaglast ) val += *diag * *p;
			z(r*stride,ic) += val*alpha;
			TVar *zp = &z((r-offset+1)*stride,ic);
			val = x.g(r*stride,ic);
			diag = fElem[r] + offset-1;
			while( diag > diaglast ) {
				*zp += alpha * *diag-- * val;
				zp += stride;
			}
		}
	}
}

void TPZMulticamadaOrthotropic::AnalyticTensor(TPZVec<REAL> &co, TPZFMatrix<REAL> &tensor) {
	
	REAL z = co[2];
	REAL x = co[0];
	REAL y = co[1];
	tensor.Redim(3,3);
	REAL height = Height();
	REAL zrel = z-(fZMax-fZMin)/2.;
	if(height <= 0.) return;
	REAL height3 = height*height*height;
	tensor(0,0) += (fNX[2]+fLinearX*fdNXdX[2]*x)/height;
	tensor(1,1) += (fNY[2]+fLinearY*fdNYdY[2]*y)/height;
	tensor(1,0) += (fNXY[2]+fLinearX*fdNXYdX[2]*x+fLinearY*fdNXYdY[2]*y)/height;
	tensor(0,1) = tensor(1,0);
	tensor(0,0) += REAL(12.)*(fMX[2]+(x*REAL(fLinearX)*fdMXdX[2]))*zrel/(height3);
	tensor(1,1) += 12.*(fMY[2]+y*REAL(fLinearY)*fdMYdY[2])*zrel/height3;
	tensor(0,1) += 12.*(fMXY[2]+REAL(fLinearX)*x*fdMXYdX[2]+REAL(fLinearY)*y*fdMXYdY[2])*zrel/height3;
	tensor(1,0) = tensor(0,1);
	tensor(0,2) += -6.*(fQX[2]+REAL(fLinearX)*x*fdQXdX[2])*(zrel*zrel-height*height/4.)/height3;
	tensor(2,0) = tensor(0,2);
	tensor(1,2) += -6.*(fQY[2]+REAL(fLinearY)*y*fdQYdY[2])*(zrel*zrel-height*height/4.)/height3;
	tensor(2,1)  = tensor(1,2);
	
}

void TPZStepSolver<TVar>::Solve(const TPZFMatrix<TVar> &F, TPZFMatrix<TVar> &result, TPZFMatrix<TVar> *residual){
	if(!this->Matrix()) {
		cout << "TPZMatrixSolver::Solve called without a matrix pointer\n";
		DebugStop();
	}
	
	TPZAutoPointer<TPZMatrix<TVar> > mat = this->Matrix();
    // update the matrix to which the preconditioner refers
    if(fPrecond)
    {
        
        fPrecond->UpdateFrom(this->Matrix());
    }
    
	if(result.Rows() != mat->Rows() || result.Cols() != F.Cols()) {
		result.Redim(mat->Rows(),F.Cols());
	}
	
	if(this->fScratch.Rows() != result.Rows() || this->fScratch.Cols() != result.Cols()) {
		this->fScratch.Redim(result.Rows(),result.Cols());
	}
	
	TVar tol = fTol;
	int numiterations = fNumIterations;
	switch(fSolver) {
		case TPZStepSolver::ENoSolver:
		default:
			cout << "TPZMatrixSolver::Solve called without initialized solver, Jacobi used\n";
			SetJacobi(1,0.,0);
		case TPZStepSolver::EJacobi:
			//    cout << "fScratch dimension " << fScratch.Rows() << ' ' << fScratch.Cols() << endl;
			mat->SolveJacobi(numiterations,F,result,residual,this->fScratch,tol,fFromCurrent);
			break;
		case TPZStepSolver::ESOR:
			mat->SolveSOR(numiterations,F,result,residual,this->fScratch,fOverRelax,tol,fFromCurrent);
			break;
		case TPZStepSolver::ESSOR:
			mat->SolveSSOR(numiterations,F,result,residual,this->fScratch,fOverRelax,tol,fFromCurrent);
			break;
		case TPZStepSolver::ECG:
			mat->SolveCG(numiterations,*fPrecond,F,result,residual,tol,fFromCurrent);
#ifdef LOG4CXX
		{
			std::stringstream sout;
			sout << "Number of equations " << mat->Rows() << std::endl;
			sout << "Number of CG iterations " << numiterations << " tol = " << tol;
			LOGPZ_DEBUG(logger,sout.str().c_str());
		}
#endif
			break;
		case TPZStepSolver::EGMRES: {
			TPZFMatrix<TVar> H(fNumVectors+1,fNumVectors+1,0.);
			mat->SolveGMRES(numiterations,*fPrecond,H,fNumVectors,F,result,residual,tol,fFromCurrent);
			if(numiterations == fNumIterations || tol >= fTol)
			{
				std::cout << "GMRes tolerance was not achieved : numiter " << numiterations <<
				" tol " << tol << endl;
			}
#ifdef LOG4CXX
			{
				std::stringstream sout;
				sout << "Number of GMRES iterations " << numiterations << " tol = " << tol;
				LOGPZ_DEBUG(logger,sout.str().c_str());
			}
#endif
		}
			break;
		case TPZStepSolver::EBICGSTAB: 
			mat->SolveBICGStab(numiterations, *fPrecond, F, result,residual,tol,fFromCurrent);
			
			if(numiterations == fNumIterations || tol >= fTol)
			{
				std::cout << "BiCGStab tolerance was not achieved : numiter " << numiterations <<
				" tol " << tol << endl;
			}
#ifdef LOG4CXX
		{
			std::stringstream sout;
			sout << "Number of BiCGStab iterations " << numiterations << " tol = " << tol;
			LOGPZ_DEBUG(logger,sout.str().c_str());
		}
#endif
			break;
		case TPZStepSolver::EDirect:
			result = F;
			mat->SolveDirect(result,fDecompose,fSingular);
			if(residual) residual->Redim(F.Rows(),F.Cols());
			break;
		case TPZStepSolver::EMultiply:
			mat->Multiply(F,result);
			if(residual) mat->Residual(result,F,*residual);
			
	}
}

//	void TPZPoint::GetSideHDivPermutation(int side, TPZVec<int64_t> &id, TPZVec<int> &permgather)
//	{
//		permgather[0] = 0;
//		return;
//	}
    void TPZPoint::ComputeDirections(int side, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, TPZVec<int> &sidevectors)
    {
        directions.Redim(3, 0);
        sidevectors.Resize(0);
    }

void Forcing(const TPZVec<REAL> &pt, TPZVec<REAL> &res,TPZFMatrix<STATE> &disp){
	disp.Redim(2,1);
//	double x = pt[0];
  //  double y = pt[1];
    res[0]= 0.;
}