From 4e3412ecd2e691bf2799eabdef15752d4fb31eb6 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Thu, 16 Oct 2008 23:00:11 +0000 Subject: [PATCH] Changed the AddSpecies interface back to use const element compositions. This entailed creating a temp vector inside the routine, that is potentially changed. --- Cantera/src/thermo/Constituents.cpp | 817 ++++++++++++++-------------- Cantera/src/thermo/Constituents.h | 4 +- 2 files changed, 409 insertions(+), 412 deletions(-) diff --git a/Cantera/src/thermo/Constituents.cpp b/Cantera/src/thermo/Constituents.cpp index 68e280e52..894536de8 100755 --- a/Cantera/src/thermo/Constituents.cpp +++ b/Cantera/src/thermo/Constituents.cpp @@ -7,8 +7,6 @@ /* $Author$ * $Date$ * $Revision$ - * - */ // Copyright 2001 California Institute of Technology @@ -24,444 +22,443 @@ using namespace std; namespace Cantera { + /* + * Constructor sets all base variable types to zero. Also, it + * sets the pointer to the Elements object for this object to the + * default value of BaseElements. If the BaseElements Elements + * object doesn't exist, it creates it. + * + * Input + * -------- + * ptr_Elements: If the Constituents object requires a different + * Elements object than the default one, input + * address here. This argument defaults to null, + * in which case the default Elements Object is + * chosen. + */ + + /* + * DGG: I have reversed the role of ptr_Elements. In this version, + * the default is that a new Elements object is created, so this + * Constituents object is independent of any other object. But if + * ptr_Elements is supplied, it will be used. This way, a class + * implementing a multi-phase mixture is responsible for + * maintaining the global elements list for the mixture, and no + * static global element list is required. + */ + Constituents::Constituents(Elements* ptr_Elements) : + m_kk(0), + m_speciesFrozen(false) , + m_Elements(ptr_Elements) + { + if (!m_Elements) m_Elements = new Elements(); + + // Register subscription to Elements object whether or not we + // created it here. + m_Elements->subscribe(); + } + + /** + * Destructor for class Constituents. + * + * Some cleanup of of the Global_Elements_List array is + * effected by unsubscribing to m_Elements. + */ + Constituents::~Constituents() + { + int ileft = m_Elements->unsubscribe(); /* - * Constructor sets all base variable types to zero. Also, it - * sets the pointer to the Elements object for this object to the - * default value of BaseElements. If the BaseElements Elements - * object doesn't exist, it creates it. - * - * Input - * -------- - * ptr_Elements: If the Constituents object requires a different - * Elements object than the default one, input - * address here. This argument defaults to null, - * in which case the default Elements Object is - * chosen. + * Here we may delete Elements Objects or not. Right now, we + * will delete them. We also delete the global pointer entry + * to keep everything consistent. */ - - /* - * DGG: I have reversed the role of ptr_Elements. In this version, - * the default is that a new Elements object is created, so this - * Constituents object is independent of any other object. But if - * ptr_Elements is supplied, it will be used. This way, a class - * implementing a multi-phase mixture is responsible for - * maintaining the global elements list for the mixture, and no - * static global element list is required. - */ - Constituents::Constituents(Elements* ptr_Elements) : - m_kk(0), - m_speciesFrozen(false) , - m_Elements(ptr_Elements) - { - if (!m_Elements) m_Elements = new Elements(); - - // Register subscription to Elements object whether or not we - // created it here. - m_Elements->subscribe(); - } - - /** - * Destructor for class Constituents. - * - * Some cleanup of of the Global_Elements_List array is - * effected by unsubscribing to m_Elements. - */ - Constituents::~Constituents() - { - int ileft = m_Elements->unsubscribe(); - /* - * Here we may delete Elements Objects or not. Right now, we - * will delete them. We also delete the global pointer entry - * to keep everything consistent. - */ - if (ileft <= 0) { - vector::iterator it; - for (it = Elements::Global_Elements_List.begin(); - it != Elements::Global_Elements_List.end(); ++it) { - if (*it == m_Elements) { - Elements::Global_Elements_List.erase(it); - break; - } + if (ileft <= 0) { + vector::iterator it; + for (it = Elements::Global_Elements_List.begin(); + it != Elements::Global_Elements_List.end(); ++it) { + if (*it == m_Elements) { + Elements::Global_Elements_List.erase(it); + break; } - delete m_Elements; } - } - - int Constituents::nElements() const { return m_Elements->nElements(); } - - - /** - * Return the Atomic weight of element m. - * units = Kg / Kmol - */ - doublereal Constituents::atomicWeight(int m) const { - return m_Elements->atomicWeight(m); + delete m_Elements; } + } + + int Constituents::nElements() const { return m_Elements->nElements(); } + + + /** + * Return the Atomic weight of element m. + * units = Kg / Kmol + */ + doublereal Constituents::atomicWeight(int m) const { + return m_Elements->atomicWeight(m); + } doublereal Constituents::entropyElement298(int m) const { return m_Elements->entropyElement298(m); } - /** - * returns a reference to the vector of atomic weights pertinent - * to this constituents object - * units = kg / Kmol - */ - const vector_fp& Constituents::atomicWeights() const { - return m_Elements->atomicWeights(); + /** + * returns a reference to the vector of atomic weights pertinent + * to this constituents object + * units = kg / Kmol + */ + const vector_fp& Constituents::atomicWeights() const { + return m_Elements->atomicWeights(); + } + + + /** + * Return the atomic number of element m. + */ + int Constituents::atomicNumber(int m) const { + return m_Elements->atomicNumber(m); + } + + + /** + * Add an element to the set. + * @param symbol symbol string + * @param weight atomic weight in kg/mol. + * + * If weight is not given, then a lookup is performed in the + * element object + * + */ + void Constituents:: + addElement(const std::string& symbol, doublereal weight) + { + m_Elements->addElement(symbol, weight); + } + + void Constituents:: + addElement(const XML_Node& e) + { + m_Elements->addElement(e); + } + + /* + * Add a unique element to the set. A check on the symbol is made + * If the symbol is already an element, then a new element is + * not created. + * + * @param symbol symbol string + * @param weight atomic weight in kg/mol. + * + * If weight is not given, then a lookup is performed in the + * element object + * + * -> Passthrough to the Element lvl. + */ + void Constituents:: + addUniqueElement(const std::string& symbol, doublereal weight, + int atomicNumber, doublereal entropy298) + { + m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298); + } + + void Constituents:: + addUniqueElement(const XML_Node& e) + { + m_Elements->addUniqueElement(e); + } + + void Constituents::addElementsFromXML(const XML_Node& phase) { + m_Elements->addElementsFromXML(phase); + } + + /* + * -> Passthrough to the Element lvl. + */ + void Constituents::freezeElements() { + m_Elements->freezeElements(); + } + + /* + * -> Passthrough to the Element lvl. + */ + bool Constituents::elementsFrozen() { + return m_Elements->elementsFrozen(); + } + + /* + * Index of element named \a name. The index is an integer + * assigned to each element in the order it was added, + * beginning with 0 for the first element. If \a name is not + * the name of an element in the set, then the value -1 is + * returned. + * + * + * -> Passthrough to the Element class. + */ + int Constituents::elementIndex(std::string name) const { + return (m_Elements->elementIndex(name)); + } + + /* + * Name of the element with index m. + * + * This is a passthrough routine to the Element object. + * @param m @{ Element index. @} + * \exception If m < 0 or m >= nElements(), the + * exception, ElementRangeError, is thrown. + */ + string Constituents::elementName(int m) const { + return (m_Elements->elementName(m)); + } + + /******************************************************************* + * + * elementNames(): + * + * Returns a read-only reference to the vector of element names. + * @code + * Constituents c; + * ... + * const vector& enames = c.elementNames(); + * int n = enames.size(); + * for (int i = 0; i < n; i++) cout << enames[i] << endl; + * @endcode + * + * + * -> Passthrough to the Element lvl. + */ + const vector& Constituents::elementNames() const { + return m_Elements->elementNames(); + } + + /********************************************************************** + * + * molecularWeight() + * + * Returns the molecular weight of a species given the species index + * + * units = kg / kmol. + */ + doublereal Constituents::molecularWeight(int k) const { + if (k < 0 || k >= nSpecies()) { + throw SpeciesRangeError("Constituents::molecularWeight", + k, nSpecies()); } + return m_weight[k]; + } + /********************************************************************** + * + * molecularWeights() + * + * Returns a const reference to the vector of molecular weights + * for all of the species defined in the object. + * + * units = kg / kmol. + */ + const array_fp& Constituents::molecularWeights() const { + return m_weight; + } - /** - * Return the atomic number of element m. - */ - int Constituents::atomicNumber(int m) const { - return m_Elements->atomicNumber(m); + /********************************************************************** + * + * charge(): + * + * Electrical charge of one species k molecule, divided by + * \f$ e = 1.602 \times 10^{-19}\f$ Coulombs. + */ + doublereal Constituents::charge(int k) const { + return m_speciesCharge[k]; + } + + /* + * addSpecies() + * + * Add a species to a Constituents object. Note, no check is made + * as to whether the species has a unique name. + * + * Input + * --------- + * name = string containing the name + * comp[] + * charge = + * weight = weight of the species. Default = 0.0. + * Note, the weight is a bit redundent and potentially + * harmful. If weight is less than or equal to zero, + * the weight is calculated from the element composition + * and it need not be supplied on the command line. + */ + void Constituents:: + addSpecies(const std::string& name, const doublereal* comp, + doublereal charge, doublereal size) { + m_Elements->freezeElements(); + m_speciesNames.push_back(name); + m_speciesCharge.push_back(charge); + m_speciesSize.push_back(size); + int m_mm = m_Elements->nElements(); + // Create a changeable copy of the element composition. We now change the charge potentially + vector_fp compNew(m_mm); + for (int m = 0; m < m_mm; m++) { + compNew[m] = comp[m]; } - - - /** - * Add an element to the set. - * @param symbol symbol string - * @param weight atomic weight in kg/mol. - * - * If weight is not given, then a lookup is performed in the - * element object - * - */ - void Constituents:: - addElement(const std::string& symbol, doublereal weight) - { - m_Elements->addElement(symbol, weight); - } - - void Constituents:: - addElement(const XML_Node& e) - { - m_Elements->addElement(e); - } - - /* - * Add a unique element to the set. A check on the symbol is made - * If the symbol is already an element, then a new element is - * not created. - * - * @param symbol symbol string - * @param weight atomic weight in kg/mol. - * - * If weight is not given, then a lookup is performed in the - * element object - * - * -> Passthrough to the Element lvl. - */ - void Constituents:: - addUniqueElement(const std::string& symbol, doublereal weight, - int atomicNumber, doublereal entropy298) - { - m_Elements->addUniqueElement(symbol, weight, atomicNumber, entropy298); - } - - void Constituents:: - addUniqueElement(const XML_Node& e) - { - m_Elements->addUniqueElement(e); - } - - void Constituents::addElementsFromXML(const XML_Node& phase) { - m_Elements->addElementsFromXML(phase); - } - - /* - * -> Passthrough to the Element lvl. - */ - void Constituents::freezeElements() { - m_Elements->freezeElements(); - } - - /* - * -> Passthrough to the Element lvl. - */ - bool Constituents::elementsFrozen() { - return m_Elements->elementsFrozen(); - } - - /* - * Index of element named \a name. The index is an integer - * assigned to each element in the order it was added, - * beginning with 0 for the first element. If \a name is not - * the name of an element in the set, then the value -1 is - * returned. - * - * - * -> Passthrough to the Element class. - */ - int Constituents::elementIndex(std::string name) const { - return (m_Elements->elementIndex(name)); - } - - /* - * Name of the element with index m. - * - * This is a passthrough routine to the Element object. - * @param m @{ Element index. @} - * \exception If m < 0 or m >= nElements(), the - * exception, ElementRangeError, is thrown. - */ - string Constituents::elementName(int m) const { - return (m_Elements->elementName(m)); - } - - /******************************************************************* - * - * elementNames(): - * - * Returns a read-only reference to the vector of element names. - * @code - * Constituents c; - * ... - * const vector& enames = c.elementNames(); - * int n = enames.size(); - * for (int i = 0; i < n; i++) cout << enames[i] << endl; - * @endcode - * - * - * -> Passthrough to the Element lvl. - */ - const vector& Constituents::elementNames() const { - return m_Elements->elementNames(); - } - - /********************************************************************** - * - * molecularWeight() - * - * Returns the molecular weight of a species given the species index - * - * units = kg / kmol. - */ - doublereal Constituents::molecularWeight(int k) const { - if (k < 0 || k >= nSpecies()) { - throw SpeciesRangeError("Constituents::molecularWeight", - k, nSpecies()); - } - return m_weight[k]; - } - - /********************************************************************** - * - * molecularWeights() - * - * Returns a const reference to the vector of molecular weights - * for all of the species defined in the object. - * - * units = kg / kmol. - */ - const array_fp& Constituents::molecularWeights() const { - return m_weight; - } - - /********************************************************************** - * - * charge(): - * - * Electrical charge of one species k molecule, divided by - * \f$ e = 1.602 \times 10^{-19}\f$ Coulombs. - */ - doublereal Constituents::charge(int k) const { - return m_speciesCharge[k]; - } - - /* - * - * addSpecies() - * - * Add a species to a Constituents object. Note, no check is made - * as to whether the species has a unique name. - * - * Input - * --------- - * name = string containing the name - * comp[] - * charge = - * weight = weight of the species. Default = 0.0. - * Note, the weight is a bit redundent and potentially - * harmful. If weight is less than or equal to zero, - * the weight is calculated from the element composition - * and it need not be supplied on the command line. - */ - void Constituents:: - addSpecies(const std::string& name, doublereal* comp, - doublereal charge, doublereal size) { - m_Elements->freezeElements(); - m_speciesNames.push_back(name); - m_speciesCharge.push_back(charge); - m_speciesSize.push_back(size); - vector_fp compNew;; - double wt = 0.0; - int m_mm = m_Elements->nElements(); - const vector_fp &aw = m_Elements->atomicWeights(); - if (charge != 0.0) { - int eindex = m_Elements->elementIndex("E"); - if (eindex >= 0) { - doublereal ecomp = comp[eindex]; - if (fabs (charge + ecomp) > 0.001) { - if (ecomp != 0.0) { - throw CanteraError("Constituents::addSpecies", - "Input charge and element E compositions differ for species " + name); - } else { - // Just fix up the element E composition based on the input species charge - comp[eindex] = -charge; - } - } - } else { - m_Elements->m_elementsFrozen = false; - addUniqueElement("E", 0.000545, 0, 0.0); - m_Elements->m_elementsFrozen = true; - m_mm = m_Elements->nElements(); - if (m_kk > 0) { - vector_fp old(m_speciesComp); - m_speciesComp.resize(m_kk*m_mm, 0.0); - for (int k = 0; k < m_kk; k++) { - int m_old = m_mm - 1; - for (int m = 0; m < m_old; m++) { - m_speciesComp[k * m_mm + m] = old[k * (m_old) + m]; - } - m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0; - } + double wt = 0.0; + const vector_fp &aw = m_Elements->atomicWeights(); + if (charge != 0.0) { + int eindex = m_Elements->elementIndex("E"); + if (eindex >= 0) { + doublereal ecomp = compNew[eindex]; + if (fabs (charge + ecomp) > 0.001) { + if (ecomp != 0.0) { + throw CanteraError("Constituents::addSpecies", + "Input charge and element E compositions differ for species " + name); + } else { + // Just fix up the element E composition based on the input species charge + compNew[eindex] = -charge; } - eindex = m_Elements->elementIndex("E"); - compNew.resize(m_mm); - for (int m = 0; m < m_mm-1; m++) { - compNew[m] = comp[m]; - } - compNew[m_mm-1] = - charge; - comp = DATA_PTR(compNew); - //comp[eindex] = -charge; - // throw CanteraError("Constituents::addSpecies", - // "Element List doesn't include E, yet this species has charge:" + name); - } - } - for (int m = 0; m < m_mm; m++) { - m_speciesComp.push_back(comp[m]); - wt += comp[m] * aw[m]; - } - m_weight.push_back(wt); - m_kk++; - } - - /* - * - * addUniqueSpecies(): - * - * Add a species to a Constituents object. This routine will - * first check to see if the species is already part of the - * phase. It does this via a string comparison with the - * existing species in the phase. - */ - void Constituents:: - addUniqueSpecies(const std::string& name, doublereal* comp, - doublereal charge, doublereal size) { - vector::const_iterator it = m_speciesNames.begin(); - for (int k = 0; k < m_kk; k++) { - if (*it == name) { - /* - * We have found a match. At this point we could do some - * compatibility checks. However, let's just return for the - * moment without specifying any error. - */ - int m_mm = m_Elements->nElements(); - for (int i = 0; i < m_mm; i++) { - if (comp[i] != m_speciesComp[m_kk * m_mm + i]) { - throw CanteraError("addUniqueSpecies", - "Duplicate species have different " - "compositions: " + *it); + } + } else { + m_Elements->m_elementsFrozen = false; + addUniqueElement("E", 0.000545, 0, 0.0); + m_Elements->m_elementsFrozen = true; + m_mm = m_Elements->nElements(); + if (m_kk > 0) { + vector_fp old(m_speciesComp); + m_speciesComp.resize(m_kk*m_mm, 0.0); + for (int k = 0; k < m_kk; k++) { + int m_old = m_mm - 1; + for (int m = 0; m < m_old; m++) { + m_speciesComp[k * m_mm + m] = old[k * (m_old) + m]; } + m_speciesComp[k * (m_mm) + (m_mm-1)] = 0.0; } - if (charge != m_speciesCharge[m_kk]) { + } + eindex = m_Elements->elementIndex("E"); + compNew.resize(m_mm); + compNew[m_mm-1] = - charge; + //comp[eindex] = -charge; + // throw CanteraError("Constituents::addSpecies", + // "Element List doesn't include E, yet this species has charge:" + name); + } + } + for (int m = 0; m < m_mm; m++) { + m_speciesComp.push_back(compNew[m]); + wt += compNew[m] * aw[m]; + } + m_weight.push_back(wt); + m_kk++; + } + + /* + * + * addUniqueSpecies(): + * + * Add a species to a Constituents object. This routine will + * first check to see if the species is already part of the + * phase. It does this via a string comparison with the + * existing species in the phase. + */ + void Constituents:: + addUniqueSpecies(const std::string& name, const doublereal* comp, + doublereal charge, doublereal size) { + vector::const_iterator it = m_speciesNames.begin(); + for (int k = 0; k < m_kk; k++) { + if (*it == name) { + /* + * We have found a match. At this point we could do some + * compatibility checks. However, let's just return for the + * moment without specifying any error. + */ + int m_mm = m_Elements->nElements(); + for (int i = 0; i < m_mm; i++) { + if (comp[i] != m_speciesComp[m_kk * m_mm + i]) { throw CanteraError("addUniqueSpecies", "Duplicate species have different " - "charges: " + *it); + "compositions: " + *it); } - if (size != m_speciesSize[m_kk]) { - throw CanteraError("addUniqueSpecies", - "Duplicate species have different " - "sizes: " + *it); - } - return; - } - ++it; + } + if (charge != m_speciesCharge[m_kk]) { + throw CanteraError("addUniqueSpecies", + "Duplicate species have different " + "charges: " + *it); + } + if (size != m_speciesSize[m_kk]) { + throw CanteraError("addUniqueSpecies", + "Duplicate species have different " + "sizes: " + *it); + } + return; } - addSpecies(name, comp, charge, size); + ++it; } + addSpecies(name, comp, charge, size); + } - /* - * - * freezeSpecies() - * Set the boolean indicating that we are no longer allowing - * species to be added to the Constituents class object. - */ - void Constituents::freezeSpecies() { - m_speciesFrozen = true; - } + /* + * + * freezeSpecies() + * Set the boolean indicating that we are no longer allowing + * species to be added to the Constituents class object. + */ + void Constituents::freezeSpecies() { + m_speciesFrozen = true; + } - /* - * - * speciesIndex() - * - * Index of species named \c name. The first species added - * will have index 0, and the last one index nSpecies() - 1. - * - * Note, the [] operator shouldn't be used for map's because it - * creates new entries. Here, we use find() to look up entries. - * - * If name isn't in the list, then a -1 is returned. - */ - int Constituents::speciesIndex(std::string name) const { - vector::const_iterator it = m_speciesNames.begin(); - for (int k = 0; k < m_kk; k++) { - if (*it == name) { - /* - * We have found a match. - */ - return k; - } - ++it; + /* + * + * speciesIndex() + * + * Index of species named \c name. The first species added + * will have index 0, and the last one index nSpecies() - 1. + * + * Note, the [] operator shouldn't be used for map's because it + * creates new entries. Here, we use find() to look up entries. + * + * If name isn't in the list, then a -1 is returned. + */ + int Constituents::speciesIndex(std::string name) const { + vector::const_iterator it = m_speciesNames.begin(); + for (int k = 0; k < m_kk; k++) { + if (*it == name) { + /* + * We have found a match. + */ + return k; } - return -1; + ++it; } + return -1; + } - /* - * - * speciesName() - * - * Name of the species with index k - */ - string Constituents::speciesName(int k) const { - if (k < 0 || k >= nSpecies()) - throw SpeciesRangeError("Constituents::speciesName", - k, nSpecies()); - return m_speciesNames[k]; - } + /* + * + * speciesName() + * + * Name of the species with index k + */ + string Constituents::speciesName(int k) const { + if (k < 0 || k >= nSpecies()) + throw SpeciesRangeError("Constituents::speciesName", + k, nSpecies()); + return m_speciesNames[k]; + } - /* - * - * speciesNames() - * - * Return a const reference to the vector of species names - */ - const vector& Constituents::speciesNames() const { - return m_speciesNames; - } + /* + * + * speciesNames() + * + * Return a const reference to the vector of species names + */ + const vector& Constituents::speciesNames() const { + return m_speciesNames; + } - /* - * - * ready(): - * True if both elements and species have been frozen - */ - bool Constituents::ready() const { - return (m_Elements->elementsFrozen() && m_speciesFrozen); - } + /* + * + * ready(): + * True if both elements and species have been frozen + */ + bool Constituents::ready() const { + return (m_Elements->elementsFrozen() && m_speciesFrozen); + } /* * Returns the number of atoms of element \c m in species \c k. diff --git a/Cantera/src/thermo/Constituents.h b/Cantera/src/thermo/Constituents.h index 65195b558..8cba4a52f 100755 --- a/Cantera/src/thermo/Constituents.h +++ b/Cantera/src/thermo/Constituents.h @@ -264,7 +264,7 @@ namespace Cantera { * They are not usually called by user programs. */ //@{ - void addSpecies(const std::string& name, doublereal* comp, + void addSpecies(const std::string& name, const doublereal* comp, doublereal charge = 0.0, doublereal size = 1.0); //! Add a species to the phase, checking for uniqueness of the name @@ -278,7 +278,7 @@ namespace Cantera { * @param charge Charge of the species. Defaults to zero. * @param size Size of the species (meters). Defaults to 1 meter. */ - void addUniqueSpecies(const std::string& name, doublereal* comp, + void addUniqueSpecies(const std::string& name, const doublereal* comp, doublereal charge = 0.0, doublereal size = 1.0);