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

void
StructuralInterfaceElementPhF :: computeTraction(FloatArray &traction, IntegrationPoint *ip, FloatArray &jump, TimeStep *tStep)
{
    // Returns the traction in global coordinate system
    FloatMatrix rotationMatGtoL, F;
    this->computeTransformationMatrixAt(ip, rotationMatGtoL);
    jump.rotatedWith(rotationMatGtoL, 'n');      // transform jump to local coord system

    double damage = this->computeDamageAt(ip, VM_Total, tStep);
    
    this->giveEngTraction(traction, ip, jump, damage, tStep);
    
    
    traction.rotatedWith(rotationMatGtoL, 't');     // transform traction to global coord system
}

void
InteractionPFEMParticle :: givePrescribedUnknownVector(FloatArray &answer, const IntArray &dofIDArry,
                                                       ValueModeType mode, TimeStep *stepN)
{
    answer.resize( dofIDArry.giveSize() );


    FloatArray velocities;
    FluidStructureProblem *fsiProblem = this->giveFluidStructureMasterProblem();
    if (fsiProblem) {
        StructuralEngngModel *structuralProblem = this->giveStructuralProblem();
        if ( structuralProblem ) {
            int j = 1;
            if (fsiProblem->giveIterationNumber() < 1) {  
                for (int dofid: dofIDArry ) {
                    answer.at(j++) = this->giveDofWithID( dofid )->giveBcValue(mode, stepN);
                }
            } else {
                DofManager *dman = structuralProblem->giveDomain(1)->giveDofManager(coupledNode);
                //dman->giveUnknownVectorOfType(velocities, VelocityVector, VM_Velocity, stepN);
                //dman->giveUnknownVector(velocities, dofIDArry, VM_Velocity, stepN);
                dman->giveCompleteUnknownVector(velocities, VM_Velocity, stepN);
                for ( int dofid: dofIDArry) {
                    answer.at(dofid) = velocities.at( dofid );
                }
            }
        }
    }

    // Transform to global c.s.
    FloatMatrix L2G;
    if (this->computeL2GTransformation(L2G, dofIDArry)) {
        answer.rotatedWith(L2G, 'n');
    }
}

void DofManager :: giveUnknownVector(FloatArray &answer, const IntArray &dofIDArry,
                                     PrimaryField &field, ValueModeType mode, TimeStep *tStep, bool padding)
{
    answer.resize( dofIDArry.giveSize() );

    int k = 0;
    for ( auto &dofid: dofIDArry ) {
        auto pos = this->findDofWithDofId( ( DofIDItem ) dofid );
        if ( pos == this->end() ) {
            if ( padding ) {
                answer.at(++k) = 0.;
                continue;
            } else {
                continue;
            }
        }
        answer.at(++k) = (*pos)->giveUnknown(field, mode, tStep);
    }
    answer.resizeWithValues(k);

    // Transform to global c.s.
    FloatMatrix L2G;
    if ( this->computeL2GTransformation(L2G, dofIDArry) ) {
        answer.rotatedWith(L2G, 'n');
    }
}

void
StructuralInterfaceElement :: computeTraction(FloatArray &traction, IntegrationPoint *ip, FloatArray &jump, TimeStep *tStep)
{
    // Returns the traction in global coordinate system
    FloatMatrix rotationMatGtoL, F;
    this->computeTransformationMatrixAt(ip, rotationMatGtoL);
    jump.rotatedWith(rotationMatGtoL, 'n');      // transform jump to local coord system

    if ( this->nlGeometry == 0 ) {
        this->giveEngTraction(traction, ip, jump, tStep);
    } else if ( this->nlGeometry == 1 ) {
        ///@todo compute F in a proper way
        F.beUnitMatrix();
        F.rotatedWith(rotationMatGtoL, 'n');
        this->giveFirstPKTraction(traction, ip, jump, F, tStep);
    }

    traction.rotatedWith(rotationMatGtoL, 't');     // transform traction to global coord system
}

void
CCTPlate3d :: computeBodyLoadVectorAt(FloatArray &answer, Load *forLoad, TimeStep *stepN, ValueModeType mode)
// Computes numerically the load vector of the receiver due to the body loads, at stepN.
// load is assumed to be in global cs.
// load vector is then transformed to coordinate system in each node.
// (should be global coordinate system, but there may be defined
//  different coordinate system in each node)
{
    double dens, dV, load;
    GaussPoint *gp = NULL;
    FloatArray force;
    FloatMatrix T;

    if ( ( forLoad->giveBCGeoType() != BodyLoadBGT ) || ( forLoad->giveBCValType() != ForceLoadBVT ) ) {
        _error("computeBodyLoadVectorAt: unknown load type");
    }

    GaussIntegrationRule irule(1, this, 1, 5);
    irule.SetUpPointsOnTriangle(1, _2dPlate);

    // note: force is assumed to be in global coordinate system.
    forLoad->computeComponentArrayAt(force, stepN, mode);

    if ( force.giveSize() ) {
        gp = irule.getIntegrationPoint(0);

        dens = this->giveMaterial()->give('d', gp);
        dV   = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness);

        answer.resize(18);
        answer.zero();

        load = force.at(1) * dens * dV / 3.0;
        answer.at(1) = load;
        answer.at(7) = load;
        answer.at(13) = load;

        load = force.at(2) * dens * dV / 3.0;
        answer.at(2) = load;
        answer.at(8) = load;
        answer.at(14) = load;

        load = force.at(3) * dens * dV / 3.0;
        answer.at(3) = load;
        answer.at(9) = load;
        answer.at(15) = load;

        // transform result from global cs to local element cs.
        if ( this->computeGtoLRotationMatrix(T) ) {
            answer.rotatedWith(T, 'n');
        }
    } else {
        answer.resize(0);          // nil resultant
    }
}

void DofManager :: givePrescribedUnknownVector(FloatArray &answer, const IntArray &dofIDArry,
                                               ValueModeType mode, TimeStep *tStep)
{
    answer.resize(dofIDArry.giveSize());

    int j = 1;
    for ( int dofid: dofIDArry ) {
        answer.at(j++) = this->giveDofWithID( dofid )->giveBcValue(mode, tStep);
    }

    // Transform to global c.s.
    FloatMatrix L2G;
    if ( this->computeL2GTransformation(L2G, dofIDArry) ) {
        answer.rotatedWith(L2G, 'n');
    }
}

void
RerShell :: computeBodyLoadVectorAt(FloatArray &answer, Load *forLoad, TimeStep *stepN, ValueModeType mode)
// Computes numerically the load vector of the receiver due to the body
// loads, at stepN. - no support for momentum bodyload
{
    double dens, dV, load;
    GaussPoint *gp = integrationRulesArray [ 0 ]->getIntegrationPoint(0);
    FloatArray f;
    FloatMatrix T;


    forLoad->computeComponentArrayAt(f, stepN, mode);
    //f.times( this->giveMaterial()->give('d') );

    if ( f.giveSize() == 0 ) {
        answer.resize(0);
        return;                                             // nil resultant
    } else {
        dens = this->giveMaterial()->give('d', gp);
        dV = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness);

        answer.resize(18);

        load = f.at(1) * dens * dV / 3.0;
        answer.at(1) = load;
        answer.at(7) = load;
        answer.at(13) = load;

        load = f.at(2) * dens * dV / 3.0;
        answer.at(2) = load;
        answer.at(8) = load;
        answer.at(14) = load;

        load = f.at(3) * dens * dV / 3.0;
        answer.at(3) = load;
        answer.at(9) = load;
        answer.at(15) = load;

        // transform result from global cs to local element  cs.
        if ( this->computeGtoLRotationMatrix(T) ) {
            answer.rotatedWith(T, 'n');
        }

        return;
    }
}

void
Beam2d :: computeBoundaryEdgeLoadVector(FloatArray &answer, BoundaryLoad *load, int edge, CharType type, ValueModeType mode, TimeStep *tStep, bool global)
{
    answer.clear();

    if ( edge != 1 ) {
        OOFEM_ERROR("Beam2D only has 1 edge (the midline) that supports loads. Attempted to apply load to edge %d", edge);
    }

    if ( type != ExternalForcesVector ) {
        return;
    }

    double l = this->computeLength();
    FloatArray coords, t;
    FloatMatrix N, T;

    answer.clear();
    for ( auto &gp : *this->giveDefaultIntegrationRulePtr() ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        this->computeNmatrixAt(lcoords, N);
        if ( load ) {
            this->computeGlobalCoordinates(coords, lcoords);
            load->computeValues(t, tStep, coords, { D_u, D_w, R_v }, mode);
        } else {
            load->computeValues(t, tStep, lcoords, { D_u, D_w, R_v }, mode);
        }

        if ( load->giveCoordSystMode() == Load :: CST_Global ) {
            if ( this->computeLoadGToLRotationMtrx(T) ) {
                t.rotatedWith(T, 'n');
            }
        }

        double dl = gp->giveWeight() * 0.5 * l;
        answer.plusProduct(N, t, dl);
    }

    if (global) {
      // Loads from sets expects global c.s.
      this->computeGtoLRotationMatrix(T);
      answer.rotatedWith(T, 't');
    }
}

void
GradDpElement :: computeLocForceLoadVector(FloatArray &answer, TimeStep *tStep, ValueModeType mode)
// computes the part of load vector, which is imposed by force loads acting
// on element volume (surface).
// When reactions forces are computed, they are computed from element::GiveRealStressVector
// in this vector a real forces are stored (temperature part is subtracted).
// so we need further sobstract part corresponding to non-nodeal loading.
{
    FloatMatrix T;
    NLStructuralElement *elem = this->giveNLStructuralElement();
    elem->computeLocalForceLoadVector(answer, tStep, mode);

    // transform result from global cs to nodal cs. if necessary
    if ( answer.isNotEmpty() ) {
        if ( elem->computeGtoLRotationMatrix(T) ) {
            // first back to global cs from element local
            answer.rotatedWith(T, 't');
        }
    } else {
        answer.resize(locSize);
        answer.zero();
    }
}

//.........这里部分代码省略.........
                residualStrength = E * e0 * exp(-1. * ( tempKappa - e0 ) / ef);
            } else {
                OOFEM_ERROR("Unknown softening type for cohesive crack model.");
            }
        }

        if ( status ) {
            status->setResidualTensileStrength(residualStrength);
        }

#endif
    } else {
        // damage grows
        tempKappa = kappa;

        FloatArray crackPlaneNormal(3);
        for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
            crackPlaneNormal.at(i) = principalDir.at(i, 1);
        }
        this->initDamaged(tempKappa, crackPlaneNormal, gp, tStep);
        this->computeDamage(omega, tempKappa, gp);
    }

    answer.zero();

    if ( omega > 0. ) {
        tempNominalStress = tempEffectiveStress;

        if ( this->isotropic ) {
            // convert effective stress to nominal stress
            tempNominalStress.times(1. - omega);
            answer.add(tempNominalStress);
        } else {
            // stress transformation matrix
            FloatMatrix Tstress;

            // compute principal nominal stresses by multiplying effective stresses by damage
            for ( int i = 1; i <= principalStress.giveSize(); i++ ) {
                if ( principalStress.at(i) > 0. ) {
                    // convert principal effective stress to nominal stress
                    principalStress.at(i) *= ( 1. - omega );
                }
            }

            if ( mode == _PlaneStress ) {
                principalStress.resizeWithValues(3);
                givePlaneStressVectorTranformationMtrx(Tstress, principalDir, true);
            } else {
                principalStress.resizeWithValues(6);
                giveStressVectorTranformationMtrx(Tstress, principalDir, true);
            }

            principalStress.rotatedWith(Tstress, 'n');


            if ( mode == _PlaneStress ) { // plane stress
                answer.add(principalStress);
            } else if ( this->giveSizeOfVoigtSymVector(mode) != this->giveSizeOfVoigtSymVector(_3dMat) ) { // mode = _PlaneStrain or axial symmetry
                StressVector redFormStress(mode);
                redFormStress.convertFromFullForm(principalStress, mode);
                answer.add(redFormStress);
            } else { // 3D
                answer.add(principalStress);
            }
        }
    } else {
        answer.add(tempEffectiveStress);
    }

#ifdef supplementary_info

    if ( ( omega == 0. ) || ( sigma1 <= 0 ) ) {
        status->setCrackWidth(0.);
    } else {
        FloatArray principalStrains;
        this->computePrincipalValues(principalStrains, totalStrain, principal_strain);

        double crackWidth;
        //case when the strain localizes into narrow band and after a while the stresses relax

        double strainWithoutTemperShrink = principalStrains.at(1);
        strainWithoutTemperShrink -=  status->giveTempThermalStrain();
        strainWithoutTemperShrink -=  status->giveTempDryingShrinkageStrain();
        strainWithoutTemperShrink -=  status->giveTempAutogenousShrinkageStrain();


        crackWidth = status->giveCharLength() * omega * strainWithoutTemperShrink;

        status->setCrackWidth(crackWidth);
    }

#endif

    // update gp
    status->letTempStrainVectorBe(totalStrain);
    status->letTempStressVectorBe(answer);
    status->letTempViscoelasticStressVectorBe(tempEffectiveStress);
    status->setTempKappa(tempKappa);
    status->setTempDamage(omega);
}

void
StructuralEngngModel :: giveInternalForces(FloatArray &answer, bool normFlag, int di, TimeStep *stepN)
{
    // Simply assembles contributions from each element in domain
    Element *element;
    IntArray loc;
    FloatArray charVec;
    FloatMatrix R;
    int nelems;
    EModelDefaultEquationNumbering dn;

    answer.resize( this->giveNumberOfEquations( EID_MomentumBalance ) );
    answer.zero();
    if ( normFlag ) internalForcesEBENorm = 0.0;

    // Update solution state counter
    stepN->incrementStateCounter();

#ifdef __PARALLEL_MODE
    if ( isParallel() ) {
        exchangeRemoteElementData( RemoteElementExchangeTag  );
    }
#endif

    nelems = this->giveDomain(di)->giveNumberOfElements();
    this->timer.resumeTimer(EngngModelTimer :: EMTT_NetComputationalStepTimer);

    for ( int i = 1; i <= nelems; i++ ) {
        element = this->giveDomain(di)->giveElement(i);
#ifdef __PARALLEL_MODE
        // Skip remote elements (these are used as mirrors of remote elements on other domains
        // when nonlocal constitutive models are used. Their introduction is necessary to
        // allow local averaging on domains without fine grain communication between domains).
        if ( element->giveParallelMode() == Element_remote ) continue;

#endif
        element->giveLocationArray( loc, EID_MomentumBalance, dn );
        element->giveCharacteristicVector( charVec, InternalForcesVector, VM_Total, stepN );
        if ( element->giveRotationMatrix(R, EID_MomentumBalance) ) {
            charVec.rotatedWith(R, 't');
        }
        answer.assemble(charVec, loc);

        // Compute element norm contribution.
        if ( normFlag ) internalForcesEBENorm += charVec.computeSquaredNorm();
    }

    this->timer.pauseTimer( EngngModelTimer :: EMTT_NetComputationalStepTimer );

#ifdef __PARALLEL_MODE
    this->updateSharedDofManagers(answer, InternalForcesExchangeTag);

    if ( normFlag ) {
        // Exchange norm contributions.
        double localEBENorm = internalForcesEBENorm;
        internalForcesEBENorm = 0.0;
        MPI_Allreduce(& localEBENorm, & internalForcesEBENorm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
    }
#endif

    // Remember last internal vars update time stamp.
    internalVarUpdateStamp = stepN->giveSolutionStateCounter();
}