C++ FloatArray类代码示例

void
BoundaryLoad :: computeValueAt(FloatArray &answer, TimeStep *tStep, FloatArray &coords, ValueModeType mode)
{
    // Evaluates the value at specific integration point
    int i, j, nSize;
    double value, factor;
    FloatArray N;

    if ( ( mode != VM_Total ) && ( mode != VM_Incremental ) ) {
        OOFEM_ERROR("unknown mode");
    }

    answer.resize(this->nDofs);

    this->computeNArray(N, coords);
    nSize = N.giveSize();

    if ( ( this->componentArray.giveSize() / nSize ) != nDofs ) {
        OOFEM_ERROR("componentArray size mismatch");
    }

    for ( i = 1; i <= nDofs; i++ ) {
        for ( value = 0., j = 1; j <= nSize; j++ ) {
            value += N.at(j) * this->componentArray.at(i + ( j - 1 ) * nDofs);
        }

        answer.at(i) = value;
    }

    // time distribution

    /*
     * factor = this -> giveTimeFunction() -> at(tStep->giveTime()) ;
     * if ((mode==VM_Incremental) && (!tStep->isTheFirstStep()))
     * //factor -= this->giveTimeFunction()->at(tStep->givePreviousStep()->giveTime()) ;
     * factor -= this->giveTimeFunction()->at(tStep->giveTime()-tStep->giveTimeIncrement());
     */
    factor = this->giveTimeFunction()->evaluate(tStep, mode);

    answer.times(factor);
}

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

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

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

    force.resize(0);
    if ( gravity.giveSize() ) {
        IntegrationRule *ir = integrationRulesArray [ 0 ]; ///@todo Other/higher integration for lumped mass matrices perhaps?
        for ( int i = 0; i < ir->getNumberOfIntegrationPoints(); ++i) {
            gp = ir->getIntegrationPoint(i);

            this->interp_lin.evalN(n, *gp->giveCoordinates(), FEIElementGeometryWrapper(this));
            dV = this->computeVolumeAround(gp) * this->giveCrossSection()->give(CS_Thickness);
            load = this->giveMaterial()->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.resize(0);
    }
}

double HeMoTKMaterial :: computeCapacityCoeff(MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
{
    if ( mode == Capacity_ww ) {
        return 1.0 * rho;
    } else if ( mode == Capacity_wh ) {
        return 0.0;
    } else if ( mode == Capacity_hw ) {
        TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
        FloatArray s;
        double w, t;

        s = status->giveTempField();
        if ( s.isEmpty() ) {
            OOFEM_ERROR("undefined state vector");
        }

        w = s.at(2);
        t = s.at(1);
        return get_b(w, t) * get_latent(w, t);
    } else if ( mode == Capacity_hh ) {
        TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
        FloatArray s;
        double w, t;

        s = status->giveTempField();
        if ( s.isEmpty() ) {
            OOFEM_ERROR("undefined state vector");
        }

        w = s.at(2);
        t = s.at(1);
        return get_ceff(w, t);
    } else {
        OOFEM_ERROR("Unknown MatResponseMode");
    }

    return 0.0; // to make compiler happy
}

double HeMoTKMaterial :: computeCapacityCoeff(MatResponseMode mode, GaussPoint *gp, TimeStep *atTime)
{
    if ( mode == Capacity_ww ) {
        return 1.0 * rho;
    } else if ( mode == Capacity_wh ) {
        return 0.0;
    } else if ( mode == Capacity_hw ) {
        TransportMaterialStatus *status = ( TransportMaterialStatus * ) this->giveStatus(gp);
        FloatArray s;
        double w, t;

        s = status->giveTempStateVector();
        if ( s.isEmpty() ) {
            _error("computeCapacityCoeff: undefined state vector");
        }

        w = s.at(2);
        t = s.at(1);
        return get_b(w, t) * get_latent(w, t);
    } else if ( mode == Capacity_hh ) {
        TransportMaterialStatus *status = ( TransportMaterialStatus * ) this->giveStatus(gp);
        FloatArray s;
        double w, t;

        s = status->giveTempStateVector();
        if ( s.isEmpty() ) {
            _error("computeCapacityCoeff: undefined state vector");
        }

        w = s.at(2);
        t = s.at(1);
        return get_ceff(w, t);
    } else {
        _error("Unknown MatResponseMode");
    }

    return 0.0; // to make compiler happy
}

void
LEPlic :: giveMaterialMixtureAt(FloatArray &answer, FloatArray &position)
{
    answer.resize(2);
    Element *elem = domain->giveSpatialLocalizer()->giveElementContainingPoint(position);
    LEPlicElementInterface *interface = ( LEPlicElementInterface * ) elem->giveInterface(LEPlicElementInterfaceType);
    if ( interface ) {
        Polygon pg;
        FloatArray n;
        interface->giveTempInterfaceNormal(n);
        interface->formVolumeInterfacePoly(pg, this, n, interface->giveTempLineConstant(), false);
        if ( pg.testPoint( position.at(1), position.at(2) ) ) {
            answer.at(1) = 1.0;
            answer.at(2) = 0.0;
        } else {
            answer.at(1) = 0.0;
            answer.at(2) = 1.0;
        }
    } else {
        answer.at(1) = 1.0;
        answer.at(2) = 0.0;
    }
}

void
LIBeam2dNL :: computeStrainVectorInLayer(FloatArray &answer, GaussPoint *masterGp,
                                         GaussPoint *slaveGp, TimeStep *tStep)
//
// returns full 3d strain vector of given layer (whose z-coordinate from center-line is
// stored in slaveGp) for given tStep
//
{
    FloatArray masterGpStrain;
    double layerZeta, layerZCoord, top, bottom;

    this->computeStrainVector(masterGpStrain, masterGp, tStep);
    top    = masterGp->giveElement()->giveCrossSection()->give(CS_TopZCoord);
    bottom = masterGp->giveElement()->giveCrossSection()->give(CS_BottomZCoord);
    layerZeta = slaveGp->giveCoordinate(3);
    layerZCoord = 0.5 * ( ( 1. - layerZeta ) * bottom + ( 1. + layerZeta ) * top );

    answer.resize(6); // {Exx,Eyy,Ezz,GMyz,GMzx,GMxy}
    answer.zero();

    answer.at(1) = masterGpStrain.at(1) + masterGpStrain.at(2) * layerZCoord;
    answer.at(5) = masterGpStrain.at(3);
}

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
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 :: computeStressSpaceHardeningVars(FloatArray &answer, GaussPoint *gp,
                                                      const FloatArray &strainSpaceHardeningVariables)
{
    // in full stress strain space
    int count = 0, size = this->giveSizeOfFullHardeningVarsVector(), isize, rSize;
    IntArray mask;

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

    answer.resize(size);
    StructuralMaterial :: giveVoigtSymVectorMask( mask, gp->giveMaterialMode() );
    isize = mask.giveSize();
    rSize = this->giveSizeOfReducedHardeningVarsVector(gp);

    /* kinematic hardening variables are first */
    if ( this->kinematicHardeningFlag ) {
        for ( int i = 1; i <= isize; i++ ) {
            // to be consistent with equivalent plastic strain formulation
            // we multiply by (sqrt(2.)*2./3.)
            answer.at( mask.at(i) ) = ( sqrt(2.) * 2. / 3. ) * this->kinematicModuli * strainSpaceHardeningVariables.at(i);
        }

        count = 6;
    }

    if ( this->isotropicHardeningFlag ) {
        answer.at(count + 1) = this->isotropicModuli *
                               strainSpaceHardeningVariables.at(rSize);
    }

    answer.negated();
}

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
LEPlic :: findCellLineConstant(double &p, FloatArray &fvgrad, int ie, bool coord_upd, bool temp_vof_flag)
{
    /* The line constatnt p is solved from the general non-linear function
     * F(p) = V(p)-V = 0,
     * where V(p) is the truncated volume resulting from the intersection between
     *            assumed line segment and the reference cell
     *       V    is the given Volume of Fluid ratio
     * To find zero of this function, Brent's method is been used.
     */
    Element *elem = domain->giveElement(ie);
    LEPlicElementInterface *interface = ( LEPlicElementInterface * ) elem->giveInterface(LEPlicElementInterfaceType);
    int ivert, nelemnodes = elem->giveNumberOfNodes();
    double ivof, fvi;
    double ivx, ivy, pp, target_vof;
    if ( temp_vof_flag ) {
        target_vof = interface->giveTempVolumeFraction();
    } else {
        target_vof = interface->giveVolumeFraction();
    }

    computeLEPLICVolumeFractionWrapper wrapper(interface, this, fvgrad, target_vof, coord_upd);
    /*
     * Initial part: find lower and uper bounds for Brent algorithm
     * Here lines with given interface normal are passed through each vertex
     * and corresponding volume fractions are computed. The two lines forming
     * truncation volumes that bound the actual material in the cell provide
     * upper and lower bounds for line constant
     */
    double upper_vof = 10.0, lower_vof = -10.0;
    double upper_p = 0.0, lower_p = 0.0;

    if ( temp_vof_flag ) {
        fvi = interface->giveTempVolumeFraction();
    } else {
        fvi = interface->giveVolumeFraction();
    }

    if ( ( fvi > LEPLIC_ZERO_VOF ) && ( fvi < 1.0 ) ) {
        // boundary cell


        for ( ivert = 1; ivert <= nelemnodes; ivert++ ) {
            if ( coord_upd ) {
                ivx = giveUpdatedXCoordinate( elem->giveNode(ivert)->giveNumber() );
                ivy = giveUpdatedYCoordinate( elem->giveNode(ivert)->giveNumber() );
            } else {
                ivx = elem->giveNode(ivert)->giveCoordinate(1);
                ivy = elem->giveNode(ivert)->giveCoordinate(2);
            }

            // determine line constant for vertex ivert
            pp = -( fvgrad.at(1) * ivx + fvgrad.at(2) * ivy );
            ivof = interface->computeLEPLICVolumeFraction(fvgrad, pp, this, coord_upd);
            if ( ( ( ivof - target_vof ) >= 0. ) && ( ivof < upper_vof ) ) {
                upper_vof = ivof;
                upper_p = pp;
            } else if ( ( ( target_vof - ivof ) >= 0. ) && ( ivof > lower_vof ) ) {
                lower_vof = ivof;
                lower_p = pp;
            }
        }

        if ( ( lower_vof >= 0. ) && ( upper_vof <= 1.00000000001 ) ) {
            // now use brent's method to find minima of V(p)-V function
            brent(lower_p, 0.5 * ( lower_p + upper_p ), upper_p,
                  mem_fun< computeLEPLICVolumeFractionWrapper >(& wrapper, & computeLEPLICVolumeFractionWrapper :: eval),
                  LEPLIC_BRENT_EPS, p);
            //interface->setTempLineConstant (p);


#ifdef __OOFEG
            /*
             * Polygon grp, cell;
             * // check here
             * //Polygon testvolpoly;
             * //interface->formMaterialVolumePoly(testvolpoly, this, fvgrad, p, true);
             * //double debug_vof = interface->truncateMatVolume(testvolpoly);
             *
             * //printf ("Cell %d-> target_vof is %e, debug val is %e\n", ie, target_vof, debug_vof);
             * interface->formMyVolumePoly (cell, this, true);
             * GraphicObj *goc = cell.draw(::gc[0], false);
             * EASValsSetLayer(OOFEG_DEBUG_LAYER);
             * EASValsSetColor(::gc[0].getBcIcColor());
             * EGWithMaskChangeAttributes(COLOR_MASK | LAYER_MASK, goc);
             * EMDrawGraphics (ESIModel(), goc);
             *
             * interface->formVolumeInterfacePoly (grp, this, fvgrad, p, true);
             * GraphicObj *go = grp.draw(::gc[0], true);
             * EASValsSetColor(::gc[0].getActiveCrackColor());
             * EGWithMaskChangeAttributes(COLOR_MASK | LAYER_MASK, go);
             * EMDrawGraphics (ESIModel(), go);
             * //ESIEventLoop (YES, "findCellLineConstant -> Press Ctrl-p to continue");
             */
#endif
        } else {
            OOFEM_ERROR("LEPlic::findCellLineConstant: finding lower and uper bounds of line constant value failed");
        }
    }
}

void
LEPlic :: doCellDLS(FloatArray &fvgrad, int ie, bool coord_upd, bool vof_temp_flag)
{
    int i, ineighbr, nneighbr;
    double fvi, fvk, wk, dx, dy;
    bool isBoundaryCell = false;
    LEPlicElementInterface *interface, *ineghbrInterface;
    FloatMatrix lhs(2, 2);
    FloatArray rhs(2), xi(2), xk(2);
    IntArray currCell(1), neighborList;
    ConnectivityTable *contable = domain->giveConnectivityTable();

    if ( ( interface = ( LEPlicElementInterface * ) ( domain->giveElement(ie)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
        if ( vof_temp_flag ) {
            fvi = interface->giveTempVolumeFraction();
        } else {
            fvi = interface->giveVolumeFraction();
        }

        if ( ( fvi > 0. ) && ( fvi <= 1.0 ) ) {
            // potentially boundary cell

            if ( ( fvi > 0. ) && ( fvi < 1.0 ) ) {
                isBoundaryCell = true;
            }

            /* DLS (Differential least square reconstruction)
             *
             * In the DLS method, volume fraction Taylor series expansion of vf (volume fraction)
             * is formed from each reference cell volume fraction vf at element center x(i) to each
             * cell neighbor at point x(k). The sum (vf(i)-vf(k))^2 over all immediate neighbors
             * is then minimized inthe least square sense.
             */
            // get list of neighbours to current cell including current cell
            currCell.at(1) = ie;
            contable->giveElementNeighbourList(neighborList, currCell);
            // loop over neighbors to assemble normal equations
            nneighbr = neighborList.giveSize();
            interface->giveElementCenter(this, xi, coord_upd);
            lhs.zero();
            rhs.zero();
            for ( i = 1; i <= nneighbr; i++ ) {
                ineighbr = neighborList.at(i);
                if ( ineighbr == ie ) {
                    continue;         // skip itself
                }

                if ( ( ineghbrInterface =
                          ( LEPlicElementInterface * ) ( domain->giveElement(ineighbr)->giveInterface(LEPlicElementInterfaceType) ) ) ) {
                    if ( vof_temp_flag ) {
                        fvk = ineghbrInterface->giveTempVolumeFraction();
                    } else {
                        fvk = ineghbrInterface->giveVolumeFraction();
                    }

                    if ( fvk < 1.0 ) {
                        isBoundaryCell = true;
                    }

                    ineghbrInterface->giveElementCenter(this, xk, coord_upd);
                    wk = xk.distance(xi);
                    dx = ( xk.at(1) - xi.at(1) ) / wk;
                    dy = ( xk.at(2) - xi.at(2) ) / wk;
                    lhs.at(1, 1) += dx * dx;
                    lhs.at(1, 2) += dx * dy;
                    lhs.at(2, 2) += dy * dy;

                    rhs.at(1) += ( fvi - fvk ) * dx / wk;
                    rhs.at(2) += ( fvi - fvk ) * dy / wk;
                }
            }

            if ( isBoundaryCell ) {
                // symmetry
                lhs.at(2, 1) = lhs.at(1, 2);

                // solve normal equation for volume fraction gradient
                lhs.solveForRhs(rhs, fvgrad);

                // compute unit normal
                fvgrad.normalize();
                fvgrad.negated();
#ifdef __OOFEG
                /*
                 * EASValsSetLayer(OOFEG_DEBUG_LAYER);
                 * WCRec p[2];
                 * double tx = -fvgrad.at(2), ty=fvgrad.at(1);
                 * p[0].x=xi.at(1)-tx*0.1;
                 * p[0].y=xi.at(2)-ty*0.1;
                 * p[1].x=xi.at(1)+tx*0.1;
                 * p[1].y=xi.at(2)+ty*0.1;
                 * p[0].z = p[1].z = 0.0;
                 * GraphicObj *go = CreateLine3D(p);
                 * EGWithMaskChangeAttributes(LAYER_MASK, go);
                 * EMAddGraphicsToModel(ESIModel(), go);
                 * ESIEventLoop (YES, "Cell DLS finished; Press Ctrl-p to continue");
                 */
#endif
            } else {
                fvgrad.zero();
//.........这里部分代码省略.........

void
B3SolidMaterial :: computeCharCoefficients(FloatArray &answer, double tStep)
{
    /*
     * If EmoduliMode == 0 then analysis of continuous retardation spectrum is used for
     * computing characteristic coefficients (moduli) of Kelvin chain
     * Else least-squares method is used
     */
    if ( this->EmoduliMode == 0 ) {
      
	int mu;
	double tau0, tauMu;

      // constant "n" is assumed to be equal to 0.1 (typical value)

      // modulus of elasticity of the first unit of Kelvin chain.
      // (aging elastic spring with retardation time = 0)
      double lambda0ToPowN = pow(lambda0, 0.1);
      tau0 = pow(2 * this->giveCharTime(1) / sqrt(10.0), 0.1);
      EspringVal = 1.e6 / ( q2 * log(1.0 + tau0 / lambda0ToPowN) - q2 * tau0 / ( 10.0 * lambda0ToPowN + 10.0 * tau0 ) );
      
      // evaluation of moduli of elasticity for the remaining units
      // (Solidifying kelvin units with retardation times tauMu)
      answer.resize(nUnits);
      answer.zero();
      for ( mu = 1; mu <= this->nUnits; mu++ ) {
	tauMu = pow(2 * this->giveCharTime(mu), 0.1);
        answer.at(mu) = 10.e6 * pow(1 + tauMu / lambda0ToPowN, 2) / ( log(10.0) * q2 * ( tauMu / lambda0ToPowN ) * ( 0.9 + tauMu / lambda0ToPowN ) );
        this->charTimes.at(mu) *= 1.35;
      }
      
      answer.at(nUnits) /= 1.2;   // modulus of the last unit is reduced

    } else {   // moduli computed using the least-squares method
        int rSize;
        double taui, tauj, tti, ttj;
        FloatArray rhs(this->nUnits);
        FloatMatrix A(this->nUnits, this->nUnits);

        const FloatArray &rTimes = this->giveDiscreteTimes();
        rSize = rTimes.giveSize();
        FloatArray discreteComplianceFunctionVal(rSize);

        // compute values of the compliance function at specified times rTimes
        // (can be done directly, since the compliance function is available)

        for ( int i = 1; i <= rSize; i++ ) {
	  discreteComplianceFunctionVal.at(i) = this->computeCreepFunction( rTimes.at(i), 0. );
        }

        // assemble the matrix of the set of linear equations
        // for computing the optimal compliances
        // !!! chartime exponents are assumed to be equal to 1 !!!
        for ( int i = 1; i <= this->nUnits; i++ ) {
            taui = this->giveCharTime(i);
            for ( int j = 1; j <= this->nUnits; j++ ) {
                tauj = this->giveCharTime(j);
                double sum = 0.;
                for ( int r = 1; r <= rSize; r++ ) {
                    tti = rTimes.at(r) / taui;
                    ttj = rTimes.at(r) / tauj;
                    sum += ( 1. - exp(-tti) ) * ( 1. - exp(-ttj) );
                }

                A.at(i, j) = sum;
            }

            // assemble rhs
            // !!! chartime exponents are assumed to be equal to 1 !!!
            double sumRhs = 0.;
            for ( int r = 1; r <= rSize; r++ ) {
                tti = rTimes.at(r) / taui;
                sumRhs += ( 1. - exp(-tti) ) * discreteComplianceFunctionVal.at(r);
            }

            rhs.at(i) = sumRhs;
        }

        // solve the linear system
        A.solveForRhs(rhs, answer);

        // convert compliances into moduli
        for ( int i = 1; i <= this->nUnits; i++ ) {
            answer.at(i) = 1.e6 / answer.at(i);
        }
    }
}

NM_Status
SubspaceIteration :: solve(SparseMtrx &a, SparseMtrx &b, FloatArray &_eigv, FloatMatrix &_r, double rtol, int nroot)
//
// this function solve the generalized eigenproblem using the Generalized
// jacobi iteration
//
{
    if ( a.giveNumberOfColumns() != b.giveNumberOfColumns() ) {
        OOFEM_ERROR("matrices size mismatch");
    }

    FloatArray temp, w, d, tt, f, rtolv, eigv;
    FloatMatrix r;
    int nn, nc1, ij = 0, is;
    double rt, art, brt, eigvt;
    FloatMatrix ar, br, vec;
    std :: unique_ptr< SparseLinearSystemNM > solver( GiveClassFactory().createSparseLinSolver(ST_Direct, domain, engngModel) );

    GJacobi mtd(domain, engngModel);
    int nc = min(2 * nroot, nroot + 8);
    nn = a.giveNumberOfColumns();
    if ( nc > nn ) {
        nc = nn;
    }

    ar.resize(nc, nc);
    ar.zero();
    br.resize(nc, nc);
    br.zero();

    //
    // creation of initial iteration vectors
    //
    nc1 = nc - 1;

    w.resize(nn);
    w.zero();
    d.resize(nc);
    d.zero();
    tt.resize(nn);
    tt.zero();
    rtolv.resize(nc);
    rtolv.zero();
    vec.resize(nc, nc);
    vec.zero();                   // eigen vectors of reduced problem

    //
    // create work arrays
    //
    r.resize(nn, nc);
    r.zero();
    eigv.resize(nc);
    eigv.zero();

    FloatArray h(nn);
    for ( int i = 1; i <= nn; i++ ) {
        h.at(i) = 1.0;
        w.at(i) = b.at(i, i) / a.at(i, i);
    }

    b.times(h, tt);
    r.setColumn(tt, 1);

    for ( int j = 2; j <= nc; j++ ) {
        rt = 0.0;
        for ( int i = 1; i <= nn; i++ ) {
            if ( fabs( w.at(i) ) >= rt ) {
                rt = fabs( w.at(i) );
                ij = i;
            }
        }

        tt.at(j) = ij;
        w.at(ij) = 0.;
        for ( int i = 1; i <= nn; i++ ) {
            if ( i == ij ) {
                h.at(i) = 1.0;
            } else {
                h.at(i) = 0.0;
            }
        }

        b.times(h, tt);
        r.setColumn(tt, j);
    } // (r = z)

# ifdef DETAILED_REPORT
    OOFEM_LOG_INFO("SubspaceIteration :: solveYourselfAt: Degrees of freedom invoked by initial vectors :\n");
    tt.printYourself();
    OOFEM_LOG_INFO("SubspaceIteration :: solveYourselfAt: initial vectors for iteration:\n");
    r.printYourself();
# endif

    //ish = 0;
    a.factorized();
    //
    // start of iteration loop
    //
    for ( int nite = 0; ; ++nite ) {               // label 100
# ifdef DETAILED_REPORT
//.........这里部分代码省略.........

void VectorAssembler :: vectorFromNodeLoad(FloatArray& vec, DofManager& dman, NodalLoad* load, TimeStep* tStep, ValueModeType mode) const { vec.clear(); }