C++ vector_fp类代码示例

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

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

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

  /*
   * Vector xpts contains a monotonic sequence of grid points, and 
   * vector fpts contains function values defined at these points.
   * The value returned is the linear interpolate at point x.
   * If x is outside the range of xpts, the value of fpts at the 
   * nearest end is returned.
   *
   * @param x value of the x coordinate
   * @param xpts value of the grid points
   * @param fpts value of the interpolant at the grid points
   *
   * @return Returned value is the value of of the interpolated
   *         function at x.
   */
  doublereal linearInterp(doublereal x, const vector_fp& xpts, 
			  const vector_fp& fpts) {
    if (x <= xpts[0]) 
      return fpts[0];
    if (x >= xpts.back()) 
      return fpts.back();
    vector_fp::const_iterator loc = 
      lower_bound(xpts.begin(), xpts.end(), x);
    int iloc = int(loc - xpts.begin()) - 1;
    doublereal ff = fpts[iloc] + 
      (x - xpts[iloc])*(fpts[iloc + 1] 
			- fpts[iloc])/(xpts[iloc + 1] - xpts[iloc]);
    return ff;
  }

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

/*!
 * @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);
      }
    }
  }

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 Phase::restoreState(const vector_fp& state)
{
    restoreState(state.size(),&state[0]);
}

doublereal Phase::mean_X(const vector_fp& Q) const
{
    return m_mmw*std::inner_product(m_ym.begin(), m_ym.end(), Q.begin(), 0.0);
}

  /**
   *
   * Read species data from THERMO section records. 
   * 
   * @param names        List of species names (input).
   * @param species      Table of species objects holding data from records
   *                     in THERMO section (output). 
   * @param temp         Devault vector of temperature region boundaries
   *                     There are one more temperatures than there are
   *                     temperature regions.
   * @param allowExtThermoData   True if 'THERMO' specified, false if 
   *                             'THERMO ALL' specified.
   *
   * @return            True, if the THERMO section exists and the species
   *                    have all been successfully processed. False, if
   *                    the THERMO section doesn't exist or there were
   *                    additional problems.
   */
  bool CKParser::readNASA9ThermoSection(std::vector<string>& names, 
					speciesTable& species, vector_fp& temp, 
					int& optionFlag, std::ostream& log) {
    // String buffer for lines
    string s;
    vector<string> toks;
    string defaultDate="";
    int nreg = 2;
    int i;

    int nsp = static_cast<int>(names.size());

    // Comment string
    string comment;


    // if "THERMO ALL" specified, or if optionFlag is set to HasTempRange,
    // then the next line must be the default temperatures for the database.
    //
    //  This line will have nreg+2 tokens on it
    //    The last token is a date.
    if (0) {
      if (optionFlag == NoThermoDatabase || optionFlag == HasTempRange) {
	getCKLine(s, comment);
	getTokens(s, static_cast<int>(s.size()), toks);
	nreg = toks.size();
	if (nreg >= 1) {
	  temp.resize(nreg+1);
	  for (i = 0; i <= nreg; i++) {
	    temp[i] = de_atof(toks[i]);
	  }
	  defaultDate = toks[nreg+1];
	}
      
	if (verbose) {
	  log.flags(ios::showpoint | ios::fixed);
	  log.precision(2);
	  log << endl << " Default # of temperature regions: " << nreg << endl;
	  log << "          ";
	  for (i = 0; i <= nreg; i++) {
	    log << temp[i] << "  ";
	  }
	  log << endl;
	}
	checkNASA9Temps(log, temp);
      }
    }
    
    // Check to see that we expect to be reading a NASA9 formatted file
    if (!m_nasa9fmt) {
      throw CK_SyntaxError(log, 
			   "In NASA9 parser. However, we expect a different file format",
			   -1);
    }

    // now read in all species records that have names in list 'names'

    bool getAllSpecies = (nsp > 0 && match(names[0], "<ALL>"));
    if (getAllSpecies) names.clear();

    // Map between the number of times a species name appears in the database
    map<string, int> dup; // used to check for duplicate THERMO records
    bool already_read;

    while (1 > 0) {
      // If we don't have any more species to read, break
      if (nsp == 0) break;
      already_read = false;

      // Read a new species record from the section
      Species spec;
      readNASA9ThermoRecord(spec);

      // we signal the end of the section by putting <END> as a
      // species name. Break if you find the end of the section.
      if (spec.name == "<END>") {
	break;
      }
        
      // check for duplicate thermo data
      if (dup[spec.name] == 2) {
	log << "Warning: more than one THERMO record for "
//.........这里部分代码省略.........

void Phase::saveState(vector_fp& state) const
{
    state.resize(nSpecies() + 2);
    saveState(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");

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