C++ FloatArray::add方法代码示例

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
TEST(ArrayTest, ReserveAndAdd)
{
    FloatArray a;
    a.reserve(5);
    a.add(1.0f);
    a.add(3.3f);
    a.add(5.5f);

    ASSERT_EQ(3u, a.size());
    ASSERT_TRUE(a.capacity() >= 5);
    ASSERT_EQ(1.0f, a[0]);
    ASSERT_EQ(3.3f, a[1]);
    ASSERT_EQ(5.5f, a[2]);

    // To test reuse of buffer
    float* before = a.ptr();

    a.reserve(3);
    ASSERT_TRUE(a.capacity() >= 5);

    float* after = a.ptr();

    // Check that no realloc has been done
    ASSERT_EQ(before, after);
    ASSERT_TRUE(a.capacity() >= 5);
}

void Tr21Stokes :: computeLoadVector(FloatArray &answer, TimeStep *tStep)
{
    int i, load_number, load_id;
    Load *load;
    bcGeomType ltype;
    FloatArray vec;

    int nLoads = this->boundaryLoadArray.giveSize() / 2;
    answer.resize(15);
    answer.zero();
    for ( i = 1; i <= nLoads; i++ ) {  // For each Neumann boundary condition
        load_number = this->boundaryLoadArray.at(2 * i - 1);
        load_id = this->boundaryLoadArray.at(2 * i);
        load = this->domain->giveLoad(load_number);
        ltype = load->giveBCGeoType();

        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeBCSubVectorAt(vec, load, load_id, tStep);
            answer.add(vec);
        }
    }

    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( i = 1; i <= nLoads; i++ ) {
        load  = domain->giveLoad( bodyLoadArray.at(i) );
        ltype = load->giveBCGeoType();
        if ( ltype == BodyLoadBGT && load->giveBCValType() == ForceLoadBVT ) {
            this->computeBodyLoadVectorAt(vec, load, tStep);
            answer.add(vec);
        }
    }
}

void
CoupledFieldsElement :: giveInternalForcesVectorGen(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord, 
    void (*Nfunc)(GaussPoint*, FloatMatrix), void (*Bfunc)(GaussPoint*, FloatMatrix, int, int), //(GaussPoint*, FloatMatrix)
    void (*NStressFunc)(GaussPoint*, FloatArray), void (*BStressFunc)(GaussPoint*, FloatArray),
    double (*dVFunc)(GaussPoint*))
{
    // General implementation of internal forces that computes
    // f = sum_gp( N^T*GenStress_N + B^T*GenStress_B ) * dV

    FloatArray NStress, BStress, vGenStress, NS, BS;
    FloatMatrix N, B;

    for ( int j = 0; j < this->giveNumberOfIntegrationRules(); j++ ) {
        for ( auto &gp: this->giveIntegrationRule(j) ) {
            double dV  = this->computeVolumeAround(gp);
            
            // compute generalized stress measures
            if ( NStressFunc && Nfunc ) {
                Nfunc(gp, N);
                NStressFunc(gp, NStress);
                NS.beTProductOf(N, NStress);
                answer.add(dV, NS);
            }

            if ( BStressFunc && Bfunc ) {
                Bfunc(gp, B, 1, 3);
                BStressFunc(gp, BStress);
                BS.beTProductOf(B, BStress);
                answer.add(dV, BS);
            }

            
        }
    }
}

void Tr21Stokes :: computeExternalForcesVector(FloatArray &answer, TimeStep *tStep)
{
    FloatArray vec;

    answer.clear();

    int nLoads = this->boundaryLoadArray.giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {  // For each Neumann boundary condition
        int load_number = this->boundaryLoadArray.at(2 * i - 1);
        int load_id = this->boundaryLoadArray.at(2 * i);
        Load *load = this->domain->giveLoad(load_number);
        bcGeomType ltype = load->giveBCGeoType();

        if ( ltype == EdgeLoadBGT ) {
            this->computeBoundarySurfaceLoadVector(vec, static_cast< BoundaryLoad * >(load), load_id, ExternalForcesVector, VM_Total, tStep);
            answer.add(vec);
        }
    }

    BodyLoad *bload;
    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load = domain->giveLoad( bodyLoadArray.at(i) );
	if ((bload = dynamic_cast<BodyLoad*>(load))) {
	  bcGeomType ltype = load->giveBCGeoType();
	  if ( ltype == BodyLoadBGT && load->giveBCValType() == ForceLoadBVT ) {
            this->computeLoadVector(vec, bload, ExternalForcesVector, VM_Total, tStep);
            answer.add(vec);
	  }
        }
    }
}

int
TR_SHELL01 :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *atTime)
{
    if ( type == IST_ShellForceMomentumTensor ) {
        FloatArray aux;
        GaussPoint *membraneGP = membrane->giveDefaultIntegrationRulePtr()->getIntegrationPoint(gp->giveNumber()-1);
        GaussPoint *plateGP = plate->giveDefaultIntegrationRulePtr()->getIntegrationPoint(gp->giveNumber()-1);
        
        plate->giveIPValue(answer, plateGP, IST_ShellForceMomentumTensor, atTime);
        membrane->giveIPValue(aux, membraneGP, IST_ShellForceMomentumTensor, atTime);
        answer.add(aux);
        return 1;
    } else if ( type == IST_ShellStrainCurvatureTensor ) {
        FloatArray aux;
        GaussPoint *membraneGP = membrane->giveDefaultIntegrationRulePtr()->getIntegrationPoint(gp->giveNumber()-1);
        GaussPoint *plateGP = plate->giveDefaultIntegrationRulePtr()->getIntegrationPoint(gp->giveNumber()-1);

        plate->giveIPValue(answer, plateGP, IST_ShellStrainCurvatureTensor, atTime);
        membrane->giveIPValue(aux, membraneGP, IST_ShellStrainCurvatureTensor, atTime);
        answer.add(aux);

        return 1;
    } else {
        return StructuralElement::giveIPValue(answer, gp, type, atTime);
    }
}

void
SUPGElement2 :: computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
{
    int nLoads;

    answer.clear();

    int rule = 0;
    IntegrationRule *iRule = this->integrationRulesArray [ rule ];
    FloatArray un, gVector, s, helpLoadVector;
    FloatMatrix b, nu;

    // add body load (gravity) termms
    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load = domain->giveLoad( bodyLoadArray.at(i) );
        bcGeomType ltype = load->giveBCGeoType();
        if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
            load->computeComponentArrayAt(gVector, tStep, VM_Total);
            if ( gVector.giveSize() ) {
                for ( GaussPoint *gp: *iRule ) {
                    this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
                    this->computeNuMatrix(nu, gp);
                    double dV  = this->computeVolumeAround(gp);
                    double rho = this->giveMaterial()->give('d', gp);
                    answer.plusProduct(b, gVector, t_supg * rho * dV);
                    answer.plusProduct(nu, gVector, rho * dV);
                }
            }
        }
    }

    // integrate tractions
    // if no traction bc applied but side marked as with traction load
    // then zero traction is assumed !!!

    // loop over boundary load array
    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {
        int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
        int id = boundaryLoadArray.at(i * 2);
        Load *load  = domain->giveLoad(n);
        bcGeomType ltype = load->giveBCGeoType();
        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else if ( ltype == SurfaceLoadBGT ) {
            this->computeSurfaceLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else {
            OOFEM_ERROR("unsupported load type class");
        }
    }
}

void
SUPGElement2 :: computeBCRhsTerm_MC(FloatArray &answer, TimeStep *tStep)
{
    int nLoads;
    FloatArray s, gVector, helpLoadVector;
    FloatMatrix g;

    int rule = 1;

    answer.clear();

    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load  = domain->giveLoad( bodyLoadArray.at(i) );
        bcGeomType ltype = load->giveBCGeoType();
        if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
            load->computeComponentArrayAt(gVector, tStep, VM_Total);
            if ( gVector.giveSize() ) {
                for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
                    this->computeGradPMatrix(g, gp);
                    double dV = this->computeVolumeAround(gp);
                    answer.plusProduct(g, gVector, t_pspg * dV);
                }
            }
        }
    }

    // integrate tractions
    // if no traction bc applied but side marked as with traction load
    // then zero traction is assumed !!!

    // loop over boundary load array
    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {
        int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
        int id = boundaryLoadArray.at(i * 2);
        Load *load = domain->giveLoad(n);
        bcGeomType ltype = load->giveBCGeoType();
        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeLoadVector_MC(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else if ( ltype == SurfaceLoadBGT ) {
            this->computeSurfaceLoadVector_MC(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else {
            OOFEM_ERROR("unsupported load type class");
        }
    }
}

void
FEI3dTrQuad :: surfaceEvalBaseVectorsAt(FloatArray &G1, FloatArray &G2, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    // Note: These are not normalized. Returns the two tangent vectors to the surface.
    FloatMatrix dNdxi;
    this->surfaceEvaldNdxi(dNdxi, lcoords);

    G1.clear();
    G2.clear();
    for ( int i = 0; i < 6; ++i ) {
        G1.add( dNdxi(i, 1), * cellgeo.giveVertexCoordinates(i) );
        G2.add( dNdxi(i, 2), * cellgeo.giveVertexCoordinates(i) );
    }
}

void
GradDpElement :: giveLocalInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    NLStructuralElement *elem = this->giveNLStructuralElement();
    int nlGeo = elem->giveGeometryMode();
    FloatArray BS, vStress;
    FloatMatrix B;
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {

        if ( nlGeo == 0 || elem->domain->giveEngngModel()->giveFormulation() == AL ) {
            elem->computeBmatrixAt(gp, B);
        } else if ( nlGeo == 1 ) {
            elem->computeBHmatrixAt(gp, B);
        }
        vStress = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() )->giveTempStressVector();

        if ( vStress.giveSize() == 0 ) {         /// @todo is this check really necessary?
            break;
        }

        // Compute nodal internal forces at nodes as f = B^T*Stress dV
        double dV  = elem->computeVolumeAround(gp);
        BS.beTProductOf(B, vStress);
        answer.add(dV, BS);
    }
}

void
GradDpElement :: giveNonlocalInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    double localCumulatedStrain = 0.;
    NLStructuralElement *elem = this->giveNLStructuralElement();
    FloatMatrix stiffKappa;
    FloatArray Nk;
    FloatArray aux, dKappa, stress;

    int size = nSecVars * nSecNodes;

    //set displacement and nonlocal location array
    this->setDisplacementLocationArray();
    this->setNonlocalLocationArray();

    answer.resize(size);
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {
        this->computeNkappaMatrixAt(gp, Nk);

        double dV = elem->computeVolumeAround(gp);
        this->computeStressVectorAndLocalCumulatedStrain(stress, localCumulatedStrain, gp, tStep);
        aux.add(-dV * localCumulatedStrain, Nk);
    }

    this->computeStiffnessMatrix_kk(stiffKappa, TangentStiffness, tStep);
    this->computeNonlocalDegreesOfFreedom(dKappa, tStep);
    answer.beProductOf(stiffKappa, dKappa);
    answer.add(aux);
}

void PolygonLine :: giveGlobalCoordinates(FloatArray &oGlobalCoord, const double &iArcPos) const
{
    double L = computeLength();
    double xSegStart = 0.0, xSegEnd = 0.0;
    double xiSegStart = 0.0, xiSegEnd = 0.0;
    size_t numSeg = mVertices.size() - 1;
    const double xiTol = 1.0e-9;
    if ( iArcPos < xiTol ) {
        oGlobalCoord = mVertices [ 0 ];
        return;
    }

    for ( size_t i = 0; i < numSeg; i++ ) {
        xSegEnd += mVertices [ i ].distance(mVertices [ i + 1 ]);

        xiSegStart = xSegStart / L;
        xiSegEnd        = xSegEnd / L;

        if ( iArcPos > xiSegStart-xiTol && iArcPos < xiSegEnd+xiTol ) {
            // Point is within the segment
            FloatArray p;
            double elXi = ( iArcPos - xiSegStart ) / ( xiSegEnd - xiSegStart );
            p.beScaled( ( 1.0 - elXi ), mVertices [ i ] );
            p.add(elXi, mVertices [ i + 1 ]);
            oGlobalCoord = p;
            return;
        }

    }
}

void IntermediateConvectionDiffusionAssembler :: vectorFromElement(FloatArray &vec, Element &element, TimeStep *tStep, ValueModeType mode) const
{
    FloatArray vec_p;
    static_cast< CBSElement & >( element ).computeConvectionTermsI(vec, tStep);
    static_cast< CBSElement & >( element ).computeDiffusionTermsI(vec_p, tStep);
    vec.add(vec_p);
}

void SimpleVitrificationMaterial :: giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp,
                                                            const FloatArray &reducedStrain, TimeStep *tStep)
{
    FloatArray strainVector;
    FloatMatrix d;
    FloatArray deltaStrain;

    StructuralMaterialStatus *status = dynamic_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );

    this->giveStressDependentPartOfStrainVector(strainVector, gp, reducedStrain, tStep, VM_Total);

    deltaStrain.beDifferenceOf( strainVector, status->giveStrainVector() );

    this->give3dMaterialStiffnessMatrix(d, TangentStiffness, gp, tStep);

    FloatArray deltaStress;
    deltaStress.beProductOf(d, deltaStrain);

    answer = status->giveStressVector();
    answer.add(deltaStress);

    // update gp
    status->letTempStrainVectorBe(reducedStrain);
    status->letTempStressVectorBe(answer);
}

void
TrabBoneNL3D :: giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp,
                                        const FloatArray &totalStrain, TimeStep *tStep)
{
    TrabBoneNL3DStatus *nlStatus = static_cast< TrabBoneNL3DStatus * >( this->giveStatus(gp) );

    this->initTempStatus(gp);

    double tempDam;
    FloatArray effStress, totalStress, densStress;

    performPlasticityReturn(gp, totalStrain, tStep);
    tempDam = computeDamage(gp, tStep);
    effStress = nlStatus->giveTempEffectiveStress();

    totalStress = ( 1 - tempDam ) * effStress;

    for ( int i = 1; i <= 6; i++ ) {
        if ( sqrt( totalStress.at(i) * totalStress.at(i) ) < 1e-8 ) {
            totalStress.at(i) = 0.;
        }
    }

    computePlasStrainEnerDensity(gp, totalStrain, totalStress);

    if ( densCrit != 0. ) {
        computeDensificationStress(densStress, gp, totalStrain, tStep);
        answer.add(densStress);
    }

    answer = totalStress;
    nlStatus->setTempDam(tempDam);
    nlStatus->letTempStrainVectorBe(totalStrain);
    nlStatus->letTempStressVectorBe(answer);
}

void Tr1Darcy :: computeInternalForcesVector(FloatArray &answer, TimeStep *atTime)
{
    FloatArray *lcoords, w, a, gradP, I;
    FloatMatrix BT;
    GaussPoint *gp;

    TransportMaterial *mat = ( TransportMaterial * ) this->domain->giveMaterial(this->material);
    IntegrationRule *iRule = integrationRulesArray [ 0 ];

    this->computeVectorOf(EID_ConservationEquation, VM_Total, atTime, a);

    answer.resize(3);
    answer.zero();

    for ( int i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        lcoords = gp->giveCoordinates();

        double detJ = this->interpolation_lin.giveTransformationJacobian( * lcoords, FEIElementGeometryWrapper(this) );
        this->interpolation_lin.evaldNdx( BT, * lcoords, FEIElementGeometryWrapper(this) );

        gradP.beTProductOf(BT, a);

        mat->giveFluxVector(w, gp, gradP, atTime);

        I.beProductOf(BT, w);
        answer.add(- gp->giveWeight() * detJ, I);
    }
}