GCC Code Coverage Report

Directory: ./
File: Solver/autoregulation.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 149 202 73.8%
Branches: 146 345 42.3%
Line Branch Exec Source
1 // ______ ______ _ _ _____ ______
2 // | ____| ____| | (_)/ ____| | ____|
3 // | |__ | |__ | | _| (___ ___| |__
4 // | __| | __| | | | |\___ \ / __| __|
5 // | | | |____| |____| |____) | (__| |____
6 // |_| |______|______|_|_____/ \___|______|
7 // Finite Elements for Life Sciences and Engineering
8 //
9 // License: LGL2.1 License
10 // FELiScE default license: LICENSE in root folder
11 //
12 // Main authors:
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/autoregulation.hpp"
21 #include "Core/felisceTools.hpp"
22
23 namespace felisce {
24
25
26 /*! \brief Initialization of the class, parameters taken from FelisceParam
27
28 \param[in] elemBoundaryType std::string with the keyword associated to the boundary element type
29 \return a std::unordered_map that associates to a group of labels defined in a std::string the corresponding labels stored in a std::vector
30 */
31 std::unordered_map <std::string,std::vector< int> >
32 8 Autoregulation::initialize(std::string elemBoundaryType) {
33
34 /*! shape parameters for the function sigmoid() : #m_slope abd #m_center are std::set from data in FelisceParam */
35
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 double k( FelisceParam::instance().q / ( 1 - FelisceParam::instance().q ) );
36
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 m_slope = 2 * std::log(k)/(mmHg2DynPerCm2 * FelisceParam::instance().up - mmHg2DynPerCm2 * FelisceParam::instance().down );
37
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 m_center = (mmHg2DynPerCm2 * FelisceParam::instance().up + mmHg2DynPerCm2 * FelisceParam::instance().down )/2;
38
39 /*! reference values and regulation parameters to handle autoregulation in windkesel are std::set */
40
3/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 4 times.
 8 if ( FelisceParam::instance().useRegulationInWindkessel ) {
41
42
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 m_RdistRef = FelisceParam::instance().lumpedModelBC_Rdist; // resistances of the different outlets
43
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 m_CRef = FelisceParam::instance().lumpedModelBC_C; // capacitances
44
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 m_regParam = FelisceParam::instance().lumpedModelBC_regParam ; // parameter for the regulation part
45
46 // initialization of current quantities
47
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_Rdist.resize(m_RdistRef.size(),0.0);
48
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 m_C.resize(m_CRef.size(),0.0);
49 }
50
51 /*! #m_nC the number of cardiac cycles is computed */
52 8 m_nC = static_cast<int>(std::min(ceil(FelisceParam::instance().timeMax/m_cycleDuration),
53
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 8 ceil(FelisceParam::instance().timeIterationMax*FelisceParam::instance().timeStep/m_cycleDuration)));
54
55 // once that m_nC is known mean inlet quantities are resized
56
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_QInflow.resize(m_nC,0.0);
57
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_PInflow.resize(m_nC,0.0);
58
59 8 std::unordered_map<std::string,std::vector<int> > group2Labels;
60
3/6
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 8 times.
✗ Branch 8 not taken.
 8 std::map<int, std::string> label2Group = FelisceParam::instance().felName2MapintRef2DescLine[ elemBoundaryType ];
61
62 //Here the std::unordered_map stored in FelisceParam is "reversed"
63
2/2
✓ Branch 4 taken 60 times.
✓ Branch 5 taken 8 times.
 68 for(auto it = label2Group.begin(); it != label2Group.end(); ++it) {
64
2/4
✓ Branch 2 taken 60 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 60 times.
✗ Branch 7 not taken.
 60 group2Labels[ it->second ].push_back( it->first );
65 }
66
67 /*! the list of labels: #m_labelInlet, #m_labelOutlet and #m_labelStucture are filled */
68
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_labelInlet = group2Labels["inflow"] ;
69
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_labelOutlet = group2Labels["outflow"] ;
70
3/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 8 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 8 times.
✗ Branch 9 not taken.
 8 m_labelStructure = group2Labels["lateral"] ;
71
72
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
 8 if ( FelisceParam::verbose() > 1 ) {
73 ✗ std::stringstream master;
74 ✗ master<<" inflow labels: ";
75 ✗ for ( std::size_t iLab(0); iLab < m_labelInlet.size(); ++iLab ) {
76 ✗ master<<m_labelInlet[iLab]<<'\t';
77 }
78 ✗ master<<std::endl;
79 ✗ master<<" outflow labels: ";
80 ✗ for ( std::size_t iLab(0); iLab < m_labelOutlet.size(); ++iLab ) {
81 ✗ master<<m_labelOutlet[iLab]<<'\t';
82 }
83 ✗ master<<std::endl;
84 ✗ master<<" structure labels: ";
85 ✗ for ( std::size_t iLab(0); iLab < m_labelStructure.size(); ++iLab ) {
86 ✗ master<<m_labelStructure[iLab]<<'\t';
87 }
88 ✗ master<<std::endl;
89 ✗ PetscPrintf(MpiInfo::petscComm(),"%s",master.str().c_str());
90 }
91 16 return group2Labels;
92 8 }
93
94 /*! \brief post-process function, displays quantities to video and possibly to a summary*/
95 void
96 44 Autoregulation::update( const double& time, double Q, double P, std::vector<double> pout, const int& rank ) {
97
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 33 times.
 44 if( rank == MASTER ) {
98 11 double n = time / m_cycleDuration;
99 11 std::size_t nf = static_cast<std::size_t>(floor(n));
100
101 11 std::cout<<"labelInlet "<<m_labelInlet[0]<<" T "<<time<<" nf "<<nf
102 11 <<" Q "<<Q<<" P [mmHg]"<<P*m_cycleDuration/mmHg2DynPerCm2<<std::endl;
103
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
 11 if( Tools::equalTol(time, (nf+1)*m_cycleDuration , FelisceParam::instance().timeStep/2) ) {
104 ✗ m_QInflow[nf] += FelisceParam::instance().timeStep*Q/2;
105 ✗ m_PInflow[nf] += FelisceParam::instance().timeStep*P/2;
106 ✗ if( nf + 1<m_QInflow.size()) {
107 ✗ m_QInflow[nf+1] += FelisceParam::instance().timeStep*Q/2;
108 ✗ m_PInflow[nf+1] += FelisceParam::instance().timeStep*P/2;
109 }
110 ✗ exportSummary(rank);
111
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
 11 } else if( Tools::equalTol(time, nf*m_cycleDuration , FelisceParam::instance().timeStep/2) ) {
112 ✗ if( nf<m_QInflow.size() ) {
113 ✗ m_QInflow[nf] += FelisceParam::instance().timeStep*Q/2;
114 ✗ m_PInflow[nf] += FelisceParam::instance().timeStep*P/2;
115 }
116 ✗ m_QInflow[nf-1] += FelisceParam::instance().timeStep*Q/2;
117 ✗ m_PInflow[nf-1] += FelisceParam::instance().timeStep*P/2;
118 ✗ exportSummary(rank);
119 } else {
120
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 if( nf<m_QInflow.size() ) {
121 11 m_QInflow[nf] += FelisceParam::instance().timeStep*Q;
122 11 m_PInflow[nf] += FelisceParam::instance().timeStep*P;
123 }
124 }
125
2/2
✓ Branch 1 taken 22 times.
✓ Branch 2 taken 11 times.
 33 for( std::size_t kp(0); kp<m_labelOutlet.size(); kp++) {
126 22 pout[kp]=pout[kp]/mmHg2DynPerCm2;
127 22 std::cout<<"pout label+"<<m_labelOutlet[kp]<<" = "<<pout[kp]<<std::endl;
128 }
129 }
130 44 }
131
132 /*! \brief it returns the pressure at the inlet at the current time step */
133 double
134 169127 Autoregulation::inletPressure(const double &t, const int& caseTest) const
135 {
136 169127 double timeInCycle = ( t - static_cast<int>( t / m_cycleDuration )*m_cycleDuration );
137
138 169127 double peakPressure(0), diastPres(0);
139
2/3
✓ Branch 0 taken 8000 times.
✓ Branch 1 taken 161127 times.
✗ Branch 2 not taken.
 169127 switch(caseTest)
140 {
141 8000 case 0:
142
2/4
✓ Branch 1 taken 8000 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8000 times.
✗ Branch 5 not taken.
 8000 meanPressure2PeakAndDiast( FelisceParam::instance().nominalPressure *mmHg2DynPerCm2, peakPressure, diastPres);
143 8000 break;
144 161127 case 3:
145 case 9:
146 case 10:
147 case 11:
148 case 26:
149
2/4
✓ Branch 1 taken 161127 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 161127 times.
✗ Branch 5 not taken.
 161127 meanPressure2PeakAndDiast( FelisceParam::instance().incomingPressure *mmHg2DynPerCm2, peakPressure, diastPres);
150 161127 break;
151 ✗ default:
152 ✗ FEL_ERROR("non existing (or removed) case required");
153 ✗ break;
154 }
155
156 169127 bool systole = timeInCycle < m_systoleDuration;
157
2/2
✓ Branch 0 taken 53703 times.
✓ Branch 1 taken 115424 times.
 169127 if ( systole )
158 53703 return - peakPressure * std::sin( PI * timeInCycle / m_systoleDuration ) - diastPres;
159 else // diastole
160 115424 return - diastPres;
161 }
162
163 /*! \brief it enforce the autoregulation in the windkessels, by calling setLumpedModel()*/
164 void
165 44 Autoregulation::handleLumpedModel(const double& t, LumpedModelBC* lumpedModel ) {
166
1/2
✓ Branch 1 taken 44 times.
✗ Branch 2 not taken.
 44 if ( FelisceParam::instance().useRegulationInWindkessel ) {
167 //computation of the regulated windkessels
168 44 setLumpedModel(t);
169 //here we change the parameter in FelisceParam ...
170 44 FelisceParam::instance().lumpedModelBC_Rdist = m_Rdist;
171 44 FelisceParam::instance().lumpedModelBC_C = m_C;
172 44 lumpedModel->setRd(m_Rdist);
173 44 lumpedModel->setC(m_C);
174 }
175 44 }
176 /*! \brief for each cardiac cycle the difference between the mean pressure and the reference one is computed*/
177 void
178 8 Autoregulation::computeExtraPressure( const int &nC ) {
179
180 8 double startingTime(0.0);
181 8 double currentP(0.0);
182
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 double nominalP = computeMeanInflowPress( FelisceParam::instance().idControlFreeCase , startingTime);
183
184
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 m_extraPressure.resize(nC);
185
2/2
✓ Branch 0 taken 160 times.
✓ Branch 1 taken 8 times.
 168 for (int k(0); k<nC; k++) {
186
2/4
✓ Branch 1 taken 160 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 160 times.
✗ Branch 5 not taken.
 160 currentP = computeMeanInflowPress( FelisceParam::instance().idCase , k * m_cycleDuration);
187 160 m_extraPressure[k] = nominalP-currentP;
188 }
189 8 }
190
191 /*! \brief it returns the active part of the fibers pre-stress*/
192 double
193 182901 Autoregulation::control( const double &t ) const {
194
1/2
✓ Branch 1 taken 182901 times.
✗ Branch 2 not taken.
 182901 double extraP( interpPressure( t ) );
195
1/2
✓ Branch 1 taken 182901 times.
✗ Branch 2 not taken.
 182901 return FelisceParam::instance().maxT * sigmoid( extraP );
196 }
197
198 /*! \brief post-process function, it writes a summary */
199 void
200 4 Autoregulation::exportSummary( const int& rank ) const {
201
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 3 times.
 4 if( rank == MASTER ) {
202
5/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
 2 std::string nameFile = FelisceParam::instance().resultDir + "../" + FelisceParam::instance().nameTest + ".out";
203
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 std::ofstream file( nameFile.c_str() );
204
205
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 1 file<<FelisceParam::instance().nameTest <<std::endl;
206
207
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 file<<"Flow at inlet [cgs]"<<std::endl;
208
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for ( int k(0); k<m_nC; k++ ) {
209
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 file<<m_QInflow[k]<<'\t';
210 }
211
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 file<<std::endl;
212
213
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 file<<"Mean pressure at inlet [cgs]"<<std::endl;
214
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 for ( int k(0); k<m_nC; k++ ) {
215
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 1 file<<m_PInflow[k]<<'\t';
216 }
217
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 file<<std::endl;
218 1 }
219 4 }
220 /*! \brief post-process function, it writes a summary */
221 void
222 8 Autoregulation::exportControlFunctions(const int& numCycle) {
223
224 8 double dt = 0.1;
225 8 unsigned int N = static_cast<unsigned int>(numCycle * m_cycleDuration / dt);
226 unsigned int i;
227
228 8 std::string filename;
229
230
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 filename = FelisceParam::instance().resultDir + "extraP.txt" ;
231
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 std::ofstream extraP( filename.c_str() );
232 //Time
233
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i=0; i < N; i++) {
234
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 1280 extraP << dt*i<<'\t';
235 }
236
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 extraP<<std::endl;
237 //ExtraPressure interpolated
238
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i=0; i < N; i++) {
239
3/6
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1280 times.
✗ Branch 8 not taken.
 1280 extraP << interpPressure(dt*i)<<'\t';
240 }
241
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 extraP<<std::endl;
242
243
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 extraP.close();
244
245
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 filename = FelisceParam::instance().resultDir + "cont.txt" ;
246
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 std::ofstream cont(filename.c_str());
247 //Time
248
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i=0; i < N; i++) {
249
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 1280 cont << dt*i<<'\t';
250 }
251
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 cont<<std::endl;
252 //Force value
253
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i=0; i < N; i++) {
254
3/6
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1280 times.
✗ Branch 8 not taken.
 1280 cont << control(dt*i)<<'\t';
255 }
256
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 cont<<std::endl;
257
258
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 cont.close();
259
260
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
 8 filename = FelisceParam::instance().resultDir + "W.txt" ;
261
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 std::ofstream Windkessel(filename.c_str());
262 //Time
263
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i=0; i < N; i++) {
264
2/4
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1280 times.
✗ Branch 5 not taken.
 1280 Windkessel << dt*i<<'\t';
265 }
266
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Windkessel<<std::endl;
267 //Distal resistances
268
2/2
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 8 times.
 16 for ( std::size_t nbOfW(0); nbOfW < m_Rdist.size(); nbOfW++) {
269
2/2
✓ Branch 0 taken 1280 times.
✓ Branch 1 taken 8 times.
 1288 for ( i = 0; i<N; i++) {
270
1/2
✓ Branch 1 taken 1280 times.
✗ Branch 2 not taken.
 1280 setLumpedModel(dt*i);
271
2/4
✓ Branch 2 taken 1280 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1280 times.
✗ Branch 6 not taken.
 1280 Windkessel << m_Rdist[nbOfW]<<'\t';
272 }
273
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Windkessel<<std::endl;
274 }
275
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 Windkessel.close();
276 8 }
277 /*! \brief given a mean P value it computes parameters for the shape of the inlet pressure*/
278 void
279 169127 Autoregulation::meanPressure2PeakAndDiast(double const& meanP, double& peakP, double& diastP ) const {
280 //the mean pressure is divided between the diastolic and the sistolic minus diastolic
281 169127 double proportion = FelisceParam::instance().propCardiacCycle ;
282 169127 peakP = proportion*meanP*2*PI*m_cycleDuration;
283 169127 diastP = ( 1. - proportion )*meanP;
284 169127 }
285
286 /*! \brief ??*/
287 double
288 168 Autoregulation::computeMeanInflowPress( const int &Case, const double &t0) const {
289 //quadrature points on a time cycle
290 168 int N = 1000;
291
292 //quadrature step
293 168 double dt = m_cycleDuration/N;
294 //first quadrature point
295 168 double time = t0 + dt/2;
296 //result
297 168 double Pmean(0.0);
298
299
2/2
✓ Branch 0 taken 168000 times.
✓ Branch 1 taken 168 times.
 168168 for (int k(0); k<N; k++){
300 //integral
301
1/2
✓ Branch 1 taken 168000 times.
✗ Branch 2 not taken.
 168000 Pmean += dt*inletPressure(time, Case);
302 168000 time+=dt;
303 }
304
305 168 return Pmean/m_cycleDuration;
306 }
307
308 /*! \brief at the current time step it returns an interpolated value for extraP
309
310 \param[in] t current time step
311 \return the interpolated value
312
313 for each cardiac cycle the mean pressure is computed and its difference with respect
314 to the reference value is also computed.\n
315 A piecewise linear interpolation is then used to obtain the value at the current time step.\n
316 The value of m_extraPressure at the i-th cycle is given to last instant of the cardiac cycle.
317 */
318 double
319 185505 Autoregulation::interpPressure( const double &t) const {
320
321 //The value of the cycle is assigned to right node, the first node is arbitrarly std::set to zero
322 185505 int N = static_cast<int>(floor( t / m_cycleDuration));
323 185505 double tLoc = t-m_cycleDuration*N;
324
2/2
✓ Branch 0 taken 181857 times.
✓ Branch 1 taken 3648 times.
 185505 if (N==0) {
325 181857 return tLoc/m_cycleDuration*m_extraPressure[0];
326 } else {
327 3648 return m_extraPressure[N-1]+tLoc/m_cycleDuration*(m_extraPressure[N]-m_extraPressure[N-1]);
328 }
329 }
330 /*! \brief it computes the current value for the windkessel parameters
331
332 \param[in] t current time step
333
334 At a given time step t, extraP is computed using interpPressure.
335 Then, for each outlet, the distal resistance is updated by using
336 \f$
337 R_{dist}= q R_{dist,ref} + R_{dist,ref} ( 1 - q ) ( 1 + \mbox{regParam}\,\mbox{sigmoid}( \delta p ) )
338 \f$
339 */
340 1324 void Autoregulation::setLumpedModel(const double& t) {
341
1/2
✓ Branch 1 taken 1324 times.
✗ Branch 2 not taken.
 1324 double extraP = interpPressure( t );
342
1/2
✓ Branch 1 taken 1324 times.
✗ Branch 2 not taken.
 1324 double q = FelisceParam::instance().windkesselProp ;
343
2/2
✓ Branch 1 taken 2648 times.
✓ Branch 2 taken 1324 times.
 3972 for( std::size_t i(0);i<m_RdistRef.size();i++) {
344 2648 m_Rdist[i] = m_RdistRef[i]*q + m_RdistRef[i]* (1-q) * ( 1 + m_regParam[i] * sigmoid(extraP));
345 2648 m_C[i] = m_CRef[i]/( q + (1-q) * m_regParam[i] * sigmoid(extraP) );
346 }
347 1324 }
348 /*! \brief the MASTER proc displays several value
349 */
350 void
351 ✗ Autoregulation::print(const int & rank) {
352
353 ✗ if (rank==MASTER) {
354 ✗ if (FelisceParam::verbose() > 3 ) {
355 ✗ std::cout<<"Shape Parameters for the control functions: "<<std::endl;
356 ✗ std::cout<<" Maximum value: "<<FelisceParam::instance().maxT <<std::endl;
357 ✗ std::cout<<" Minimum value: "<<FelisceParam::instance().maxT *sigmoid(-1.e5)<<std::endl;
358 ✗ std::cout<<" Center point: "<<m_center<<std::endl;
359 ✗ std::cout<<" Slope: "<<m_slope<<std::endl;
360 ✗ std::cout<<std::endl;
361 ✗ std::cout<<std::endl;
362
363 ✗ std::cout<<"Computed extraPressure"<<std::endl;
364 ✗ for ( std::size_t i(0); i<m_extraPressure.size(); i++)
365 ✗ std::cout<<i+1<<'\t'<<m_extraPressure[i]<<std::endl;
366 ✗ std::cout<<std::endl;
367 ✗ std::cout<<std::endl;
368
369 ✗ Tools::printVector(m_labelInlet,"Inflow labels");
370 ✗ Tools::printVector(m_labelOutlet,"Outflow labels");
371 ✗ Tools::printVector(m_labelStructure,"Structural labels");
372
373 ✗ if ( FelisceParam::instance().useRegulationInWindkessel ) {
374 ✗ std::cout<<"Windkessel parameters"<<std::endl;
375 ✗ Tools::printVector(m_regParam, "Regulation Parameter");
376 ✗ Tools::printVector(m_RdistRef, "Reference Distal Resistance value");
377 ✗ Tools::printVector(m_CRef, "Reference capacitance value");
378 }
379 else
380 ✗ std::cout<<"Regulation in windkessel switched off"<<std::endl;
381 }
382 }
383 }
384
385 }
386