Added ability for solver to account for leakage in the vessel, as well as relevant test cases for validation.

Code/Source/.vscode/settings.json

@@ -0,0 +1,6 @@
+{
+    "files.associations": {
+        "iostream": "cpp",
+        "ostream": "cpp"
+    }
+}

Code/Source/cvOneDBFSolver.cxx

@@ -1119,6 +1119,8 @@ void cvOneDBFSolver::QuerryModelInformation(void)
       jointList.push_back(feajoint);
     }
 
+    cout << "Line 1122" << endl;
+
     // Loop on Segments
     long temp = 0;
     for (i=0; i<is; i++){
@@ -1468,10 +1470,11 @@ void cvOneDBFSolver::GenerateSolution(void){
 
       // Check Newton-Raphson Convergence
       if((currentTime != deltaTime || (currentTime == deltaTime && iter != 0)) && normf < convCriteria && norms < convCriteria){
-        cout << "    iter: " << std::to_string(iter) << " ";
-        cout << "normf: " << normf << " ";
-        cout << "norms: " << norms << " ";
-        cout << "time: " << ((float)(tend_iter-tstart_iter))/CLOCKS_PER_SEC << endl;
+        // !!!!! REMEMBER TO UNCOMMENT THIS !!!!!
+        // cout << "    iter: " << std::to_string(iter) << " ";
+        // cout << "normf: " << normf << " ";
+        // cout << "norms: " << norms << " ";
+        // cout << "time: " << ((float)(tend_iter-tstart_iter))/CLOCKS_PER_SEC << endl;
         break;
       }
 
@@ -1512,10 +1515,10 @@ void cvOneDBFSolver::GenerateSolution(void){
          }
         }
 
-        if(negArea==1) {
+      if(negArea==1) {
         postprocess_Text();
         assert(0);
-        }
+      }
 
       if(cvOneDGlobal::debugMode){
         printf("(Debug) Printing Solution...\n");

Code/Source/cvOneDMaterial.h

@@ -1,34 +1,34 @@
-/* Copyright (c) Stanford University, The Regents of the University of
- *               California, and others.
- *
- * All Rights Reserved.
- *
- * See Copyright-SimVascular.txt for additional details.
- *
- * Permission is hereby granted, free of charge, to any person obtaining
- * a copy of this software and associated documentation files (the
- * "Software"), to deal in the Software without restriction, including
- * without limitation the rights to use, copy, modify, merge, publish,
- * distribute, sublicense, and/or sell copies of the Software, and to
- * permit persons to whom the Software is furnished to do so, subject
- * to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included
- * in all copies or substantial portions of the Software.
- *
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
- * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
- * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
- * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
- * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
- */
-
+/* Copyright (c) Stanford University, The Regents of the University of
+ *               California, and others.
+ *
+ * All Rights Reserved.
+ *
+ * See Copyright-SimVascular.txt for additional details.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject
+ * to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+ * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+ * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+ * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+ * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+ * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+ * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+ * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+ * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+ * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
 #ifndef CVONEDMATERIAL_H
 #define CVONEDMATERIAL_H
 
@@ -82,10 +82,9 @@ class cvOneDMaterial{
     virtual double GetProperty(char* what) const = 0;
     virtual double GetIntegralpS (double area, double z) const = 0;
 
-
-    virtual double GetN(double S) const = 0;//not really dependent on S actually IV 080703
+    virtual double GetN(double S) const = 0;//not really dependent on S actually IV 080703
     virtual double GetEHR(double z) const = 0;
-
+
     virtual double GetOutflowFunction(double pressure, double z) const = 0;
     virtual double GetDpDz(double area, double z) const = 0;
     virtual double GetDOutflowDp(double pressure, double z) const = 0;

Code/Source/cvOneDMaterialLinear.cxx

@@ -83,12 +83,16 @@ cvOneDMaterialLinear& cvOneDMaterialLinear::operator=(const cvOneDMaterialLinear
     ehr = that.ehr;
     p1_=that.PP1_; //impose P refrence here otherwise p1_ is not the value set in the input file.
     printf("this that set ehr =%f p1_=%f \n",ehr, p1_);
+    P_ambient = that.P_ambient;
+    L_P = that.L_P;
   }
   return *this;
 }
 
 void cvOneDMaterialLinear::SetEHR(double ehr_val, double pref_val){
-  ehr = ehr_val;
+  ehr = 4.0/3.0*ehr_val; // JR 15/11/23: multiplied EHR by the correct factor (since downstream analysis using EHR
+  // in the segmentModel.cxx file does not multiply this by the value, and this is not specified in the documentation that
+  // the user should pre-multiply the E/h/r value by our 4/3 constant.
   PP1_= pref_val;  // add additional commend to set P refrence otherwise p1_ is not the value set in the input file.
 }
 
@@ -163,7 +167,7 @@ double cvOneDMaterialLinear::GetArea(double pressure, double z)const{
   double EHR  = GetEHR(z);//*4/3
 
   // This is the area computation using the "pressure-strain" modulus, EHR.
-  double area = So_*pow(1.0+(pres-p1_)/EHR,2.0);
+  double area = So_/pow(1.0-(pres-p1_)/EHR,2.0);
   if(cvOneDGlobal::debugMode){
     printf("So_: %e\n",So_);
     printf("pres: %e\n",pres);
@@ -181,7 +185,7 @@ double cvOneDMaterialLinear::GetPressure(double S, double z)const{
   // Then we impliment Olufsen's constitutive law...
   double So_   = GetS1(z);
   double EHR   = GetEHR(z);  // From Olufsen's paper
-  double press = p1_ + EHR*(sqrt(S/So_)-1.0);// for linear model dynes/cm^2
+  double press = p1_ + EHR*(1.0-sqrt(So_/S));// for linear model dynes/cm^2
   return press;
 }
 
@@ -190,7 +194,7 @@ double cvOneDMaterialLinear::GetDpDS(double S, double z)const{
   double EHR = GetEHR(z);
   double So_ = GetS1(z);
   double ro  = Getr1(z);
-  double dpds=0.5* EHR/sqrt(So_*S) ;// for linear model
+  double dpds=0.5* EHR* sqrt(So_/S)/S ;// for linear model
   return dpds;
 }
 
@@ -201,7 +205,7 @@ double cvOneDMaterialLinear::GetD2pDS2( double area, double z) const{
 }
 
 double cvOneDMaterialLinear::GetOutflowFunction(double pressure, double z)const{
-  return 0.; // This is not used in our model
+  return L_P * (pressure - P_ambient); // JR: 13-11-23: added outflow model
 }
 
 double cvOneDMaterialLinear::GetDOutflowDp(double pressure, double z)const{

Code/Source/cvOneDMaterialLinear.h

@@ -54,6 +54,9 @@ class cvOneDMaterialLinear:public cvOneDMaterial{
     void   SetStop(double S){Stop = S;}
     void   SetSbottom(double S){Sbot = S;}
     void   SetLength(double length){len = length;}
+    void   SetHydraulicConductivity(double value) {L_P = value;};
+    void   SetStarlingAmbientPressure(double value) {P_ambient = value;};
+    double GetStarlingAmbientPressure() {return P_ambient;}
     double GetProperty( char* what) const;
     double GetArea( double pressure, double z) const;
     double GetPressure( double S, double z) const;
@@ -72,7 +75,7 @@ class cvOneDMaterialLinear:public cvOneDMaterial{
     double GetMette2(double area,double z) const;
     double GetLinCompliance(double z) const;
     double GetnonLinCompliance(double area,double z) const;
-    double GetN(double S) const{return 0.0;};
+    double GetN(double S) const{return N;}; // JR 15/11/23: not sure why this returned 0.0 instead of just returning N. This has now bee fixed.
     void   SetPeriod(double period){};
 
   private:
@@ -81,6 +84,9 @@ class cvOneDMaterialLinear:public cvOneDMaterial{
     double Sbot;
     double len;
 
+    double L_P;
+    double P_ambient;
+
     double ehr;
     double PP1_;
 

Code/Source/cvOneDMaterialManager.cxx

@@ -57,27 +57,31 @@ int cvOneDMaterialManager::AddNewMaterial(MaterialType type, cvOneDMaterial* mat
 }
 
 int cvOneDMaterialManager::AddNewMaterialOlufsen(double density, double dynamicViscosity,
-                                                 double profile_exponent, double pRef,
-                                                 double *params){
+                                                 double profile_exponent, double pRef, double L_P,
+                                                 double P_ambient, double *params){
   cvOneDMaterialOlufsen* olfmat = new cvOneDMaterialOlufsen();
   olfmat->SetDensity(density);
   olfmat->SetDynamicViscosity(dynamicViscosity);
   olfmat->SetProfileExponent(profile_exponent);
   olfmat->SetReferencePressure(pRef);
+  olfmat->SetHydraulicConductivity(L_P);
+  olfmat->SetStarlingAmbientPressure(P_ambient);
   olfmat->SetMaterialType(params,pRef);
   printf("new cvOneMaterialOlufsen called check pRef %f \n", olfmat->GetReferencePressure());
   return cvOneDGlobal::gMaterialManager->AddNewMaterial(MaterialType_MATERIAL_OLUFSEN,(cvOneDMaterial*)olfmat);
 }
 
 int cvOneDMaterialManager::AddNewMaterialLinear(double density, double dynamicViscosity,
-                                                double profile_exponent, double pRef,
-                                                double EHR){
+                                                double profile_exponent, double pRef, double L_P,
+                                                double P_ambient, double EHR){
   cvOneDMaterialLinear* linearmat = new cvOneDMaterialLinear();
   linearmat->SetDensity(density);
   linearmat->SetDynamicViscosity(dynamicViscosity);
   linearmat->SetProfileExponent(profile_exponent);
   linearmat->SetReferencePressure(pRef);
   linearmat->SetEHR(EHR,pRef);
+  linearmat->SetHydraulicConductivity(L_P);
+  linearmat->SetStarlingAmbientPressure(P_ambient);
   return cvOneDGlobal::gMaterialManager->AddNewMaterial(MaterialType_MATERIAL_LINEAR,(cvOneDMaterial*)linearmat);
 }
 
@@ -88,6 +92,7 @@ cvOneDMaterial* cvOneDMaterialManager::GetNewInstance(int matID){
   //  printf("In GetNewInstance cvOneDMaterialOlufsen is called  matID=%i \n",matID);
     *olfmat = *((cvOneDMaterialOlufsen*)(materials[matID]));
   //  printf("In GetNewInstance cvOneDMaterialOlufsen* materials is called \n");
+    olfmat->GetStarlingAmbientPressure();
     return (cvOneDMaterial*)olfmat;
   }else if (types[matID] == MaterialType_MATERIAL_LINEAR) {
     cvOneDMaterialLinear* linearmat = new cvOneDMaterialLinear();

Code/Source/cvOneDMaterialManager.h

@@ -45,10 +45,12 @@ class cvOneDMaterialManager{
     int AddNewMaterial(MaterialType type, cvOneDMaterial* mat);
 
     int AddNewMaterialOlufsen(double density, double dynamicViscosity,
-                              double profile_exponent, double pRef, double *params);
+                              double profile_exponent, double pRef, 
+                              double L_P, double P_ambient, double *params);
 
     int AddNewMaterialLinear(double density, double dynamicViscosity,
-                             double profile_exponent, double pRef, double EHR);
+                             double profile_exponent, double pRef,
+                             double L_P, double P_ambient, double EHR);
 
     // caller must deallocate material instance to avoid memory leak
     cvOneDMaterial* GetNewInstance(int matID);

Code/Source/cvOneDMaterialOlufsen.cxx

@@ -63,35 +63,37 @@ cvOneDMaterialOlufsen::cvOneDMaterialOlufsen(){
  // printf("call cvOneMaterialsOlufsen p1_=%f K3_=%f \n",p1_,K3_);
 }
 
-cvOneDMaterialOlufsen::~cvOneDMaterialOlufsen(){
-}
+cvOneDMaterialOlufsen::~cvOneDMaterialOlufsen(){}
 
 cvOneDMaterialOlufsen::cvOneDMaterialOlufsen (const cvOneDMaterialOlufsen &rhs){
  // this seems not used
   cvOneDMaterial::operator=(rhs);
   double k1 = 0,k2 = 0,k3 = 0;
-  double pref=0;
-  rhs.GetParams( &k1, &k2, &k3, &pref);
+  double pref=0; double P_amb_ref=0; double L_P_ref = 0;
+  rhs.GetParams( &k1, &k2, &k3, &pref, &P_amb_ref, &L_P_ref);
   K1_ = k1;
   K2_ = k2;
   K3_ = k3;
   p1_= pref;
+  P_ambient = P_amb_ref;
+  L_P = L_P_ref;
  // printf("call cvOneDmaterialOlufsen &rhs\n");
 
 }
 
-
 cvOneDMaterialOlufsen& cvOneDMaterialOlufsen::operator= (const cvOneDMaterialOlufsen &that){
   if (this != &that) {
     cvOneDMaterial::operator=(that);
     double k1 = 0,k2 = 0,k3 = 0;
 
-    double pref=0;
-    that.GetParams( &k1, &k2, &k3, &pref);
+    double pref=0; double P_amb_ref=0; double L_P_ref = 0;
+    that.GetParams( &k1, &k2, &k3, &pref, &P_amb_ref, &L_P_ref);
     K1_ = k1;
     K2_ = k2;
     K3_ = k3;
     p1_= pref;
+    P_ambient = P_amb_ref;
+    L_P = L_P_ref;
   //  printf("call cvOneDMaterialOlufsen that this K3_=%f p1_=%f \n",K3_,p1_ );
   }
   return *this;
@@ -102,8 +104,8 @@ void cvOneDMaterialOlufsen::SetMaterialType(double *mType,double Pref){
   K2_ = mType[1];
   K3_ = mType[2];
   PP1_=Pref;
-  cout<< "Setting material K's "<< K1_ <<" "<< K2_<<" "<< K3_<< " ..." << endl;
-  cout<< "Setting reference Pressure "<< PP1_<<endl;
+  // cout<< "Setting material K's "<< K1_ <<" "<< K2_<<" "<< K3_<< " ..." << endl;
+  // cout<< "Setting reference Pressure "<< PP1_<<endl;
  //  printf("call SetMaterialType K3_ %f \n",K3_);
 }
 
@@ -112,7 +114,6 @@ void cvOneDMaterialOlufsen::SetPeriod(double per){
 }
 
 double cvOneDMaterialOlufsen::GetProperty(char* what)const{
-
   // Nothing to change here...
   if( strcmp( what, "density") == 0)
     return density;
@@ -132,6 +133,18 @@ double cvOneDMaterialOlufsen::GetProperty(char* what)const{
   }
 }
 
+void cvOneDMaterialOlufsen::SetHydraulicConductivity(double value) {
+  L_P = value;
+}
+
+void cvOneDMaterialOlufsen::SetStarlingAmbientPressure(double value) {
+  P_ambient = value; 
+}
+
+double cvOneDMaterialOlufsen::GetStarlingAmbientPressure() {
+  return P_ambient;
+}
+
 double cvOneDMaterialOlufsen::GetEHR(double z)const{
   double ro = Getr1(z);
   double ans =4./3.*( K1_*exp(K2_*ro) + K3_);//dyne/cm^2=g/cm/s^2
@@ -184,7 +197,6 @@ double cvOneDMaterialOlufsen::GetArea(double pressure, double z)const{
   //         This property comes from the subdomain
   //
   //       o Po is the the zero transmural pressure
-
   double pres = pressure;
   double So_  = GetS1(z);
   double EHR  = GetEHR(z);
@@ -207,7 +219,6 @@ double cvOneDMaterialOlufsen::GetDpDS(double S, double z)const{
   double So_ = GetS1(z);
   double ro=Getr1(z);
   double dpds=0.5* EHR * sqrt(So_/S)/S ;
-
   return dpds;
 }
 
@@ -218,7 +229,7 @@ double cvOneDMaterialOlufsen::GetD2pDS2(double area, double z)const{
 }
 
 double cvOneDMaterialOlufsen::GetOutflowFunction(double pressure, double z)const{
-  return 0.0; // This is not used in our model
+  return L_P * (pressure - P_ambient);
 }
 
 double cvOneDMaterialOlufsen::GetDOutflowDp(double pressure, double z)const{