C++ vector_fp::size方法代码示例

std::string vec2str(const vector_fp& v, const std::string& fmt,
                    const std::string& sep)
{
    char buf[64];
    std::stringstream o;
    for (size_t i = 0; i < v.size(); i++) {
        SNPRINTF(buf, 63, fmt.c_str(), v[i]);
        o << v[i];
        if (i != v.size() - 1) {
            o << sep;
        }
    }
    return o.str();
}

void Phase::init(const vector_fp& mw)
{
    m_kk = mw.size();
    m_rmolwts.resize(m_kk);
    m_y.resize(m_kk, 0.0);
    m_ym.resize(m_kk, 0.0);
    copy(mw.begin(), mw.end(), m_molwts.begin());
    for (size_t k = 0; k < m_kk; k++) {
        if (m_molwts[k] < 0.0) {
            throw CanteraError("Phase::init",
                               "negative molecular weight for species number "
                               + int2str(k));
        }

        // Some surface phases may define species representing empty sites
        // that have zero molecular weight. Give them a very small molecular
        // weight to avoid dividing by zero.
        if (m_molwts[k] < Tiny) {
            m_molwts[k] = Tiny;
        }
        m_rmolwts[k] = 1.0/m_molwts[k];
    }

    // Now that we have resized the State object, let's fill it with a valid
    // mass fraction vector that sums to one. The Phase object should never
    // have a mass fraction vector that doesn't sum to one. We will assume that
    // species 0 has a mass fraction of 1.0 and mass fraction of all other
    // species is 0.0.
    m_y[0] = 1.0;
    m_ym[0] = m_y[0] * m_rmolwts[0];
    m_mmw = 1.0 / m_ym[0];
}

void Phase::getMolecularWeights(vector_fp& weights) const
{
    const vector_fp& mw = molecularWeights();
    if (weights.size() < mw.size()) {
        weights.resize(mw.size());
    }
    copy(mw.begin(), mw.end(), weights.begin());
}

/*!
 * @param rxn  Reaction index of the current reaction. This is used
 *             as an index into vectors which have length n_total_rxn.
 * @param k    This is a vector of integer values specifying the
 *             species indices. The length of this vector species
 *             the number of different species in the description.
 *             The value of the entries are the species indices.
 *             These are used as indexes into vectors which have
 *             length n_total_species.
 *  @param order This is a vector of the same length as vector k.
 *         The order is used for the routine power(), which produces
 *         a power law expression involving the species vector.
 *  @param stoich  This is used to handle fractional stoichiometric coefficients
 *                 on the product side of irreversible reactions.
 */
void StoichManagerN::add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
                         const vector_fp& stoich) {
    //printf ("add called\n");
    if (order.size() != k.size()) {
        throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ");
    }
    if (stoich.size() != k.size()) {
        throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ");
    }
    bool frac = false;
    for (size_t n = 0; n < stoich.size(); n++) {
        if (fmod(stoich[n], 1.0) || fmod(order[n], 1.0)) {
            frac = true;
            break;
        }
    }
    if (frac || k.size() > 3) {
        m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
    } else {
        // Try to express the reaction with unity stoichiometric
        // coefficients (by repeating species when necessary) so that the
        // simpler 'multiply' function can be used to compute the rate
        // instead of 'power'.
        std::vector<size_t> kRep;
        for (size_t n = 0; n < k.size(); n++) {
            for (size_t i = 0; i < stoich[n]; i++)
                kRep.push_back(k[n]);
        }

        switch (kRep.size()) {
        case 1:
            m_c1_list.push_back(C1(rxn, kRep[0]));
            break;
        case 2:
            m_c2_list.push_back(C2(rxn, kRep[0], kRep[1]));
            break;
        case 3:
            m_c3_list.push_back(C3(rxn, kRep[0], kRep[1], kRep[2]));
            break;
        default:
            m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
        }
    }
}

  /**
   *  Check validity of the temperatures defining the 
   *  temperature ranges for the NASA9 polynomial species thermodynamic
   *  property fits. 
   *  @param log     log file output stream
   *  @param temp    Vector of temperatures
   */
  static void checkNASA9Temps(std::ostream& log, vector_fp &temp) {
    int i;
    for (i = 1; i <= (int) temp.size(); i++) {
      double tlow = temp[i-1];
      double thigh = temp[i];
      if (thigh <= tlow) {
	string sss =  "error reading temperature";
	throw CK_SyntaxError(log, sss);
      }
    }
  }

void ThermoPhase::setElementPotentials(const vector_fp& lambda)
{
    size_t mm = nElements();
    if (lambda.size() < mm) {
        throw CanteraError("setElementPotentials", "lambda too small");
    }
    if (!m_hasElementPotentials) {
        m_lambdaRRT.resize(mm);
    }
    scale(lambda.begin(), lambda.end(), m_lambdaRRT.begin(), 1.0/RT());
    m_hasElementPotentials = true;
}

static void printProgress(const vector<string> &spName,
                          const vector_fp &soln,
                          const vector_fp &ff)
{
    double sum = 0.0;
    plogf(" --- Summary of current progress:\n");
    plogf(" ---                   Name           Moles  -       SSGibbs \n");
    plogf(" -------------------------------------------------------------------------------------\n");
    for (size_t k = 0; k < soln.size(); k++) {
        plogf(" ---      %20s %12.4g  - %12.4g\n", spName[k], soln[k], ff[k]);
        sum += soln[k] * ff[k];
    }
    plogf(" ---  Total sum to be minimized = %g\n", sum);
}

    /// print the falloff parameters for a pressure-dependent reaction
    bool writeFalloff(int type, const vector_fp& c, std::ostream& log) {

        log.precision(6);
        log.width(0);
        log.flags(ios::uppercase);

        //    bool ok = true;
        switch (type) {

        case Lindemann:
            log << "   Lindemann falloff function" << endl;
            return true;

        case Troe:
            log << "   Troe falloff function: " << endl;
            if (c.size() == 3) {
                log << "      alpha, T***, T* = (" << c[0] << ", " << c[1] 
                    << ", " << c[2] << ")" << endl;
            }
            else if (c.size() == 4) {
                log << "      alpha, T***, T*, T** = (" << c[0] << ", " << c[1] 
                    << ", " << c[2] << ", " << c[3] << ")" << endl;
            }
            else {
				for (size_t n = 0; n < c.size(); n++) {
					log << c[n] << ", "; log << endl;
				}
                log << "###### ERROR #####   incorrect number of parameters" << endl;
                return false;
            }
            return true;

        case SRI:
            log << "   SRI falloff function: " << endl;
            if (c.size() == 3) {
                log << "      a, b, c = (" << c[0] << ", " << c[1] 
                    << ", " << c[2] << ")" << endl;
            }
            else if (c.size() == 5) {
                log << "      a, b, c, d, e = (" << c[0] << ", " << c[1] 
                    << ", " << c[2] << ", " << c[3] << ", " << c[4] 
                    << ")" << endl;
            }
            else {
				for (size_t n = 0; n < c.size(); n++) {
					log << c[n] << ", "; log << endl;
				}
                log << "##### ERROR #####  incorrect number of parameters" << endl;
                return false;
            }
            return true;

        default:
            log << "unknown falloff type: " << type << endl;
            return false;
        }
    }

void IonFlow::setElectronTransport(vector_fp& tfix, vector_fp& diff_e,
                                   vector_fp& mobi_e)
{
    m_import_electron_transport = true;
    size_t degree = 5;
    size_t n = tfix.size();
    vector_fp tlog;
    for (size_t i = 0; i < n; i++) {
        tlog.push_back(log(tfix[i]));
    }
    vector_fp w(n, -1.0);
    m_diff_e_fix.resize(degree + 1);
    m_mobi_e_fix.resize(degree + 1);
    polyfit(n, degree, tlog.data(), diff_e.data(), w.data(), m_diff_e_fix.data());
    polyfit(n, degree, tlog.data(), mobi_e.data(), w.data(), m_mobi_e_fix.data());
}

void Phase::saveState(vector_fp& state) const
{
    state.resize(nSpecies() + 2);
    saveState(state.size(), &state[0]);
}

void Phase::restoreState(const vector_fp& state)
{
    restoreState(state.size(),&state[0]);
}

/*
 * Choose the optimum basis for the calculations. This is done by
 * choosing the species with the largest mole fraction
 * not currently a linear combination of the previous components.
 * Then, calculate the stoichiometric coefficient matrix for that
 * basis.
 *
 * Calculates the identity of the component species in the mechanism.
 * Rearranges the solution data to put the component data at the
 * front of the species list.
 *
 * Then, calculates SC(J,I) the formation reactions for all noncomponent
 * species in the mechanism.
 *
 * Input
 * ---------
 * mphase          Pointer to the multiphase object. Contains the
 *                 species mole fractions, which are used to pick the
 *                 current optimal species component basis.
 * orderVectorElement
 *                 Order vector for the elements. The element rows
 *                 in the formula matrix are
 *                 rearranged according to this vector.
 * orderVectorSpecies
 *                 Order vector for the species. The species are
 *                 rearranged according to this formula. The first
 *                 nCompoments of this vector contain the calculated
 *                 species components on exit.
 * doFormRxn       If true, the routine calculates the formation
 *                 reaction matrix based on the calculated
 *                 component species. If false, this step is skipped.
 *
 * Output
 * ---------
 * usedZeroedSpecies = If true, then a species with a zero concentration
 *                     was used as a component. The problem may be
 *                     converged.
 * formRxnMatrix
 *
 * Return
 * --------------
 * returns the number of components.
 *
 *
 */
size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
                              MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies,
                              std::vector<size_t>& orderVectorElements,
                              vector_fp& formRxnMatrix)
{

    size_t  j, jj, k=0, kk, l, i, jl, ml;
    bool lindep;
    std::string ename;
    std::string sname;
    /*
     * Get the total number of elements defined in the multiphase object
     */
    size_t ne = mphase->nElements();
    /*
     * Get the total number of species in the multiphase object
     */
    size_t nspecies = mphase->nSpecies();
    doublereal tmp;
    doublereal const USEDBEFORE = -1;

    /*
     * Perhaps, initialize the element ordering
     */
    if (orderVectorElements.size() < ne) {
        orderVectorElements.resize(ne);
        for (j = 0; j < ne; j++) {
            orderVectorElements[j] = j;
        }
    }

    /*
     * Perhaps, initialize the species ordering
     */
    if (orderVectorSpecies.size() != nspecies) {
        orderVectorSpecies.resize(nspecies);
        for (k = 0; k < nspecies; k++) {
            orderVectorSpecies[k] = k;
        }
    }

#ifdef DEBUG_MODE
    double molSave = 0.0;
    if (BasisOptimize_print_lvl >= 1) {
        writelog("   ");
        for (i=0; i<77; i++) {
            writelog("-");
        }
        writelog("\n");
        writelog("   --- Subroutine BASOPT called to ");
        writelog("calculate the number of components and ");
        writelog("evaluate the formation matrix\n");
        if (BasisOptimize_print_lvl > 0) {
            writelog("   ---\n");

//.........这里部分代码省略.........

void Phase::restoreState(const vector_fp& state)
{
    restoreState(state.size(),&state[0]);
    compositionChanged();
}

void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
                   MultiPhase* mphase,
                   std::vector<size_t>& orderVectorSpecies,
                   std::vector<size_t>& orderVectorElements)
{
    size_t nelements = mphase->nElements();
    // Get the total number of species in the multiphase object
    size_t nspecies = mphase->nSpecies();

    if (BasisOptimize_print_lvl > 0) {
        writelog("   ");
        writeline('-', 77);
        writelog("   --- Subroutine ElemRearrange() called to ");
        writelog("check stoich. coefficient matrix\n");
        writelog("   ---    and to rearrange the element ordering once\n");
    }

    // Perhaps, initialize the element ordering
    if (orderVectorElements.size() < nelements) {
        orderVectorElements.resize(nelements);
        for (size_t j = 0; j < nelements; j++) {
            orderVectorElements[j] = j;
        }
    }

    // Perhaps, initialize the species ordering. However, this is dangerous, as
    // this ordering is assumed to yield the component species for the problem
    if (orderVectorSpecies.size() != nspecies) {
        orderVectorSpecies.resize(nspecies);
        for (size_t k = 0; k < nspecies; k++) {
            orderVectorSpecies[k] = k;
        }
    }

    // If the elementAbundances aren't input, just create a fake one based on
    // summing the column of the stoich matrix. This will force elements with
    // zero species to the end of the element ordering.
    vector_fp eAbund(nelements,0.0);
    if (elementAbundances.size() != nelements) {
        for (size_t j = 0; j < nelements; j++) {
            eAbund[j] = 0.0;
            for (size_t k = 0; k < nspecies; k++) {
                eAbund[j] += fabs(mphase->nAtoms(k, j));
            }
        }
    } else {
        copy(elementAbundances.begin(), elementAbundances.end(),
             eAbund.begin());
    }

    vector_fp sa(nelements,0.0);
    vector_fp ss(nelements,0.0);
    vector_fp sm(nelements*nelements,0.0);

    // Top of a loop of some sort based on the index JR. JR is the current
    // number independent elements found.
    size_t jr = 0;
    while (jr < nComponents) {
        // Top of another loop point based on finding a linearly independent
        // element
        size_t k = nelements;
        while (true) {
            // Search the element vector. We first locate elements that are
            // present in any amount. Then, we locate elements that are not
            // present in any amount. Return its identity in K.
            size_t kk;
            for (size_t ielem = jr; ielem < nelements; ielem++) {
                kk = orderVectorElements[ielem];
                if (eAbund[kk] != USEDBEFORE && eAbund[kk] > 0.0) {
                    k = ielem;
                    break;
                }
            }
            for (size_t ielem = jr; ielem < nelements; ielem++) {
                kk = orderVectorElements[ielem];
                if (eAbund[kk] != USEDBEFORE) {
                    k = ielem;
                    break;
                }
            }

            if (k == nelements) {
                // When we are here, there is an error usually.
                // We haven't found the number of elements necessary.
                if (BasisOptimize_print_lvl > 0) {
                    writelogf("Error exit: returning with nComponents = %d\n", jr);
                }
                throw CanteraError("ElemRearrange", "Required number of elements not found.");
            }

            // Assign a large negative number to the element that we have
            // just found, in order to take it out of further consideration.
            eAbund[kk] = USEDBEFORE;

            // CHECK LINEAR INDEPENDENCE OF CURRENT FORMULA MATRIX
            // LINE WITH PREVIOUS LINES OF THE FORMULA MATRIX

            // Modified Gram-Schmidt Method, p. 202 Dalquist
            // QR factorization of a matrix without row pivoting.
            size_t jl = jr;
//.........这里部分代码省略.........

size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
                     std::vector<size_t>& orderVectorSpecies,
                     std::vector<size_t>& orderVectorElements,
                     vector_fp& formRxnMatrix)
{
    // Get the total number of elements defined in the multiphase object
    size_t ne = mphase->nElements();

    // Get the total number of species in the multiphase object
    size_t nspecies = mphase->nSpecies();

    // Perhaps, initialize the element ordering
    if (orderVectorElements.size() < ne) {
        orderVectorElements.resize(ne);
        iota(orderVectorElements.begin(), orderVectorElements.end(), 0);
    }

    // Perhaps, initialize the species ordering
    if (orderVectorSpecies.size() != nspecies) {
        orderVectorSpecies.resize(nspecies);
        iota(orderVectorSpecies.begin(), orderVectorSpecies.end(), 0);
    }

    if (BasisOptimize_print_lvl >= 1) {
        writelog("   ");
        writeline('-', 77);
        writelog("   --- Subroutine BASOPT called to ");
        writelog("calculate the number of components and ");
        writelog("evaluate the formation matrix\n");
        if (BasisOptimize_print_lvl > 0) {
            writelog("   ---\n");

            writelog("   ---      Formula Matrix used in BASOPT calculation\n");
            writelog("   ---      Species | Order | ");
            for (size_t j = 0; j < ne; j++) {
                size_t jj = orderVectorElements[j];
                writelog(" {:>4.4s}({:1d})", mphase->elementName(jj), j);
            }
            writelog("\n");
            for (size_t k = 0; k < nspecies; k++) {
                size_t kk = orderVectorSpecies[k];
                writelog("   --- {:>11.11s} |   {:4d} |",
                         mphase->speciesName(kk), k);
                for (size_t j = 0; j < ne; j++) {
                    size_t jj = orderVectorElements[j];
                    double num = mphase->nAtoms(kk,jj);
                    writelogf("%6.1g  ", num);
                }
                writelog("\n");
            }
            writelog("   --- \n");
        }
    }

    // Calculate the maximum value of the number of components possible. It's
    // equal to the minimum of the number of elements and the number of total
    // species.
    size_t nComponents = std::min(ne, nspecies);
    size_t nNonComponents = nspecies - nComponents;

    // Set this return variable to false
    *usedZeroedSpecies = false;

    // Create an array of mole numbers
    vector_fp molNum(nspecies,0.0);
    mphase->getMoles(molNum.data());

    // Other workspace
    vector_fp sm(ne*ne, 0.0);
    vector_fp ss(ne, 0.0);
    vector_fp sa(ne, 0.0);
    if (formRxnMatrix.size() < nspecies*ne) {
        formRxnMatrix.resize(nspecies*ne, 0.0);
    }

    // For debugging purposes keep an unmodified copy of the array.
    vector_fp molNumBase = molNum;
    double molSave = 0.0;
    size_t jr = 0;

    // Top of a loop of some sort based on the index JR. JR is the current
    // number of component species found.
    while (jr < nComponents) {
        // Top of another loop point based on finding a linearly independent
        // species
        size_t k = npos;
        while (true) {
            // Search the remaining part of the mole number vector, molNum for
            // the largest remaining species. Return its identity. kk is the raw
            // number. k is the orderVectorSpecies index.
            size_t kk = max_element(molNum.begin(), molNum.end()) - molNum.begin();
            size_t j;
            for (j = 0; j < nspecies; j++) {
                if (orderVectorSpecies[j] == kk) {
                    k = j;
                    break;
                }
            }
            if (j == nspecies) {
                throw CanteraError("BasisOptimize", "orderVectorSpecies contains an error");
//.........这里部分代码省略.........