C++ VectorType::resize方法代码示例

 void BeamPointPressureCondition::InitializeSystemMatrices( MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector) {
     
       rLeftHandSideMatrix.resize( 6, 6, false );
       noalias( rLeftHandSideMatrix ) = ZeroMatrix( 6, 6 ); //resetting LHS
       rRightHandSideVector.resize( 6, false );
       rRightHandSideVector = ZeroVector( 6 ); //resetting RHS
 }

//************************************************************************************
//************************************************************************************
void ThermalFace2D::CalculateAll(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo, bool CalculateStiffnessMatrixFlag, bool CalculateResidualVectorFlag)
{
    KRATOS_TRY

    unsigned int number_of_nodes = GetGeometry().size();

    //resizing as needed the LHS
    unsigned int MatSize=number_of_nodes;

    ConvectionDiffusionSettings::Pointer my_settings = rCurrentProcessInfo.GetValue(CONVECTION_DIFFUSION_SETTINGS);

    const Variable<double>& rUnknownVar = my_settings->GetUnknownVariable();

    const Variable<double>& rSurfaceSourceVar = my_settings->GetSurfaceSourceVariable();

    //calculate lenght
    double x21 = GetGeometry()[1].X() - GetGeometry()[0].X();
    double y21 = GetGeometry()[1].Y() - GetGeometry()[0].Y();
    double lenght = x21*x21 + y21*y21;
    lenght = sqrt(lenght);

    const Properties& ConstProp = GetProperties();
    const double& ambient_temperature = ConstProp[AMBIENT_TEMPERATURE];
    double StefenBoltzmann = 5.67e-8;
    double emissivity = ConstProp[EMISSIVITY];
    double convection_coefficient = ConstProp[CONVECTION_COEFFICIENT];

    const double& T0 = GetGeometry()[0].FastGetSolutionStepValue(rUnknownVar);
    const double& T1 = GetGeometry()[1].FastGetSolutionStepValue(rUnknownVar);

    const double& q0 =GetGeometry()[0].FastGetSolutionStepValue(rSurfaceSourceVar);
    const double& q1 =GetGeometry()[1].FastGetSolutionStepValue(rSurfaceSourceVar);

    if (CalculateStiffnessMatrixFlag == true) //calculation of the matrix is required
    {
        if(rLeftHandSideMatrix.size1() != MatSize )
            rLeftHandSideMatrix.resize(MatSize,MatSize,false);
        noalias(rLeftHandSideMatrix) = ZeroMatrix(MatSize,MatSize);

        rLeftHandSideMatrix(0,0) = ( convection_coefficient + emissivity*StefenBoltzmann*4.0*pow(T0,3)  )* 0.5 * lenght;
        rLeftHandSideMatrix(1,1) = ( convection_coefficient + emissivity*StefenBoltzmann*4.0*pow(T1,3)  )* 0.5 * lenght;
    }

    //resizing as needed the RHS
    double aux = pow(ambient_temperature,4);
    if (CalculateResidualVectorFlag == true) //calculation of the matrix is required
    {
        if(rRightHandSideVector.size() != MatSize )
            rRightHandSideVector.resize(MatSize,false);

        rRightHandSideVector[0] =  q0 - emissivity*StefenBoltzmann*(pow(T0,4) - aux)  -  convection_coefficient * ( T0 - ambient_temperature);

        rRightHandSideVector[1] =  q1   - emissivity*StefenBoltzmann*(pow(T1,4) - aux) -  convection_coefficient * ( T1 - ambient_temperature);

        rRightHandSideVector *= 0.5*lenght;

    }

    KRATOS_CATCH("")
}

/**
 * calculates this contact element's local contributions
 */
void SlaveContactFace3DNewmark::CalculateLocalSystem( MatrixType& rLeftHandSideMatrix,
        VectorType& rRightHandSideVector,
        ProcessInfo& rCurrentProcessInfo)
{
    rRightHandSideVector.resize(0,false);
    rLeftHandSideMatrix(0,0);
}

/**
 * calculates only the RHS vector (certainly to be removed due to contact algorithm)
 */
void MasterContactPoint2D::CalculateRightHandSide( VectorType& rRightHandSideVector,
        ProcessInfo& rCurrentProcessInfo)
{
    unsigned int ndof = GetGeometry().size()*2;
    if( rRightHandSideVector.size() != ndof )
        rRightHandSideVector.resize(ndof,false);
    rRightHandSideVector = ZeroVector(ndof);
}

void SolidFace3D::CalculateRightHandSide(VectorType& rRightHandSideVector,
                                         ProcessInfo& r_process_info)
{
    const unsigned int number_of_nodes = GetGeometry().size();
    unsigned int MatSize = number_of_nodes * 3;
    
    if (rRightHandSideVector.size() != MatSize)
    {
        rRightHandSideVector.resize(MatSize, false);
    }
    rRightHandSideVector = ZeroVector(MatSize);
    
    std::vector<SphericParticle*>& rNeighbours = this->mNeighbourSphericParticles;
    
    for (unsigned int i=0; i<rNeighbours.size(); i++)
    {
        if(rNeighbours[i]->Is(BLOCKED)) continue; //Inlet Generator Spheres are ignored when integrating forces.
        
        std::vector<DEMWall*>& rRFnei = rNeighbours[i]->mNeighbourRigidFaces;
                
        for (unsigned int i_nei = 0; i_nei < rRFnei.size(); i_nei++)
        {
            int Contact_Type = rNeighbours[i]->mContactConditionContactTypes[i_nei];
            
            if ( ( rRFnei[i_nei]->Id() == this->Id() ) && (Contact_Type > 0 ) )
            {
                
                array_1d<double, 4> weights_vector = rNeighbours[i]->mContactConditionWeights[i_nei];
                double weight = 0.0;
                
                double ContactForce[3] = {0.0};

                const array_1d<double, 3>& neighbour_rigid_faces_contact_force = rNeighbours[i]->mNeighbourRigidFacesTotalContactForce[i_nei];

                ContactForce[0] = neighbour_rigid_faces_contact_force[0];
                ContactForce[1] = neighbour_rigid_faces_contact_force[1];
                ContactForce[2] = neighbour_rigid_faces_contact_force[2];

                for (unsigned int k=0; k< number_of_nodes; k++)
                {
                    weight = weights_vector[k];
  
                    unsigned int w =  k * 3;

                    rRightHandSideVector[w + 0] += -ContactForce[0] * weight;
                    rRightHandSideVector[w + 1] += -ContactForce[1] * weight;
                    rRightHandSideVector[w + 2] += -ContactForce[2] * weight;
                }
                
            }//if the condition neighbour of my sphere neighbour is myself.
        }//Loop spheres neighbours (condition)
    }//Loop condition neighbours (spheres)
}//CalculateRightHandSide

inline
void StackContextBase<TSuperClass>::setSlotVariable(const VariableSlotID slot,
                                                    const UnitType &newValue,
                                                    VectorType &container) const
{
    if(slot < container.size())
        container.replace(slot, newValue);
    else
    {
        container.resize(slot + 1);
        container.replace(slot, newValue);
    }
}

void AbsoluteEulerAngleDecoder::Encode(array_view<const DirectX::Quaternion> rots, VectorType & x)
{
	int n = rots.size();

	Eigen::Vector4f qs;
	XMVECTOR q;
	qs.setZero();
	x.resize(n * 3);
	for (int i = 0; i < n; i++)
	{
		q = XMLoad(rots[i]);
		q = XMQuaternionEulerAngleYawPitchRoll(q); // Decompsoe in to euler angle
		XMStoreFloat4(qs.data(), q);
		x.segment<3>(i * 3) = qs.head<3>();
	}
}

void PeriodicConditionLM2D2N::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
    if(rLeftHandSideMatrix.size1() != 6 || rLeftHandSideMatrix.size2() != 6) rLeftHandSideMatrix.resize(6, 6,false);
    noalias(rLeftHandSideMatrix) = ZeroMatrix(6, 6);

    if(rRightHandSideVector.size() != 6) rRightHandSideVector.resize(6, false);
	noalias( rRightHandSideVector ) = ZeroVector(6);

	// Nodal IDs = [slave ID, master ID]
	
	GeometryType& geom = GetGeometry();

	// get current values and form the system matrix

	Vector currentValues(6);
	currentValues(0) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_X);
	currentValues(1) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_Y);
	currentValues(2) = geom[1].FastGetSolutionStepValue(DISPLACEMENT_X);
	currentValues(3) = geom[1].FastGetSolutionStepValue(DISPLACEMENT_Y);
	currentValues(4) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_LAGRANGE_X);
	currentValues(5) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_LAGRANGE_Y);
	
	//KRATOS_WATCH(currentValues);

	rLeftHandSideMatrix(4,0) =  1.0;
	rLeftHandSideMatrix(4,2) = -1.0;
	rLeftHandSideMatrix(0,4) =  1.0;
	rLeftHandSideMatrix(2,4) = -1.0;
	rLeftHandSideMatrix(5,1) =  1.0;
	rLeftHandSideMatrix(5,3) = -1.0;
	rLeftHandSideMatrix(1,5) =  1.0;
	rLeftHandSideMatrix(3,5) = -1.0;

	// form residual

	noalias(rRightHandSideVector) -= prod( rLeftHandSideMatrix, currentValues );
}

//************************************************************************************
//************************************************************************************
void Monolithic2DNeumann::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{

    if(rLeftHandSideMatrix.size1() != 4)
    {
        rLeftHandSideMatrix.resize(4,4,false);
        rRightHandSideVector.resize(4,false);

    }

    noalias(rLeftHandSideMatrix) = ZeroMatrix(4,4);

    //calculate normal to element.(normal follows the cross rule)
    array_1d<double,2> An,edge;
    edge[0] = GetGeometry()[1].X() - GetGeometry()[0].X();
    edge[1] = GetGeometry()[1].Y() - GetGeometry()[0].Y();


    double norm = edge[0]*edge[0] + edge[1]*edge[1];
    norm = pow(norm,0.5);

    An[0] = -edge[1];
    An[1] = edge[0];
    //An /= norm; this is then simplified by length of element in integration so is not divided.

    double mean_ex_p = 0.0;

    for(unsigned int i = 0; i<2 ; i++)
        mean_ex_p += 0.5*GetGeometry()[i].FastGetSolutionStepValue(EXTERNAL_PRESSURE);

    double p0 = GetGeometry()[0].FastGetSolutionStepValue(EXTERNAL_PRESSURE);
    rRightHandSideVector[0] = -0.5*An[0]*p0;
    rRightHandSideVector[1] = -0.5*An[1]*p0;

    double p1 = GetGeometry()[1].FastGetSolutionStepValue(EXTERNAL_PRESSURE);
    rRightHandSideVector[2] = -0.5*An[0]*p1;
    rRightHandSideVector[3] = -0.5*An[1]*p1;

    //	if(mean_ex_p !=GetGeometry()[0].FastGetSolutionStepValue(EXTERNAL_PRESSURE))
    //		mean_ex_p = 0.0;
    //KRATOS_WATCH(mean_ex_p);

    /*			for(unsigned int ii = 0; ii< 2; ++ii)
    				{

    					int id = (2 + 1)*(ii);
    					rRightHandSideVector[id] = mean_ex_p * An[0]* 0.5;
    					rRightHandSideVector[id + 1] = mean_ex_p * An[1]* 0.5;
    					rRightHandSideVector[id + 2] = 0.0;
    				//KRATOS_WATCH(An);
    				}*/
// 			KRATOS_WATCH(An);
//KRATOS_WATCH(p0);
//KRATOS_WATCH(p1);
//KRATOS_WATCH("TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT");
    //	KRATOS_WATCH(rRightHandSideVector);




}

 void resize(VectorType & vec, size_t new_size)
 {
   vec.resize(new_size); 
 }

//***********************************************************************************
//***********************************************************************************
void LineForce::CalculateRightHandSide( VectorType& rRightHandSideVector,
        ProcessInfo& rCurrentProcessInfo )
{
    KRATOS_TRY

    const unsigned int number_of_nodes = GetGeometry().size();
    const unsigned int dim = GetGeometry().WorkingSpaceDimension();
    unsigned int MatSize = number_of_nodes * dim;

    //resizing as needed the RHS
    if ( rRightHandSideVector.size() != MatSize )
        rRightHandSideVector.resize( MatSize, false );
    rRightHandSideVector = ZeroVector( MatSize ); //resetting RHS

    //reading integration points and local gradients
    const GeometryType::IntegrationPointsArrayType& integration_points =
        GetGeometry().IntegrationPoints();

    // DN_DeContainer is the array of shape function gradients at each integration points
    const GeometryType::ShapeFunctionsGradientsType& DN_DeContainer =
        GetGeometry().ShapeFunctionsLocalGradients();

    // Ncontainer is the array of shape function values at each integration points
    const Matrix& Ncontainer = GetGeometry().ShapeFunctionsValues();

    //loop over integration points
    Vector Load( dim );
    Vector LoadOnNode( dim );
    for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); ++PointNumber )
    {
        noalias( Load ) = ZeroVector( dim );

        for ( unsigned int n = 0; n < GetGeometry().size(); n++ )
        {
            noalias( LoadOnNode ) = ( GetGeometry()[n] ).GetSolutionStepValue( FACE_LOAD );

            for ( unsigned int i = 0; i < dim; i++ )
            {
                Load( i ) += LoadOnNode( i ) * Ncontainer( PointNumber, n );
            }
        }
        
        double IntegrationWeight = GetGeometry().IntegrationPoints()[PointNumber].Weight();

//        if(dim == 2) IntegrationWeight *= GetProperties()[THICKNESS]; // TODO: check

        Vector t = ZeroVector( dim );//tangential vector
        for ( unsigned int n = 0; n < GetGeometry().size(); ++n )
        {
            t[0] += GetGeometry().GetPoint( n ).X0() * DN_DeContainer[PointNumber]( n, 0 );
            t[1] += GetGeometry().GetPoint( n ).Y0() * DN_DeContainer[PointNumber]( n, 0 );
            if(dim == 3)
                t[2] += GetGeometry().GetPoint( n ).Z0() * DN_DeContainer[PointNumber]( n, 0 );
        }

        // calculating length
        double dL = norm_2(t);
        
        // RIGHT HAND SIDE VECTOR
        for ( unsigned int prim = 0; prim < GetGeometry().size(); ++prim )
            for ( unsigned int i = 0; i < dim; ++i )
                rRightHandSideVector( prim * dim + i ) +=
                    Ncontainer( PointNumber, prim ) * Load( i ) * IntegrationWeight * dL;
    }

    KRATOS_CATCH( "" )
}

//************************************************************************************
//************************************************************************************
void Electrostatic2D::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
    KRATOS_TRY

    const unsigned int number_of_points = GetGeometry().size();

    if(rLeftHandSideMatrix.size1() != number_of_points)
        rLeftHandSideMatrix.resize(number_of_points,number_of_points,false);

    if(rRightHandSideVector.size() != number_of_points)
        rRightHandSideVector.resize(number_of_points,false);

    //getting data for the given geometry
    double Area;
    GeometryUtils::CalculateGeometryData(GetGeometry(), msDN_DX, msN, Area);

    //reading properties and conditions
    array_1d<double,3> permittivity = GetProperties()[ELECTRICAL_PERMITTIVITY];
    msD(0,0)=permittivity[0];
    msD(1,1)=permittivity[1];
    msD(1,0)=0.0;
    msD(0,1)=0.0;
    //point_sources[0] = GetGeometry()[0].FastGetSolutionStepValue(ELECTROSTATIC_POINT_CHARGE);
    //point_sources[1] = GetGeometry()[1].FastGetSolutionStepValue(ELECTROSTATIC_POINT_CHARGE);
    //point_sources[2] = GetGeometry()[2].FastGetSolutionStepValue(ELECTROSTATIC_POINT_CHARGE);

    //double surface_sources = (this)->GetValue(ELECTROSTATIC_SURFACE_CHARGE);
    surface_sources[0] = (this)->GetValue(ELECTROSTATIC_SURFACE_CHARGE);
    surface_sources[1] = (this)->GetValue(ELECTROSTATIC_SURFACE_CHARGE);
    surface_sources[2] = (this)->GetValue(ELECTROSTATIC_SURFACE_CHARGE);
    //surface_sources[0] = GetGeometry()[0].FastGetSolutionStepValue(ELECTROSTATIC_SURFACE_CHARGE);
    //surface_sources[1] = GetGeometry()[1].FastGetSolutionStepValue(ELECTROSTATIC_SURFACE_CHARGE);
    //surface_sources[2] = GetGeometry()[2].FastGetSolutionStepValue(ELECTROSTATIC_SURFACE_CHARGE);

    // main loop
    const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints();

    noalias(rLeftHandSideMatrix) = prod(msDN_DX,Matrix(prod(msD,trans(msDN_DX))));

    /*		for(unsigned int k = 1; k<integration_points.size(); k++)	//integration points
    		{
    			double w_detj = integration_points[k].Weight()*mDetJo[k];
    */
    rLeftHandSideMatrix *= Area;
//		}


    //Point charge contribution
    //noalias(rRightHandSideVector) = point_sources;
    //noalias(rRightHandSideVector) = point_sources/3.0;
    //+surface_sources*Area/3.0;
    noalias(rRightHandSideVector) = surface_sources*Area/3.0;
    //subtracting the dirichlet term
    // RHS -= LHS*ELECTROSTATIC_POTENTIALs
    for(unsigned int iii = 0; iii<number_of_points; iii++)
        ms_temp[iii] = GetGeometry()[iii].FastGetSolutionStepValue(ELECTROSTATIC_POTENTIAL);
    noalias(rRightHandSideVector) -= prod(rLeftHandSideMatrix,ms_temp);

    //multiplying by area, rho and density
    //rRightHandSideVector *= (Area * permittivity);
    //rLeftHandSideMatrix *= (Area * permittivity);
    KRATOS_CATCH("");
}

//----------------------
//-----  PRIVATE  ------
//----------------------
//***********************************************************************************
void FaceHeatConvection::CalculateAll( MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, const ProcessInfo& rCurrentProcessInfo, bool CalculateStiffnessMatrixFlag, bool CalculateResidualVectorFlag )
{
    KRATOS_TRY

    const unsigned int number_of_nodes = GetGeometry().size();
    unsigned int MatSize = number_of_nodes;

    //resizing as needed the LHS
    if ( CalculateStiffnessMatrixFlag == true ) //calculation of the matrix is required
    {
        if ( rLeftHandSideMatrix.size1() != MatSize )
            rLeftHandSideMatrix.resize( MatSize, MatSize, false );
        noalias( rLeftHandSideMatrix ) = ZeroMatrix( MatSize, MatSize ); //resetting LHS
    }

    //resizing as needed the RHS
    if ( CalculateResidualVectorFlag == true ) //calculation of the matrix is required
    {
        if ( rRightHandSideVector.size() != MatSize )
            rRightHandSideVector.resize( MatSize );
        rRightHandSideVector = ZeroVector( MatSize ); //resetting RHS
    }

    //reading integration points and local gradients
    const GeometryType::IntegrationPointsArrayType& integration_points = GetGeometry().IntegrationPoints();
    const GeometryType::ShapeFunctionsGradientsType& DN_DeContainer = GetGeometry().ShapeFunctionsLocalGradients();
    const Matrix& Ncontainer = GetGeometry().ShapeFunctionsValues();

    //calculating actual jacobian
    GeometryType::JacobiansType J;
    J = GetGeometry().Jacobian( J );

    //auxiliary terms
    for ( unsigned int PointNumber = 0; PointNumber < integration_points.size(); PointNumber++ )
    {
        double convection_coefficient = 0.0;
        double T_ambient = 0.0;
        double T = 0.0;


        for ( unsigned int n = 0; n < GetGeometry().size(); n++ )
        {
            convection_coefficient += ( GetGeometry()[n] ).GetSolutionStepValue( CONVECTION_COEFFICIENT ) * Ncontainer( PointNumber, n );
            T_ambient += ( GetGeometry()[n] ).GetSolutionStepValue( AMBIENT_TEMPERATURE ) * Ncontainer( PointNumber, n );
            T += ( GetGeometry()[n] ).GetSolutionStepValue( TEMPERATURE ) * Ncontainer( PointNumber, n );
        }

//         if ( PointNumber == 1 )
//             std::cout << "CONDITION ### HeatConvection:  h= " << convection_coefficient << ",\t T_ambient= " << T_ambient << ",\t T_surface= " << T << std::endl;

        double IntegrationWeight = GetGeometry().IntegrationPoints()[PointNumber].Weight();
        Vector t1 = ZeroVector( 3 );//first tangential vector
        Vector t2 = ZeroVector( 3 );//second tangential vector

        for ( unsigned int n = 0; n < number_of_nodes; n++ )
        {
            t1[0] += GetGeometry().GetPoint( n ).X0() * DN_DeContainer[PointNumber]( n, 0 );
            t1[1] += GetGeometry().GetPoint( n ).Y0() * DN_DeContainer[PointNumber]( n, 0 );
            t1[2] += GetGeometry().GetPoint( n ).Z0() * DN_DeContainer[PointNumber]( n, 0 );
            t2[0] += GetGeometry().GetPoint( n ).X0() * DN_DeContainer[PointNumber]( n, 1 );
            t2[1] += GetGeometry().GetPoint( n ).Y0() * DN_DeContainer[PointNumber]( n, 1 );
            t2[2] += GetGeometry().GetPoint( n ).Z0() * DN_DeContainer[PointNumber]( n, 1 );
        }

        //calculating normal
        Vector v3 = ZeroVector( 3 );
        v3[0] = t1[1] * t2[2] - t1[2] * t2[1];
        v3[1] = t1[2] * t2[0] - t1[0] * t2[2];
        v3[2] = t1[0] * t2[1] - t1[1] * t2[0];

        double	dA = sqrt( v3[0] * v3[0] + v3[1] * v3[1] + v3[2] * v3[2] );

        if ( CalculateResidualVectorFlag == true ) //calculation of the matrix is required
        {
            for ( unsigned int n = 0; n < number_of_nodes; n++ )
                rRightHandSideVector( n ) -= Ncontainer( PointNumber, n )
                                             * convection_coefficient * ( T - T_ambient )
                                             * IntegrationWeight * dA; // W/(m^2*°C) = kg*s^-3*°C°^-1
        }

        if ( CalculateStiffnessMatrixFlag == true ) //calculation of the matrix is required
        {
            for ( unsigned int n = 0; n < number_of_nodes; n++ )
                rLeftHandSideMatrix( n, n ) += Ncontainer( PointNumber, n )
                                               * convection_coefficient
                                               * Ncontainer( PointNumber, n ) * IntegrationWeight * dA; // W/(m^2*°C) = kg*s^-3*°C°^-1
        }

    }
    KRATOS_CATCH( "" )
}

//************************************************************************************
//************************************************************************************
//calculation by component of the fractional step velocity corresponding to the first stage
void NDFluid2DCrankNicolson::Stage1(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector,
                                    ProcessInfo& rCurrentProcessInfo, unsigned int ComponentIndex)
{
    KRATOS_TRY;

    const unsigned int number_of_points = 3;

    if(rLeftHandSideMatrix.size1() != number_of_points)
        rLeftHandSideMatrix.resize(number_of_points,number_of_points,false); //false says not to preserve existing storage!!

    if(rRightHandSideVector.size() != number_of_points)
        rRightHandSideVector.resize(number_of_points,false); //false says not to preserve existing storage!!

    //getting data for the given geometry
    double Area;
    GeometryUtils::CalculateGeometryData(GetGeometry(),msDN_DX,msN,Area);

    //getting the velocity vector on the nodes

    //getting the velocity on the nodes
    const array_1d<double,3>& fv0 = GetGeometry()[0].FastGetSolutionStepValue(FRACT_VEL,0);
    const array_1d<double,3>& fv0_old = GetGeometry()[0].FastGetSolutionStepValue(VELOCITY,1);
    const array_1d<double,3>& w0 = GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY);
    const array_1d<double,3>& w0_old = GetGeometry()[0].FastGetSolutionStepValue(MESH_VELOCITY,1);
    const array_1d<double,3>& proj0 = GetGeometry()[0].FastGetSolutionStepValue(CONV_PROJ);
    const array_1d<double,3>& proj0_old = GetGeometry()[0].FastGetSolutionStepValue(CONV_PROJ,1);
    double p0old = GetGeometry()[0].FastGetSolutionStepValue(PRESSURE_OLD_IT);
    const double nu0 = GetGeometry()[0].FastGetSolutionStepValue(VISCOSITY);
    const double rho0 = GetGeometry()[0].FastGetSolutionStepValue(DENSITY);
    const double fcomp0 = GetGeometry()[0].FastGetSolutionStepValue(BODY_FORCE)[ComponentIndex];
    const double eps0 = GetGeometry()[0].FastGetSolutionStepValue(POROSITY);
    const double dp0 = GetGeometry()[0].FastGetSolutionStepValue(DIAMETER);


    const array_1d<double,3>& fv1 = GetGeometry()[1].FastGetSolutionStepValue(FRACT_VEL);
    const array_1d<double,3>& fv1_old = GetGeometry()[1].FastGetSolutionStepValue(VELOCITY,1);
    const array_1d<double,3>& w1 = GetGeometry()[1].FastGetSolutionStepValue(MESH_VELOCITY);
    const array_1d<double,3>& w1_old = GetGeometry()[1].FastGetSolutionStepValue(MESH_VELOCITY,1);
    const array_1d<double,3>& proj1 = GetGeometry()[1].FastGetSolutionStepValue(CONV_PROJ);
    double p1old = GetGeometry()[1].FastGetSolutionStepValue(PRESSURE_OLD_IT);
    const array_1d<double,3>& proj1_old = GetGeometry()[1].FastGetSolutionStepValue(CONV_PROJ,1);
    const double nu1 = GetGeometry()[1].FastGetSolutionStepValue(VISCOSITY);
    const double rho1 = GetGeometry()[1].FastGetSolutionStepValue(DENSITY);
    const double fcomp1 = GetGeometry()[1].FastGetSolutionStepValue(BODY_FORCE)[ComponentIndex];
    const double eps1 = GetGeometry()[1].FastGetSolutionStepValue(POROSITY);
    const double dp1 = GetGeometry()[1].FastGetSolutionStepValue(DIAMETER);


    const array_1d<double,3>& fv2 = GetGeometry()[2].FastGetSolutionStepValue(FRACT_VEL);
    const array_1d<double,3>& fv2_old = GetGeometry()[2].FastGetSolutionStepValue(VELOCITY,1);
    const array_1d<double,3>& w2 = GetGeometry()[2].FastGetSolutionStepValue(MESH_VELOCITY);
    const array_1d<double,3>& w2_old = GetGeometry()[2].FastGetSolutionStepValue(MESH_VELOCITY,1);
    const array_1d<double,3>& proj2 = GetGeometry()[2].FastGetSolutionStepValue(CONV_PROJ);
    const array_1d<double,3>& proj2_old = GetGeometry()[2].FastGetSolutionStepValue(CONV_PROJ,1);
    double p2old = GetGeometry()[2].FastGetSolutionStepValue(PRESSURE_OLD_IT);
    const double nu2 = GetGeometry()[2].FastGetSolutionStepValue(VISCOSITY);
    const double rho2 = GetGeometry()[2].FastGetSolutionStepValue(DENSITY);
    const double fcomp2 = GetGeometry()[2].FastGetSolutionStepValue(BODY_FORCE)[ComponentIndex];
    const double eps2 = GetGeometry()[2].FastGetSolutionStepValue(POROSITY);
    const double dp2 = GetGeometry()[2].FastGetSolutionStepValue(DIAMETER);


    //
    // vel_gauss = sum( N[i]*(vel[i]-mesh_vel[i]), i=0, number_of_points) (note that the fractional step vel is used)
    ms_vel_gauss[0] =  msN[0]*(fv0[0]-w0[0]) + msN[1]*(fv1[0]-w1[0]) +  msN[2]*(fv2[0]-w2[0]);
    ms_vel_gauss[1] =  msN[0]*(fv0[1]-w0[1]) + msN[1]*(fv1[1]-w1[1]) +  msN[2]*(fv2[1]-w2[1]);

    //vel_gauss = sum( N[i]*(vel[i]-mesh_vel[i]), i=0, number_of_points) (note that the fractional step vel is used)
    ms_vel_gauss_old[0] =  msN[0]*(fv0_old[0]-w0_old[0]) + msN[1]*(fv1_old[0]-w1_old[0]) +  msN[2]*(fv2_old[0]-w2_old[0]);
    ms_vel_gauss_old[1] =  msN[0]*(fv0_old[1]-w0_old[1]) + msN[1]*(fv1_old[1]-w1_old[1]) +  msN[2]*(fv2_old[1]-w2_old[1]);


    //ms_vel_gauss = v at (n+1)/2;
    ms_vel_gauss[0] += ms_vel_gauss_old[0];
    ms_vel_gauss[0] *= 0.5;
    ms_vel_gauss[1] += ms_vel_gauss_old[1];
    ms_vel_gauss[1] *= 0.5;

    //calculating viscosity
    double nu = 0.333333333333333333333333*(nu0 + nu1 + nu2 );
    double density = 0.3333333333333333333333*(rho0 + rho1 + rho2 );

    //DIAMETER of the element
    double dp = 0.3333333333333333333333*(dp0 + dp1 + dp2);
    //POROSITY of the element: average value of the porosity
    double eps = 0.3333333333333333333333*(eps0 + eps1 + eps2 );

    //1/PERMEABILITY of the element: average value of the porosity
    double kinv = 0.0;

    //Calculate the elemental Kinv in function of the nodal K of each element.

    //Version 1: we can calculate the elemental kinv from the nodal kinv;
    //THERE IS AN ERROR: IN THE INTERPHASE ELEMENTS A WATER NODE HAS TO BE ''MORE IMPORTANT'' THAN A POROUS ONE!!!
// 		if(kinv0 != 0.0 || kinv1 != 0.0 || kinv2 != 0.0) //if it is a fluid element
// 		{	double k0 = 0.0;
// 			double k1 = 0.0;
//.........这里部分代码省略.........

void TwoStepPeriodicCondition<TDim>::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
    rLeftHandSideMatrix.resize(0,0,false);
    rRightHandSideVector.resize(0,false);
}