fixed bug that leads to 1D linear material not behaving properly in simulations

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,10 @@ 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

@@ -88,7 +88,9 @@ cvOneDMaterialLinear& cvOneDMaterialLinear::operator=(const cvOneDMaterialLinear
 }
 
 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 +165,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 +183,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 +192,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;
 }
 

Code/Source/cvOneDMaterialLinear.h

@@ -72,7 +72,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;};
     void   SetPeriod(double period){};
 
   private:

tests/cases/bifurcation_R_linear.in

@@ -0,0 +1,107 @@
+# Modelled from aortic bifurcation test case from 
+# Xiao, N., Alastruey, J., Figueroa, C.A. A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models. Int J Numer Meth Bio, 2014; 30:204-231
+
+MODEL results_bifurcation_R_linear_
+
+NODE 0 0.0 0.0 0.0
+NODE 1 0.0 0.0 -8.6
+NODE 2 0.0 -3.25280917510326 -7.85297602634594
+NODE 3 0.0 3.25280917510326 -7.85297602634594
+
+JOINT JOINT1 1 INSEGS OUTSEGS
+JOINTINLET INSEGS 1 0
+JOINTOUTLET OUTSEGS 2 1 2
+
+SEGMENT seg0 0 8.6 50 0 1 2.32352192659501 2.32352192659501 0.0 MAT1 NONE 0.0 0 0 NOBOUND NONE
+SEGMENT seg1 1 8.5 50 1 3 1.13097335529233 1.13097335529233 0.0 MAT1 NONE 0.0 0 0 RESISTANCE R_VALS
+SEGMENT seg2 2 8.5 50 1 2 1.13097335529233 1.13097335529233 0.0 MAT1 NONE 0.0 0 0 RESISTANCE R_VALS
+
+DATATABLE R_VALS LIST
+0.0 991.36
+ENDDATATABLE
+
+DATATABLE STEADY_FLOW LIST
+0.0 7.985
+1.0 7.985
+ENDDATATABLE
+
+DATATABLE PULS_FLOW LIST
+0.0 0.0
+0.019668108360095 -4.11549971450822
+0.055247073448669 -7.16517105402019
+0.089913757381125 -1.16130916560675
+0.113633067440174 9.27967911020849
+0.133703252874754 20.4428655623545
+0.150124313684865 32.8883708080991
+0.162896249870507 43.9606301469767
+0.173843623743914 54.6495650354764
+0.186615559929556 67.3593157640345
+0.204861183051901 79.6320314281343
+0.234784004972547 86.0097422945109
+0.26872086398011 78.7468811589053
+0.294264736351393 64.7319679994526
+0.314334921785973 52.023318573387
+0.32893142028385 41.4040199668672
+0.34535248109396 29.7287601931208
+0.363598104216306 17.3938305987427
+0.380019165026417 6.42196693062957
+0.396440225836527 -5.11194051566898
+0.422389556499419 -18.901261909892
+0.455347523345814 -23.8602212809087
+0.496182965572016 -18.2411583475954
+0.533485128399922 -10.4792705793446
+0.575247332435512 -4.00466780439001
+0.640931575675956 -2.79835885098868
+0.682440368279291 -5.30400645008189
+0.726077816913567 -8.52809746163143
+0.762721110017611 -9.03919130533637
+0.802405340308712 -7.30944473123762
+0.846498929521047 -5.55578067364513
+0.885944229033165 -5.87317544184486
+0.925563296384544 -7.14939949266443
+0.968258054490832 -6.10169723438909
+1.00949316274733 -4.10721983893183
+1.05237037708484 -3.31554531122748
+1.087 -1.77083891076532
+ENDDATATABLE
+
+
+# Ref Pressure 1333.22*85 where 85 is in mmHg
+# Rigid for now, but can be elastic with the following parameters:
+# h_0 = 1.032mm, E_0 = 500 kPa, h_1 = h_2 = 0.72 mm, E_1 = E_2 = 700 kPa
+MATERIAL MAT1 LINEAR 1.06 0.04 0 2.0 1.8e10 
+
+SOLVEROPTIONS 0.001087 50 1000 2 STEADY_FLOW FLOW 1.0e-6 1 1 
+
+OUTPUT TEXT
+
+
+# analytical solution
+
+# parameters
+# viscosity		mu	0.04
+# seg0 length		L_0	8.6
+# seg0 radius		r_0	0.86
+# seg1 length		L_1	8.5 (same as seg2 length)
+# seg1 radius		r_1	0.6 (same as seg2 radius)
+# resistance BC		R	991.36
+# steady Inlet Flow	Q_0	7.985
+# Distal pressure	Pd	0
+
+# reference solution
+# assumes no pressure losses at bifurcation and purely parallel resistances
+# total 1D model resistance	Rtot = (R_0 + 0.5*(R_2 + R)
+# Pressure at junction		P_0_out = P_1_in = P_2_in
+
+# Results to be checked
+# P_0_in	3997.46433118937
+# P_0_out	3984.67700709878
+# P_1_in	3984.67700709878 (same as P_2_in)
+# P_1_out	3958.0048 (same as P_2_out because symmetric)
+# Q_1_in	3.9925 (same as Q_1_out)
+
+
+
+
+
+

tests/cases/tube_pressure_linear.in

@@ -0,0 +1,37 @@
+MODEL results_tube_pressure_linear_
+
+NODE 0 0.0 0.0 0.0
+NODE 1 0.0 0.0 10.0
+
+SEGMENT seg0 0 10.0 50 0 1 1.0 1.0 0.0 MAT1 NONE 0.0 0 0 PRESSURE OUTLETDATA
+
+DATATABLE INLETDATA LIST
+0.0 100.0
+10.0 100.0
+ENDDATATABLE
+
+DATATABLE OUTLETDATA LIST
+0.0 10000.0
+ENDDATATABLE
+
+SOLVEROPTIONS 0.001 1000 1000 2 INLETDATA FLOW 1.0e-8 1 1
+
+MATERIAL MAT1 LINEAR 1.06 0.04 0 2.0 1.8e10
+
+OUTPUT TEXT
+
+# analytical solution
+
+# parameters
+# viscosity		mu	0.04
+# vessel length		L	10
+# vessel cross-section	A	1
+# vessel radius		r	0.5641895835
+# pressure outlet       Pout	10000
+# prescribed inflow	Q	100
+
+# reference solution
+# vessel resistance	R1 = 8*mu*L*PI/A^2 = 10.05309649
+
+# results to be checked
+# pressure inlet	Pin = Pout + Q*R1 = 11005.30965

tests/cases/tube_r_stab_linear.in

@@ -0,0 +1,38 @@
+MODEL results_tube_r_stab_linear_
+
+NODE 0 0.0 0.0 0.0
+NODE 1 0.0 0.0 10.0
+
+SEGMENT seg0 0 10.0 50 0 1 1.0 1.0 0.0 MAT1 NONE 0.0 0 0 RESISTANCE OUTLETDATA
+
+DATATABLE INLETDATA LIST
+0.0 100.0
+10.0 100.0
+ENDDATATABLE
+
+DATATABLE OUTLETDATA LIST
+0.0 100.0
+ENDDATATABLE
+
+SOLVEROPTIONS 0.001 500 500 2 INLETDATA FLOW 1.0e-8 1 0
+
+MATERIAL MAT1 LINEAR 1.06 0.04 0 2.0 1.0e11
+
+OUTPUT TEXT
+
+# analytical solution
+
+# parameters
+# viscosity             mu      0.04
+# vessel length         L       10
+# vessel cross-section  A       1
+# vessel radius         r       0.5641895835
+# resistance BC         R2      100
+# prescribed inflow     Q       100
+
+# reference solution
+# vessel resistance     R1 = 8*mu*L*PI/A^2 = 10.05309649
+
+# results to be checked
+# pressure inlet	Pin  = Q*(R1+R2) =	11005.30965
+# pressure outlet	Pout = Q*R2 =		10000

tests/test_integration.py

@@ -103,6 +103,14 @@ def test_tube_pressure(tmpdir):
     assert np.isclose(get_result(results, 'area', 0, -1, -1, 'point'), 1.0, rtol=1.0e-5)
 
 
+def test_tube_pressure_linear(tmpdir):
+    results = run_test_case_by_name('tube_pressure_linear', tmpdir)
+    assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 11005.30965, rtol=1.0e-6)
+    assert np.isclose(get_result(results, 'pressure', 0, -1, -1, 'point'), 10000.0, rtol=1.0e-8)
+    assert np.isclose(get_result(results, 'flow', 0, -1, -1, 'point'), 100.0, rtol=1.0e-16)
+    assert np.isclose(get_result(results, 'area', 0, -1, -1, 'point'), 1.0, rtol=1.0e-5)
+
+
 def test_tube_pressure_wave(tmpdir):
     results = run_test_case_by_name('tube_pressure_wave', tmpdir)
     assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 10000.0 , rtol=1.0e-8)
@@ -151,6 +159,14 @@ def test_tube_r_stab(tmpdir):
     assert np.isclose(get_result(results, 'area', 0, -1, -1, 'point'), 1.0, rtol=1.0e-5)
 
 
+def test_tube_r_stab_linear(tmpdir):
+    results = run_test_case_by_name('tube_r_stab_linear', tmpdir)
+    assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 11005.30965, rtol=1.0e-6)
+    assert np.isclose(get_result(results, 'pressure', 0, -1, -1, 'point'), 10000.0, rtol=1.0e-8)
+    assert np.isclose(get_result(results, 'flow', 0, -1, -1, 'point'), 100.0, rtol=1.0e-16)
+    assert np.isclose(get_result(results, 'area', 0, -1, -1, 'point'), 1.0, rtol=1.0e-5)
+
+
 def test_tube_stenosis_r(tmpdir):
     results = run_test_case_by_name('tube_stenosis_r', tmpdir)
     assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 10150.68211, rtol=1.0e-6)
@@ -183,6 +199,18 @@ def test_bifurcation_R(tmpdir):
     assert np.isclose(get_result(results, 'flow', 2, -1, -1, 'point'), 3.9925, rtol=1e-5)
 
 
+def test_bifurcation_R_linear(tmpdir):
+    results = run_test_case_by_name('bifurcation_R_linear', tmpdir)
+    assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 3997.46433118937, rtol=1e-5)
+    assert np.isclose(get_result(results, 'pressure', 0, -1, -1, 'point'), 3984.67700709878, rtol=1e-5)
+    assert np.isclose(get_result(results, 'pressure', 1, 0, -1, 'point'), 3984.67700709878, rtol=1e-5)
+    assert np.isclose(get_result(results, 'pressure', 2, 0, -1, 'point'), 3984.67700709878, rtol=1e-5)
+    assert np.isclose(get_result(results, 'pressure', 1, -1, -1, 'point'), 3958.0048, rtol=1e-5)
+    assert np.isclose(get_result(results, 'pressure', 2, -1, -1, 'point'), 3958.0048, rtol=1e-5)
+    assert np.isclose(get_result(results, 'flow', 1, -1, -1, 'point'), 3.9925, rtol=1e-5)
+    assert np.isclose(get_result(results, 'flow', 2, -1, -1, 'point'), 3.9925, rtol=1e-5)
+
+
 def test_bifurcation_R_stab(tmpdir):
     results = run_test_case_by_name('bifurcation_R_stab', tmpdir)
     assert np.isclose(get_result(results, 'pressure', 0, 0, -1, 'point'), 3997.46433118937, rtol=1e-6)