C++ TPZVec::size方法代码示例

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

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

//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 RefinamentoUniforme(TPZAutoPointer<TPZGeoMesh> gmesh, int nref,TPZVec<int> dims)
{
    
    int ir, iel, k;
    int nel=0, dim=0;
    int ndims = dims.size();
	for(ir = 0; ir < nref; ir++ )
    {
		TPZVec<TPZGeoEl *> filhos;
        nel = gmesh->NElements();
        
		for (iel = 0; iel < nel; iel++ )
        {
			TPZGeoEl * gel = gmesh->ElementVec()[iel];
            if(!gel) DebugStop();
            
            dim = gel->Dimension();
            
            for(k = 0; k<ndims; k++)
            {
                if(dim == dims[k])
                {
                    gel->Divide (filhos);
                    break;
                }
            }
		}
	}
    
}

void TPZMaterial::Solution(TPZVec<TPZMaterialData> &datavec, int var, TPZVec<REAL> &Solout){
	if (datavec.size()==1) {
		this->Solution(datavec[0], var, Solout);
	}
	else {
		this->Solution(datavec, var, Solout);
	}
}

void TPZTracerFlow::FillBoundaryConditionDataRequirement(int type,TPZVec<TPZMaterialData > &datavec){
    int nref = datavec.size();
	for(int i = 0; i<nref; i++)
	{
        datavec[i].fNeedsSol = true;
		datavec[i].fNeedsNormal = false;
	}
}

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

int TPZMaterial::IntegrationRuleOrder(TPZVec<int> elPMaxOrder) const
{
	int pmax = 0;
	for (int ip=0;  ip<elPMaxOrder.size(); ip++) 
	{
		if(elPMaxOrder[ip] > pmax) pmax = elPMaxOrder[ip];  
	}
	
	return  2*pmax;
}

void TPZMaterial::FillDataRequirements(TPZVec<TPZMaterialData > &datavec)
{
	int nref = datavec.size();
	for(int i = 0; i<nref; i++ )
	{
		datavec[i].SetAllRequirements(true);
		datavec[i].fNeedsNeighborSol = false;
		datavec[i].fNeedsNeighborCenter = false;
		datavec[i].fNeedsNormal = false;
	}
	
}

void TPZMatPoissonD3::ContributeBCInterface(TPZMaterialData &data, TPZVec<TPZMaterialData> &dataleft, REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
{
    
    
#ifdef PZDEBUG
	int nref =  dataleft.size();
	if (nref != 2 ) {
        std::cout << " Error. This implementation needs only two computational meshes. \n";
		DebugStop();
	}
#endif
    
#ifdef PZDEBUG
	int bref =  bc.Val2().Rows();
	if (bref != 2 ) {
        std::cout << " Erro. The size of the datavec is different from 2 \n";
		DebugStop();
	}
#endif
    
    REAL Qn = bc.Val2()(0,0); // cuidado para, na hora de passar os valores de cond contorno, seguir essa ordem
    REAL Pd = 0.0; // = bc.Val2()(1,0); // fluxo normal na primeira casa e pressao na segunda
    
	TPZManVector<REAL,3> &normal = data.normal;
	//REAL n1 = normal[0];
	//REAL n2 = normal[1];
    
    //REAL v2;
    if(bc.HasForcingFunction())
    {
		TPZManVector<STATE> res(3);
		bc.ForcingFunction()->Execute(dataleft[0].x,res);
		Pd = res[0];
        Qn = res[0];
	}else
    {
        Pd = bc.Val2()(1,0);
    }
    

    // Setting the phis
    TPZFMatrix<REAL>  &phiQ =  dataleft[0].phi;
    TPZFMatrix<REAL>  &phip =  dataleft[1].phi;
	//TPZFMatrix<REAL> &dphiQ = datavec[0].dphix;
    //TPZFMatrix<REAL> &dphip = datavec[1].dphix;

    int phrq, phrp;
    phrp = phip.Rows();
    phrq = dataleft[0].fVecShapeIndex.NElements();

	//Calculate the matrix contribution for boundary conditions
    for (int iq = 0; iq<phrq; iq++)
    {
        int ivecind = dataleft[0].fVecShapeIndex[iq].first;
        int ishapeind = dataleft[0].fVecShapeIndex[iq].second;
        TPZFNMatrix<3> ivec(3,1);
        ivec(0,0) = dataleft[0].fNormalVec(0,ivecind);
        ivec(1,0) = dataleft[0].fNormalVec(1,ivecind);
        ivec(2,0) = dataleft[0].fNormalVec(2,ivecind);
        ivec *= phiQ(ishapeind,0);
        
        
        REAL NormalProjectioni = 0.;
        for(int iloc=0; iloc<fDim; iloc++)
        {
            NormalProjectioni += ivec(iloc,0)*normal[iloc];
        }
        
        for (int jp=0; jp<phrp; jp++)
        {
            
            REAL integration = weight*NormalProjectioni*phip(jp,0);
            
            //para a equacao do fluxo - 1o conjunto da formulacao
            ek(iq, phrq+jp) += (-1.0)*integration;
            
            // para a equacao da pressao - 2o conjunto da formulacao
            ek(phrq+jp, iq) += (-1.0)*integration;
            
        }
    }
    

    //if (bc.Type()==0){std::cout << "...." << std::endl;}

    switch (bc.Type())
    {  
        case 0:  // Dirichlet
        {
            //REAL InvK = 1./fK;
            
                        //termo fonte referente a equacao do fluxo
            for (int iq = 0; iq<phrq; iq++)
            {
                int ivecind = dataleft[0].fVecShapeIndex[iq].first;
                int ishapeind = dataleft[0].fVecShapeIndex[iq].second;
                TPZFNMatrix<3> ivec(3,1);
                ivec(0,0) = dataleft[0].fNormalVec(0,ivecind);
                ivec(1,0) = dataleft[0].fNormalVec(1,ivecind);
                ivec(2,0) = dataleft[0].fNormalVec(2,ivecind);
//.........这里部分代码省略.........

void TPZTracerFlow::ContributeBC(TPZVec<TPZMaterialData> &datavec,REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc){
    
    if (fPressureEquationFilter == false)
    {
        return;
    }
#ifdef PZDEBUG
    int nref =  datavec.size();
	if (nref != 3 ) {
        std::cout << " Erro.!! datavec tem que ser de tamanho 2 \n";
		DebugStop();
	}
#endif
	
	TPZFMatrix<REAL>  &phiQ = datavec[1].phi;
	int phrQ = phiQ.Rows();//datavec[1].fVecShapeIndex.NElements();
    int phrS = datavec[0].phi.Rows();
    
	REAL v2;
	v2 = bc.Val2()(1,0);
    
    STATE BigNum = gBigNumber;
	
	switch (bc.Type()) {
		case 0 :		// Dirichlet condition
			//primeira equacao
			for(int iq=0; iq<phrQ; iq++)
            {
                //the contribution of the Dirichlet boundary condition appears in the flow equation
                ef(iq+phrS,0) += (-1.)*v2*phiQ(iq,0)*weight;
            }
            break;
			
		case 1 :			// Neumann condition
			//primeira equacao
            if(IsZero(v2)) BigNum = 1.e10;
			for(int iq=0; iq<phrQ; iq++)
            {
                ef(iq+phrS,0)+= BigNum*v2*phiQ(iq,0)*weight;
                for (int jq=0; jq<phrQ; jq++) {
                    
                    ek(iq+phrS,jq+phrS)+= BigNum*phiQ(iq,0)*phiQ(jq,0)*weight;
                }
            }
			break;
            
        case 2 :			// mixed condition
            for(int iq = 0; iq < phrQ; iq++) {
                
				ef(iq+phrS,0) += v2*phiQ(iq,0)*weight;
				for (int jq = 0; jq < phrQ; jq++) {
					ek(iq+phrS,jq+phrS) += weight*bc.Val1()(0,0)*phiQ(iq,0)*phiQ(jq,0);
				}
			}
            
            break;
            
        case 5 :        // Neumann(pressure)-Inflow(saturation)
			//primeira equacao
			for(int iq=0; iq<phrQ; iq++)
            {
                ef(iq+phrS,0)+= BigNum*v2*phiQ(iq,0)*weight;
                for (int jq=0; jq<phrQ; jq++) {
                    
                    ek(iq+phrS,jq+phrS)+= BigNum*phiQ(iq,0)*phiQ(jq,0)*weight;
                }
            }
			break;
            
        case 6 :		// Dirichlet(pressure)-Outflow(saturation)
			//primeira equacao
			for(int iq=0; iq<phrQ; iq++)
            {
                //the contribution of the Dirichlet boundary condition appears in the flow equation
                ef(iq+phrS,0) += (-1.)*v2*phiQ(iq,0)*weight;
            }
            break;
            
        case 7 :		// Dirichlet(pressure)-Inflow(saturation)
			//primeira equacao
			for(int iq=0; iq<phrQ; iq++)
            {
                //the contribution of the Dirichlet boundary condition appears in the flow equation
                ef(iq+phrS,0) += (-1.)*v2*phiQ(iq,0)*weight;
            }
            break;
	}
}

void TPZTracerFlow::Contribute(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef){
    
#ifdef PZDEBUG
	int nref =  datavec.size();
	if (nref != 3 ) {
        std::cout << " Erro. The size of the datavec is different from 3 \n";
		DebugStop();
	}
#endif
    
    // Setting the phis
    TPZFMatrix<REAL>  &phiQ =  datavec[1].phi;
    TPZFMatrix<REAL> &dphiQ = datavec[1].dphix;
    TPZFMatrix<REAL>  &phiP =  datavec[2].phi;
    TPZFMatrix<REAL>  &phiS =  datavec[0].phi;
    TPZFMatrix<REAL>  &dphiS =  datavec[0].dphix;
    
	
	TPZFMatrix<REAL> &axesS = datavec[0].axes;
	int phrQ = datavec[1].fVecShapeIndex.NElements();//phiQ.Rows();
    int phrP = phiP.Rows();
    int phrS = phiS.Rows();
    
    //current state n+1: stiffness matrix
    if(gState == ECurrentState)
    {
        if (fPressureEquationFilter == true)
        {
            STATE permeability = fk;
            if(fForcingFunction) {
                TPZManVector<STATE> res(1);
                fForcingFunction->Execute(datavec[1].x,res);
                permeability = res[0];
            }
            
            //Calculate the matrix contribution for flux. Matrix A
            REAL ratiok = fVisc/permeability;
            for(int iq=0; iq<phrQ; iq++)
            {
                ef(iq+phrS, 0) += 0.;
                
                int ivecind = datavec[1].fVecShapeIndex[iq].first;
                int ishapeind = datavec[1].fVecShapeIndex[iq].second;
                for (int jq=0; jq<phrQ; jq++)
                {
                    int jvecind = datavec[1].fVecShapeIndex[jq].first;
                    int jshapeind = datavec[1].fVecShapeIndex[jq].second;
                    REAL prod = datavec[1].fNormalVec(0,ivecind)*datavec[1].fNormalVec(0,jvecind)+
                    datavec[1].fNormalVec(1,ivecind)*datavec[1].fNormalVec(1,jvecind)+
                    datavec[1].fNormalVec(2,ivecind)*datavec[1].fNormalVec(2,jvecind);//dot product between u and v
                    ek(iq+phrS,jq+phrS) += ratiok*weight*phiQ(ishapeind,0)*phiQ(jshapeind,0)*prod;
                }
            }
            
            // Coupling terms between flux and pressure. Matrix B
            for(int iq=0; iq<phrQ; iq++)
            {
                int ivecind = datavec[1].fVecShapeIndex[iq].first;
                int ishapeind = datavec[1].fVecShapeIndex[iq].second;
                
                TPZFNMatrix<3> ivec(3,1);
                ivec(0,0) = datavec[1].fNormalVec(0,ivecind);
                ivec(1,0) = datavec[1].fNormalVec(1,ivecind);
                ivec(2,0) = datavec[1].fNormalVec(2,ivecind);
                TPZFNMatrix<3> axesvec(3,1);
                datavec[1].axes.Multiply(ivec,axesvec);
                
                REAL divwq = 0.;
                for(int iloc=0; iloc<fDim; iloc++)
                {
                    divwq += axesvec(iloc,0)*dphiQ(iloc,ishapeind);
                }
                for (int jp=0; jp<phrP; jp++) {
                    
                    REAL fact = (-1.)*weight*phiP(jp,0)*divwq;
                    // Matrix B
                    ek(iq+phrS, phrS+phrQ+jp) += fact;
                    
                    // Matrix B^T
                    ek(phrS+phrQ+jp,iq+phrS) += fact;
                }
            }
            
            //Right side of the pressure equation
            REAL fXfLocP = fxfPQ;
//            if(fForcingFunction) {
//                TPZManVector<STATE> res(1);
//                fForcingFunction->Execute(datavec[2].x,res);
//                fXfLocP = res[0];
//            }
            for(int ip=0; ip<phrP; ip++){
                ef(phrS+phrQ+ip,0) += (-1.)*weight*fXfLocP*phiP(ip,0);
            }
        }
        else
        {
            STATE DeltaT = fTimeStep;
            STATE forceSat = fxfS;
            //Calculate the matrix contribution for saturation.
            for(int in = 0; in < phrS; in++ ) {
//.........这里部分代码省略.........