vcs_VolPhase upgrade: making more data private

Cantera/src/equil/vcs_Exception.h

@@ -38,7 +38,13 @@ protected:
  *              where the assertion failed
  * @ingroup errorhandling
  */
-#define AssertThrowVCS(expr, proc)  ((expr) ? (void) 0 : throw vcsError(proc, std::string("failed assert: ") + #expr,-1))
+#define AssertThrowVCS(expr, proc)  ((expr) ? (void) 0 : throw vcsError(proc, std::string("failed Assert: ") + #expr,-1))
+
+#ifdef DEBUG_HKM
+#define DebugAssertThrowVCS(expr, proc)  ((expr) ? (void) 0 : throw vcsError(proc, std::string("failed debugAssert: ") + #expr,-1))
+#else
+#define DebugAssertThrowVCS(expr, proc)  
+#endif
 
 }
 

Cantera/src/equil/vcs_MultiPhaseEquil.cpp

@@ -1007,7 +1007,8 @@ namespace VCSnonideal {
 
     VolPhase->m_speciesUnknownType.resize(ns, VCS_SPECIES_TYPE_MOLNUM);
 
-    VolPhase->ElGlobalIndex.resize(ne, -1);
+    VolPhase->elemResize(ne);
+    //ElGlobalIndex.resize(ne, -1);
 
     
     e = 0;
@@ -1115,6 +1116,7 @@ namespace VCSnonideal {
        */
       iSurPhase = -1;
       tPhase = &(mphase->phase(iphase));
+      int nelem = tPhase->nElements();
     
       /*
        * Query Cantera for the equation of state type of the
@@ -1144,7 +1146,7 @@ namespace VCSnonideal {
        */
 
       vcs_VolPhase *VolPhase = vprob->VPhaseList[iphase];
-      VolPhase->resize(iphase, nSpPhase, phaseName.c_str(), 0.0);
+      VolPhase->resize(iphase, nSpPhase, nelem, phaseName.c_str(), 0.0);
       VolPhase->m_gasPhase = gasPhase;
       /*
        * Tell the vcs_VolPhase pointer about cantera

Cantera/src/equil/vcs_VolPhase.cpp

@@ -38,10 +38,9 @@ namespace VCSnonideal {
     m_eqnState(VCS_EOS_CONSTANT),
     nElemConstraints(0),
     ChargeNeutralityElement(-1),
-    ElGlobalIndex(0),
+    m_elemGlobalIndex(0),
     NVolSpecies(0),
     m_totalMolesInert(0.0),
-    m_molarVolInert(1000.),
     m_activityConvention(0),
     m_isIdealSoln(false),
     m_existence(0),
@@ -174,7 +173,7 @@ namespace VCSnonideal {
       }
 
       m_speciesUnknownType = b.m_speciesUnknownType;
-      ElGlobalIndex = b.ElGlobalIndex;
+      m_elemGlobalIndex    = b.m_elemGlobalIndex;
       NVolSpecies         = b.NVolSpecies;
       PhaseName           = b.PhaseName;
       m_totalMolesInert   = b.m_totalMolesInert;
@@ -243,7 +242,7 @@ namespace VCSnonideal {
   }
   /***************************************************************************/
 
-  void vcs_VolPhase::resize(int phaseNum, int nspecies, const char *phaseName,
+  void vcs_VolPhase::resize(int phaseNum, int nspecies, int numElem, const char *phaseName,
 			    double molesInert) {
     if (nspecies <= 0) {
       plogf("nspecies Error\n");
@@ -280,7 +279,7 @@ namespace VCSnonideal {
       m_singleSpecies = true;
     }
 
-    if (NVolSpecies == nspecies) {
+    if (NVolSpecies == nspecies && numElem == nElemConstraints) {
       return;
     }
  
@@ -289,6 +288,7 @@ namespace VCSnonideal {
       m_singleSpecies = false;
     }
 
+
     IndSpecies.resize(nspecies,-1);
 
     if ((int) ListSpeciesPtr.size() >= NVolSpecies) {
@@ -325,6 +325,29 @@ namespace VCSnonideal {
     m_UpToDate_VolPM      = false;
     m_UpToDate_GStar      = false;
     m_UpToDate_G0         = false;
+
+
+    elemResize(numElem);
+
+  }
+  /***************************************************************************/
+
+  void vcs_VolPhase::elemResize(const int numElemConstraints) {
+
+    ElName.resize(numElemConstraints);
+
+    ElActive.resize(numElemConstraints, 1);
+    m_elType.resize(numElemConstraints, VCS_ELEM_TYPE_ABSPOS);
+
+
+    FormulaMatrix.resize(numElemConstraints, NVolSpecies, 0.0);
+
+
+    m_elemGlobalIndex.resize(numElemConstraints, -1);
+
+    
+
+    nElemConstraints = numElemConstraints;
   }
   /***************************************************************************/
 
@@ -1009,11 +1032,12 @@ namespace VCSnonideal {
       m_VCS_UnitsFormat = VCS_UNITS_MKS;
       m_phi = TP_ptr->electricPotential();
       int nsp = TP_ptr->nSpecies();
+      int nelem = TP_ptr->nElements();
       if (nsp !=  NVolSpecies) {
 	if (NVolSpecies != 0) {
 	  plogf("Warning Nsp != NVolSpeces: %d %d \n", nsp, NVolSpecies);
 	}
-	resize(VP_ID, nsp, PhaseName.c_str());
+	resize(VP_ID, nsp, nelem, PhaseName.c_str());
       }
       TP_ptr->getMoleFractions(VCS_DATA_PTR(Xmol));
       _updateMoleFractionDependencies();
@@ -1289,5 +1313,19 @@ namespace VCSnonideal {
     return m_totalMolesInert;
   }
   /**********************************************************************/
+
+   //! Returns the global index of the local element index for the phase
+  int vcs_VolPhase::elemGlobalIndex(const int e) const {
+    DebugAssertThrowVCS(e >= 0, " vcs_VolPhase::elemGlobalIndex") ;
+    DebugAssertThrowVCS(e < nElemConstraints, " vcs_VolPhase::elemGlobalIndex") ;
+    return m_elemGlobalIndex[e];
+  }
+
+   //! Returns the global index of the local element index for the phase
+  void vcs_VolPhase::setElemGlobalIndex(const int eLocal, const int eGlobal) {
+    DebugAssertThrowVCS(eLocal >= 0, "vcs_VolPhase::setElemGlobalIndex");
+    DebugAssertThrowVCS(eLocal < nElemConstraints, "vcs_VolPhase::setElemGlobalIndex");
+    m_elemGlobalIndex[eLocal] = eGlobal;
+  }
 }
 

Cantera/src/equil/vcs_VolPhase.h

@@ -143,9 +143,11 @@ namespace VCSnonideal {
      *  @param phaseName   String name for the phase
      *  @param molesInert  kmoles of inert in the phase (defaults to zero)
      */
-    void resize(int phaseNum, int numSpecies, const char *phaseName,
+    void resize(int phaseNum, int numSpecies, int numElem, const char *phaseName,
 		double molesInert = 0.0);
 
+    void elemResize(const int numElemConstraints);
+
     //! Evaluate activity coefficients and return the kspec coefficient
     /*!
      *   We carry out a calculation whenever UpTODate_AC is false. Specifically
@@ -465,6 +467,16 @@ namespace VCSnonideal {
      */
     double totalMolesInert() const;
 
+    //! Returns the global index of the local element index for the phase
+    int elemGlobalIndex(const int e) const;
+
+    //! sets a local phase element to a global index value
+    /*!
+     * @param eLocal Local phase element index
+     * @param eGlobal Global phase element index
+     */
+    void setElemGlobalIndex(const int eLocal, const int eGlobal);
+
   private:
 
     //! Evaluate the activity coefficients at the current conditions
@@ -635,10 +647,12 @@ namespace VCSnonideal {
      */
     std::vector<int> m_speciesUnknownType;
 
+  private:
     //! Index of the element number in the global list of elements
     //!  storred in VCS_PROB or VCS_SOLVE       
-    std::vector<int> ElGlobalIndex;
+    std::vector<int> m_elemGlobalIndex;
 
+  public:
     //! Number of species in the phase       
     int NVolSpecies;
 
@@ -648,13 +662,8 @@ namespace VCSnonideal {
   private:
     //!  Total moles of inert in the phase 
     double m_totalMolesInert;
-  public:
-    //! molar volume of the inert species if present
-    /*!
-     *  units m**3 / kg
-     */
-    double m_molarVolInert;
 
+  public:
     //! Convention for the activity formulation
     /*!
      *  0 = molar based activities (default)

Cantera/src/equil/vcs_elem_rearrange.cpp

@@ -232,17 +232,17 @@ namespace VCSnonideal {
     }
 #endif  
     /*
-     * Change the element Global Index list in each phase object
+     * Change the element Global Index list in each vcs_VolPhase object
      * to reflect the switch in the element positions.
      */
     for (int iph = 0; iph < m_numPhases; iph++) {
       volPhase = m_VolPhaseList[iph];
       for (int e = 0; e < volPhase->nElemConstraints; e++) {
-	if (volPhase->ElGlobalIndex[e] == ipos) {
-          volPhase->ElGlobalIndex[e] = jpos; 
+	if (volPhase->elemGlobalIndex(e) == ipos) {
+	  volPhase->setElemGlobalIndex(e, jpos);
 	}
-	if (volPhase->ElGlobalIndex[e] == jpos) {
-          volPhase->ElGlobalIndex[e] =ipos; 
+	if (volPhase->elemGlobalIndex(e) == jpos) {
+	  volPhase->setElemGlobalIndex(e, ipos);
 	}
       } 
     }

Cantera/src/equil/vcs_prob.cpp

@@ -378,7 +378,7 @@ namespace VCSnonideal {
       for (e = 0; e < ne; e++) {
 	en = ElName[e];
 	if (!strcmp(enVP.c_str(), en.c_str())) {
-	  volPhase->ElGlobalIndex[eVP] = e;
+	  volPhase->setElemGlobalIndex(eVP, e);
 	  foundPos = e;
 	}
       }
@@ -386,7 +386,7 @@ namespace VCSnonideal {
 	int elType = volPhase->m_elType[eVP];
 	int elactive = volPhase->ElActive[eVP];
 	e = addElement(enVP.c_str(), elType, elactive);
-	volPhase->ElGlobalIndex[eVP] = e;
+	volPhase->setElemGlobalIndex(eVP, e);
       }
     }
   }
@@ -444,10 +444,12 @@ namespace VCSnonideal {
     }
     double *const *const fm = volPhase->FormulaMatrix.baseDataAddr();
     for (eVP = 0; eVP < volPhase->nElemConstraints; eVP++) {
-      e = volPhase->ElGlobalIndex[eVP];
+      e = volPhase->elemGlobalIndex(eVP);
+#ifdef DEBUG_MODE
       if (e < 0) {
-	exit(-1);
+	std::exit(-1);
       }
+#endif
       FormulaMatrix[e][kT] = fm[eVP][k];
     }
     /*