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

double HeMoKunzelMaterial :: computeCapacityCoeff(MatResponseMode mode, GaussPoint *gp, TimeStep *atTime)
{
    //     if  (gp->giveElement()->giveNumber() == 4)
    //       double bzzz = 20;


    if ( mode == Capacity_ww ) {
        TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
        FloatArray s;
        double h;
        double dw_dh;

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

        h = s.at(2);
        dw_dh = this->sorptionIsothermDerivative(h);

        return dw_dh;

        // CONSTANT
        //return 10.;
    } else if ( mode == Capacity_wh ) {
        return 0.0;
    } else if ( mode == Capacity_hw ) {
        return 0.0;
    } else if ( mode == Capacity_hh ) {
        TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
        FloatArray s;
        double h, w;
        double dHs_dT, dHw_dT;

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

        h = s.at(2);

        w = this->sorptionIsotherm(h);

        dHs_dT = cs * give('d', NULL);
        dHw_dT = cw * w;

        return ( dHs_dT + dHw_dT );

        // CONSTANT	return 1.7e6;
    } else {
        OOFEM_ERROR("Unknown MatResponseMode");
    }

    return 0.0; // to make compiler happy
}

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
HeMoKunzelMaterial :: matcond1d(FloatMatrix &d, GaussPoint *gp, MatResponseMode mode, TimeStep *atTime)
{
    double k = 0.0, h = 0.0, t = 0.0;
    TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    FloatArray s;

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

    h = s.at(2);
    t = s.at(1);

    if ( mode == Conductivity_ww ) {
        k = perm_mm(h, t);
    } else if ( mode == Conductivity_wh ) {
        k = perm_mh(h, t);
    } else if ( mode == Conductivity_hw ) {
        k = perm_hm(h, t);
    } else if ( mode == Conductivity_hh ) {
        k = perm_hh(h, t);
    } else {
        OOFEM_ERROR("Unknown MatResponseMode");
    }

    d.resize(1, 1);
    d.at(1, 1) = k;
}

void
TR_SHELL02 :: giveCharacteristicVector(FloatArray &answer, CharType mtrx, ValueModeType mode, TimeStep *tStep)
//
// returns characteristics vector of receiver accordind to mtrx
//
{
    FloatArray aux;

    answer.resize(18);
    answer.zero();

    plate->giveCharacteristicVector(aux, mtrx, mode, tStep);
    if ( !aux.isEmpty() ) answer.assemble(aux, loc_plate);

    membrane->giveCharacteristicVector(aux, mtrx, mode, tStep);
    if ( !aux.isEmpty() ) answer.assemble(aux, loc_membrane);
}

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
}

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
}

void
HydratingHeMoMaterial :: updateInternalState(const FloatArray &vec, GaussPoint *gp, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    FloatArray aux;
    if ( ms ) {
        ms->letTempStateVectorBe(vec);
        if ( hydration ) {
            /* OBSOLETE
             * FloatArray s = ms->giveStateVector ();
             * if (vec.isEmpty()) OOFEM_ERROR("empty new state vector");
             * aux.resize(2);
             * aux.at(1) = vec.at(1);
             * if (s.isEmpty()||(tStep->giveTime()<=0)) aux.at(2) = initialHydrationDegree; // apply initial conditions
             * else {
             *  aux.at(2) = s.at(2);
             *  if (!castAt || (tStep->giveTime()>=castAt)) aux.at(2) += hydrationModel->dksi (s.at(2), vec.at(1), tStep->giveTimeIncrement()); // compute hydration degree increment
             * }
             */
            // it is necessary to convert the passed state vector to relative humidity expected by the hydration model
            //!!! might be cleaner to choose wc / h in hydration model, but it must be defined which one is passed anyway; so relative humidity was chosen
            //!!! also, the humidity vector might be evaluated by a function (ensure 2 elements and set humidity)
            FloatArray vech = vec;
            if ( vech.giveSize() >= 2 ) {
                vech.at(2) = inverse_sorption_isotherm( vec.at(2) );           // compute relative humidity
            } else {
                vech.resize(2);
                vech.at(2) = 1.; // saturated if undefined
            }

            HydrationModelInterface :: updateInternalState(vech, gp, tStep);

            // additional file output !!!
            if ( teplotaOut && ( gp->giveNumber() == 1 ) && giveStatus(gp) ) {
                FILE *vyst = fopen("teplota.out", "a");
                computeInternalSourceVector(aux, gp, tStep, VM_Incremental);
                if ( aux.isEmpty() ) {
                    aux.resize(1);
                    aux.zero();
                }

                aux.times( 1. / give('d', gp) );
                fprintf( vyst, "Elem %.3d krok %.2d: t= %.0f, dt=%.0f, %ld. it, ksi= %.12f, T= %.8f, heat=%.8f\n", gp->giveElement()->giveNumber(), tStep->giveNumber(),
                        tStep->giveTargetTime(), tStep->giveTimeIncrement(), tStep->giveSolutionStateCounter(),
                        giveHydrationDegree(gp, tStep, VM_Total), vec.at(1), aux.at(1) * tStep->giveTimeIncrement() );
                fclose(vyst);
            }
        }
    }
}

void
MazarsMaterial :: computeEquivalentStrain(double &kappa, const FloatArray &strain, GaussPoint *gp, TimeStep *tStep)
{
    double posNorm = 0.0;
    FloatArray principalStrains, strainb;

    if ( strain.isEmpty() ) {
        kappa = 0.;
        return;
    }

    StructuralMaterial :: giveFullSymVectorForm( strainb, strain, gp->giveMaterialMode() );
    // if plane stress mode -> compute strain in z-direction from condition of zero stress in corresponding direction
    int ndim = giveNumberOfSpatialDimensions(gp);
    if ( ndim == 2 ) {
        strainb.at(3) = -nu * ( strainb.at(1) + strainb.at(2) ) / ( 1. - nu );
    } else if ( ndim == 1 ) {
        strainb.at(2) = -nu *strainb.at(1);
        strainb.at(3) = -nu *strainb.at(1);
    }

    if ( ndim == 1 ) {
        principalStrains.resize(3);
        for ( int i = 1; i <= 3; i++ ) {
            principalStrains.at(i) = strainb.at(i);
        }
    } else {
        this->computePrincipalValues(principalStrains, strainb, principal_strain);
    }

    /*
     * // simple check
     * double i1 = strainb.at(1)+strainb.at(2)+strainb.at(3) - principalStrains.at(1)-principalStrains.at(2)-principalStrains.at(3);
     * double i2 = strainb.at(1)*strainb.at(3)+strainb.at(2)*strainb.at(3)+strainb.at(1)*strainb.at(2) -
     * 0.25*(strainb.at(4)*strainb.at(4)+strainb.at(5)*strainb.at(5)+strainb.at(6)*strainb.at(6)) -
     * principalStrains.at(1)*principalStrains.at(3)+principalStrains.at(2)*principalStrains.at(3)+principalStrains.at(1)*principalStrains.at(2);
     * if ((fabs(i1) > 1.e-6) || (fabs(i2) > 1.e-6)) {
     * printf("v");
     * }
     * // end simple check
     */
    for ( int i = 1; i <= 3; i++ ) {
        if ( principalStrains.at(i) > 0.0 ) {
            posNorm += principalStrains.at(i) * principalStrains.at(i);
        }
    }

    kappa = sqrt(posNorm);
}

void
HeMoTKMaterial :: matcond3d(FloatMatrix &d, GaussPoint *gp, MatResponseMode mode, TimeStep *atTime)
//  function creates conductivity matrix of the
//  isotropic heat material for 3D problems
//
//  d - conductivity matrix of the material
//  25.9.2001
{
    double k = 0.0, w = 0.0, t = 0.0;
    TransportMaterialStatus *status = ( TransportMaterialStatus * ) this->giveStatus(gp);
    FloatArray s;


    //  w = Tm->ip[ipp].av[0];
    //  t = Tm->ip[ipp].av[1];
    s = status->giveTempStateVector();
    if ( s.isEmpty() ) {
        _error("matcond3d: undefined state vector");
    }

    w = s.at(2);
    t = s.at(1);

    if ( mode == Conductivity_ww ) {
        k = perm_ww(w, t);
    } else if ( mode == Conductivity_wh ) {
        k = perm_wt(w, t);
    } else if ( mode == Conductivity_hw ) {
        k = perm_ww(w, t) * get_latent(w, t);
    } else if ( mode == Conductivity_hh ) {
        k = get_chi(w, t) + get_latent(w, t) * perm_wt(w, t);
    } else {
        _error("Unknown MatResponseMode");
    }

    d.resize(3, 3);
    d.at(1, 1) = k;
    d.at(1, 2) = 0.0;
    d.at(1, 3) = 0.0;
    d.at(2, 1) = 0.0;
    d.at(2, 2) = k;
    d.at(2, 3) = 0.0;
    d.at(3, 1) = 0.0;
    d.at(3, 2) = 0.0;
    d.at(3, 3) = k;
}

double
J2Mat :: computeJ2InvariantAt(const FloatArray &stressVector)
{
    double answer;
    double v1, v2, v3;

    if ( stressVector.isEmpty() ) {
        return 0.0;
    }

    v1 = ( ( stressVector.at(1) - stressVector.at(2) ) * ( stressVector.at(1) - stressVector.at(2) ) );
    v2 = ( ( stressVector.at(2) - stressVector.at(3) ) * ( stressVector.at(2) - stressVector.at(3) ) );
    v3 = ( ( stressVector.at(3) - stressVector.at(1) ) * ( stressVector.at(3) - stressVector.at(1) ) );

    answer = ( 1. / 6. ) * ( v1 + v2 + v3 ) + stressVector.at(4) * stressVector.at(4) +
             stressVector.at(5) * stressVector.at(5) + stressVector.at(6) * stressVector.at(6);

    return answer;
}

double
HydrationModel :: giveCharacteristicValue(const FloatArray &vec, MatResponseMode rmode, GaussPoint *gp, TimeStep *tStep)
// Transport status needs to be obtained from master status, it's better to pass as a parameter
// to enable usage from structural material
{
    double answer = 0;

    if ( ( rmode == IntSource ) || ( rmode == IntSource_hh ) || ( rmode == IntSource_ww ) || ( rmode == IntSource_wh ) || ( rmode == IntSource_hw ) ) {
        if ( vec.isEmpty() ) {
            OOFEM_ERROR("undefined state vector.");
        }

        answer = computeIntSource(vec, gp, tStep, rmode);
    } else {
        OOFEM_ERROR("wrong MatResponseMode.");
    }

    return answer;
}

void
HydratingIsoHeatMaterial :: updateInternalState(const FloatArray &vec, GaussPoint *gp, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    FloatArray aux;
    if ( ms ) {
        ms->letTempStateVectorBe(vec);
        if ( hydration ) {
            /* OBSOLETE
             * FloatArray s = ms->giveStateVector ();
             * if (vec.isEmpty()) OOFEM_ERROR("empty new state vector");
             * aux.resize(2);
             * aux.at(1) = vec.at(1);
             * if (s.isEmpty()||(tStep->giveTime()<=0)) aux.at(2) = initialHydrationDegree; // apply initial conditions
             * else {
             *  aux.at(2) = s.at(2);
             *  if (!castAt || (tStep->giveTime()>=castAt)) aux.at(2) += hydrationModel->dksi (s.at(2), vec.at(1), tStep->giveTimeIncrement()); // compute hydration degree increment
             * }
             */
            HydrationModelInterface :: updateInternalState(vec, gp, tStep);

            // additional file output !!!
            if ( ( gp->giveNumber() == 1 ) && giveStatus(gp) ) {
                FILE *vyst = fopen("teplota.out", "a");
                computeInternalSourceVector(aux, gp, tStep, VM_Incremental);
                if ( aux.isEmpty() ) {
                    aux.resize(1);
                    aux.zero();
                }

                aux.times( 1. / give('d', gp, tStep) );
                fprintf( vyst, "Elem %.3d krok %.2d: t= %.0f, dt=%.0f, %ld. it, ksi= %.12f, T= %.8f, heat=%.8f\n", gp->giveElement()->giveNumber(), tStep->giveNumber(),
                        tStep->giveTargetTime(), tStep->giveTimeIncrement(), tStep->giveSolutionStateCounter(),
                        giveHydrationDegree(gp, tStep, VM_Total), vec.at(1), aux.at(1) * tStep->giveTimeIncrement() );
                fclose(vyst);
            }
        }
    }
}

void
HeMoTKMaterial :: matcond1d(FloatMatrix &d, GaussPoint *gp, MatResponseMode mode, TimeStep *tStep)
//  function creates conductivity matrix of the
//  isotropic heat material for 1D problems
//
//  d - conductivity matrix of the material
//  25.9.2001
{
    double k = 0.0, w = 0.0, t = 0.0;
    TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    FloatArray s;


    //  w = Tm->ip[ipp].av[0];
    //  t = Tm->ip[ipp].av[1];
    s = status->giveTempField();
    if ( s.isEmpty() ) {
        OOFEM_ERROR("undefined state vector");
    }

    w = s.at(2);
    t = s.at(1);

    if ( mode == Conductivity_ww ) {
        k = perm_ww(w, t);
    } else if ( mode == Conductivity_wh ) {
        k = perm_wt(w, t);
    } else if ( mode == Conductivity_hw ) {
        k = perm_ww(w, t) * get_latent(w, t);
    } else if ( mode == Conductivity_hh ) {
        k = get_chi(w, t) + get_latent(w, t) * perm_wt(w, t);
    } else {
        OOFEM_ERROR("Unknown MatResponseMode");
    }

    d.resize(1, 1);
    d.at(1, 1) = k;
}

void
HeMoKunzelMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &grad, const FloatArray &field, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );

    FloatArray s;
    //     s = ms->giveTempStateVector();
    s = ms->giveTempField();
    if ( s.isEmpty() ) {
        OOFEM_ERROR("matcond1d: undefined state vector");
    }
    double h = s.at(2);
    double t = s.at(1);

    FloatArray ans_w, ans_t;
    FloatArray grad_w, grad_t;
    int size = grad.giveSize() / 2;
    for ( int i = 1; i <= size; ++i ) {
        grad_w.at(i) = grad.at(i);
    }
    for ( int i = size + 1; i <= size * 2; ++i ) {
        grad_t.at(i) = grad.at(i);
    }

    ans_w.beScaled(perm_mm(h, t), grad_w);
    ans_w.beScaled(perm_mh(h, t), grad_t);
    ans_t.beScaled(perm_hm(h, t), grad_w);
    ans_t.beScaled(perm_hh(h, t), grad_t);

    answer.resize(size * 2);
    answer.zero();
    answer.addSubVector(ans_w, 1);
    answer.addSubVector(ans_t, size + 1);

    ms->setTempField(field);
    ms->setTempGradient(grad);
    ms->setTempFlux(answer);
}

void
MazarsMaterial :: computeDamageParam(double &omega, double kappa, const FloatArray &strain, GaussPoint *gp)
{
    // positive_flag = 0, negat_count = 0;
    FloatMatrix de, ce;
    FloatArray sigp, epsti, epsi;
    double gt, gc, alpha_t, alpha_c, alpha, eqStrain2;

    if ( kappa <= this->e0 ) { // strain below damage threshold
        omega = 0.0;
        return;
    }

    // strain above damage threshold

    int ndim = giveNumberOfSpatialDimensions(gp);
    if ( !strain.isEmpty() ) {
        this->computePrincipalValues(epsi, strain, principal_strain);
    } else {
        epsi.resize(ndim);
        epsi.zero();
    }

    // construct the normal block of elastic compliance matrix
    giveNormalBlockOfElasticCompliance(ce, gp);
    // compute the normal block of elastic stiffness matrix
    de.beInverseOf(ce);

    // compute principal stresses
    sigp.beProductOf(de, epsi);
    // take positive part
    for ( int i = 1; i <= ndim; i++ ) {
        if ( sigp.at(i) < 0. ) {
            sigp.at(i) = 0.;
        }
    }

    // compute principal strains due to positive stresses
    epsti.beProductOf(ce, sigp);

    // extend strains to 3 dimensions
    if ( ndim == 2 ) {
        epsi.resize(3);
        epsti.resize(3);
        epsi.at(3) = -nu * ( epsi.at(1) + epsi.at(2) ) / ( 1. - nu );
        epsti.at(3) = -nu * ( epsti.at(1) + epsti.at(2) ) / ( 1. - nu );
    } else if ( ndim == 1 ) {
        epsi.resize(3);
        epsti.resize(3);
        epsi.at(2) = epsi.at(3) = -nu *epsi.at(1);
        epsti.at(2) = epsti.at(3) = -nu *epsti.at(1);
    }

    /* the following section "improves" biaxial compression
     * // but it may lead to convergence problems, probably due to the condition > 1.e-6
     * // therefore it was commented out by Milan Jirasek on 22 Nov 2009
     *
     * positive_flag = negat_count = 0;
     *    for ( i = 1; i <= 3; i++ ) {
     *        if ( sigp.at(i) > 1.e-6 ) {
     *            positive_flag = 1;
     *            break;
     *        } else if ( sigp.at(i) < 1.e-6 ) {
     *            negat_count++;
     *        }
     *    }
     *
     *    if ( ( positive_flag == 0 ) && ( negat_count > 1 ) ) {
     *        // adjust equivalent strain to improve biaxial compression
     *        double f = 0.0, g = 0.0;
     *        for ( i = 1; i <= 3; i++ ) {
     *            if ( sigp.at(i) < 0 ) {
     *                f += sigp.at(i) * sigp.at(i);
     *                g += fabs( sigp.at(i) );
     *            }
     *        }
     *
     *        f = sqrt(f) / g;
     *        kappa *= f;
     *    }
     *
     *    // test adjusted kappa
     *    if ( kappa < this->e0 ) {
     *        omega = 0.0;
     *        return;
     *    }
     *    // end of section that improves biaxial compression
     */

    // evaluation of damage functions
    gt = computeGt(kappa, gp);
    gc = computeGc(kappa, gp);

    // evaluation of factors alpha_t and alpha_c
    alpha = 0.0;
    eqStrain2 = 0.0;
    for ( int i = 1; i <= 3; i++ ) {
        if ( epsi.at(i) > 0.0 ) {
            eqStrain2 += epsi.at(i) * epsi.at(i);
            alpha += epsti.at(i) * epsi.at(i);
//.........这里部分代码省略.........