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

void
TrabBoneNL3D :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                        FloatArray &rcontrib, TimeStep *tStep)
{
    TrabBoneNL3DStatus *nlStatus = static_cast< TrabBoneNL3DStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );

    FloatMatrix b;

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);

    double kappa = nlStatus->giveKappa();
    double tempKappa = nlStatus->giveTempKappa();
    double dKappa = tempKappa - kappa;
    if ( dKappa < 10.e-9 ) {
        dKappa = 0;
    }

    if ( dKappa > 0.0 ) {
        FloatArray remoteNu, prodTensor;
        const FloatArray &plasFlowDirec = nlStatus->givePlasFlowDirec();
        const FloatMatrix &SSaTensor = nlStatus->giveSSaTensor();
        double beta = nlStatus->giveBeta();

        prodTensor.beTProductOf(SSaTensor, plasFlowDirec);
        remoteNu = 1 / beta * prodTensor;
        rcontrib.beTProductOf(b, remoteNu);
    } else {
        rcontrib.resize(b.giveNumberOfColumns());
        rcontrib.zero();
    }
}

void Line2SurfaceTension :: computeLoadVector(FloatArray &answer, ValueModeType mode, TimeStep *tStep)
{
    ///@todo Support axisymm.
    //domainType dt = this->giveDomain()->giveDomainType();
    IntegrationRule *iRule = this->integrationRulesArray [ 0 ];
    double t = 1, gamma_s;
    ///@todo Should i use this? Not used in FM module (but perhaps it should?) / Mikael.
    //t = this->giveDomain()->giveCrossSection(1)->give(CS_Thickness);
    gamma_s = this->giveMaterial()->give('g', NULL);

    FloatMatrix xy(2, 3);
    Node *node;
    for ( int i = 1; i <= 3; i++ ) {
        node = giveNode(i);
        xy.at(1, i) = node->giveCoordinate(1);
        xy.at(2, i) = node->giveCoordinate(2);
    }

    FloatArray A;
    FloatArray dNdxi(3);
    FloatArray es(2); // tangent vector to curve
    FloatMatrix BJ(2, 6);
    BJ.zero();

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

    for ( int k = 0; k < iRule->getNumberOfIntegrationPoints(); k++ ) {
        GaussPoint *gp = iRule->getIntegrationPoint(k);
        //interpolation.evaldNdx(dN, domain, dofManArray, * gp->giveCoordinates(), 0.0);
        double xi = gp->giveCoordinate(1);

        // Some simplifications can be performed, since the mapping J is a scalar.
        dNdxi.at(1) = -0.5 + xi;
        dNdxi.at(2) =  0.5 + xi;
        dNdxi.at(3) = -2.0 * xi;

        es.beProductOf(xy, dNdxi);
        double J = es.computeNorm();
        es.times(1 / J); //es.normalize();

        // dNds = dNdxi/J
        // B.at(1,1) = dNds.at(1); and so on.

        BJ.at(1, 1) = BJ.at(2, 2) = dNdxi.at(1);
        BJ.at(1, 3) = BJ.at(2, 4) = dNdxi.at(2);
        BJ.at(1, 5) = BJ.at(2, 6) = dNdxi.at(3);

        A.beTProductOf(BJ, es);
        answer.add( - gamma_s * t * gp->giveWeight(), A); // Note! Negative sign!
    }
}

void FloatArray :: rotatedWith(FloatMatrix &r, char mode)
// Returns the receiver 'a' rotated according the change-of-base matrix r.
// If mode = 't', the method performs the operation  a = r(transp) * a .
// If mode = 'n', the method performs the operation  a = r * a .
{
    FloatArray rta;

    if ( mode == 't' ) {
        rta.beTProductOf(r, * this);
    } else if ( mode == 'n' ) {
        rta.beProductOf(r, * this);
    } else {
        OOFEM_ERROR("FloatArray :: rotatedWith: unsupported mode");
    }

    * this = rta;
}

void
LIBeam3dNL :: computeTempCurv(FloatArray &answer, TimeStep *tStep)
{
    Material *mat = this->giveMaterial();
    IntegrationRule *iRule = integrationRulesArray [ giveDefaultIntegrationRule() ];
    GaussPoint *gp = iRule->getIntegrationPoint(0);
    ;
    FloatArray ui(3), xd(3), curv(3), ac(3), PrevEpsilon;
    FloatMatrix sc(3, 3), tmid(3, 3);

    answer.resize(3);

    // update curvature at midpoint
    // first, compute Tmid
    // ask increments
    this->computeVectorOf(EID_MomentumBalance, VM_Incremental, tStep, ui);

    ac.at(1) = 0.5 * ( ui.at(10) - ui.at(4) );
    ac.at(2) = 0.5 * ( ui.at(11) - ui.at(5) );
    ac.at(3) = 0.5 * ( ui.at(12) - ui.at(6) );
    this->computeSMtrx(sc, ac);
    sc.times(1. / 2.);
    // compute I+sc
    sc.at(1, 1) += 1.0;
    sc.at(2, 2) += 1.0;
    sc.at(3, 3) += 1.0;
    tmid.beProductOf(sc, this->tc);

    // update curvature at centre
    ac.at(1) = ( ui.at(10) - ui.at(4) );
    ac.at(2) = ( ui.at(11) - ui.at(5) );
    ac.at(3) = ( ui.at(12) - ui.at(6) );

    answer.beTProductOf(tmid, ac);
    answer.times(1 / this->l0);

    // ask for previous kappa
    PrevEpsilon = ( ( StructuralMaterialStatus * ) mat->giveStatus(gp) )->giveStrainVector();
    if ( PrevEpsilon.giveSize() ) {
        answer.at(1) += PrevEpsilon.at(4);
        answer.at(2) += PrevEpsilon.at(5);
        answer.at(3) += PrevEpsilon.at(6);
    }
}

void PrescribedGradientBCWeak :: computeIntForceGPContrib(FloatArray &oContrib_disp, IntArray &oDisp_loc_array, FloatArray &oContrib_trac, IntArray &oTrac_loc_array,TracSegArray &iEl, GaussPoint &iGP, int iDim, TimeStep *tStep, const FloatArray &iBndCoord, const double &iScaleFac, ValueModeType mode, CharType type, const UnknownNumberingScheme &s)
{

    SpatialLocalizer *localizer = domain->giveSpatialLocalizer();

	FloatMatrix contrib;
	assembleTangentGPContributionNew(contrib, iEl, iGP, iScaleFac, iBndCoord);

    // Compute vector of traction unknowns
    FloatArray tracUnknowns;
    iEl.mFirstNode->giveUnknownVector(tracUnknowns, giveTracDofIDs(), mode, tStep);

    iEl.giveTractionLocationArray(oTrac_loc_array, type, s);

    FloatArray dispElLocCoord, closestPoint;
    Element *dispEl = localizer->giveElementClosestToPoint(dispElLocCoord, closestPoint, iBndCoord );

    // Compute vector of displacement unknowns
    FloatArray dispUnknowns;
    int numDMan = dispEl->giveNumberOfDofManagers();
    for(int i = 1; i <= numDMan; i++) {
    	FloatArray nodeUnknowns;
    	DofManager *dMan = dispEl->giveDofManager(i);

    	IntArray dispIDs = giveRegularDispDofIDs();
        if(domain->hasXfemManager()) {
        	XfemManager *xMan = domain->giveXfemManager();
        	dispIDs.followedBy(xMan->giveEnrichedDofIDs(*dMan));
        }

        dMan->giveUnknownVector(nodeUnknowns, dispIDs,mode, tStep);
        dispUnknowns.append(nodeUnknowns);

    }

    dispEl->giveLocationArray(oDisp_loc_array, s);


    oContrib_disp.beTProductOf(contrib, tracUnknowns);
    oContrib_disp.negated();

    oContrib_trac.beProductOf(contrib, dispUnknowns);
    oContrib_trac.negated();
}