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

int
LinQuad3DPlaneStress :: giveIPValue(FloatArray &answer, GaussPoint *gp, InternalStateType type, TimeStep *tStep)
{
    FloatMatrix globTensor;
    CharTensor cht;

    answer.resize(6);

    if ( type == IST_ShellForceTensor || type == IST_ShellStrainTensor ) {
        double c = 1.0;
        if ( type == IST_ShellForceTensor ) {
            cht = GlobalForceTensor;
        } else {
            cht = GlobalStrainTensor;
            c = 2.0;
        }

        this->giveCharacteristicTensor(globTensor, cht, gp, tStep);

        answer.at(1) = globTensor.at(1, 1); //xx
        answer.at(2) = globTensor.at(2, 2); //yy
        answer.at(3) = globTensor.at(3, 3); //zz
        answer.at(4) = c * globTensor.at(2, 3); //yz
        answer.at(5) = c * globTensor.at(1, 3); //xz
        answer.at(6) = c * globTensor.at(2, 3); //yz

        return 1;
    } else if ( type == IST_ShellMomentTensor || type == IST_ShellCurvatureTensor ) {
        answer.clear();
        return 1;
    } else {
        answer.clear();
        return 0;
    }
}

void
FEI3dTrQuad :: surfaceEvalBaseVectorsAt(FloatArray &G1, FloatArray &G2, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    // Note: These are not normalized. Returns the two tangent vectors to the surface.
    FloatMatrix dNdxi;
    this->surfaceEvaldNdxi(dNdxi, lcoords);

    G1.clear();
    G2.clear();
    for ( int i = 0; i < 6; ++i ) {
        G1.add( dNdxi(i, 1), * cellgeo.giveVertexCoordinates(i) );
        G2.add( dNdxi(i, 2), * cellgeo.giveVertexCoordinates(i) );
    }
}

void
RankinePlasticMaterial :: computeStressSpaceHardeningVarsReducedGradient(FloatArray &answer, functType ftype, int isurf, GaussPoint *gp,
                                                                         const FloatArray &stressVector,
                                                                         const FloatArray &stressSpaceHardeningVars)
{
    answer.clear();
}

void
POIExportModule :: exportIntVarAs(InternalStateType valID, FILE *stream, TimeStep *tStep)
{
    int i, region;
    IntArray toMap(1);
    Domain *d = emodel->giveDomain(1);
    FloatArray poiCoords(3);
    FloatArray val;

    toMap.at(1) = ( int ) valID;

    // loop over POIs
    for ( auto &poi: POIList ) {
        poiCoords.at(1) = poi.x;
        poiCoords.at(2) = poi.y;
        poiCoords.at(3) = poi.z;
        region = poi.region;

        this->giveMapper()->__init(d, toMap, poiCoords, * d->giveSet(region), tStep);
        if ( !this->giveMapper()->__mapVariable(val, poiCoords, valID, tStep) ) {
            OOFEM_WARNING("Failed to map variable");
            val.clear();
        }
        fprintf(stream, "%10d ", poi.id);
        for ( i = 1; i <= val.giveSize(); i++ ) {
            fprintf( stream, " %15e", val.at(i) );
        }

        fprintf(stream, "\n");
    }
}

void
SimpleCrossSection :: giveTemperatureVector(FloatArray &answer, GaussPoint *gp, TimeStep *tStep)
{
    Element *elem = gp->giveElement();
    answer.clear();
    //sum up all prescribed temperatures over an element
    StructuralElement *selem = dynamic_cast< StructuralElement * >(elem);
    selem->computeResultingIPTemperatureAt(answer, tStep, gp, VM_Total);

    /* add external source, if provided */
    FieldManager *fm = this->domain->giveEngngModel()->giveContext()->giveFieldManager();
    FM_FieldPtr tf;

    if ( ( tf = fm->giveField(FT_Temperature) ) ) {
        // temperature field registered
        FloatArray gcoords, et2;
        int err;
        elem->computeGlobalCoordinates( gcoords, gp->giveNaturalCoordinates() );
        if ( ( err = tf->evaluateAt(et2, gcoords, VM_Total, tStep) ) ) {
            OOFEM_ERROR("tf->evaluateAt failed, element %d, error code %d", elem->giveNumber(), err);
        }

        if ( et2.isNotEmpty() ) {
            if ( answer.isEmpty() ) {
                answer = et2;
            } else {
                answer.at(1) += et2.at(1);
            }
        }
    }
}

void
SUPGElement2 :: computeAdvectionTerm_MB(FloatArray &answer, TimeStep *tStep)
{
    FloatMatrix n, b, bn;
    FloatArray u, v;

    answer.clear();

    this->computeVectorOfVelocities(VM_Total, tStep, u);

    int rule = 2;
    /* consistent part + supg stabilization term */
    for ( GaussPoint *gp: *this->integrationRulesArray [ rule ] ) {
        this->computeNuMatrix(n, gp);
        this->computeUDotGradUMatrix( bn, gp, tStep->givePreviousStep() );
        this->computeUDotGradUMatrix(b, gp, tStep);
        v.beProductOf(b, u);
        double dV = this->computeVolumeAround(gp);
        double rho = this->giveMaterial()->give('d', gp);
        /* consistent part */
        answer.plusProduct(n, v, rho * dV);

        /* supg stabilization */
        answer.plusProduct(bn, v, t_supg * rho * dV);
    }
}

void MixedGradientPressureWeakPeriodic :: integrateTractionDev(FloatArray &answer, Element *el, int boundary, const FloatMatrix &ddev)
{
    // Computes the integral: int dt . dx dA
    FloatMatrix mMatrix;
    FloatArray normal, coords, vM_dev;

    FEInterpolation *interp = el->giveInterpolation(); // Geometry interpolation. The displacements or velocities must have the same interpolation scheme (on the boundary at least).

    int maxorder = this->order + interp->giveInterpolationOrder() * 3;
    std :: unique_ptr< IntegrationRule >ir( interp->giveBoundaryIntegrationRule(maxorder, boundary) );
    answer.clear();

    for ( GaussPoint *gp: *ir ) {
        const FloatArray &lcoords = gp->giveNaturalCoordinates();
        FEIElementGeometryWrapper cellgeo(el);

        double detJ = interp->boundaryEvalNormal(normal, boundary, lcoords, cellgeo);
        // Compute v_m = d_dev . x
        interp->boundaryLocal2Global(coords, boundary, lcoords, cellgeo);
        vM_dev.beProductOf(ddev, coords);

        this->constructMMatrix(mMatrix, coords, normal);

        answer.plusProduct(mMatrix, vM_dev, detJ * gp->giveWeight());
    }
}

void Tr21Stokes :: NodalAveragingRecoveryMI_computeNodalValue(FloatArray &answer, int node, InternalStateType type, TimeStep *tStep)
{
    if ( type == IST_Pressure ) {
        answer.resize(1);
        if ( node == 1 || node == 2 || node == 3 ) {
            answer.at(1) = this->giveNode(node)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
        } else {
            double a, b;
            if ( node == 4 ) {
                a = this->giveNode(1)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
                b = this->giveNode(2)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
            } else if ( node == 5 ) {
                a = this->giveNode(2)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
                b = this->giveNode(3)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
            } else { /*if ( node == 6 )*/
                a = this->giveNode(3)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
                b = this->giveNode(1)->giveDofWithID(P_f)->giveUnknown(VM_Total, tStep);
            }

            answer.at(1) = ( a + b ) / 2;
        }
    } else {
        answer.clear();
    }
}

void
MisesMatNl :: giveRemoteNonlocalStiffnessContribution(GaussPoint *gp, IntArray &rloc, const UnknownNumberingScheme &s,
                                                      FloatArray &rcontrib, TimeStep *tStep)
{
    double kappa, tempKappa;
    MisesMatNlStatus *status = static_cast< MisesMatNlStatus * >( this->giveStatus(gp) );
    StructuralElement *elem = static_cast< StructuralElement * >( gp->giveElement() );
    FloatMatrix b;
    LinearElasticMaterial *lmat = this->giveLinearElasticMaterial();
    double E = lmat->give('E', gp);

    elem->giveLocationArray(rloc, s);
    elem->computeBmatrixAt(gp, b);
    kappa = status->giveCumulativePlasticStrain();
    tempKappa = status->giveTempCumulativePlasticStrain();

    rcontrib.clear();
    if ( ( tempKappa - kappa ) > 0 ) {
        const FloatArray &stress = status->giveTempEffectiveStress();
        if ( gp->giveMaterialMode() == _1dMat ) {
            double coeff = sgn( stress.at(1) ) * E / ( E + H );
            rcontrib.plusProduct(b, stress, coeff);
            return;
        }
    }
    rcontrib.resize(b.giveNumberOfColumns());
    rcontrib.zero();
}

void
StructuralInterfaceElementPhF :: giveInternalForcesVector_u(FloatArray &answer, TimeStep *tStep, int useUpdatedGpRecord)
{

    
    FloatMatrix N;
    FloatArray u, traction, jump;
    
    IntArray dofIdArray;
    this->giveDofManDofIDMask_u(dofIdArray);
    this->computeVectorOf(dofIdArray, VM_Total, tStep, u);

    // subtract initial displacements, if defined
    if ( initialDisplacements.giveSize() ) {
        u.subtract(initialDisplacements);
    }
    
    // zero answer will resize accordingly when adding first contribution
    answer.clear();
    
    for ( auto &ip: *this->giveDefaultIntegrationRulePtr() ) {
        
        this->computeNmatrixAt(ip, N);
        
        jump.beProductOf(N, u);
        this->computeTraction(traction, ip, jump, tStep);
        //traction.resize(2);
        
        // compute internal cohesive forces as f = N^T*traction dA
        double dA = this->computeAreaAround(ip);
        answer.plusProduct(N, traction, dA);
    }

}

void Tr21Stokes :: computeLoadVector(FloatArray &answer, BodyLoad *load, CharType type, ValueModeType mode, TimeStep *tStep)
{
    if ( type != ExternalForcesVector ) {
        answer.clear();
        return;
    }

    FluidDynamicMaterial *mat = static_cast< FluidCrossSection * >( this->giveCrossSection() )->giveFluidMaterial();
    FloatArray N, gVector, temparray(12);

    load->computeComponentArrayAt(gVector, tStep, VM_Total);
    temparray.zero();
    if ( gVector.giveSize() ) {
        for ( GaussPoint *gp: *integrationRulesArray [ 0 ] ) {
            const FloatArray &lcoords = gp->giveNaturalCoordinates();

            double rho = mat->give('d', gp);
            double detJ = fabs( this->interpolation_quad.giveTransformationJacobian( lcoords, FEIElementGeometryWrapper(this) ) );
            double dA = detJ * gp->giveWeight();

            this->interpolation_quad.evalN( N, lcoords, FEIElementGeometryWrapper(this) );
            for ( int j = 0; j < 6; j++ ) {
                temparray(2 * j)     += N(j) * rho * gVector(0) * dA;
                temparray(2 * j + 1) += N(j) * rho * gVector(1) * dA;
            }
        }
    }

    answer.resize(15);
    answer.zero();
    answer.assemble(temparray, this->momentum_ordering);
}

void Tr21Stokes :: computeExternalForcesVector(FloatArray &answer, TimeStep *tStep)
{
    FloatArray vec;

    answer.clear();

    int nLoads = this->boundaryLoadArray.giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {  // For each Neumann boundary condition
        int load_number = this->boundaryLoadArray.at(2 * i - 1);
        int load_id = this->boundaryLoadArray.at(2 * i);
        Load *load = this->domain->giveLoad(load_number);
        bcGeomType ltype = load->giveBCGeoType();

        if ( ltype == EdgeLoadBGT ) {
            this->computeBoundarySurfaceLoadVector(vec, static_cast< BoundaryLoad * >(load), load_id, ExternalForcesVector, VM_Total, tStep);
            answer.add(vec);
        }
    }

    BodyLoad *bload;
    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load = domain->giveLoad( bodyLoadArray.at(i) );
	if ((bload = dynamic_cast<BodyLoad*>(load))) {
	  bcGeomType ltype = load->giveBCGeoType();
	  if ( ltype == BodyLoadBGT && load->giveBCValType() == ForceLoadBVT ) {
            this->computeLoadVector(vec, bload, ExternalForcesVector, VM_Total, tStep);
            answer.add(vec);
	  }
        }
    }
}

void
SUPGElement2 :: computeBCRhsTerm_MB(FloatArray &answer, TimeStep *tStep)
{
    int nLoads;

    answer.clear();

    int rule = 0;
    IntegrationRule *iRule = this->integrationRulesArray [ rule ];
    FloatArray un, gVector, s, helpLoadVector;
    FloatMatrix b, nu;

    // add body load (gravity) termms
    nLoads = this->giveBodyLoadArray()->giveSize();
    for ( int i = 1; i <= nLoads; i++ ) {
        Load *load = domain->giveLoad( bodyLoadArray.at(i) );
        bcGeomType ltype = load->giveBCGeoType();
        if ( ( ltype == BodyLoadBGT ) && ( load->giveBCValType() == ForceLoadBVT ) ) {
            load->computeComponentArrayAt(gVector, tStep, VM_Total);
            if ( gVector.giveSize() ) {
                for ( GaussPoint *gp: *iRule ) {
                    this->computeUDotGradUMatrix( b, gp, tStep->givePreviousStep() );
                    this->computeNuMatrix(nu, gp);
                    double dV  = this->computeVolumeAround(gp);
                    double rho = this->giveMaterial()->give('d', gp);
                    answer.plusProduct(b, gVector, t_supg * rho * dV);
                    answer.plusProduct(nu, gVector, rho * dV);
                }
            }
        }
    }

    // integrate tractions
    // if no traction bc applied but side marked as with traction load
    // then zero traction is assumed !!!

    // loop over boundary load array
    nLoads = this->giveBoundaryLoadArray()->giveSize() / 2;
    for ( int i = 1; i <= nLoads; i++ ) {
        int n = boundaryLoadArray.at(1 + ( i - 1 ) * 2);
        int id = boundaryLoadArray.at(i * 2);
        Load *load  = domain->giveLoad(n);
        bcGeomType ltype = load->giveBCGeoType();
        if ( ltype == EdgeLoadBGT ) {
            this->computeEdgeLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else if ( ltype == SurfaceLoadBGT ) {
            this->computeSurfaceLoadVector_MB(helpLoadVector, load, id, tStep);
            if ( helpLoadVector.giveSize() ) {
                answer.add(helpLoadVector);
            }
        } else {
            OOFEM_ERROR("unsupported load type class");
        }
    }
}

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
FEI3dTrQuad :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
{
    FloatArray n;
    this->evalN(n, lcoords, cellgeo);
    answer.clear();
    for ( int i = 1; i <= 6; ++i ) {
        answer.add( n.at(i), * cellgeo.giveVertexCoordinates(i) );
    }
}