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

void TPZStructMatrixCS::Assemble(TPZMatrix<STATE> & stiffness, TPZFMatrix<STATE> & rhs,TPZAutoPointer<TPZGuiInterface> guiInterface){

    ass_stiff.start();

    if (fEquationFilter.IsActive()) {
        int64_t neqcondense = fEquationFilter.NActiveEquations();
#ifdef PZDEBUG
        if (stiffness.Rows() != neqcondense) {
            DebugStop();
        }
#endif
        TPZFMatrix<STATE> rhsloc(neqcondense,rhs.Cols(),0.);
        if(this->fNumThreads){
            this->MultiThread_Assemble(stiffness,rhsloc,guiInterface);
        }
        else{
            this->Serial_Assemble(stiffness,rhsloc,guiInterface);
        }
        
        fEquationFilter.Scatter(rhsloc, rhs);
    }
    else
    {
        if(this->fNumThreads){
            this->MultiThread_Assemble(stiffness,rhs,guiInterface);
        }
        else{
            this->Serial_Assemble(stiffness,rhs,guiInterface);
        }
    }
    ass_stiff.stop();
}

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

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())
    }

void TPZArtDiff::ContributeFastestImplDiff_dim(TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
{
    REAL delta = Delta(deltaX, sol);
    REAL constant = /*-*/ weight * delta * timeStep;
    REAL buff;
	
    TPZVec<TPZVec<STATE> > TauDiv;
    TPZVec<TPZDiffMatrix<STATE> > dTauDiv;
	
    PrepareFastestDiff<dim>( jacinv, sol, dsol, phi, dphi, TauDiv, dTauDiv);
	
    int i, j, k, l;
    int nshape = dphi.Cols();
    int nstate = dim + 2;
    int neq = nstate * nshape;
	
    // ODotProduct speeded up
	
    for(l=0;l<nshape;l++)
		for(i=0;i<nstate;i++)
			for(k=0;k<dim;k++)
			{
				buff = dphi(k,l) * constant;
				ef(i+l*nstate,0) += buff * TauDiv[k][i];
				for(j=0;j<neq;j++)
					ek(i+l*nstate,j) -= buff * dTauDiv[k](i,j);
			}
}

void TPZErrorIndicator::Sort(TPZFMatrix &error, TPZFMatrix &perm) {
  int i,j,k;
  int imin = 0;
  int imax = error.Rows();
  perm.Resize(imax,error.Cols());
  for (i=0;i<imax;i++)
    for (j=0;j<error.Cols();j++) perm(i,j) = i;
  	for(i=imin; i<imax; i++) {
    	for(j=i+1; j<imax; j++) {
      	for (k=0;k<error.Cols();k++){
					if(error((int)perm(i,k)) < error((int)perm(j,k))) {
	  				int kp = (int) perm(i,k);
	  				perm(i,k) = perm(j,k);
	  				perm(j,k) = kp;
					}
      	}
    	}
  }
}

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
	 */
}

void TPZDohrMatrix<TVar,TSubStruct>::MultAddTBB(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,
                                                const TVar alpha,const TVar beta,const int opt) const
{

#ifdef USING_TBB
    
	if ((!opt && this->Cols() != x.Rows()) || this->Rows() != x.Rows())
		this->Error( "Operator* <matrixs with incompatible dimensions>" );
	if(x.Cols() != y.Cols() || x.Cols() != z.Cols() || x.Rows() != y.Rows() || x.Rows() != z.Rows()) {
		this->Error ("TPZFMatrix::MultiplyAdd incompatible dimensions\n");
	}
	this->PrepareZ(y,z,beta,opt);
	
    
    unsigned int nglob = fGlobal.size();
    TPZAutoPointer<TPZDohrAssembleList<TVar> > assemblelist = new TPZDohrAssembleList<TVar>(nglob,z,this->fAssembly);
    
    ParallelAssembleTaskMatrix<TVar,TSubStruct> multwork(x,alpha,fAssembly,assemblelist);
    typename std::list<TPZAutoPointer<TSubStruct> >::const_iterator iter;
    int isub=0;
    for (iter=fGlobal.begin(); iter!=fGlobal.end(); iter++,isub++) {
        TPZDohrThreadMultData<TSubStruct> data(isub,*iter);
        
        multwork.addWorkItem(data);
    }
    TPZVec<pthread_t> AllThreads(1);
    
    multwork.run_parallel_for(pzenviroment.fSubstructurePartitioner);
    
    PZ_PTHREAD_CREATE(&AllThreads[0], 0, TPZDohrAssembleList<TVar>::Assemble, 
                      assemblelist.operator->(), __FUNCTION__);
    
    void *result;
    PZ_PTHREAD_JOIN(AllThreads[0], &result, __FUNCTION__);
#endif    
    
}

/** This function creates a Data Explorer file that allow to visualization of the value of a matrix passed as parameter */
void VisualMatrixDX(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 << "# Graphical Visualization of Matrix." << endl;
	out << "# Positions as the indexes of the matrix, beginning by column." << endl;
	out << "# The number of elements in x direction correspond to the number of the columns of the matrix." << endl;
	out << "# The number of elements in y direction correspond to the number of the rows of the matrix." << endl;

	out  << "object 1 class gridpositions counts " << nelx+1 << " " << nely +1 << endl;
	out << "origin 0. 0." << endl;
	out << "delta 1. 0." << endl;
	out << "delta 0. 1." << endl;
	out << "attribute \"dep\" string \"positions\"" << endl;
	out << endl;


	out << "object 2 class gridconnections counts " << nelx+1 << " " << nely +1 << endl;

 	out << "attribute \"element type\" string \"quads\"" << endl;
	out << "attribute \"ref\" string \"positions\"" << endl;

	out.precision(5);
	out  << "object 3 class array type float rank 0 items " << neltotal << " data follows" << endl;
	for (i = 0; i < nelx; i++) {
		for(j=0; j< nely ; j++) out << matrix(i,j) << endl;
	}
	out << "attribute \"dep\" string \"connections\" " << endl;
	out << endl;

	out << "object 4 class field" << endl;
	out << "component \"data\" value 3" << endl;
	out << "component \"positions\" value 1" << endl;
	out << "component \"connections\" value 2" << endl;
	out << "attribute \"name\" string \"Matrix\"" << endl;

	out << endl;
	out << "end" << endl;

	out.close();

	cout << "Data Explorer file " << outfilename << " was created with success!\n";
}

void TPZArtDiff::ContributeExplDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,  TPZFMatrix<REAL> &dphix, TPZFMatrix<STATE> &ef, REAL weight, REAL timeStep, REAL deltaX)
{
    REAL delta = Delta(deltaX, sol);
    REAL constant = /*-*/ weight * delta * timeStep;
	
    TPZVec<TPZVec<STATE> > TauDiv;
	
    PrepareFastDiff(dim, jacinv, sol, dsol, dphix, TauDiv, NULL);
	
    int i, k, l;
    int nshape = dphix.Cols();
    int nstate = dim + 2;
	
    // ODotProduct speeded up
    for(l=0;l<nshape;l++)
		for(i=0;i<nstate;i++)
			for(k=0;k<dim;k++)
				ef(i+l*nstate,0) += dphix(k,l) * TauDiv[k][i] * constant;
}

void TPZArtDiff::Divergent(TPZFMatrix<STATE> &dsol,
						   TPZFMatrix<REAL> & dphi,
						   TPZVec<TPZDiffMatrix<STATE> > & Ai,
						   TPZVec<STATE> & Div,
						   TPZDiffMatrix<STATE> * dDiv)
{
	int nstate = Ai[0].Cols();
	int dim = nstate - 2;
	int nshape = dphi.Cols();
	Div.Resize(nstate);
	Div.Fill(0.);
	
	int i, j, k;
	
	// computing the divergent:
	// A.du/dx + B.du/dy + C.du/dz
	for(k=0;k<dim;k++)
		for(i=0;i<nstate; i++)
			for(j=0;j<nstate;j++)
			{
				Div[i]+=Ai[k](i,j)*dsol(k,j);
			}
	
	if(!dDiv)return;
	// computing an approximation to the divergent derivative:
	
	// dDiv/dUj ~= A.d2U/dUidx + B.d2U/dUidy + C.d2U/dUidz
	dDiv->Redim(nstate, nstate * nshape);
	int l;
	REAL buff;
	for(l=0;l<nshape;l++)
		for(j=0;j<nstate;j++)
			for(i=0;i<nstate; i++)
			{
				buff =0.;
				for(k=0;k<dim;k++)
				{
					buff+=Ai[k](i,j)*dphi(k,l);
				}
				dDiv->operator()(i,j+l*nstate)=buff;
			}
	
}

 /** @brief Computing operator for the parallel for. */
 void operator()(const blocked_range<size_t>& range) const
 {
     
     for(size_t i=range.begin(); i!=range.end(); ++i )
     {
         TPZDohrThreadMultData<TSubStruct> runner = mWorkItems[i];
         TPZFMatrix<TVar> xlocal;
         fAssembly->Extract(runner.fisub,*(fInput),xlocal);
         TPZAutoPointer<TPZDohrAssembleItem<TVar> > assembleItem = new TPZDohrAssembleItem<TVar>(runner.fisub,xlocal.Rows(),xlocal.Cols());
         runner.fSub->ContributeKULocal(fAlpha,xlocal,assembleItem->fAssembleData);
         fAssemblyStructure->AddItem(assembleItem);
     }
 }

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 TPZGeoTriangle::VecHdiv(TPZFMatrix<REAL> & coord, TPZFMatrix<REAL> & fNormalVec,TPZVec<int> &fVectorSide) {
    if(coord.Rows()!=3)
    {
        cout<< "Erro na dimens�o das linhas de coord"<< endl;
    }
    if(coord.Cols()!=3)
    {
        cout<< "Erro na dimens�o das colunas de coord"<< endl;
    }
    TPZVec<REAL> p1(3), p2(3), p3(3),result(3);
    for(int j=0; j<3; j++)
    {
        p1[j]=coord(j,0);
        p2[j]=coord(j,1);
        p3[j]=coord(j,2);
    }
    fNormalVec.Resize(14, 3);
    fVectorSide.Resize(14);
    int count=0;

    //primeira face
    for(int j=0; j<3; j++) //v0
    {
        fNormalVec(0,j) = coord(j,0)- coord(j,2);
    }
    fVectorSide[count]=0;
    count++;
    for(int j=0; j<3; j++) //v1
    {
        fNormalVec(1,j) = coord(j,1)- coord(j,2);
    }
    fVectorSide[count]=1;
    count++;
    //v2
    ComputeNormal(p1,p2,p3,result);
    fNormalVec(2,0) = -result[0];
    fNormalVec(2,1) = -result[1];
    fNormalVec(2,2) = -result[2];
    fVectorSide[count]=3;
    count++;
    //segunda face
    for(int j=0; j<3; j++) //v3
    {
        fNormalVec(3,j) = coord(j,1)- coord(j,0);
    }
    fVectorSide[count]=1;
    count++;
    for(int j=0; j<3; j++) //v4
    {
        fNormalVec(4,j) = coord(j,2)- coord(j,0);
    }
    fVectorSide[count]=2;
    count++;
    //v5
    ComputeNormal(p2,p3,p1,result);
    fNormalVec(5,0) = -result[0];
    fNormalVec(5,1) = -result[1];
    fNormalVec(5,2) = -result[2];
    fVectorSide[count]=4;
    count++;
    //terceira face
    for(int j=0; j<3; j++) //v6
    {
        fNormalVec(6,j) = coord(j,2)- coord(j,1);
    }
    fVectorSide[count]=2;
    count++;
    for(int j=0; j<3; j++) //v7
    {
        fNormalVec(7,j) = coord(j,0)- coord(j,1);
    }
    fVectorSide[count]=0;
    count++;
    //v8
    ComputeNormal(p3,p1,p2,result);
    fNormalVec(8,0) = -result[0];
    fNormalVec(8,1) = -result[1];
    fNormalVec(8,2) = -result[2];
    fVectorSide[count]=5;
    count++;
    // internos tangentes
    for(int j=0; j<3; j++) //v9
    {
        fNormalVec(9,j) = coord(j,1)- coord(j,0);
    }
    fVectorSide[count]=3;
    count++;
    for(int j=0; j<3; j++) //v10
    {
        fNormalVec(10,j) = coord(j,2)- coord(j,1);
    }
    fVectorSide[count]=4;
    count++;

    for(int j=0; j<3; j++) //v11
    {
        fNormalVec(11,j) = coord(j,0)- coord(j,2);
    }
    fVectorSide[count]=5;
    count++;
//.........这里部分代码省略.........

void TPZTransform::SetMatrix(TPZFMatrix<REAL> &mult, TPZFMatrix<REAL> &sum) {
	fRow = mult.Rows();
	fCol = mult.Cols();
	fMult = mult;
	fSum = sum;
}