TPZMatrix<STATE> * TPZSpStructMatrix::Create(){
    int64_t neq = fEquationFilter.NActiveEquations();
	/*    if(fMesh->FatherMesh()) {
	 TPZSubCompMesh *smesh = (TPZSubCompMesh *) fMesh;
	 neq = smesh->NumInternalEquations();
	 }*/
    TPZFYsmpMatrix<STATE> * mat = new TPZFYsmpMatrix<STATE>(neq,neq);
	
    /**
     *Longhin implementation
	 */
    TPZStack<int64_t> elgraph;
    TPZVec<int64_t> elgraphindex;
    //    int nnodes = 0;
    fMesh->ComputeElGraph(elgraph,elgraphindex);
    /**Creates a element graph*/
    TPZMetis metis;
    metis.SetElementsNodes(elgraphindex.NElements() -1 ,fMesh->NIndependentConnects());
    metis.SetElementGraph(elgraph,elgraphindex);
	
    TPZManVector<int64_t> nodegraph;
    TPZManVector<int64_t> nodegraphindex;
    /**
     *converts an element graph structure into a node graph structure
     *those vectors have size ZERO !!!
     */
    metis.ConvertGraph(elgraph,elgraphindex,nodegraph,nodegraphindex);
    
#ifdef LOG4CXX2
    if(logger->isDebugEnabled()){
        std::stringstream sout;
        sout << "Node graph \n";
        metis.TPZRenumbering::Print(nodegraph, nodegraphindex);
        LOGPZ_DEBUG(logger, sout.str())
    }

void TPZArtDiff::ContributeImplDiff(int dim, TPZFMatrix<REAL> &jacinv, TPZVec<FADREAL> &sol, TPZVec<FADREAL> &dsol, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, REAL weight,  REAL timeStep, REAL deltaX)
{
    TPZVec<STATE> solReal(sol.NElements());
    int i;
    for(i = 0; i < sol.NElements(); i++)
		solReal[i] = sol[i].val();
	
    REAL delta = Delta(deltaX, solReal);
    REAL constant = /*-*/ delta * weight * timeStep;
	
    TPZVec<TPZVec<FADREAL> > TauDiv;
	
    PrepareFastDiff(dim, jacinv, sol, dsol, TauDiv);
	
    TPZVec<FADREAL> Diff;
    TPZVec<REAL> gradv(dim);
	
    int j, k, l;
    int nstate = dim + 2;
    int neq = sol[0].size();
    int nshape = neq/nstate;
	
    for(l=0;l<nshape;l++)
	{
		for(k=0;k<dim;k++)
			gradv[k] = dsol[k].dx/*fastAccessDx*/(/*k+*/l*nstate);// always retrieving this information from the first state variable...
		ODotOperator(gradv, TauDiv, Diff);
		for(i=0;i<nstate;i++)
		{
			ef(i+l*nstate,0) += constant * Diff[i].val();
			for(j=0;j<neq;j++)
				ek(i+l*nstate, j) -= constant * Diff[i].dx/*fastAccessDx*/(j);
		}
	}
}

void ReadSolution(ifstream &arq, TPZVec<REAL> &sol, TPZCompMesh *cmesh, int &nstate, TPZVec<int> &dimstate){

  int i,j,totaldim=0;
  for(i=0;i<nstate;i++) totaldim += dimstate[i];
  TPZVec<REAL> pt(3,0.);
  TPZVec<REAL> coord(3,0.);
  TPZVec<REAL> auxsol(totaldim,0.);

  int iter = 0;
  int nel = cmesh->NElements();
  int solsize = totaldim * nel;
  sol.Resize(solsize);
  sol.Fill(0.);

  for(i=0; i<nel; i++){
    TPZCompEl *el = cmesh->ElementVec()[i];
    if (!el) continue;
    el->Reference()->CenterPoint(el->Reference()->NSides()-1,pt);
    el->Reference()->X(pt,coord);
    EvaluateSolution(coord,auxsol);
    for (j=0;j<totaldim;j++){
      sol[iter] = auxsol[j];
      iter++;
    }
  }
}

void TPZEulerConsLawDEP::Flux(TPZVec<REAL> &x,TPZVec<REAL> &Sol,TPZFMatrix &DSol,
						   TPZFMatrix &axes,TPZVec<REAL> &flux) {
	TPZVec<REAL> Fx,Fy,Fz;
	Flux(Sol,Fx,Fy,Fz);
	int cap = Sol.NElements();
	int nstate = NStateVariables(),i;
	if(cap != nstate){
		PZError << "\nTPZEulerConsLawDEP::Flux data size error\n";
		flux.Resize(0);
		return;
	}
	if(nstate == 3){
		flux.Resize(3);
		for(i=0;i<3;i++) flux[i] = Fx[i];
		return;
	} else
		if(nstate == 4){
			flux.Resize(8);
			for(i=0;i<4;i++) flux[i] = Fx[i];
			for(i=4;i<8;i++) flux[i] = Fy[i];
			return;
		} else
			if(nstate == 5){
				flux.Resize(15);
				for(i=00;i<05;i++) flux[i] = Fx[i];
				for(i=05;i<10;i++) flux[i] = Fy[i];
				for(i=10;i<15;i++) flux[i] = Fz[i];
			}
}

void RSNAMeshPoints(TPZVec< TPZVec<REAL> > & pt, TPZVec< TPZVec< int64_t> > &elms)
{
   REAL x1 = 0.,
        x2 = 4.12791,
	y1 = 0.,
	y2 = 1.;

   pt.Resize(6);
   TPZVec<REAL> coord(3);

   coord[0] = x1;
   coord[1] = y1;
   coord[2] = 0.;
   pt[0] = coord;

   coord[0] = x2/2.;
   coord[1] = y1;
   coord[2] = 0.;
   pt[1] = coord;

   coord[0] = x2;
   coord[1] = y1;
   coord[2] = 0.;
   pt[2] = coord;

   coord[0] = x1;
   coord[1] = y2;
   coord[2] = 0.;
   pt[3] = coord;

   coord[0] = x2/2;
   coord[1] = y2;
   coord[2] = 0.;
   pt[4] = coord;

   coord[0] = x2;
   coord[1] = y2;
   coord[2] = 0.;
   pt[5] = coord;

// quadrilateral data

   TPZVec< int64_t > nodes(4);

   elms.Resize(2);

   nodes[0] = 0;
   nodes[1] = 1;
   nodes[2] = 4;
   nodes[3] = 3;
   elms[0] = nodes;

   nodes[0] = 1;
   nodes[1] = 2;
   nodes[2] = 5;
   nodes[3] = 4;
   elms[1] = nodes;

}

/**
 * @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 TPZArtDiff::SUPG(int dim, TPZVec<T> & sol, TPZVec<TPZDiffMatrix<T> > & Ai, TPZVec<TPZDiffMatrix<T> > & Tau){
	
#ifdef FASTESTDIFF
	
	TPZDiffMatrix<T> RTM, RMi,
	X, Xi,
	Temp, INVA2B2,
	LambdaSUPG;
	T us, c;
	
	TPZEulerConsLaw::uRes(sol, us);
	TPZEulerConsLaw::cSpeed(sol, 1.4, c);
	
	RMMatrix(sol, us, fGamma, RTM, RMi);
	
	EigenSystemSUPG(sol, us, c, fGamma, X, Xi, LambdaSUPG);
	
	
	RTM.     Multiply(X, Temp);
	Temp.   Multiply(LambdaSUPG, INVA2B2);
	INVA2B2.Multiply(Xi, Temp);
	Temp.   Multiply(RMi, INVA2B2);
	
	for(int i = 0; i < Ai.NElements();i++)
	{
		Ai[i].Multiply(INVA2B2, Tau[i]);
	}
#else
	
	TPZDiffMatrix<T> Rot, RotT,
	X, Xi,
	M, Mi,
	Temp, INVA2B2,
	LambdaSUPG;
	T us, c;
	
	TPZEulerConsLaw::uRes(sol, us);
	TPZEulerConsLaw::cSpeed(sol, 1.4, c);
	
	RotMatrix(sol, us, Rot, RotT);
	MMatrix(sol, us, fGamma, M, Mi);
	EigenSystemSUPG(sol, us, c, fGamma, X, Xi, LambdaSUPG);
	
	RotT.   Multiply(M, Temp);
	Temp.   Multiply(X, INVA2B2);
	INVA2B2.Multiply(LambdaSUPG, Temp);
	Temp.   Multiply(Xi, INVA2B2);
	INVA2B2.Multiply(Mi, Temp);
	Temp.   Multiply(Rot, INVA2B2);
	
	for(int i = 0; i < Ai.NElements();i++)
	{
		Ai[i].Multiply(INVA2B2, Tau[i]);
	}
	
#endif
}

void TPZErrorIndicator::SetError(TPZVec<REAL> &maxerror, TPZVec<REAL> &minerror, TPZVec<int> &erantype){
  if  (maxerror.NElements() != fState.NElements() || minerror.NElements() != fState.NElements() || erantype.NElements() != fState.NElements()){
    cout << "TPZErrorIndicator::SetError - Error : The error vector must have the dimension equal to the number of state variables\n";
    cout << "Try to set the number of state variables (by SetSolution) first...\n";
    exit (-1);
  }
  fMaxError = maxerror;
  fMinError = minerror;
  fErAnType = erantype;
}

void TPZGeoTriangle::VectorialProduct(TPZVec<REAL> &v1, TPZVec<REAL> &v2,TPZVec<REAL> &result) {
    if(v1.NElements()!=3||v2.NElements()!=3)
    {
        cout << " o tamanho do vetores eh diferente de 3"<< endl;
    }
    REAL x1=v1[0], y1=v1[1],z1=v1[2];
    REAL x2=v2[0], y2=v2[1],z2=v2[2];
    result.Resize(v1.NElements());
    result[0]=y1*z2-z1*y2;
    result[1]=z1*x2-x1*z2;
    result[2]=x1*y2-y1*x2;
}

void TPZCheckGeom::CreateMesh() {
	
	if(fMesh) delete fMesh;
	fMesh = new TPZGeoMesh();
	int noind[12];
	int no;
	for(no=0; no<12; no++) {
		noind[no] = fMesh->NodeVec().AllocateNewElement();
		TPZVec<REAL> coord(3);
		coord[0] = nodeco[no][0];
		coord[1] = nodeco[no][1];
		coord[2] = nodeco[no][2];
		fMesh->NodeVec()[noind[no]].Initialize(coord,*fMesh);
	}
	int matid = 1;
	TPZVec<int> nodeindex;
	int nel;
	for(nel=0; nel<7; nel++) {
		int in;
		nodeindex.Resize(numnos[nel]);
		for(in=0; in<numnos[nel]; in++) {
			nodeindex[in] = nodind[nel][in];
		}
		int index;  
		switch(nel) {
			case 0:
				fMesh->CreateGeoElement(ECube, nodeindex, matid, index);
				break;
			case 1:
				fMesh->CreateGeoElement(EPiramide, nodeindex,matid, index);
				break;
			case 2:
				fMesh->CreateGeoElement(ETetraedro, nodeindex,matid, index);
				break;
			case 3:
				fMesh->CreateGeoElement(EPrisma, nodeindex,matid, index);
				break;
			case 4:
				fMesh->CreateGeoElement(EOned, nodeindex,matid, index);
				break;
			case 5:
				fMesh->CreateGeoElement(EQuadrilateral, nodeindex,matid, index);
				break;
			case 6:
				fMesh->CreateGeoElement(ETriangle, nodeindex,matid, index);
				break;
			default:
				break;
		}
	}
	fMesh->BuildConnectivity();
}

void TPZArtDiff::ODotOperator(TPZVec<REAL> &dphi, TPZVec<TPZDiffMatrix<T> > &M, TPZDiffMatrix<T> &Result){
	
	int dim = M.NElements();
	int size = dphi.NElements();
	if(size<1 || size>3){
		PZError << "TPZArtDiff::PointOperator: error data size";
	}
	
	Result.Redim(M[0].fRows, M[0].fCols());
	
	int i;
	for (i=0;i<dim;i++)Result.Add(M[i], dphi[i]);
}

void TPZEuler::Solution(TPZVec<REAL> &Sol,TPZFMatrix<REAL> &DSol,TPZFMatrix<REAL> &axes,int var,
						TPZVec<REAL> &Solout){
	if(var == 1) {
		Solout.Resize(1);
		Solout[0] = gEul.Pressure(Sol);
	} else if(var ==2) {
		Solout.Resize(1);
		Solout[0] = Sol[0];
	} else if(var == 3) {
		Solout.Resize(2);
		Solout[0] = Sol[1];
		Solout[1] = Sol[2];
	} else TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
}

void TPZSkylMatrix<TVar>::InitializeElem(const TPZVec<int> &skyline, TPZManVector<REAL> &storage, TPZVec<REAL *> &point) {
	int dim = skyline.NElements();
	int nel = NumElements(skyline);
	storage.Resize(nel);
	storage.Fill(0.);
	int i;
	point.Resize(dim+1);
	if(dim) {
		point[0] = &storage[0];
		point[dim] = &storage[0]+nel;
	} else {
		point[0] = 0;
	}
	for(i=1; i<dim+1; i++) point[i] = point[i-1]+(i-1)-skyline[i-1]+1;
}

void TPZErrorIndicator::SetSolution(TPZVec<REAL> &sol, int nstate, TPZVec<int> &dimstate, TPZVec<int> &statetoanalyse){
  if (nstate <= 0){
    cout << "TPZErrorIndicator::SetSolution - Error : The number of state variables must be greater than zero\n";
    cout << "The adaptive modulus is going down...\n";
    exit (-1);
  }
  
  fNState = nstate;
  fState = statetoanalyse;
  fDim = dimstate;
  
  int i,dims  = 0;
  for (i=0;i<fNState;i++) dims += fDim[i];
  
  fNDataEl = dims;

  
  if (sol.NElements() % dims != 0){
    cout << "TPZErrorIndicator::SetSolution - Warning : The number of elements in solution vector is not multiple of state variables seted.\n";
  }
  
  fSolution = sol;
  fNElements = fSolution.NElements() / dims;
  if (fMesh && fMesh->ElementVec().NElements() != fNElements){
    cout << "TPZErrorIndicator::SetSolution - Warning : The number of elements in solution vector is not equal to the number of elements in the mesh.\n";
  }
  
}

// metodo para computar erros
void TPZMatPoissonD3::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout){
    
    Solout.Resize( this->NSolutionVariables(var));

    
    if(var == 1){ //function (state variable Q)
        for (int ip = 0; ip<fDim; ip++)
        {
            Solout[ip] = data.sol[0][ip];
        }
        
        return;
    }
    
    if(var == 2){ //function (state variable p)
        
        TPZVec<STATE> SolP;
        SolP = data.sol[0];
        
        Solout[0] = SolP[0];
        return;
    }
    
    
}

void TPZMixedDarcyFlow::FillBoundaryConditionDataRequirement(int type, TPZVec<TPZMaterialData> &datavec){
    int ndata = datavec.size();
    for (int idata=0; idata < ndata ; idata++) {
        datavec[idata].SetAllRequirements(false);
        datavec[idata].fNeedsSol = true;
    }
}

//Contribution of skeletal elements.
void TPZLagrangeMultiplier::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
{
    int nmesh = datavec.size();
    if (nmesh!=2) DebugStop();

    TPZFMatrix<REAL>  &phiQ = datavec[0].phi;
    TPZFMatrix<REAL> &phiP = datavec[1].phi;
    int phrq = phiQ.Rows();
    int phrp = phiP.Rows();
    
//------- Block of matrix B ------
    int iq, jp;
	for(iq = 0; iq<phrq; iq++) {
		for(jp=0; jp<phrp; jp++) {
            ek(iq, phrq+jp) += fMultiplier*weight*phiQ(iq,0)*phiP(jp,0);
		}
	}
    
    
//------- Block of matrix B^T ------
    int ip, jq;
	for(ip=0; ip<phrp; ip++) {
		for(jq=0; jq<phrq; jq++) {
			ek(ip + phrq,jq) += fMultiplier*weight*phiP(ip,0)*phiQ(jq,0);
		}
	}
}

//----
void TPZBndCond::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<REAL> &ek, TPZFMatrix<REAL> &ef) {

    //this->UpdataBCValues(datavec);

    int typetmp = fType;
    if (fType == 50) {
//		int i;
#ifdef DEBUG2
        {
            for(int iref=0; iref < datavec.size(); iref++) {
                std::stringstream sout;
                sout << __PRETTY_FUNCTION__ << datavec[iref].sol << " " << datavec[iref].x;
                LOGPZ_DEBUG(logger,sout.str().c_str());
            }
        }
#endif
        //for (i = 0; i <data.sol.NElements(); i++){
//			fBCVal2(i,0) = gBigNumber*data.sol[i];
//			fBCVal1(i,i) = gBigNumber;
//		}
//		fType = 2;
    }

    this->fMaterial->ContributeBC(datavec,weight,ek,ef,*this);
    fType = typetmp;
}

	void TPZPrism::CenterPoint(int side, TPZVec<REAL> &center) {
		center.Resize(Dimension);
		int i;
		for(i=0; i<Dimension; i++) {
			center[i] = MidSideNode[side][i];
		}
	}

void TPZMaterial::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<REAL> &ek, 
							   TPZFMatrix<REAL> &ef, TPZBndCond &bc){
	int nref=datavec.size();
	if (nref== 1) {
		this->ContributeBC(datavec[0], weight, ek,ef,bc);
	}
}