C++ TransportMaterialStatus类代码示例

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
HeMoTKMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &grad, const FloatArray &field, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );

    double w = field.at(2);
    double t = field.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_ww(w, t), grad_w);
    ans_w.beScaled(perm_wt(w, t), grad_t);
    ans_t.beScaled(perm_ww(w, t) * get_latent(w, t), grad_w);
    ans_t.beScaled(get_chi(w, t) + get_latent(w, t) * perm_wt(w, 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
TransportMaterial :: updateInternalState(const FloatArray &stateVec, GaussPoint *gp, TimeStep *)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    if ( ms ) {
        ms->letTempStateVectorBe(stateVec);
    }
}

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
IsotropicMoistureTransferMaterial :: giveFluxVector(FloatArray &answer, GaussPoint *gp, const FloatArray &grad, const FloatArray &field, TimeStep *tStep)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );

    ///@todo Shouldn't the permeability typically depend on the primary field and/or its gradient?
    answer.beScaled(-this->givePermeability(gp, tStep), grad);

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

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

    answer.beProductOf(k, grad);
    answer.negated();

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

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

    double gradPNorm = grad.computeNorm();
    answer.beScaled( -(1. + C * pow(gradPNorm, alpha)), grad);

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

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);
            }
        }
    }
}

int
NonlinearMassTransferMaterial :: giveIPValue(FloatArray &answer, GaussPoint *aGaussPoint, InternalStateType type, TimeStep *atTime)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(aGaussPoint) );
 
    switch ( type ) {
    case IST_Velocity:
        answer = ms->giveFlux();
        break;
    case IST_PressureGradient:
        answer = ms->giveGradient();
        break;
    default:
      return TransportMaterial :: giveIPValue(answer, aGaussPoint, type, atTime);
    }
    return 1;
}

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;
}

void
NonlinearMassTransferMaterial :: giveCharacteristicMatrix(FloatMatrix &answer,
                                                          MatResponseMode mode,
                                                          GaussPoint *gp,
                                                          TimeStep *atTime)
{
    MaterialMode mMode = gp->giveMaterialMode();
    TransportMaterialStatus *status = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    FloatArray eps = status->giveTempGradient();
    double gradPNorm;
    FloatMatrix t1, t2;

    gradPNorm = eps.computeNorm();

    t1.beDyadicProductOf(eps, eps);
    if ( gradPNorm != 0.0 ) {
        t1.times( C * alpha * pow(gradPNorm, alpha - 2) );
    }

    switch  ( mMode ) {
    case _1dHeat:
        t2.resize(1, 1);
        t2.at(1, 1) = 1;
        break;
    case _2dHeat:
        t2.resize(2, 2);
        t2.at(1, 1) = t2.at(2, 2) = 1;
        break;
    case _3dHeat:
        t2.resize(3, 3);
        t2.at(1, 1) = t2.at(2, 2) = t2.at(3, 3) = 1;
        break;
    default:
        _error2( "giveCharacteristicMatrix : unknown mode (%s)", __MaterialModeToString(mMode) );
    }

    answer.beEmptyMtrx();
    answer.add(t1);
    answer.add(1 + C * pow(gradPNorm, alpha), t2);

}

double
HeMoTKMaterial :: giveHumidity(GaussPoint *gp, ValueModeType mode)
{
    TransportMaterialStatus *ms = static_cast< TransportMaterialStatus * >( this->giveStatus(gp) );
    const FloatArray &tempState = ms->giveTempField();
    if ( tempState.giveSize() < 2 ) {
        OOFEM_ERROR("undefined moisture status!");
    }

    FloatArray state = ms->giveField();

    if ( mode == VM_Total ) {
        return inverse_sorption_isotherm( tempState.at(2) );
    } else if ( mode == VM_Incremental ) {
        return inverse_sorption_isotherm( tempState.at(2) ) - inverse_sorption_isotherm( state.at(2) );
    } else if ( mode == VM_Velocity ) { // VM_Previous
        return inverse_sorption_isotherm( state.at(2) );
    }

    return 1.;
}

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);
            }
        }
    }
}

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
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);
}