Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
wants to merge 13 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions Code/Source/.vscode/settings.json
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
{
"files.associations": {
"iostream": "cpp",
"ostream": "cpp"
}
}
15 changes: 9 additions & 6 deletions Code/Source/cvOneDBFSolver.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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++){
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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");
Expand Down
67 changes: 33 additions & 34 deletions Code/Source/cvOneDMaterial.h
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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;
Expand Down
14 changes: 9 additions & 5 deletions Code/Source/cvOneDMaterialLinear.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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.
}

Expand Down Expand Up @@ -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);
Expand All @@ -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;
}

Expand All @@ -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;
}

Expand All @@ -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{
Expand Down
8 changes: 7 additions & 1 deletion Code/Source/cvOneDMaterialLinear.h
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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:
Expand All @@ -81,6 +84,9 @@ class cvOneDMaterialLinear:public cvOneDMaterial{
double Sbot;
double len;

double L_P;
double P_ambient;

double ehr;
double PP1_;

Expand Down
13 changes: 9 additions & 4 deletions Code/Source/cvOneDMaterialManager.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand All @@ -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();
Expand Down
6 changes: 4 additions & 2 deletions Code/Source/cvOneDMaterialManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
37 changes: 24 additions & 13 deletions Code/Source/cvOneDMaterialOlufsen.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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_);
}

Expand All @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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);
Expand All @@ -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;
}

Expand All @@ -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{
Expand Down
Loading
Loading