diff --git a/Code/Source/cvOneDMaterial.h b/Code/Source/cvOneDMaterial.h index c4e958f..789da9b 100644 --- a/Code/Source/cvOneDMaterial.h +++ b/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; diff --git a/Code/Source/cvOneDMaterialLinear.cxx b/Code/Source/cvOneDMaterialLinear.cxx index c00b1d7..31cf939 100644 --- a/Code/Source/cvOneDMaterialLinear.cxx +++ b/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; } diff --git a/Code/Source/cvOneDMaterialLinear.h b/Code/Source/cvOneDMaterialLinear.h index 83bbfb7..54977f7 100644 --- a/Code/Source/cvOneDMaterialLinear.h +++ b/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: diff --git a/svOneDSolver b/svOneDSolver new file mode 160000 index 0000000..c4b8aac --- /dev/null +++ b/svOneDSolver @@ -0,0 +1 @@ +Subproject commit c4b8aac1156764d9faa17b0c9eb1d6a68ca65c85 diff --git a/tests/cases/bifurcation_R_linear.in b/tests/cases/bifurcation_R_linear.in new file mode 100644 index 0000000..0a0ecb1 --- /dev/null +++ b/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) + + + + + + diff --git a/tests/cases/tube_pressure_linear.in b/tests/cases/tube_pressure_linear.in new file mode 100644 index 0000000..36ee659 --- /dev/null +++ b/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 diff --git a/tests/cases/tube_r_stab_linear.in b/tests/cases/tube_r_stab_linear.in new file mode 100644 index 0000000..57bc47c --- /dev/null +++ b/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 diff --git a/tests/test_integration.py b/tests/test_integration.py index aa957a4..da144ef 100644 --- a/tests/test_integration.py +++ b/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)