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

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

/**
 * @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 TPZMaterialTest::ContributeBC(TPZMaterialData &data,
                                   REAL weight,
                                   TPZFMatrix<STATE> &ek,
                                   TPZFMatrix<STATE> &ef,
                                   TPZBndCond &bc) {
	TPZFMatrix<REAL> &phi = data.phi;
	
	if(bc.Material() != this) {
		PZError << "TPZMat1dLin.apply_bc warning : this material didn't create the boundary condition!\n";
	}
	
	if(bc.Type() < 0 && bc.Type() > 2){
		PZError << "TPZMat1dLin.aplybc, unknown boundary condition type :"  <<
		bc.Type() << " boundary condition ignored\n";
	}
	
	
	int numdof = NStateVariables();
	int numnod = ek.Rows()/numdof;
	int r = numdof;
	
	int idf,jdf,in,jn;
	switch(bc.Type()) {
			
		case 0:
			for(in=0 ; in<numnod ; ++in){
				for(idf = 0;idf<r;idf++) {
					(ef)(in*r+idf,0) += gBigNumber*phi(in,0)*bc.Val2()(idf,0)*weight;
				}
				for(jn=0 ; jn<numnod ; ++jn) {
					for(idf = 0;idf<r;idf++) {
						ek(in*r+idf,jn*r+idf) += gBigNumber*phi(in,0)*phi(jn,0)*weight;
					}
				}
			}
			break;
			
		case 1:
			for(in=0 ; in<numnod ; ++in){
				for(idf = 0;idf<r;idf++) {
					(ef)(in*r+idf,0) += phi(in,0)*bc.Val2()(idf,0)*weight;
				}
			}
			break;
			
		case 2:
			for(in=0 ; in<numnod ; ++in){
				for(idf = 0;idf<r;idf++) {
					(ef)(in*r+idf,0) += phi(in,0)*bc.Val2()(idf,0)*weight;
				}
				for(jn=0 ; jn<numnod ; ++jn) {
					for(idf = 0;idf<r;idf++) {
						for(jdf = 0;jdf<r;jdf++) {
							ek(in*r+idf,jn*r+jdf) += bc.Val1()(idf,jdf)*phi(in,0)*phi(jn,0)*weight;
						}
					}
				}
			}//fim switch
	}
}

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

void TPZMaterialTest3D::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
    else if(var == 2)
    {
        Solout[0] = DSol(0,0);//derivate
        Solout[1] = DSol(1,0);//derivate
        if(DSol.Rows()>2) Solout[2] = DSol(2,0);//derivate
    }
    else TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
}

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 TPZGeoTriangle::X(TPZFMatrix<REAL> & coord, TPZVec<REAL> & loc,TPZVec<REAL> &result) {

    REAL spacephi[3],spacedphi[6];
    TPZFMatrix<REAL> phi(3,1,spacephi,3);
    TPZFMatrix<REAL> dphi(2,3,spacedphi,6);
    Shape(loc,phi,dphi);
    int space = coord.Rows();

    for(int i = 0; i < space; i++) {
        result[i] = 0.0;
        for(int j = 0; j < 3; j++) result[i] += phi(j,0)*coord(i,j);
    }
}

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    
    
}

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

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