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

void
NonLinearDynamic :: assembleIncrementalReferenceLoadVectors(FloatArray &_incrementalLoadVector,
                                                            FloatArray &_incrementalLoadVectorOfPrescribed,
                                                            SparseNonLinearSystemNM :: referenceLoadInputModeType _refMode,
                                                            Domain *sourceDomain, EquationID ut, TimeStep *tStep)
{
    EModelDefaultEquationNumbering en;

    _incrementalLoadVector.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations(sourceDomain->giveNumber(), EModelDefaultEquationNumbering()) );
    _incrementalLoadVector.zero();
    _incrementalLoadVectorOfPrescribed.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations(sourceDomain->giveNumber(), EModelDefaultPrescribedEquationNumbering()) );
    _incrementalLoadVectorOfPrescribed.zero();

    if ( _refMode == SparseNonLinearSystemNM :: rlm_incremental ) {
        this->assembleVector(_incrementalLoadVector, tStep, ut, ExternalForcesVector,
                             VM_Incremental, EModelDefaultEquationNumbering(), sourceDomain);

        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ut, ExternalForcesVector,
                             VM_Incremental, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    } else {
        this->assembleVector(_incrementalLoadVector, tStep, ut, ExternalForcesVector,
                             VM_Total, EModelDefaultEquationNumbering(), sourceDomain);
        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ut, ExternalForcesVector,
                             VM_Total, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    }

#ifdef __PARALLEL_MODE
    this->updateSharedDofManagers(_incrementalLoadVector, LoadExchangeTag);
#endif
}

void
NonLinearStatic :: assembleIncrementalReferenceLoadVectors(FloatArray &_incrementalLoadVector,
                                                           FloatArray &_incrementalLoadVectorOfPrescribed,
                                                           SparseNonLinearSystemNM :: referenceLoadInputModeType _refMode,
                                                           Domain *sourceDomain, TimeStep *tStep)
{
    _incrementalLoadVector.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations( sourceDomain->giveNumber(), EModelDefaultEquationNumbering() ) );
    _incrementalLoadVector.zero();
    _incrementalLoadVectorOfPrescribed.resize( sourceDomain->giveEngngModel()->giveNumberOfDomainEquations( sourceDomain->giveNumber(), EModelDefaultPrescribedEquationNumbering() ) );
    _incrementalLoadVectorOfPrescribed.zero();

    if ( _refMode == SparseNonLinearSystemNM :: rlm_incremental ) {
        ///@todo This was almost definitely wrong before. It never seems to be used. Is this code even relevant?
        this->assembleVector(_incrementalLoadVector, tStep, ExternalForceAssembler(),
                             VM_Incremental, EModelDefaultEquationNumbering(), sourceDomain);

        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ExternalForceAssembler(),
                             VM_Incremental, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    } else {
        this->assembleVector(_incrementalLoadVector, tStep, ExternalForceAssembler(),
                             VM_Total, EModelDefaultEquationNumbering(), sourceDomain);

        this->assembleVector(_incrementalLoadVectorOfPrescribed, tStep, ExternalForceAssembler(),
                             VM_Total, EModelDefaultPrescribedEquationNumbering(), sourceDomain);
    }

    this->updateSharedDofManagers(_incrementalLoadVector, EModelDefaultEquationNumbering(), LoadExchangeTag);
}

int
IntMatIsoDamage :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    IntMatIsoDamageStatus *status = static_cast< IntMatIsoDamageStatus * >( this->giveStatus(gp) );
    if ( type == IST_DamageTensor ) {
        answer.resize(6);
        answer.zero();
        answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
        return 1;
    } else if ( type == IST_DamageTensorTemp ) {
        answer.resize(6);
        answer.zero();
        answer.at(1) = answer.at(2) = answer.at(3) = status->giveTempDamage();
        return 1;
    } else if ( type == IST_PrincipalDamageTensor ) {
        answer.resize(3);
        answer.at(1) = answer.at(2) = answer.at(3) = status->giveDamage();
        return 1;
    } else if ( type == IST_PrincipalDamageTempTensor ) {
        answer.resize(3);
        answer.at(1) = answer.at(2) = answer.at(3) = status->giveTempDamage();
        return 1;
    } else if ( type == IST_MaxEquivalentStrainLevel ) {
        answer.resize(1);
        answer.at(1) = status->giveKappa();
        return 1;
    } else {
        return StructuralInterfaceMaterial :: giveIPValue(answer, gp, type, tStep);
    }
}

void
M1Material :: giveRealStressVector_3d(FloatArray &answer,
                                      GaussPoint *gp,
                                      const FloatArray &totalStrain,
                                      TimeStep *tStep)
{
    answer.resize(6);
    answer.zero();

    // get the status at the beginning
    M1MaterialStatus *status = static_cast< M1MaterialStatus * >( this->giveStatus(gp) );
    // prepare status at the end
    this->initTempStatus(gp);
    // get the initial values of plastic strains on microplanes (set to zero in the first step)
    FloatArray epspN = status->giveNormalMplanePlasticStrains();
    if ( epspN.giveSize() < numberOfMicroplanes ) {
        epspN.resize(numberOfMicroplanes);
        epspN.zero();
    }

    // loop over microplanes
    FloatArray sigN(numberOfMicroplanes);
    IntArray plState(numberOfMicroplanes);
    for ( int imp = 1; imp <= numberOfMicroplanes; imp++ ) {
        Microplane *mPlane = this->giveMicroplane(imp - 1, gp);
        //IntegrationPointStatus *mPlaneStatus =  this->giveMicroplaneStatus(mPlane);
        double epsN = computeNormalStrainComponent(mPlane, totalStrain);
        // evaluate trial stress on the microplane
        double sigTrial = EN * ( epsN - epspN.at(imp) );
        // evaluate the yield stress (from total microplane strain, not from its plastic part)
        double sigYield = EN * ( s0 + HN * epsN ) / ( EN + HN );
        if ( sigYield < 0. ) {
            sigYield = 0.;
        }
        // check whether the yield stress is exceeded and set the microplane stress
        if ( sigTrial > sigYield ) { //plastic yielding
            sigN.at(imp) = sigYield;
            epspN.at(imp) = epsN - sigYield / EN;
            plState.at(imp) = 1;
        } else {
            sigN.at(imp) = sigTrial;
            plState.at(imp) = 0;
        }
        // add the contribution of the microplane to macroscopic stresses
        for ( int i = 1; i <= 6; i++ ) {
            answer.at(i) += N [ imp - 1 ] [ i - 1 ] * sigN.at(imp) * microplaneWeights [ imp - 1 ];
        }
    }
    // multiply the integral over unit hemisphere by 6
    answer.times(6);

    // update status
    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(answer);
    status->letTempNormalMplaneStressesBe(sigN);
    status->letTempNormalMplanePlasticStrainsBe(epspN);
    status->letPlasticStateIndicatorsBe(plState);
}

void
SimpleCrossSection :: giveReducedCharacteristicVector(FloatArray &answer, GaussPoint *gp,
                                                      const FloatArray &charVector3d)
//
// returns reduced stressVector or strainVector from full 3d vector reduced
// to vector required by gp->giveStressStrainMode()
//
{
    MaterialMode mode = gp->giveMaterialMode();
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( gp->giveElement()->giveMaterial() );
    IntArray indx;
    int size = charVector3d.giveSize();
    int i, j;

    if ( ( mode == _3dShell ) || ( mode == _3dBeam ) || ( mode == _2dPlate ) || ( mode == _2dBeam ) ) {
        if ( size != 12 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        mat->giveStressStrainMask( indx, ReducedForm, gp->giveMaterialMode() );
        answer.resize( indx.giveSize() );
        answer.zero();

        for ( i = 1; i <= indx.giveSize(); i++ ) {
            if ( ( j = indx.at(i) ) ) {
                answer.at(i) = charVector3d.at(j);
            }
        }

        return;
    } else if ( mode == _3dBeam ) {
        if ( size != 6 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        answer =  charVector3d;
        return;
    } else if ( mode == _PlaneStressRot ) {
        if ( size != 7 ) {
            OOFEM_ERROR("SimpleCrossSection :: giveReducedCharacteristicVector - charVector3d size mismatch");
        }

        mat->giveStressStrainMask( indx, ReducedForm, gp->giveMaterialMode() );
        answer.resize( indx.giveSize() );
        answer.zero();

        for ( i = 1; i <= indx.giveSize(); i++ ) {
            if ( ( j = indx.at(i) ) ) {
                answer.at(i) = charVector3d.at(j);
            }
        }

        return;
    } else {
        StructuralCrossSection :: giveReducedCharacteristicVector(answer, gp, charVector3d);
        return;
    }
}

int
FEI3dHexaQuad :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)
{
    FloatArray res, delta, guess;
    FloatMatrix jac;
    double convergence_limit, error;

    // find a suitable convergence limit
    convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);

    // setup initial guess
    answer.resize(gcoords.giveSize());
    answer.zero();

    // apply Newton-Raphson to solve the problem
    for (int nite = 0; nite < 10; nite++) {
        // compute the residual
        this->local2global(guess, answer, cellgeo);
        res.beDifferenceOf(gcoords, guess);

        // check for convergence
        error = res.computeNorm();
        if ( error < convergence_limit ) {
            break;
        }

        // compute the corrections
        this->giveJacobianMatrixAt(jac, answer, cellgeo);
        jac.solveForRhs(res, delta, true);

        // update guess
        answer.add(delta);
    }
    if ( error > convergence_limit) { // Imperfect, could give false negatives.
        //OOFEM_ERROR ("global2local: no convergence after 10 iterations");
        answer.zero();
        return false;
    }

    // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
    bool inside = true;
    for ( int i = 1; i <= answer.giveSize(); i++ ) {
        if ( answer.at(i) < ( -1. - POINT_TOL ) ) {
            answer.at(i) = -1.;
            inside = false;
        } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
            answer.at(i) = 1.;
            inside = false;
        }
    }

    return inside;
}

void
StructuralEngngModel :: computeExternalLoadReactionContribution(FloatArray &reactions, TimeStep *tStep, int di)
{
    int numRestrDofs = this->giveNumberOfPrescribedDomainEquations(di, EID_MomentumBalance);
    reactions.resize(numRestrDofs);
    reactions.zero();
    FloatArray contribution(numRestrDofs);

    reactions.resize( this->giveNumberOfPrescribedDomainEquations(di, EID_MomentumBalance) );
    reactions.zero();
    this->assembleVector( reactions, tStep, EID_MomentumBalance, ExternalForcesVector, VM_Total,
                          EModelDefaultPrescribedEquationNumbering(), this->giveDomain(di) );
}

void RampFunction :: evaluateEnrFuncDerivAt(FloatArray &oEnrFuncDeriv, const FloatArray &iPos, const double &iLevelSet, const FloatArray &iGradLevelSet) const
{
    oEnrFuncDeriv.resize(2);
    oEnrFuncDeriv.zero();
    oEnrFuncDeriv.at(1) = iGradLevelSet.at(1) * sgn(iLevelSet);
    oEnrFuncDeriv.at(2) = iGradLevelSet.at(2) * sgn(iLevelSet);
}

void Tr21Stokes :: computeBodyLoadVectorAt(FloatArray &answer, Load *load, TimeStep *tStep)
{
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    GaussPoint *gp;
    FloatArray N, gVector, *lcoords, temparray(15);
    double dA, detJ, rho;

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( int k = 0; k < iRule->getNumberOfIntegrationPoints(); k++ ) {
            gp = iRule->getIntegrationPoint(k);
            lcoords = gp->giveCoordinates();

            rho = this->giveMaterial()->giveCharacteristicValue(MRM_Density, gp, tStep);
            detJ = fabs( this->interpolation_quad.giveTransformationJacobian(* lcoords, FEIElementGeometryWrapper(this)) );
            dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN(N, * lcoords, FEIElementGeometryWrapper(this));
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble( temparray, this->ordering );
}

void
Node2NodeContactL :: computeContactForces(FloatArray &answer, TimeStep *tStep, CharType type, ValueModeType mode,
                                const UnknownNumberingScheme &s, Domain *domain, FloatArray *eNorms)
{
  
    //Loop through all the master objects and let them do their thing
    FloatArray gap, C, Fc;
    this->computeGap(gap, tStep);
    
    GaussPoint *gp = this->integrationRule->getIntegrationPoint(0);
    FloatArray t;
    this->computeContactTractionAt(gp, t ,gap, tStep);
    this->computeCmatrixAt(gp, C, tStep);
    
    // compute load vector
    // for Lagrange: fc = traction * C^T * A (traction = lambda)
    FloatArray temp = t.at(1) *this->area * C;
    
    answer.resize( C.giveSize() + 1);
    answer.zero();
    if( gap.at(1) < 0.0 ) {
        answer.addSubVector(temp,1);
        answer.at( C.giveSize() + 1 ) = -gap.at(1);
    }
}

void
SolutionbasedShapeFunction :: setLoads(EngngModel *myEngngModel, int d)
{
    DynamicInputRecord ir;
    FloatArray gradP;

    gradP.resize( this->giveDomain()->giveNumberOfSpatialDimensions() );
    gradP.zero();
    gradP.at(d) = 1.0;

    ir.setRecordKeywordField("deadweight", 1);
    ir.setField(gradP, _IFT_Load_components);
    ir.setField(1, _IFT_GeneralBoundaryCondition_timeFunct);

    int bcID = myEngngModel->giveDomain(1)->giveNumberOfBoundaryConditions() + 1;
    GeneralBoundaryCondition *myBodyLoad;
    myBodyLoad = classFactory.createBoundaryCondition( "deadweight", bcID, myEngngModel->giveDomain(1) );
    myBodyLoad->initializeFrom(& ir);
    myEngngModel->giveDomain(1)->setBoundaryCondition(bcID, myBodyLoad);

    for ( int i = 1; i <= myEngngModel->giveDomain(1)->giveNumberOfElements(); i++ ) {
        IntArray *blArray;
        blArray = myEngngModel->giveDomain(1)->giveElement(i)->giveBodyLoadArray();
        blArray->resizeWithValues(blArray->giveSize() + 1);
        blArray->at( blArray->giveSize() ) = bcID;
    }
}

void
MisesMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                      FloatArray &rcontrib, TimeStep *tStep)
{
    double kappa, tempKappa;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp);

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);
    kappa = status->giveCumulativePlasticStrain();
    tempKappa = status->giveTempCumulativePlasticStrain();

    rcontrib.clear();
    if ( ( tempKappa - kappa ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        if ( gp->giveMaterialMode() == _1dMat ) {
            double coeff = sgn( stress.at(1) ) * E / ( E + H );
            rcontrib.plusProduct(b, stress, coeff);
            return;
        }
    }
    rcontrib.resize(b.giveNumberOfColumns());
    rcontrib.zero();
}

void Tr21Stokes :: computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
{
    if ( type != ExternalForcesVector ) {
        answer.clear();
        return;
    }

    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
    FloatArray N, gVector, temparray(12);

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();

            double rho = mat->give('d', gp);
            double detJ = fabs( this->interpolation_quad.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) );
            double dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN( N, lcoords, FEIElementGeometryWrapper(this) );
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble(temparray, this->momentum_ordering);
}

void
ConstantEdgeLoad :: computeValueAt(FloatArray &answer, TimeStep *stepN, FloatArray &coords, ValueModeType mode)
{
    // we overload general implementation on the boundary load level due
    // to implementation efficiency

    double factor;

    if ( ( mode != VM_Total ) && ( mode != VM_Incremental ) ) {
        _error("computeValueAt: mode not supported");
    }

    // ask time distribution
    /*
     * factor = this -> giveLoadTimeFunction() -> at(stepN->giveTime()) ;
     * if ((mode==VM_Incremental) && (!stepN->isTheFirstStep()))
     * //factor -= this->giveLoadTimeFunction()->at(stepN->givePreviousStep()->giveTime()) ;
     * factor -= this->giveLoadTimeFunction()->at(stepN->giveTime()-stepN->giveTimeIncrement()) ;
     */
    factor = this->giveLoadTimeFunction()->evaluate(stepN, mode);
    answer = componentArray;
    answer.times(factor);
    
    if ( !isImposed(stepN) ){
        answer.zero();
    }
}

void DynCompCol :: timesT(const FloatArray &x, FloatArray &answer) const
{
    // Check for compatible dimensions:
    if ( x.giveSize() != nRows ) {
        OOFEM_ERROR("Error in CompCol -- incompatible dimensions");
    }
    answer.resize(nColumns);
    answer.zero();

#ifndef DynCompCol_USE_STL_SETS
    int i, t;
    double r;

    for ( i = 0; i < nColumns; i++ ) {
        r = 0.0;
        for ( t = 1; t <= columns_ [ i ]->giveSize(); t++ ) {
            r += columns_ [ i ]->at(t) * x( rowind_ [ i ]->at(t) );
        }

        answer(i) = r;
    }

#else
    double r;
    for ( int i = 0; i < nColumns; i++ ) {
        double r = 0.0;
        for ( auto &val: columns [ i ] ) {
            r += val.second * x(val.first);
        }

        answer(i) = r;
    }

#endif
}