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

void
MatlabExportModule :: doOutputHomogenizeDofIDs(TimeStep *tStep,    FILE *FID)
{

    std :: vector <FloatArray*> HomQuantities;
    double Vol = 0.0;

    // Initialize vector of arrays constaining homogenized quantities
    HomQuantities.resize(internalVarsToExport.giveSize());

    for (int j=0; j<internalVarsToExport.giveSize(); j++) {
        HomQuantities.at(j) = new FloatArray;
    }

    int nelem = this->elList.giveSize();
    for (int i = 1; i<=nelem; i++) {
        Element *e = this->emodel->giveDomain(1)->giveElement(elList.at(i));
        FEInterpolation *Interpolation = e->giveInterpolation();

        Vol = Vol + e->computeVolumeAreaOrLength();

        for ( GaussPoint *gp: *e->giveDefaultIntegrationRulePtr() ) {

            for (int j=0; j<internalVarsToExport.giveSize(); j++) {
                FloatArray elementValues;
                e->giveIPValue(elementValues, gp, (InternalStateType) internalVarsToExport(j), tStep);
                double detJ=fabs(Interpolation->giveTransformationJacobian( gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e)));

                elementValues.times(gp->giveWeight()*detJ);
                if (HomQuantities.at(j)->giveSize() == 0) {
                    HomQuantities.at(j)->resize(elementValues.giveSize());
                    HomQuantities.at(j)->zero();
                };
                HomQuantities.at(j)->add(elementValues);
            }
        }
    }


    if (noscaling) Vol=1.0;

    for ( std :: size_t i = 0; i < HomQuantities.size(); i ++) {
        FloatArray *thisIS;
        thisIS = HomQuantities.at(i);
        thisIS->times(1.0/Vol);
        fprintf(FID, "\tspecials.%s = [", __InternalStateTypeToString ( InternalStateType (internalVarsToExport(i)) ) );

        for (int j = 0; j<thisIS->giveSize(); j++) {
            fprintf(FID, "%e", thisIS->at(j+1));
            if (j!=(thisIS->giveSize()-1) ) {
                fprintf(FID, ", ");
            }
        }
        fprintf(FID, "];\n");
        delete HomQuantities.at(i);
    }

}

void
MisesMat :: giveRealStressVector_1d(FloatArray &answer,
                                 GaussPoint *gp,
                                 const FloatArray &totalStrain,
                                 TimeStep *tStep)
{
    /// @note: One should obtain the same answer using the iterations in the default implementation (this is verified for this model).
#if 1
    MisesMatStatus *status = static_cast< MisesMatStatus * >( this->giveStatus(gp) );

    this->performPlasticityReturn(gp, totalStrain);
    double omega = computeDamage(gp, tStep);
    answer = status->giveTempEffectiveStress();
    answer.times(1 - omega);

    // Compute the other components of the strain:
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp), nu = lmat->give('n', gp);

    FloatArray strain = status->getTempPlasticStrain();
    strain[0] = totalStrain[0];
    strain[1] -= nu / E * status->giveTempEffectiveStress()[0];
    strain[2] -= nu / E * status->giveTempEffectiveStress()[0];

    status->letTempStrainVectorBe(strain);
    status->setTempDamage(omega);
    status->letTempStressVectorBe(answer);
#else
    StructuralMaterial :: giveRealStressVector_1d(answer, gp, totalStrain, tStep);
#endif
}

void FE2FluidMaterial :: giveDeviatoricPressureStiffness(FloatArray &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    FE2FluidMaterialStatus *ms = static_cast<FE2FluidMaterialStatus*> (this->giveStatus(gp));
    ms->computeTangents(tStep);
    if ( mode == TangentStiffness ) {
        answer = ms->giveDeviatoricPressureTangent();
#ifdef DEBUG_TANGENT
        // Numerical ATS for debugging
        FloatArray strain(3); strain.zero();
        FloatArray sig, sigh;
        double epspvol, pressure = 0.0;
        double h = 1.00; // Linear problem, size of this doesn't matter.
        computeDeviatoricStressVector (sig, epspvol, gp, strain, pressure, tStep);
        computeDeviatoricStressVector (sigh, epspvol, gp, strain, pressure+h, tStep);

        FloatArray dsigh; dsigh.beDifferenceOf(sigh,sig); dsigh.times(1/h);

        printf("Analytical deviatoric pressure tangent = "); answer.printYourself();
        printf("Numerical deviatoric pressure tangent = "); dsigh.printYourself();
        dsigh.subtract(answer);
        double norm = dsigh.computeNorm();
        if (norm > answer.computeNorm()*DEBUG_ERR && norm > 0.0) {
            OOFEM_ERROR("Error in deviatoric pressure tangent");
        }
#endif
    } else {
        OOFEM_ERROR("Mode not implemented");
    }
}

void
NonlinearFluidMaterial :: computeDeviatoricStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &eps, TimeStep *tStep)
{
    NonlinearFluidMaterialStatus *status = static_cast< NonlinearFluidMaterialStatus * >( this->giveStatus(gp) );

    double normeps2;

    answer = eps;
    if ( eps.giveSize() == 3 ) {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + 0.5 * ( eps.at(3) * eps.at(3) );
        answer.at(3) *= 0.5;
    } else if ( eps.giveSize() == 4 ) {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + eps.at(3) * eps.at(3) + 0.5 * ( eps.at(4) * eps.at(4) );
        answer.at(4) *= 0.5;
    } else {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + eps.at(3) * eps.at(3) + 0.5 * ( eps.at(4) * eps.at(4) + eps.at(5) * eps.at(5) +  eps.at(6) * eps.at(6) );
        answer.at(4) *= 0.5;
        answer.at(5) *= 0.5;
        answer.at(6) *= 0.5;
    }

    answer.times( 2.0 * viscosity * ( 1.0 + c * pow(normeps2, alpha * 0.5) ) );

    status->letTempDeviatoricStressVectorBe(answer);
    status->letTempDeviatoricStrainVectorBe(eps);
    status->letTempStrainNorm2Be(normeps2);
}

void
Node2NodeContactL :: computeContactTangent(FloatMatrix &answer, CharType type, TimeStep *tStep)
{
    answer.resize(7,7);
    answer.zero();
    
    FloatArray gap;
    this->computeGap(gap, tStep);
    
    if( gap.at(1) < 0.0 ) {
      
        GaussPoint *gp = this->integrationRule->getIntegrationPoint(0);
        
        FloatArray C;
        this->computeCmatrixAt(gp, C, tStep);
        int sz = C.giveSize();
        C.times(this->area);
        
        answer.addSubVectorCol(C, 1, sz + 1);
        answer.addSubVectorRow(C, sz + 1, 1);
    }
    
    //TODO need to add a small number for the solver
    for ( int i = 1; i <= 7; i++ ) {
        answer.at(i,i) += 1.0e-8;
    }

}

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
PhaseFieldElement :: computeStiffnessMatrix_ud(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    FloatMatrix B, DB, N_d, DN, S(3,1);
    FloatArray stress;
    NLStructuralElement *el = this->giveElement();
    StructuralCrossSection *cs = dynamic_cast<StructuralCrossSection* > (el->giveCrossSection() );

    answer.clear();

    for ( auto &gp: el->giveIntegrationRule(0) ) {
        StructuralMaterialStatus *matStat = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() );

        double dV = el->computeVolumeAround(gp);
        // compute int_V ( B^t * D_B * B )dV

        this->computeNd_matrixAt(*gp->giveNaturalCoordinates(), N_d);

        // stress   
        FloatArray reducedStrain, a_u;
        FloatMatrix B_u;
        el->computeBmatrixAt( gp, B_u );
        stress = matStat->giveTempStressVector();
        stress.times( this->computeGPrim( gp, VM_Total, tStep ) );
        
        S.setColumn(stress,1);
        DN.beProductOf(S, N_d);
        answer.plusProductUnsym(B_u, DN, dV);
    }

}

void
LIBeam2d :: computeBodyLoadVectorAt(FloatArray &answer, Load *load, TimeStep *tStep, ValueModeType mode)
{
    FloatArray lc(1);
    StructuralElement :: computeBodyLoadVectorAt(answer, load, tStep, mode);
    answer.times( this->giveCrossSection()->give(CS_Area, lc, this) );
}

void 
IntMatCoulombContact :: computeEngTraction(double &normalStress, FloatArray &shearStress,
                                         FloatArray &tempShearStressShift, double normalJump, const FloatArray &shearJump )
{

    double maxShearStress = 0.0;
    double shift = -this->kn * this->stiffCoeff * normalClearance;

    if ( normalJump + normalClearance <= 0. ) {
        normalStress = this->kn * ( normalJump + normalClearance ) + shift; //in compression and after the clearance gap closed
        maxShearStress = fabs( normalStress ) * this->frictCoeff;
    } else {
        normalStress = this->kn * this->stiffCoeff * ( normalJump + normalClearance ) + shift;
        maxShearStress = 0.;
    }

    if ( shearJump.giveSize() != 0 ) {

        shearStress = this->kn * shearJump - tempShearStressShift;
        double dp = shearStress.computeNorm();
        double eps = 1.0e-15; // small number
        if ( dp > maxShearStress ) {
            shearStress.times( maxShearStress / ( dp + eps ) );
        }
        tempShearStressShift = this->kn * shearJump - shearStress;

    } else { // 1d -> no shear stresses
        return;
    }
}

void EnrichmentItem :: calcPolarCoord(double &oR, double &oTheta, const FloatArray &iOrigin, const FloatArray &iPos, const FloatArray &iN, const FloatArray &iT, const EfInput &iEfInput, bool iFlipTangent)
{
    FloatArray q = {
        iPos.at(1) - iOrigin.at(1), iPos.at(2) - iOrigin.at(2)
    };

    const double tol = 1.0e-20;

    // Compute polar coordinates
    oR = iOrigin.distance(iPos);

    if ( oR > tol ) {
        q.times(1.0 / oR);
    }

    const double pi = M_PI;

    //    if( q.dotProduct(iT) > 0.0 ) {
    //        oTheta = asin( q.dotProduct(iN) );
    //    } else {
    //        if ( q.dotProduct(iN) > 0.0 ) {
    //            oTheta = pi - asin( q.dotProduct(iN) );
    //        } else {
    //            oTheta = -pi - asin( q.dotProduct(iN) );
    //        }
    //    }


    const double tol_q = 1.0e-3;
    double phi = iEfInput.mLevelSet;

    if ( iFlipTangent ) {
        phi *= -1.0;
    }

    double phi_r = 0.0;
    if ( oR > tol ) {
        phi_r = fabs(phi / oR);
    }

    if ( phi_r > 1.0 - XfemTolerances :: giveApproxZero() ) {
        phi_r = 1.0 - XfemTolerances :: giveApproxZero();
    }

    if ( iEfInput.mArcPos < tol_q || iEfInput.mArcPos > ( 1.0 - tol_q ) ) {
        double q_dot_n = q.dotProduct(iN);
        if ( q_dot_n > 1.0 - XfemTolerances :: giveApproxZero() ) {
            q_dot_n = 1.0 - XfemTolerances :: giveApproxZero();
        }

        oTheta = asin(q_dot_n);
    } else {
        if ( phi > 0.0 ) {
            oTheta = pi - asin( fabs(phi_r) );
        } else {
            oTheta = -pi + asin( fabs(phi_r) );
        }
    }
}

int
IntElPoint :: computeGlobalCoordinates(FloatArray &answer, const FloatArray &lcoords)
{
    answer = *this->giveNode(1)->giveCoordinates();
    answer.add(*this->giveNode(2)->giveCoordinates());
    answer.times(0.5);
    return 1;
}

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

double LineSurfaceTension :: SpatialLocalizerI_giveDistanceFromParametricCenter(const FloatArray &coords)
{
    FloatArray c;
    c = *this->giveNode(1)->giveCoordinates();
    c.add(*this->giveNode(2)->giveCoordinates());
    c.times(0.5);
    c.subtract(coords);
    return c.computeNorm();
}

double
FEI3dTrQuad :: surfaceEvalNormal(FloatArray &answer, int isurf, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    FloatArray G1, G2; // local curvilinear base vectors
    this->surfaceEvalBaseVectorsAt(G1, G2, lcoords, cellgeo);
    answer.beVectorProductOf(G1, G2);
    double J = answer.computeNorm();
    answer.times(1 / J);
    return J;
}

void
NonLinearStatic :: computeExternalLoadReactionContribution(FloatArray &reactions, TimeStep *tStep, int di)
{
    if ( ( di == 1 ) && ( tStep == this->giveCurrentStep() ) ) {
        reactions = incrementalLoadVectorOfPrescribed;
        reactions.times(loadLevel);
        reactions.add(initialLoadVectorOfPrescribed);
    } else {
        _error("computeExternalLoadReactionContribution: unable to respond due to invalid solution step or domain");
    }
}