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

void
LIBeam3dNL ::  computeRotMtrx(FloatMatrix &answer, FloatArray &psi)
{
    FloatMatrix S(3, 3), SS(3, 3);
    double psiSize;

    if ( psi.giveSize() != 3 ) {
        _error("computeSMtrx: psi param size mismatch");
    }

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

    psiSize = psi.computeNorm();
    answer.at(1, 1) = answer.at(2, 2) = answer.at(3, 3) = 1.;

    if ( psiSize <= 1.e-40 ) {
        return;
    }

    this->computeSMtrx(S, psi);
    SS.beProductOf(S, S);
    S.times(sin(psiSize) / psiSize);
    SS.times( ( 1. - cos(psiSize) ) / ( psiSize * psiSize ) );

    answer.add(S);
    answer.add(SS);
}

// returns the consistent (algorithmic) tangent stiffness matrix
void
MisesMat :: give3dSSMaterialStiffnessMatrix(FloatMatrix &answer,
                                            MatResponseMode mode,
                                            GaussPoint *gp,
                                            TimeStep *atTime)
{
    // start from the elastic stiffness
    this->giveLinearElasticMaterial()->give3dMaterialStiffnessMatrix(answer, mode, gp, atTime);
    if ( mode != TangentStiffness ) {
        return;
    }

    MisesMatStatus *status = static_cast< MisesMatStatus * >( this->giveStatus(gp) );
    double kappa = status->giveCumulativePlasticStrain();
    double tempKappa = status->giveTempCumulativePlasticStrain();
    // increment of cumulative plastic strain as an indicator of plastic loading
    double dKappa = tempKappa - kappa;

    if ( dKappa <= 0.0 ) { // elastic loading - elastic stiffness plays the role of tangent stiffness
        return;
    }

    // === plastic loading ===

    // yield stress at the beginning of the step
    double sigmaY = sig0 + H * kappa;

    // trial deviatoric stress and its norm
    StressVector trialStressDev(_3dMat);
    double trialStressVol;
    status->giveTrialStressVol(trialStressVol);
    status->giveTrialStressDev(trialStressDev);
    double trialS = trialStressDev.computeStressNorm();

    // one correction term
    FloatMatrix stiffnessCorrection(6, 6);
    stiffnessCorrection.beDyadicProductOf(trialStressDev, trialStressDev);
    double factor = -2. * sqrt(6.) * G * G / trialS;
    double factor1 = factor * sigmaY / ( ( H + 3. * G ) * trialS * trialS );
    stiffnessCorrection.times(factor1);
    answer.add(stiffnessCorrection);

    // another correction term
    stiffnessCorrection.bePinvID();
    double factor2 = factor * dKappa;
    stiffnessCorrection.times(factor2);
    answer.add(stiffnessCorrection);

    //influence of damage
    //    double omega = computeDamageParam(tempKappa);
    double omega = status->giveTempDamage();
    answer.times(1. - omega);
    FloatArray effStress;
    status->giveTempEffectiveStress(effStress);
    double omegaPrime = computeDamageParamPrime(tempKappa);
    double scalar = -omegaPrime *sqrt(6.) * G / ( 3. * G + H ) / trialS;
    stiffnessCorrection.beDyadicProductOf(effStress, trialStressDev);
    stiffnessCorrection.times(scalar);
    answer.add(stiffnessCorrection);
}

void
MisesMatGrad :: givePlaneStrainStiffMtrx(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    this->giveLinearElasticMaterial()->giveStiffnessMatrix(answer, mode, gp, tStep);
    if ( mode != TangentStiffness ) {
        return;
    }

    MisesMatGradStatus *status = static_cast< MisesMatGradStatus * >( this->giveStatus(gp) );
    double tempKappa = status->giveTempCumulativePlasticStrain();
    double kappa = status->giveCumulativePlasticStrain();
    double dKappa = tempKappa - kappa;
    double tempDamage = status->giveTempDamage();
    double damage = status->giveDamage();

    if ( dKappa <= 0.0 ) { // elastic loading - elastic stiffness plays the role of tangent stiffness
        return;
    }

    // === plastic loading ===
    // yield stress at the beginning of the step
    double sigmaY = sig0 + H * kappa;
    // trial deviatoric stress and its norm
    StressVector trialStressDev(status->giveTrialStressDev(), _PlaneStrain);
    double trialS = trialStressDev.computeStressNorm();
    // volumetric stress
    //double trialStressVol = status->giveTrialStressVol();
    // one correction term
    FloatMatrix stiffnessCorrection(4, 4);
    stiffnessCorrection.beDyadicProductOf(trialStressDev, trialStressDev);
    double factor = -2. * sqrt(6.) * G * G / trialS;
    double factor1 = factor * sigmaY / ( ( H + 3. * G ) * trialS * trialS );
    stiffnessCorrection.times(factor1);
    answer.add(stiffnessCorrection);
    // another correction term
    stiffnessCorrection.zero();
    stiffnessCorrection.at(1, 1) = stiffnessCorrection.at(2, 2) = stiffnessCorrection.at(3, 3) = 2. / 3.;
    stiffnessCorrection.at(1, 2) = stiffnessCorrection.at(1, 3) = stiffnessCorrection.at(2, 1) = -1. / 3.;
    stiffnessCorrection.at(2, 3) = stiffnessCorrection.at(3, 1) = stiffnessCorrection.at(3, 2) = -1. / 3.;
    stiffnessCorrection.at(4, 4) = 0.5;
    double factor2 = factor * dKappa;
    stiffnessCorrection.times(factor2);
    answer.add(stiffnessCorrection);
    //influence of damage
    answer.times(1 - tempDamage);
    if ( tempDamage > damage ) {
        const FloatArray &effStress = status->giveTempEffectiveStress();
        double nlKappa = status->giveNonlocalCumulatedStrain();
        kappa = mParam * nlKappa + ( 1. - mParam ) * tempKappa;
        double omegaPrime = computeDamageParamPrime(kappa);
        double scalar = -omegaPrime *sqrt(6.) * G / ( 3. * G + H ) / trialS;
        stiffnessCorrection.beDyadicProductOf(effStress, trialStressDev);
        stiffnessCorrection.times( scalar * ( 1. - mParam ) );
        answer.add(stiffnessCorrection);
    }
}

void Tr21Stokes :: giveIntegratedVelocity(FloatMatrix &answer, TimeStep *tStep )
{
    /*
    * Integrate velocity over element
    */

    IntegrationRule *iRule = integrationRulesArray [ 0 ];
    FloatMatrix v, v_gamma, ThisAnswer, boundaryV, Nmatrix;
    double detJ;
    FloatArray *lcoords, N;
    int i, j, k=0;
    Dof *d;
    GaussPoint *gp;

    v.resize(12,1);
    v.zero();
    boundaryV.resize(2,1);


    for (i=1; i<=this->giveNumberOfDofManagers(); i++) {
        for (j=1; j<=this->giveDofManager(i)->giveNumberOfDofs(); j++) {
            d = this->giveDofManager(i)->giveDof(j);
            if ((d->giveDofID()==V_u) || (d->giveDofID()==V_v)) {
                k=k+1;
                v.at(k,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);
            /*} else if (d->giveDofID()==A_x) {
                boundaryV.at(1,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);
            } else if (d->giveDofID()==A_y) {
                boundaryV.at(2,1)=d->giveUnknown(EID_ConservationEquation, VM_Total, tStep);*/
            }
        }
    }

    answer.resize(2,1);
    answer.zero();

    Nmatrix.resize(2,12);

    for (i=0; i<iRule->getNumberOfIntegrationPoints(); i++) {

        gp = iRule->getIntegrationPoint(i);

        lcoords = gp->giveCoordinates();

        this->interpolation_quad.evalN(N, *lcoords, FEIElementGeometryWrapper(this));
        detJ = this->interpolation_quad.giveTransformationJacobian(*lcoords, FEIElementGeometryWrapper(this));

        N.times(detJ*gp->giveWeight());

        for (j=1; j<=6;j++) {
            Nmatrix.at(1,j*2-1)=N.at(j);
            Nmatrix.at(2,j*2)=N.at(j);
        }

        ThisAnswer.beProductOf(Nmatrix,v);
        answer.add(ThisAnswer);

    }

}

void
LIBeam2dNL :: computeInitialStressMatrix(FloatMatrix &answer, TimeStep *tStep)
{
    int i, j, n;
    double dV;
    GaussPoint *gp;
    IntegrationRule *iRule;
    FloatArray stress;
    FloatMatrix A;
    Material *mat = this->giveMaterial();

    answer.resize(6, 6);
    answer.zero();

    iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
    // assemble initial stress matrix
    for ( i = 0; i < iRule->getNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        dV = this->computeVolumeAround(gp);
        stress = ( ( StructuralMaterialStatus * ) mat->giveStatus(gp) )->giveStressVector();
        n = stress.giveSize();
        if ( n ) {
            for ( j = 1; j <= n; j++ ) {
                // loop over each component of strain vector
                this->computeNLBMatrixAt(A, gp, j);
                if ( A.isNotEmpty() ) {
                    A.times(stress.at(j) * dV);
                    answer.add(A);
                }
            }
        }
    }
}

void
NonStationaryTransportProblem :: assembleDirichletBcRhsVector(FloatArray &answer, TimeStep *tStep,
                                                              ValueModeType mode,
                                                              const UnknownNumberingScheme &ns, Domain *d)
{
    IntArray loc, dofids;
    FloatArray rp, charVec;
    FloatMatrix s;
    FloatMatrix capacity;

    int nelem = d->giveNumberOfElements();

    for ( int ielem = 1; ielem <= nelem; ielem++ ) {
        Element *element = d->giveElement(ielem);

        element->giveElementDofIDMask(dofids);
        element->computeVectorOfPrescribed(dofids, mode, tStep, rp);
        if ( rp.containsOnlyZeroes() ) {
            continue;
        } else {
            element->giveCharacteristicMatrix(s, TangentStiffnessMatrix, tStep);
            element->giveCharacteristicMatrix(capacity, lumpedCapacityStab ? LumpedMassMatrix : MassMatrix, tStep);
            s.times(this->alpha);
            s.add(1. / tStep->giveTimeIncrement(), capacity);

            charVec.beProductOf(s, rp);
            charVec.negated();

            element->giveLocationArray(loc, ns);
            answer.assemble(charVec, loc);
        }
    } // end element loop
}

void
SUPGElement :: computeBCLhsPressureTerm_MB(FloatMatrix &answer, TimeStep *tStep)
{
    bcType boundarytype;
    int nLoads = 0;
    //bcType loadtype;
    FloatMatrix helpMatrix;
    // loop over boundary load array
    answer.clear();

    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;

    if ( nLoads ) {
        for ( int i = 1; i <= nLoads; i++ ) {
            int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
            int side = boundaryLoadArray.at(i * 2);
            Load *load = domain->giveLoad(n);
            boundarytype = load->giveType();
            if ( boundarytype == OutFlowBC ) {
                this->computeOutFlowBCTerm_MB(helpMatrix, side, tStep);
                answer.add(helpMatrix);
            } else {
                //_warning("computeForceLoadVector : unsupported load type class");
            }
        }
    }
}

void
NonLinearDynamic :: giveElementCharacteristicMatrix(FloatMatrix &answer, int num,
                                                    CharType type, TimeStep *tStep, Domain *domain)
{
    // We don't directly call element ->GiveCharacteristicMatrix() function, because some
    // engngm classes may require special modification of base types supported on
    // element class level.

    if ( type == EffectiveStiffnessMatrix ) {
        Element *element;
        FloatMatrix charMtrx;

        element = domain->giveElement(num);
        element->giveCharacteristicMatrix(answer, TangentStiffnessMatrix, tStep);
        answer.times(1 + this->delta * a1);

        element->giveCharacteristicMatrix(charMtrx, MassMatrix, tStep);
        charMtrx.times(this->a0 + this->eta * this->a1);

        answer.add(charMtrx);

        return;
    } else {
        StructuralEngngModel :: giveElementCharacteristicMatrix(answer, num, type, tStep, domain);
    }
}

void MidpointLhsAssembler :: matrixFromElement(FloatMatrix &answer, Element &el, TimeStep *tStep) const
{
    FloatMatrix capacity;
    el.giveCharacteristicMatrix(answer, TangentStiffnessMatrix, tStep);
    el.giveCharacteristicMatrix(capacity, this->lumped ? LumpedMassMatrix : MassMatrix, tStep);
    answer.times(this->alpha);
    answer.add(1. / tStep->giveTimeIncrement(), capacity);
}

void
SUPGElement :: computeBCLhsTerm_MB(FloatMatrix &answer, TimeStep *tStep)
{
    bcType boundarytype;
    int nLoads = 0;
    //bcType loadtype;
    FloatMatrix helpMatrix;
    // loop over boundary load array

    answer.clear();

    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    if ( nLoads ) {
        for ( int i = 1; i <= nLoads; i++ ) {
            int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
            int side = boundaryLoadArray.at(i * 2);
            Load *load = domain->giveLoad(n);
            boundarytype = load->giveType();
            if ( boundarytype == SlipWithFriction ) {
                this->computeSlipWithFrictionBCTerm_MB(helpMatrix, load, side, tStep);
                answer.add(helpMatrix);
            } else if ( boundarytype == PenetrationWithResistance ) {
                this->computePenetrationWithResistanceBCTerm_MB(helpMatrix, load, side, tStep);
                answer.add(helpMatrix);
            } else {
                // OOFEM_ERROR("unsupported load type class");
            }
        }
    }

    nLoads = this->giveBodyLoadArray()->giveSize();

    if ( nLoads ) {
        bcGeomType ltype;
        for ( int i = 1; i <= nLoads; i++ ) {
            Load *load = domain->giveLoad( bodyLoadArray.at(i) );
            ltype = load->giveBCGeoType();
            if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ReinforceBVT ) ) {
                this->computeHomogenizedReinforceTerm_MB(helpMatrix, load, tStep);
                answer.add(helpMatrix);
            }
        }
    }
}

void
TrPlanestressRotAllman :: computeStiffnessMatrix(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    // compute standard stiffness matrix
    TrPlaneStress2d::computeStiffnessMatrix (answer, rMode, tStep);
    // add zero energy mode stabilization
    FloatMatrix ks;
    this->computeStiffnessMatrixZeroEnergyStabilization(ks, rMode, tStep);
    answer.add(ks);
}

void
PerfectlyPlasticMaterial :: giveMaterialStiffnessMatrix(FloatMatrix &answer, MatResponseMode mode,
                                                        GaussPoint *gp,
                                                        TimeStep *tStep)
//
//
//
// computes full constitutive matrix for case of gp stress-strain state.
// it returns elasto-plastic stiffness material matrix.
// if strainIncrement == NULL a loading is assumed
// for detailed description see (W.F.Chen: Plasticity in Reinforced Concrete, McGraw-Hill, 1982,
// chapter 6.)
//
// if derived material would like to implement failure behaviour
// it must redefine basic Give3dMaterialStiffnessMatrix function
// in order to take possible failure (tension cracking) into account
//
// in this function answer is Material stiffness matrix respecting
// current stress-strain mode in gp. This is reached by using
// impose constraints functions
{
    FloatArray statusFullStressVector, statusFullPlasticVector, plasticStrainVector;
    double lambda;
    PerfectlyPlasticMaterialStatus *status = static_cast< PerfectlyPlasticMaterialStatus * >( this->giveStatus(gp) );

    // double f = domain->giveYieldCriteria(yieldCriteria)->
    //  computeValueAt(gp->giveStressVector(), gp->givePlasticStrainVector(),
    //       gp-> giveHardeningParam());

    this->giveEffectiveMaterialStiffnessMatrix(answer, mode, gp, tStep);
    // if (f < YIELD_BOUNDARY)  // linear elastic behaviour
    //  return de;

    if ( status->giveYieldFlag() == 0 ) {
        return;
    }

    // if secant stiffness requested assume initial elastic matrix
    if ( mode == SecantStiffness ) {
        return;
    }

    plasticStrainVector = status->givePlasticStrainVector();
    StructuralMaterial :: giveFullSymVectorForm( statusFullPlasticVector, plasticStrainVector, gp->giveMaterialMode() );
    StructuralMaterial :: giveFullSymVectorForm( statusFullStressVector, status->giveStressVector(), gp->giveMaterialMode() );

    // yield condition satisfied
    FloatMatrix dp;
    this->computePlasticStiffnessAt(dp, gp, & statusFullStressVector,
                                    & statusFullPlasticVector,
                                    NULL,
                                    tStep,
                                    lambda);
    answer.add(dp);
}

void
NonlinearMassTransferMaterial :: giveCharacteristicMatrix(FloatMatrix &answer,
                                                          MatResponseForm form,
                                                          MatResponseMode mode,
                                                          GaussPoint *gp,
                                                          TimeStep *atTime)
{
    MaterialMode mMode = gp->giveMaterialMode();
    AnisotropicMassTransferMaterialStatus *status = ( ( AnisotropicMassTransferMaterialStatus * ) this->giveStatus(gp) );
    FloatArray eps = status->giveGradP();
    double gradPNorm;
    FloatMatrix t1, t2;

    gradPNorm = eps.computeNorm();

    t1.beDyadicProductOf(eps, eps);
    if ( gradPNorm != 0.0 ) {
        t1.times( C * alpha * pow(gradPNorm, alpha - 2) );
    }

    switch  ( mMode ) {
    case _1dHeat:
        t2.resize(1, 1);
        t2.at(1, 1) = 1;
        break;
    case _2dHeat:
        t2.resize(2, 2);
        t2.at(1, 1) = t2.at(2, 2) = 1;
        break;
    case _3dHeat:
        t2.resize(3, 3);
        t2.at(1, 1) = t2.at(2, 2) = t2.at(3, 3) = 1;
        break;
    default:
        _error2( "giveCharacteristicMatrix : unknown mode (%s)", __MaterialModeToString(mMode) );
    }
    
    answer.beEmptyMtrx();
    answer.add(t1);
    answer.add(1 + C * pow(gradPNorm, alpha), t2);

}

void
NonlinearMassTransferMaterial :: giveCharacteristicMatrix(FloatMatrix &answer,
                                                          MatResponseMode mode,
                                                          GaussPoint *gp,
                                                          TimeStep *tStep)
{
    MaterialMode mMode = gp->giveMaterialMode();
    TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    const FloatArray &eps = status->giveTempGradient();
    double gradPNorm;
    FloatMatrix t1, t2;

    gradPNorm = eps.computeNorm();

    t1.beDyadicProductOf(eps, eps);
    if ( gradPNorm != 0.0 ) {
        t1.times( C * alpha * pow(gradPNorm, alpha - 2) );
    }

    switch  ( mMode ) {
    case _1dHeat:
        t2.resize(1, 1);
        break;
    case _2dHeat:
        t2.resize(2, 2);
        break;
    case _3dHeat:
        t2.resize(3, 3);
        break;
    default:
        OOFEM_ERROR("unknown mode (%s)", __MaterialModeToString(mMode));
    }
    t2.beUnitMatrix();

    answer.clear();
    answer.add(t1);
    answer.add(1 + C * pow(gradPNorm, alpha), t2);
}

void
GradDpElement :: computeStiffnessMatrix_uu(FloatMatrix &answer, MatResponseMode rMode, TimeStep *tStep)
{
    NLStructuralElement *elem = this->giveNLStructuralElement();
    StructuralCrossSection *cs = elem->giveStructuralCrossSection();
    FloatMatrix B, D, DB;


    int nlGeo = elem->giveGeometryMode();
    bool matStiffSymmFlag = elem->giveCrossSection()->isCharacteristicMtrxSymmetric(rMode);
    answer.clear();
    for ( GaussPoint *gp: *elem->giveIntegrationRule(0) ) {

        GradDpMaterialExtensionInterface *dpmat = dynamic_cast< GradDpMaterialExtensionInterface * >(
                    cs->giveMaterialInterface(GradDpMaterialExtensionInterfaceType, gp) );
        if ( !dpmat ) {
            OOFEM_ERROR("Material doesn't implement the required DpGrad interface!");
        }
        if ( nlGeo == 0 ) {
            elem->computeBmatrixAt(gp, B);
        } else if ( nlGeo == 1 ) {
            if ( elem->domain->giveEngngModel()->giveFormulation() == AL ) {
                elem->computeBmatrixAt(gp, B);
            } else {
                elem->computeBHmatrixAt(gp, B);
            }
        }

        dpmat->givePDGradMatrix_uu(D, rMode, gp, tStep);
        double dV = elem->computeVolumeAround(gp);
        DB.beProductOf(D, B);
        if ( matStiffSymmFlag ) {
            answer.plusProductSymmUpper(B, DB, dV);
        } else {
            answer.plusProductUnsym(B, DB, dV);
        }
    }

    if ( elem->domain->giveEngngModel()->giveFormulation() == AL ) {
        FloatMatrix initialStressMatrix;
        elem->computeInitialStressMatrix(initialStressMatrix, tStep);
        answer.add(initialStressMatrix);
    }

    if ( matStiffSymmFlag ) {
        answer.symmetrized();
    }
}