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

void
FEI2dTrQuad :: edgeEvaldNds(FloatArray &answer, int iedge,
                            const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    // I think it at least should return both dNds and J. Both are almost always needed.
    // In fact, dxdxi is also needed sometimes (surface tension)
#if 0
    IntArray edgeNodes;
    FloatArray dNdxi(3);
    FloatArray dxdxi(2);
    double xi = lcoords.at(1);
    this->computeLocalEdgeMapping(edgeNodes, iedge);
    dNdxi.at(1) = xi - 0.5;
    dNdxi.at(2) = xi + 0.5;
    dNdxi.at(3) = -2 * xi;

    dxdxi.at(1) = dNdxi.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(xind) +
    dNdxi.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(xind) +
    dNdxi.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(xind);
    dxdxi.at(2) = dNdxi.at(1) * cellgeo.giveVertexCoordinates( edgeNodes.at(1) )->at(yind) +
    dNdxi.at(2) * cellgeo.giveVertexCoordinates( edgeNodes.at(2) )->at(yind) +
    dNdxi.at(3) * cellgeo.giveVertexCoordinates( edgeNodes.at(3) )->at(yind);

    double J = dxdxi.computeNorm();
    answer = dNdxi;
    answer.times(1 / J);
    return J;

#endif
    double xi = lcoords.at(1);
    double J = edgeGiveTransformationJacobian(iedge, lcoords, cellgeo);
    answer = {
        ( xi - 0.5 ) / J,
        ( xi + 0.5 ) / J,
        -2 * xi / J
    };
}

void
FEI3dHexaLin :: giveLocalDerivative(FloatMatrix &dN, const FloatArray &lcoords)
{
    double u, v, w;
    u = lcoords.at(1);
    v = lcoords.at(2);
    w = lcoords.at(3);
    
    dN.resize(8, 3);
    
    dN.at(1, 1) = -0.125 * ( 1. - v ) * ( 1. + w );
    dN.at(2, 1) = -0.125 * ( 1. + v ) * ( 1. + w );
    dN.at(3, 1) =  0.125 * ( 1. + v ) * ( 1. + w );
    dN.at(4, 1) =  0.125 * ( 1. - v ) * ( 1. + w );
    dN.at(5, 1) = -0.125 * ( 1. - v ) * ( 1. - w );
    dN.at(6, 1) = -0.125 * ( 1. + v ) * ( 1. - w );
    dN.at(7, 1) =  0.125 * ( 1. + v ) * ( 1. - w );
    dN.at(8, 1) =  0.125 * ( 1. - v ) * ( 1. - w );

    dN.at(1, 2) = -0.125 * ( 1. - u ) * ( 1. + w );
    dN.at(2, 2) =  0.125 * ( 1. - u ) * ( 1. + w );
    dN.at(3, 2) =  0.125 * ( 1. + u ) * ( 1. + w );
    dN.at(4, 2) = -0.125 * ( 1. + u ) * ( 1. + w );
    dN.at(5, 2) = -0.125 * ( 1. - u ) * ( 1. - w );
    dN.at(6, 2) =  0.125 * ( 1. - u ) * ( 1. - w );
    dN.at(7, 2) =  0.125 * ( 1. + u ) * ( 1. - w );
    dN.at(8, 2) = -0.125 * ( 1. + u ) * ( 1. - w );

    dN.at(1, 3) =  0.125 * ( 1. - u ) * ( 1. - v );
    dN.at(2, 3) =  0.125 * ( 1. - u ) * ( 1. + v );
    dN.at(3, 3) =  0.125 * ( 1. + u ) * ( 1. + v );
    dN.at(4, 3) =  0.125 * ( 1. + u ) * ( 1. - v );
    dN.at(5, 3) = -0.125 * ( 1. - u ) * ( 1. - v );
    dN.at(6, 3) = -0.125 * ( 1. - u ) * ( 1. + v );
    dN.at(7, 3) = -0.125 * ( 1. + u ) * ( 1. + v );
    dN.at(8, 3) = -0.125 * ( 1. + u ) * ( 1. - v );
}

void
AxisymElement :: computeBHmatrixAt(GaussPoint *gp, FloatMatrix &answer)
// Returns the [ 9 x (nno*2) ] displacement gradient matrix {BH} of the receiver,
// evaluated at gp.
// BH matrix  -  9 rows : du/dx, dv/dy, dw/dz = u/r, 0, 0, du/dy,  0, 0, dv/dx
///@todo not checked if correct, is dw/dz = u/r for large deformations? /JB
{
    FloatArray n;
    FloatMatrix dnx;

    static_cast< FEInterpolation2d* > ( this->giveInterpolation() )->
        evaldNdx( dnx, gp->giveNaturalCoordinates(), *this->giveCellGeometryWrapper() );
    int nRows = dnx.giveNumberOfRows();
    answer.resize(9, nRows*2);
    answer.zero();

    double r = 0., x;
    for ( int i = 1; i <= this->giveNumberOfDofManagers(); i++ ) {
        x  = this->giveNode(i)->giveCoordinate(1);
        r += x * n.at(i);
    }


    // mode is _3dMat !!!!!! answer.at(4,*), answer.at(5,*), answer.at(7,*), and answer.at(8,*) is zero
    for ( int i = 1; i <= nRows*2; i++ ) {
        answer.at(1, 3 * i - 2) = dnx.at(i, 1);     // du/dx
        answer.at(2, 3 * i - 1) = dnx.at(i, 2);     // dv/dy
        answer.at(6, 3 * i - 2) = dnx.at(i, 2);     // du/dy
        answer.at(9, 3 * i - 1) = dnx.at(i, 1);     // dv/dx
    }

    
    for ( int i = 0; i < this->giveNumberOfDofManagers(); i++ ) {
        answer.at(3, 2*i + 1) = n.at(i+1) / r;
    }            

}

void
Quad1Mindlin :: computeBodyLoadVectorAt(FloatArray &answer, Load *forLoad, TimeStep *tStep, ValueModeType mode)
{
    // Only gravity load
    double dV, load;
    FloatArray force, gravity, n;

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

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

    force.clear();
    if ( gravity.giveSize() ) {
        ///@todo Other/higher integration for lumped mass matrices perhaps?
        for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {

            this->interp_lin.evalN( n, * gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(this) );
            dV = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness, gp);
            load = this->giveStructuralCrossSection()->give('d', gp) * gravity.at(3) * dV;

            force.add(load, n);
        }

        answer.resize(12);
        answer.zero();

        answer.at(1)  = force.at(1);
        answer.at(4)  = force.at(2);
        answer.at(7)  = force.at(3);
        answer.at(10) = force.at(4);
    } else {
        answer.clear();
    }
}

void
SolutionbasedShapeFunction :: computeDofTransformation(ActiveDof *dof, FloatArray &masterContribs)
{
    if ( !isLoaded ) {
        loadProblem();
    }

    FloatArray values, masterContribs2, d, values2;

    masterContribs.resize( this->giveDomain()->giveNumberOfSpatialDimensions() );
    masterContribs2.resize( this->giveDomain()->giveNumberOfSpatialDimensions() );

    IntArray dofIDs = {dof->giveDofID()};

    bool isPlus, isMinus, isZero, found;
    whichBoundary(* dof->giveDofManager()->giveCoordinates(), isPlus, isMinus, isZero);

    for ( int i = 1; i <= this->giveDomain()->giveNumberOfSpatialDimensions(); i++ ) {
        double factor = 1.0;
        found = false;

        modeStruct *ms = modes.at(i - 1);
        for ( size_t j = 0; j < ms->SurfaceData.size(); j++ ) {
            SurfaceDataStruct *sd = ms->SurfaceData.at(j);
            if ( sd->DofMan->giveNumber() == dof->giveDofManager()->giveNumber() ) {
                if ( sd->DofID == dof->giveDofID() ) {
                    values.resize(1);
                    values.at(1) = sd->value;
                    found = true;
                    break;
                }
            }
        }

        if ( !found ) {
            printf( "%u\n", dof->giveDofManager()->giveNumber() );
            OOFEM_ERROR("Node not found");
        }

        //giveValueAtPoint(values2, *dof->giveDofManager()->giveCoordinates(), dofIDs, *modes.at(i-1)->myEngngModel);
        //printf ("Mode %u, DofManager: %u, DofIDItem %u, value %10.10f\n", i, dof->giveDofManager()->giveNumber(), dof->giveDofID(), values.at(1));

        factor = isPlus  ? modes.at(i - 1)->ap : factor;
        factor = isMinus ? modes.at(i - 1)->am : factor;
        factor = isZero  ? 1.0 : factor;

        masterContribs.at(i) = factor * values.at(1);
    }
}

void
FEI2dQuadLin :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    double ksi, eta;

    ksi = lcoords.at(1);
    eta = lcoords.at(2);

    answer = {
        ( 1. + ksi ) * ( 1. + eta ) * 0.25,
        ( 1. - ksi ) * ( 1. + eta ) * 0.25,
        ( 1. - ksi ) * ( 1. - eta ) * 0.25,
        ( 1. + ksi ) * ( 1. - eta ) * 0.25
    };
}

void
IntMatCoulombContact :: giveEngTraction_3d( FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
{
    IntMatCoulombContactStatus *status = static_cast< IntMatCoulombContactStatus * >( this->giveStatus( gp ) );

    double normalJump = jump.at( 3 );
    FloatArray shearJump = { jump.at(1), jump.at(2) };

    double normalStress = 0.0;
    FloatArray shearStress, tempShearStressShift = status->giveShearStressShift();
    this->computeEngTraction( normalStress, shearStress, tempShearStressShift,
                              normalJump, shearJump );
    
    // Set stress components in the traction vector
    answer.resize( 3 );
    answer.at( 1 ) = shearStress.at( 1 );
    answer.at( 2 ) = shearStress.at( 2 );
    answer.at( 3 ) = normalStress;

    // Update gp
    status->setTempShearStressShift( tempShearStressShift );
    status->letTempJumpBe( jump );
    status->letTempTractionBe( answer );
}

void
IntMatCoulombContact :: giveEngTraction_1d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
{
    // Returns the (engineering) traction vector (normal stress only) in 1d based on the 
    // spatial jump. The shear stress is not relevant in this case.

    IntMatCoulombContactStatus *status = static_cast< IntMatCoulombContactStatus * >( this->giveStatus( gp ) );

    double normalJump = jump.at( 1 );
    FloatArray shearJump(0);

    double normalStress = 0.0;
    FloatArray shearStress, tempShearStressShift(0);
    this->computeEngTraction( normalStress, shearStress, tempShearStressShift,
                              normalJump, shearJump );

    // Set stress components in the traction vector
    answer.resize( 1 );
    answer.at( 1 ) = normalStress;

    // Update gp
    status->letTempJumpBe( jump );
    status->letTempTractionBe( answer );
}

void
Quad1MindlinShell3D :: computeSurfaceLoadVectorAt(FloatArray &answer, Load *load,
                                                int iSurf, TimeStep *tStep, ValueModeType mode)
{
    BoundaryLoad *surfLoad = static_cast< BoundaryLoad * >(load);
    if ( dynamic_cast< ConstantPressureLoad * >(surfLoad) ) { // Just checking the type of b.c.
        // EXPERIMENTAL CODE:
        IntegrationRule *iRule;
        FloatArray n, gcoords, pressure;

        answer.resize( 24 );
        answer.zero();

        //int approxOrder = surfLoad->giveApproxOrder() + this->giveApproxOrder();

        iRule = this->integrationRulesArray[ 0 ];
        for ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
            GaussPoint *gp = iRule->getIntegrationPoint(i);
            double dV = this->computeVolumeAround(gp);
            this->interp.evalN(n, *gp->giveCoordinates(), FEIVoidCellGeometry());
            this->interp.local2global(gcoords, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
            surfLoad->computeValueAt(pressure, tStep, gcoords, mode);

            answer.at( 3) += n.at(1) * pressure.at(1) * dV;
            answer.at( 9) += n.at(2) * pressure.at(1) * dV;
            answer.at(15) += n.at(3) * pressure.at(1) * dV;
            answer.at(21) += n.at(4) * pressure.at(1) * dV;
        }
        // Second surface is the outside;
        if ( iSurf == 2 ) {
            answer.negated();
        }
    } else {
        OOFEM_ERROR("Quad1MindlinShell3D only supports constant pressure boundary load.");
    }
}

int
HydratingIsoHeatMaterial :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    // printf ("IP %d::giveIPValue, IST %d", giveNumber(), type);
    if ( type == IST_HydrationDegree ) {
        //TransportMaterialStatus* status = (TransportMaterialStatus*) this -> giveStatus (gp);
        answer.resize(1);
        //if (hydration)
        answer.at(1) = giveHydrationDegree(gp, tStep, VM_Total);
        //else answer.at(1) = 0;
        return 1;
    } else {
        return TransportMaterial :: giveIPValue(answer, gp, type, tStep);
    }
}

void StructuralMaterialStatus :: printOutputAt(FILE *File, TimeStep *tNow)
// Prints the strains and stresses on the data file.
{
    FloatArray helpVec;
    int n;

    MaterialStatus :: printOutputAt(File, tNow);

    fprintf(File, "  strains ");
    StructuralMaterial :: giveFullSymVectorForm( helpVec, strainVector, gp->giveMaterialMode() );
    n = helpVec.giveSize();
    for ( int i = 1; i <= n; i++ ) {
        fprintf( File, " % .4e", helpVec.at(i) );
    }

    fprintf(File, "\n              stresses");
    StructuralMaterial :: giveFullSymVectorForm( helpVec, stressVector, gp->giveMaterialMode() );

    n = helpVec.giveSize();
    for ( int i = 1; i <= n; i++ ) {
        fprintf( File, " % .4e", helpVec.at(i) );
    }
    fprintf(File, "\n");
}

void PolygonLine :: giveInputRecord(DynamicInputRecord &input)
{
    input.setRecordKeywordField( "PolygonLine", 1 );

    FloatArray points;
    int nVert = mVertices.size();
    points.resize(nVert * 2);

    for ( int i = 0; i < nVert; i++ ) {
        points.at(2 * i + 1) = mVertices [ i ].at(1);
        points.at(2 * i + 2) = mVertices [ i ].at(2);
    }

    input.setField(points, _IFT_PolygonLine_points);
}

void
CohesiveInterfaceMaterial :: giveEngTraction_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &jump, TimeStep *tStep)
{
    StructuralInterfaceMaterialStatus *status = static_cast< StructuralInterfaceMaterialStatus * >( this->giveStatus(gp) );

    answer.resize(3);

    double x = jump.at(1) + transitionOpening;
    
    if (stiffCoeffKn == 1.){//tension stiffness = compression stiffness
        answer.at(1) = kn*x;
    } else {
        //composed function from two atan's to have smooth intersection between tension and compression
        answer.at(1) = (M_PI/2. + atan(smoothMag*x))/M_PI*kn*stiffCoeffKn*x + (M_PI/2.-atan(smoothMag*x))/M_PI*kn*x;
    }
    
    // shear part of elastic stress-strain law
    answer.at(2) = ks * jump.at(2);
    answer.at(3) = ks * jump.at(3);
    
    // update gp
    status->letTempJumpBe(jump);
    status->letTempTractionBe(answer);
}

void
J2MPlasticMaterial :: computeStressSpaceHardeningVarsReducedGradient(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp,
                                                                     const FloatArray &stressVector,
                                                                     const FloatArray &stressSpaceHardeningVars)
{
    /* computes stress space hardening gradient in reduced stress-strain space */

    int kcount = 0, size = this->giveSizeOfReducedHardeningVarsVector(gp);
    //double f,ax,ay,az,sx,sy,sz;
    FloatArray fullKinematicGradient, reducedKinematicGrad;

    if ( !hasHardening() ) {
        answer.clear();
        return;
    }

    answer.resize(size);

    /* kinematic hardening variables first */
    if ( this->kinematicHardeningFlag ) {
        this->computeStressGradientVector(fullKinematicGradient, ftype, isurf, gp, stressVector, stressSpaceHardeningVars);
        StructuralMaterial :: giveReducedSymVectorForm( reducedKinematicGrad, fullKinematicGradient, gp->giveMaterialMode() );

        kcount = reducedKinematicGrad.giveSize();
    }

    if ( this->kinematicHardeningFlag ) {
        for ( int i = 1; i <= kcount; i++ ) {
            answer.at(i) = reducedKinematicGrad.at(i);
        }
    }

    if ( this->isotropicHardeningFlag ) {
        answer.at(size) = sqrt(1. / 3.);
    }
}

int
Quad10_2D_SUPG :: giveIPValue(FloatArray &answer, GaussPoint *aGaussPoint, InternalStateType type, TimeStep *atTime)
{
    if ( type == IST_VOFFraction ) {
        MaterialInterface *mi = domain->giveEngngModel()->giveMaterialInterface( domain->giveNumber() );
        if ( mi ) {
            FloatArray val;
            mi->giveElementMaterialMixture( val, aGaussPoint->giveElement()->giveNumber() );
            answer.resize(1);
            answer.at(1) = val.at(1);
            return 1;
        } else {
            answer.resize(1);
            answer.at(1) = 1.0;
            return 1;
        }
    } else if ( type == IST_Density ) {
        answer.resize(1);
        answer.at(1) = this->giveMaterial()->give('d', aGaussPoint);
        return 1;
    } else {
        return SUPGElement :: giveIPValue(answer, aGaussPoint, type, atTime);
    }
}