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

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
GradDpElement :: computeDeformationGradientVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    // Computes the deformation gradient in the Voigt format at the Gauss point gp of
    // the receiver at time step tStep.
    // Order of components: 11, 22, 33, 23, 13, 12, 32, 31, 21 in the 3D.

    // Obtain the current displacement vector of the element and subtract initial displacements (if present)
    FloatArray u;
    FloatMatrix B;
    NLStructuralElement *elem = this->giveNLStructuralElement();

    this->computeDisplacementDegreesOfFreedom(u, tStep);
    // Displacement gradient H = du/dX
    elem->computeBHmatrixAt(gp, B);
    answer.beProductOf(B, u);

    // Deformation gradient F = H + I
    MaterialMode matMode = gp->giveMaterialMode();
    if ( matMode == _3dMat || matMode == _PlaneStrain ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
        answer.at(3) += 1.0;
    } else if ( matMode == _PlaneStress ) {
        answer.at(1) += 1.0;
        answer.at(2) += 1.0;
    } else if ( matMode == _1dMat ) {
        answer.at(1) += 1.0;
    } else {
        OOFEM_ERROR( "MaterialMode is not supported yet (%s)", __MaterialModeToString(matMode) );
    }
}

void
SimpleCrossSection :: giveGeneralizedStress_Beam2d(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
  /**Note: (by bp): This assumes that the behaviour is elastic
     there exist a nuumber of nonlinear integral material models for beams defined directly in terms of integral forces and moments and corresponding 
     deformations and curvatures. This would require to implement support at material model level.
     Mikael: That would not be a continuum material model, but it would highly depend on the cross-section shape, thus, it should be a special cross-section model instead.
     This cross-section assumes you can split the response into inertia moments and pure material response. This is only possible for a constant, elastic response (i.e. elastic).
     Therefore, this cross-section may only be allowed to give the elastic response.
  */
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    FloatArray elasticStrain, et, e0;
    FloatMatrix tangent;
    elasticStrain = strain;
    this->giveTemperatureVector(et, gp, tStep);
    if ( et.giveSize() > 0 ) {
        mat->giveThermalDilatationVector(e0, gp, tStep);
        double thick = this->give(CS_Thickness, gp);
        elasticStrain.at(1) -= e0.at(1) * ( et.at(1) - mat->giveReferenceTemperature() );
        if ( et.giveSize() > 1 ) {
            elasticStrain.at(2) -= e0.at(1) * et.at(2) / thick;     // kappa_x
        }
    }
    this->give2dBeamStiffMtrx(tangent, ElasticStiffness, gp, tStep);
    answer.beProductOf(tangent, elasticStrain);
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

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
LIBeam3dNL :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    GaussPoint *gp = this->giveDefaultIntegrationRulePtr()->getIntegrationPoint(0);
    FloatArray nm(6), stress, strain;
    FloatMatrix x;
    double s1, s2;

    // update temp triad
    this->updateTempTriad(tStep);

    if ( useUpdatedGpRecord == 1 ) {
        stress = static_cast< StructuralMaterialStatus * >( gp->giveMaterialStatus() )->giveStrainVector();
    } else {
        this->computeStrainVector(strain, gp, tStep);
        this->computeStressVector(stress, strain, gp, tStep);
    }

    for ( int i = 1; i <= 3; i++ ) {
        s1 = s2 = 0.0;
        for ( int j = 1; j <= 3; j++ ) {
            s1 += tempTc.at(i, j) * stress.at(j);
            s2 += tempTc.at(i, j) * stress.at(j + 3);
        }

        nm.at(i)   = s1;
        nm.at(i + 3) = s2;
    }

    this->computeXMtrx(x, tStep);
    answer.beProductOf(x, nm);
}

void
TrabBoneEmbed :: giveRealStressVector(FloatArray &answer, MatResponseForm form, GaussPoint *gp,
                                      const FloatArray &totalStrain,
                                      TimeStep *atTime)
{
    double tempDam, tempTSED;
    FloatArray newTotalDef, plasDef;
    FloatArray totalStress;
    FloatMatrix compliance, elasticity;

    this->constructIsoComplTensor(compliance, eps0, nu0);
    elasticity.beInverseOf(compliance);

    TrabBoneEmbedStatus *status = ( TrabBoneEmbedStatus * ) this->giveStatus(gp);

    this->initGpForNewStep(gp);

    performPlasticityReturn(gp, totalStrain);

    tempDam = computeDamage(gp, atTime);

    plasDef.resize(6);

    totalStress.beProductOf(elasticity, totalStrain);

    tempTSED = 0.5 * totalStrain.dotProduct(totalStress);

    answer.resize(6);
    answer = totalStress;
    status->setTempDam(tempDam);
    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(answer);
    status->setTempTSED(tempTSED);
}

void
SimpleCrossSection :: giveGeneralizedStress_Shell(FloatArray &answer, GaussPoint *gp, const FloatArray &strain, TimeStep *tStep)
{
  /**Note: (by bp): This assumes that the behaviour is elastic
     there exist a nuumber of nonlinear integral material models for beams/plates/shells
     defined directly in terms of integral forces and moments and corresponding 
     deformations and curvatures. This would require to implement support at material model level.
     Mikael: See earlier response to comment
  */
    StructuralMaterial *mat = static_cast< StructuralMaterial * >( this->giveMaterial(gp) );
    FloatArray elasticStrain, et, e0;
    FloatMatrix tangent;
    elasticStrain = strain;
    this->giveTemperatureVector(et, gp, tStep);
    if ( et.giveSize() ) {
        double thick = this->give(CS_Thickness, gp);
        mat->giveThermalDilatationVector(e0, gp, tStep);
        elasticStrain.at(1) -= e0.at(1) * ( et.at(1) - mat->giveReferenceTemperature() );
        elasticStrain.at(2) -= e0.at(2) * ( et.at(1) - mat->giveReferenceTemperature() );
        if ( et.giveSize() > 1 ) {
            elasticStrain.at(4) -= e0.at(1) * et.at(2) / thick;     // kappa_x
            elasticStrain.at(5) -= e0.at(2) * et.at(2) / thick;     // kappa_y
        }
    }
    this->give3dShellStiffMtrx(tangent, ElasticStiffness, gp, tStep);
    answer.beProductOf(tangent, elasticStrain);
    StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) );
    status->letTempStrainVectorBe(strain);
    status->letTempStressVectorBe(answer);
}

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

    aux.zero();
    int size = nSecVars * nSecNodes;

    //set displacement and nonlocal location array
    this->setDisplacementLocationArray(locU, nPrimNodes, nPrimVars, nSecNodes, nSecVars);
    this->setNonlocalLocationArray(locK, nPrimNodes, nPrimVars, nSecNodes, nSecVars);

    answer.resize(size);
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {
        this->computeNkappaMatrixAt(gp, Nk);
        for ( int j = 1; j <= nlSize; j++ ) {
            fKappa.at(j) = Nk.at(1, j);
        }

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

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

void
LIBeam3dNL :: giveInternalForcesVector(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{
    int i, j;
    Material *mat = this->giveMaterial();
    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
    GaussPoint *gp = iRule->getIntegrationPoint(0);
    FloatArray nm(6), TotalStressVector(6);
    FloatMatrix x;
    double s1, s2;

    // update temp triad
    this->updateTempTriad(tStep);

    if ( useUpdatedGpRecord == 1 ) {
        TotalStressVector = ( ( StructuralMaterialStatus * ) mat->giveStatus(gp) )
                            ->giveStressVector();
    } else {
        this->computeStressVector(TotalStressVector, gp, tStep);
    }

    for ( i = 1; i <= 3; i++ ) {
        s1 = s2 = 0.0;
        for ( j = 1; j <= 3; j++ ) {
            s1 += tempTc.at(i, j) * TotalStressVector.at(j);
            s2 += tempTc.at(i, j) * TotalStressVector.at(j + 3);
        }

        nm.at(i)   = s1;
        nm.at(i + 3) = s2;
    }

    this->computeXMtrx(x, tStep);
    answer.beProductOf(x, nm);
}

void PolygonLine :: computeIntersectionPoints(const FloatArray &iXStart, const FloatArray &iXEnd, std :: vector< FloatArray > &oIntersectionPoints) const
{
    const double detTol = 1.0e-15;

    int numSeg = this->giveNrVertices() - 1;
    for(int segIndex = 1; segIndex <= numSeg; segIndex++) {

        const FloatArray &xStart = this->giveVertex(segIndex);
        const FloatArray &xEnd = this->giveVertex(segIndex+1);

        const FloatArray t1 = {xEnd(0) - xStart(0), xEnd(1) - xStart(1)};
        const FloatArray t2 = {iXEnd(0) - iXStart(0), iXEnd(1) - iXStart(1)};

        double xi1 = 0.0, xi2 = 0.0;
        int maxIter = 1;

        for(int iter = 0; iter < maxIter; iter++) {
            FloatArray temp = {iXStart(0) + xi2*t2(0) - xStart(0) - xi1*t1(0), iXStart(1) + xi2*t2(1) - xStart(1) - xi1*t1(1)};
            FloatArray res = {-t1.dotProduct(temp), t2.dotProduct(temp)};

            //printf("iter: %d res: %e\n", iter, res.computeNorm() );

            FloatMatrix K(2,2);
            K(0,0) = t1.dotProduct(t1);
            K(0,1) = -t1.dotProduct(t2);
            K(1,0) = -t1.dotProduct(t2);
            K(1,1) = t2.dotProduct(t2);

            double detK = K.giveDeterminant();

            if(detK < detTol) {
                return;
            }

            FloatMatrix KInv;
            KInv.beInverseOf(K);

            FloatArray dxi;
            dxi.beProductOf(KInv, res);

            xi1 -= dxi(0);
            xi2 -= dxi(1);
        }

//        printf("xi1: %e xi2: %e\n", xi1, xi2);


        if(xi1 >= 0.0 && xi1 <= 1.0 && xi2 >= 0.0 && xi2 <= 1.0) {
            FloatArray pos = xStart;
            pos.add(xi1, t1);
            oIntersectionPoints.push_back(pos);
        }

    }

}

void
GradDpElement :: computeNonlocalGradient(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    FloatMatrix Bk;
    FloatArray u;

    this->computeBkappaMatrixAt(gp, Bk);
    this->computeNonlocalDegreesOfFreedom(u, tStep);
    answer.beProductOf(Bk, u);
}

void
CBSElement :: computePrescribedTermsI(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    FloatMatrix mass;
    FloatArray usp;
    this->computeConsistentMassMtrx(mass, tStep);
    this->computeVectorOf(EID_MomentumBalance, mode, tStep, usp);
    answer.beProductOf(mass, usp);
    answer.negated();
}

void
CBSElement :: computePrescribedTermsI(FloatArray &answer, TimeStep *tStep)
{
    FloatMatrix mass;
    FloatArray usp;
    this->computeConsistentMassMtrx(mass, tStep);
    this->computeVectorOfVelocities(VM_Incremental, tStep, usp);
    answer.beProductOf(mass, usp);
    answer.negated();
}

void
CBSElement :: computePrescribedTermsII(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    FloatMatrix lhs;
    FloatArray usp;
    this->computePressureLhs(lhs, tStep);
    this->computeVectorOfPressures(mode, tStep, usp);
    answer.beProductOf(lhs, usp);
    answer.negated();
}

void
SUPGElement2 :: computeDeviatoricStrain(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    FloatArray u;
    FloatMatrix b;
    this->computeVectorOfVelocities(VM_Total, tStep, u);

    this->computeBMatrix(b, gp);
    answer.beProductOf(b, u);
}