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

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

    MaterialStatus :: printOutputAt(File, tNow);

    fprintf(File, "  strains ");
    ( ( StructuralCrossSection * )
     gp->giveCrossSection() )->giveFullCharacteristicVector(helpVec, gp, strainVector);
    n = helpVec.giveSize();
    for ( i = 1; i <= n; i++ ) {
        fprintf( File, " % .4e", helpVec.at(i) );
    }

    fprintf(File, "\n              stresses");
    ( ( StructuralCrossSection * )
     gp->giveCrossSection() )->giveFullCharacteristicVector(helpVec, gp, stressVector);

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

void
LIBeam2dNL :: computeInitialStressMatrix(FloatMatrix &answer, TimeStep *tStep)
{
    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 ( int i = 0; i < iRule->giveNumberOfIntegrationPoints(); i++ ) {
        gp = iRule->getIntegrationPoint(i);
        dV = this->computeVolumeAround(gp);
        stress = static_cast< StructuralMaterialStatus * >( mat->giveStatus(gp) )->giveStressVector();
        if ( stress.giveSize() ) {
            for ( int j = 1; j <= stress.giveSize(); 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);
                }
            }
        }
    }
}

IRResultType
AnisotropicLinearElasticMaterial :: initializeFrom(InputRecord *ir)
{
    IRResultType result;                // Required by IR_GIVE_FIELD macro

    // read the stiffness coefficients arranged by rows from the diagonal to the right (21 values)
    FloatArray stiffness;
    IR_GIVE_FIELD(ir, stiffness, _IFT_AnisotropicLinearElasticMaterial_stiff);
    if ( stiffness.giveSize() != 21 ) {
        OOFEM_ERROR( "Incorrect size of stiff - should be 21, is %d\n", stiffness.giveSize() );
    }

    // put the stiffness coefficients into a 6x6 matrix
    int k = 1;
    for ( int i = 1; i <= 6; i++ ) {
        stiffmat.at(i, i) = stiffness.at(k++);
        for ( int j = i + 1; j <= 6; j++ ) {
            stiffmat.at(i, j) = stiffmat.at(j, i) = stiffness.at(k++);
        }
    }

    // read the thermal expansion coefficients (3 values)
    IR_GIVE_FIELD(ir, alpha, _IFT_AnisotropicLinearElasticMaterial_talpha);
    if ( alpha.giveSize() == 0 ) {
        alpha.resize(6);
        alpha.zero();
    } else if ( alpha.giveSize() != 6 )     {
        OOFEM_ERROR( "Incorrect size of talpha - should be 0 or 6, is %d\n", alpha.giveSize() );
    }

    return LinearElasticMaterial :: initializeFrom(ir);
}

void
StructuralFE2Material :: giveRealStressVector_3d(FloatArray &answer, GaussPoint *gp,
                                 const FloatArray &totalStrain, TimeStep *tStep)
{
    FloatArray stress;
    StructuralFE2MaterialStatus *ms = static_cast< StructuralFE2MaterialStatus * >( this->giveStatus(gp) );

#if 0
    XfemStructureManager *xMan = dynamic_cast<XfemStructureManager*>( ms->giveRVE()->giveDomain(1)->giveXfemManager() );
    if(xMan) {
        printf("Total crack length in RVE: %e\n", xMan->computeTotalCrackLength() );
    }
#endif

    ms->setTimeStep(tStep);
    // Set input
    ms->giveBC()->setPrescribedGradientVoigt(totalStrain);
    // Solve subscale problem
    ms->giveRVE()->solveYourselfAt(tStep);
    // Post-process the stress
    ms->giveBC()->computeField(stress, tStep);

    if ( stress.giveSize() == 6 ) {
        answer = stress;
    } if ( stress.giveSize() == 9 ) {
        answer = {stress[0], stress[1], stress[2], 0.5*(stress[3]+stress[6]), 0.5*(stress[4]+stress[7]), 0.5*(stress[5]+stress[8])};
    } else {
        StructuralMaterial::giveFullSymVectorForm(answer, stress, gp->giveMaterialMode() );
    }

    // Update the material status variables
    ms->letTempStressVectorBe(answer);
    ms->letTempStrainVectorBe(totalStrain);
    ms->markOldTangent(); // Mark this so that tangent is reevaluated if they are needed.
}

void
NonlinearFluidMaterial :: computeDeviatoricStressVector(FloatArray &answer, GaussPoint *gp, const FloatArray &eps, TimeStep *tStep)
{
    NonlinearFluidMaterialStatus *status = static_cast< NonlinearFluidMaterialStatus * >( this->giveStatus(gp) );

    double normeps2;

    answer = eps;
    if ( eps.giveSize() == 3 ) {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + 0.5 * ( eps.at(3) * eps.at(3) );
        answer.at(3) *= 0.5;
    } else if ( eps.giveSize() == 4 ) {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + eps.at(3) * eps.at(3) + 0.5 * ( eps.at(4) * eps.at(4) );
        answer.at(4) *= 0.5;
    } else {
        normeps2 = eps.at(1) * eps.at(1) + eps.at(2) * eps.at(2) + eps.at(3) * eps.at(3) + 0.5 * ( eps.at(4) * eps.at(4) + eps.at(5) * eps.at(5) +  eps.at(6) * eps.at(6) );
        answer.at(4) *= 0.5;
        answer.at(5) *= 0.5;
        answer.at(6) *= 0.5;
    }

    answer.times( 2.0 * viscosity * ( 1.0 + c * pow(normeps2, alpha * 0.5) ) );

    status->letTempDeviatoricStressVectorBe(answer);
    status->letTempDeviatoricStrainVectorBe(eps);
    status->letTempStrainNorm2Be(normeps2);
}

void
PetscSparseMtrx :: times(const FloatArray &x, FloatArray &answer) const
{
    if ( this->giveNumberOfColumns() != x.giveSize() ) {
        OOFEM_ERROR("Dimension mismatch");
    }

#ifdef __PARALLEL_MODE
    if ( emodel->isParallel() ) {
        OOFEM_ERROR("PetscSparseMtrx :: times - Not implemented");
    }
#endif
    Vec globX, globY;
    VecCreateSeqWithArray(PETSC_COMM_SELF, 1, x.giveSize(), x.givePointer(), & globX);
    VecCreate(PETSC_COMM_SELF, & globY);
    VecSetType(globY, VECSEQ);
    VecSetSizes(globY, PETSC_DECIDE, this->nRows);

    MatMult(this->mtrx, globX, globY);
    double *ptr;
    VecGetArray(globY, & ptr);
    answer.resize(this->nRows);
    for ( int i = 0; i < this->nRows; i++ ) {
        answer(i) = ptr [ i ];
    }

    VecRestoreArray(globY, & ptr);
    VecDestroy(&globX);
    VecDestroy(&globY);
}

void
MatlabExportModule :: doOutputHomogenizeDofIDs(TimeStep *tStep,    FILE *FID)
{

    std :: vector <FloatArray*> HomQuantities;
    double Vol = 0.0;

    // Initialize vector of arrays constaining homogenized quantities
    HomQuantities.resize(internalVarsToExport.giveSize());

    for (int j=0; j<internalVarsToExport.giveSize(); j++) {
        HomQuantities.at(j) = new FloatArray;
    }

    int nelem = this->elList.giveSize();
    for (int i = 1; i<=nelem; i++) {
        Element *e = this->emodel->giveDomain(1)->giveElement(elList.at(i));
        FEInterpolation *Interpolation = e->giveInterpolation();

        Vol = Vol + e->computeVolumeAreaOrLength();

        for ( GaussPoint *gp: *e->giveDefaultIntegrationRulePtr() ) {

            for (int j=0; j<internalVarsToExport.giveSize(); j++) {
                FloatArray elementValues;
                e->giveIPValue(elementValues, gp, (InternalStateType) internalVarsToExport(j), tStep);
                double detJ=fabs(Interpolation->giveTransformationJacobian( gp->giveNaturalCoordinates(), FEIElementGeometryWrapper(e)));

                elementValues.times(gp->giveWeight()*detJ);
                if (HomQuantities.at(j)->giveSize() == 0) {
                    HomQuantities.at(j)->resize(elementValues.giveSize());
                    HomQuantities.at(j)->zero();
                };
                HomQuantities.at(j)->add(elementValues);
            }
        }
    }


    if (noscaling) Vol=1.0;

    for ( std :: size_t i = 0; i < HomQuantities.size(); i ++) {
        FloatArray *thisIS;
        thisIS = HomQuantities.at(i);
        thisIS->times(1.0/Vol);
        fprintf(FID, "\tspecials.%s = [", __InternalStateTypeToString ( InternalStateType (internalVarsToExport(i)) ) );

        for (int j = 0; j<thisIS->giveSize(); j++) {
            fprintf(FID, "%e", thisIS->at(j+1));
            if (j!=(thisIS->giveSize()-1) ) {
                fprintf(FID, ", ");
            }
        }
        fprintf(FID, "];\n");
        delete HomQuantities.at(i);
    }

}

int
FEI3dWedgeLin :: global2local(FloatArray &answer, const FloatArray &gcoords, const FEICellGeometry &cellgeo)
{
    FloatArray res, delta, guess;
    FloatMatrix jac;
    double convergence_limit, error = 0.0;

    // find a suitable convergence limit
    convergence_limit = 1e-6 * this->giveCharacteristicLength(cellgeo);

    // setup initial guess
    answer.resize( gcoords.giveSize() );
    answer.zero();

    // apply Newton-Raphson to solve the problem
    for ( int nite = 0; nite < 10; nite++ ) {
        // compute the residual
        this->local2global(guess, answer, cellgeo);
        res.beDifferenceOf(gcoords, guess);

        // check for convergence
        error = res.computeNorm();
        if ( error < convergence_limit ) {
            break;
        }

        // compute the corrections
        this->giveJacobianMatrixAt(jac, answer, cellgeo);
        jac.solveForRhs(res, delta, true);

        // update guess
        answer.add(delta);
    }
    if ( error > convergence_limit ) { // Imperfect, could give false negatives.
        //OOFEM_ERROR("no convergence after 10 iterations");
        answer = {1. / 3., 1. / 3., 1. / 3.};
        return false;
    }

    // check limits for each local coordinate [-1,1] for quadrilaterals. (different for other elements, typically [0,1]).
    bool inside = true;
    for ( int i = 1; i <= answer.giveSize(); i++ ) {
        if ( answer.at(i) < ( 0. - POINT_TOL ) ) {
            answer.at(i) = 0.;
            inside = false;
        } else if ( answer.at(i) > ( 1. + POINT_TOL ) ) {
            answer.at(i) = 1.;
            inside = false;
        }
    }

    return inside;
}

NM_Status
IMLSolver :: solve(SparseMtrx &A, FloatArray &b, FloatArray &x)
{
    int result;

    if ( x.giveSize() != b.giveSize() ) {
        OOFEM_ERROR("size mismatch");
    }


    // check preconditioner
    if ( M ) {
        if ( ( precondInit ) || ( Lhs != &A ) || ( this->lhsVersion != A.giveVersion() ) ) {
            M->init(A);
        }
    } else {
        OOFEM_ERROR("preconditioner creation error");
    }

    Lhs = &A;
    this->lhsVersion = A.giveVersion();

#ifdef TIME_REPORT
    Timer timer;
    timer.startTimer();
#endif


    if ( solverType == IML_ST_CG ) {
        int mi = this->maxite;
        double t = this->tol;
        result = CG(* Lhs, x, b, * M, mi, t);
        OOFEM_LOG_INFO("CG(%s): flag=%d, nite %d, achieved tol. %g\n", M->giveClassName(), result, mi, t);
    } else if ( solverType == IML_ST_GMRES ) {
        int mi = this->maxite, restart = 100;
        double t = this->tol;
        FloatMatrix H(restart + 1, restart); // storage for upper Hesenberg
        result = GMRES(* Lhs, x, b, * M, H, restart, mi, t);
        OOFEM_LOG_INFO("GMRES(%s): flag=%d, nite %d, achieved tol. %g\n", M->giveClassName(), result, mi, t);
    } else {
        OOFEM_ERROR("unknown lsover type");
    }

#ifdef TIME_REPORT
    timer.stopTimer();
    OOFEM_LOG_INFO( "IMLSolver info: user time consumed by solution: %.2fs\n", timer.getUtime() );
#endif


    //solved = 1;
    return NM_Success;
}

bool
PolylineNonlocalBarrier :: isActivated(const FloatArray &c1, const FloatArray &c2)
{
    int size = vertexNodes.giveSize();
    int indx;
    double xc1, xc2, xa, xb, yc1, yc2, ya, yb;
    double a11, a12, a21, a22, b1, b2, det, t, s;
    Node *A, *B;

    int mci = max(localXCoordIndx, localYCoordIndx);
    if ( ( c1.giveSize() > mci ) || ( c2.giveSize() > mci ) ) {
        _error("PolylineNonlocalBarrier::isActivated: local coordinate index size violation");
    }

    xc1 = c1.at(localXCoordIndx);
    yc1 = c1.at(localYCoordIndx);
    xc2 = c2.at(localXCoordIndx);
    yc2 = c2.at(localYCoordIndx);

    for ( indx = 1; indx < size; indx++ ) {
        A = domain->giveNode( vertexNodes.at(indx) );
        B = domain->giveNode( vertexNodes.at(indx + 1) );

        xa = A->giveCoordinate(localXCoordIndx);
        ya = A->giveCoordinate(localYCoordIndx);
        xb = B->giveCoordinate(localXCoordIndx);
        yb = B->giveCoordinate(localYCoordIndx);

        a11 = xc2 - xc1;
        a12 = xa - xb;
        a21 = yc2 - yc1;
        a22 = ya - yb;
        b1  = xa - xc1;
        b2  = ya - yc1;
        det = a11 * a22 - a21 * a12;
        if ( det == 0 ) {
            continue;
        }

        t = ( b1 * a22 - b2 * a12 ) / det;
        if ( t < 0. || t > 1. ) {
            continue;
        }

        s = ( -b1 * a21 + b2 * a11 ) / det;
        if ( s >= 0. && s <= 1. ) {
            return true;
        }
    }

    return false;
}

void FloatArray :: beVectorProductOf(const FloatArray &v1, const FloatArray &v2)
{
    // check proper bounds
    if ( ( v1.giveSize() != 3 ) || ( v2.giveSize() != 3 ) ) {
        OOFEM_ERROR(" FloatArray::VectorProduct : size mismatch, size is not equal to 3");
    }

    this->resize(3);

    this->at(1) = v1.at(2) * v2.at(3) - v1.at(3) * v2.at(2);
    this->at(2) = v1.at(3) * v2.at(1) - v1.at(1) * v2.at(3);
    this->at(3) = v1.at(1) * v2.at(2) - v1.at(2) * v2.at(1);
}

double dot(const FloatArray &x, const FloatArray &y)
{
    //  Check for compatible dimensions:
    if ( x.giveSize() != y.giveSize() ) {
        OOFEM_ERROR("dot : incompatible dimensions");
    }

    double temp =  0;
    for ( int i = 0; i < x.giveSize(); i++ ) {
        temp += x(i) * y(i);
    }

    return temp;
}

FloatArray operator - ( const FloatArray & x, const FloatArray & y )
{
    int N = x.giveSize();
    if ( N != y.giveSize() ) {
        OOFEM_ERROR("FloatArray operator- : incompatible vector lengths");
    }

    FloatArray result(N);
    for ( int i = 0; i < N; i++ ) {
        result(i) = x(i) - y(i);
    }

    return result;
}

void 
StructuralInterfaceMaterial::giveReducedJump(FloatArray &answer, FloatArray &jump, const int size )
{
    if(size == jump.giveSize() ) {
        answer = jump;
    } else if( (size == 2) && (jump.giveSize() == 3) ) {
        answer = { jump.at(1), jump.at(3) };
    }
    else if(size == 1 && ( jump.giveSize( ) == 3 )) {
        answer = FloatArray{ jump.at(3) };
    }else{
        OOFEM_ERROR( "size must be 1, 2 or 3" );
    }

}

void
J2Mat :: computeReducedHardeningVarsLamGradient(FloatMatrix &answer, GaussPoint *gp, int actSurf,
                                                const IntArray &activeConditionMap,
                                                const FloatArray &fullStressVector,
                                                const FloatArray &strainSpaceHardeningVars,
                                                const FloatArray &gamma)
{
    int size = this->giveSizeOfReducedHardeningVarsVector(gp);
    answer.resize(size, 1);

    if ( this->kinematicHardeningFlag ) {
        int i, rsize;
        FloatArray loadGradSigVec;
        this->computeReducedStressGradientVector(loadGradSigVec, loadFunction, 1, gp, fullStressVector,
                                                 strainSpaceHardeningVars);
        rsize = loadGradSigVec.giveSize();
        for ( i = 1; i <= rsize; i++ ) {
            answer.at(i, 1) = loadGradSigVec.at(i);
        }

        answer.times( sqrt(2.) * ( 2. / 3. ) );
    }

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