GCC Code Coverage Report

Directory: ./
File: Solver/RISModel.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 293 0.0%
Branches: 0 340 0.0%
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: S.Smaldone, A.This and L.Boilevin-Kayl
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Core/felisceTools.hpp"
21 #include "Solver/RISModel.hpp"
22
23 namespace felisce {
24
25 ✗ RISModel::RISModel(FelisceTransient::Pointer fstransient , std::vector<LinearProblem*> linearProblem):
26 ✗ m_fstransient(fstransient),
27 ✗ m_linearProblem(linearProblem) {
28
29 ✗ if (!FelisceParam::instance().FusionDof) {
30 ✗ FEL_ERROR("There are no fusioned labels whereas it is required for the RIS model!");
31 }
32
33 ✗ variables_initialization(); // initialisation of attributes for the general RIS model
34
35 ✗ for(int rs = 0; rs <m_RISModelNum; rs++)
36 ✗ RIS_models_initialization(rs); // initialisation of the resistance values of all open, closed and fake surfaces
37
38 ✗ print(); // print status and resistance values of each closed and open surfaces
39
40 }
41
42
43
44 ✗ void RISModel::variables_initialization() {
45
46 ✗ m_verboseRIS = FelisceParam::instance().verboseRIS; // verbose for the model
47
48 ✗ m_R_active_surf = FelisceParam::instance().R_active_surf; // resistance of closed surface
49 ✗ m_R_inactive_surf = FelisceParam::instance().R_inactive_surf; // resistance of open surface
50
51 ✗ m_R_activeClosed_TwoValves_Intermediate = FelisceParam::instance().R_activeClosed_TwoValves_Intermediate; // resistance of closed valve just after the closing - only used for two valves
52 ✗ m_R_activeClosed_TwoValves_Final = FelisceParam::instance().R_activeClosed_TwoValves_Final; // final resistance of closed valve after several iterations after the closing - only used for two valves
53
54 ✗ m_R_FirstValve_ClosingSpeed = FelisceParam::instance().R_FirstValve_ClosingSpeed; // rate of speed for the evolution of the resistance of the closed surface of the first valve from R_Intermediate to R_Final - only used for two valves
55 ✗ m_R_SecondValve_ClosingSpeed = FelisceParam::instance().R_SecondValve_ClosingSpeed; // rate of speed for the evolution of the resistance of the closed surface of the second valve from R_Intermediate to R_Final - only used for two valves
56
57 ✗ m_ModelAdmissibleStatus = 1; // status of the model initially admissible
58 ✗ m_cycl_RIS = 0;
59
60 ✗ m_nbTimeSteps_refractory_time = FelisceParam::instance().nbTimeSteps_refractory_time; // number of time steps used for defining the refractory time
61 ✗ m_refractory_time = m_nbTimeSteps_refractory_time * m_fstransient->timeStep;
62
63 ✗ m_linear_evolution_valves = false;
64 ✗ m_rate_linear_transition = (m_R_active_surf - m_R_inactive_surf)/ std::abs(FelisceParam::instance().nbTimeSteps_linear_transition * m_fstransient->timeStep); // rate of linear evolution of the valve
65 ✗ m_rate_linear_transition_firstValve = (m_R_active_surf - m_R_inactive_surf)/ std::abs(FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[0] * m_fstransient->timeStep); // rate of the linear evolution of the first valve
66 ✗ m_rate_linear_transition_secondValve = (m_R_active_surf - m_R_inactive_surf)/ std::abs(FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[1] * m_fstransient->timeStep); // rate of the linear evolution of the second valve
67
68 ✗ if (FelisceParam::instance().nbTimeSteps_linear_transition > 0 || FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[0] > 0 || FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[1] > 0 ){
69 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n **** Linear transition in the value of the resistances of open and closed valves **** \n");
70 ✗ m_linear_evolution_valves = true;
71 }
72
73 ✗ felInt NumTotLabel = 100; // should be equal to the maximum label among open, closed and fake surfaces to spare some space in the array m_resistances
74 ✗ m_resistances.resize(NumTotLabel, 0.);
75
76 ✗ m_num_closed_surfs = FelisceParam::instance().closed_surf.size();
77 ✗ m_num_open_surfs = FelisceParam::instance().open_surf.size();
78 ✗ m_num_fake_surfs = FelisceParam::instance().fake_surf.size();
79
80 ✗ m_RISModelNum = m_num_closed_surfs / 2; // the number of RIS model is equal to the number of coupled closed surfaces
81
82 ✗ m_closed_surf.resize(m_num_closed_surfs, -1);
83 ✗ m_open_surf.resize(m_num_open_surfs, -1);
84 ✗ m_fake_surf.resize(m_num_fake_surfs, -1);
85
86 ✗ m_ValveAdmissibleStatus.resize(m_RISModelNum, 1); // status of all valves initially open
87 ✗ m_valveClosedFlag.resize(m_RISModelNum, 0); // all valves are initially considered open
88 ✗ m_closingTime.resize(m_RISModelNum, 0.);
89 ✗ m_openingTime.resize(m_RISModelNum, 0.);
90
91 ✗ m_flow_ref.resize(m_RISModelNum, FelisceParam::instance().flow_ref_oneValve);
92
93 ✗ if(m_RISModelNum == 2) {
94 ✗ m_flow_ref[0] = FelisceParam::instance().flow_ref_twoValves[0];
95 ✗ m_flow_ref[1] = FelisceParam::instance().flow_ref_twoValves[1];
96 }
97
98 ✗ if (m_RISModelNum > 2) {
99 ✗ FEL_ERROR("RIS model cannot handle more than two valves at the moment!");
100 }
101
102 }
103
104 ✗ void RISModel::RIS_models_initialization(int rs) {
105
106 ✗ if(rs == FelisceParam::instance().initiallyOpenedValveIndex){
107 ✗ m_valveClosedFlag[rs] = 0; // the valve of index initiallyOpenedValveIndex is initially open, meaning its open labels are active and its closed labels are inactive
108 }
109 else{
110 ✗ m_valveClosedFlag[rs] = 1; // the others valve are initially closed, meaning their closed labels are active and their open labels are inactive
111 }
112
113 // ---------------- closed surfaces initialization ------------------
114
115 ✗ double R_value = m_R_inactive_surf * (1 - m_valveClosedFlag[rs]) + m_valveClosedFlag[rs] * m_R_active_surf; // either R_inactive if valve is initially open or R_active if valve initally closed
116
117 ✗ m_closed_surf[2*rs] = FelisceParam::instance().closed_surf[2*rs];
118 ✗ m_closed_surf[2*rs+1] = FelisceParam::instance().closed_surf[2*rs+1];
119
120 ✗ m_resistances[ m_closed_surf[2*rs] ] = R_value;
121 ✗ m_resistances[ m_closed_surf[2*rs+1] ] = R_value;
122
123 // ---------------- open surfaces initialization ------------------
124
125 ✗ for(int io = rs; io < m_num_open_surfs/2; io++) { // to consider the case where there are more open surfaces than closed surfaces
126 ✗ R_value = m_R_inactive_surf * (m_valveClosedFlag[rs]) + (1 - m_valveClosedFlag[rs]) * m_R_active_surf; // either R_inactive if valve is initially closed or R_active if valve initially open
127
128 ✗ m_open_surf[2*io] = FelisceParam::instance().open_surf[2*io];
129 ✗ m_open_surf[2*io+1] = FelisceParam::instance().open_surf[2*io+1];
130
131 ✗ m_resistances[ m_open_surf[2*io] ] = R_value;
132 ✗ m_resistances[ m_open_surf[2*io+1] ] = R_value;
133 }
134
135 // ---------------- fake surfaces initialization ------------------
136
137 ✗ for(int is = rs; is < m_num_fake_surfs/2; is++) { // to consider the case where there are more fake surfaces than closed surfaces
138 ✗ R_value = 0.; // fake surfaces do not have any influence, so resistance is std::set to 0
139
140 ✗ m_fake_surf[2*is] = FelisceParam::instance().fake_surf[2*is];
141 ✗ m_fake_surf[2*is+1] = FelisceParam::instance().fake_surf[2*is+1];
142
143 ✗ m_resistances[ m_fake_surf[2*is] ] = 0.;
144 ✗ m_resistances[ m_fake_surf[2*is+1] ] = 0.;
145 }
146
147 }
148
149 // ******************** Update resistances ********************************
150
151 ✗ void RISModel::UpdateResistances() {
152
153 ✗ for(int rs = 0; rs <m_RISModelNum; rs++) {
154 ✗ update_open_surf(rs); // std::set resistance of open surfaces to R_active if valve is open, to R_inactive if valve is closed
155 ✗ update_closed_surf(rs); // std::set resistance of closed surfaces to R_inactive if valve is open, to R_active if valve is closed
156 }
157
158 ✗ print();
159
160 }
161
162 ✗ void RISModel::update_open_surf(int rs) {
163
164 ✗ double R_value = 0.;
165
166 ✗ for(int io = rs; io < m_num_open_surfs/2; io++) {
167 ✗ R_value = m_R_inactive_surf * (m_valveClosedFlag[rs]) + (1 - m_valveClosedFlag[rs]) * m_R_active_surf; // either R_active if open or R_inactive if closed
168
169 ✗ if (m_linear_evolution_valves) { //Linear evolution for open valve/s
170 ✗ m_tps = m_fstransient->time; // update current time
171 ✗ double rate = 0.;
172
173 ✗ if (m_RISModelNum == 2) { // set decreasing rate for the first and the second valve
174 ✗ if (rs == 0)
175 ✗ rate = m_rate_linear_transition_firstValve;
176 ✗ if (rs == 1)
177 ✗ rate = m_rate_linear_transition_secondValve;
178 } else {
179 ✗ rate = m_rate_linear_transition; // set decreasing rate for the valve
180 }
181
182 ✗ if ( m_openingTime[rs] > 0. && (m_tps - m_openingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition * m_fstransient -> timeStep) )
183 ✗ R_value = rate * (m_tps - m_openingTime[rs]) + m_R_inactive_surf;
184 ✗ else if ( m_closingTime[rs] > 0. && (m_tps - m_closingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition * m_fstransient -> timeStep) )
185 ✗ R_value = - rate * (m_tps - m_closingTime[rs]) + m_R_active_surf;
186 }
187
188 ✗ m_resistances[ m_open_surf[2*io] ] = R_value;
189 ✗ m_resistances[ m_open_surf[2*io+1] ] = R_value;
190 }
191
192 }
193
194 ✗ void RISModel::update_closed_surf(int rs) {
195
196 ✗ std::vector<double> flux_surf; // values of mean flux flow across each couple of closed surfaces (used when they are inactive)
197 ✗ std::vector<double> absolute_flux_surf;
198 ✗ std::vector<double> mesh_flux_surf;
199
200 ✗ std::vector<double> pressure_surf; // values of mean pressure on each surface of closed surfaces (used when they are active)
201
202 ✗ m_tps = m_fstransient->time;
203 ✗ m_ValveAdmissibleStatus[rs] = 1;
204 ✗ m_ModelAdmissibleStatus = 1;
205
206 ✗ m_closedLabels.resize(2,0); // consider only a couple of closed surfaces
207 ✗ m_closedLabels[0] = m_closed_surf[2*rs];
208 ✗ m_closedLabels[1] = m_closed_surf[2*rs + 1];
209
210 ✗ pressure_surf.resize(2, 0.); // size 2 for distal and proximal
211 ✗ flux_surf.resize(2, 0.); // size 2 for distal and proximal
212 ✗ absolute_flux_surf.resize(2,0.);
213 ✗ mesh_flux_surf.resize(2,0.);
214
215
216 ✗ if (m_RISModelNum == 2) {
217 ✗ two_valve_case(rs);
218 }
219
220 ✗ m_linearProblem[0]->computeMeanQuantity(pressure, m_closedLabels, pressure_surf);
221 ✗ m_linearProblem[0]->computeMeanQuantity(velocity, m_closedLabels, absolute_flux_surf);
222 ✗ if( m_fstransient->iteration > 1 && FelisceParam::instance().useALEformulation ){
223 ✗ m_linearProblem[0]->computeMeanExternalQuantity(m_linearProblem[0]->externalVec(0), m_linearProblem[0]->externalAO(0), velocity, m_closedLabels, mesh_flux_surf);
224 }
225 ✗ for(int i = 0; i < 2; i++){
226 ✗ flux_surf[i] = absolute_flux_surf[i] - (-mesh_flux_surf[i]);
227 }
228
229
230
231
232 ✗ if(m_verboseRIS>0) {
233 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n######## Valve %d (P_0 = %e on label %d, P_1 = %e on label %d)\n",
234 ✗ rs+1, pressure_surf[0], m_closedLabels[0], pressure_surf[1], m_closedLabels[1]);
235
236 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n######## Valve %d (flow = %e (std::abs %e, mesh %e) on label %d)\n",
237 ✗ rs+1, flux_surf[0],absolute_flux_surf[0],mesh_flux_surf[0], m_closedLabels[0]);
238 }
239
240 ✗ if(m_valveClosedFlag[rs] == 1) { // if valve rs is closed check if it should open
241 ✗ check_opening_condition(rs, pressure_surf); // on the pressure difference
242 }
243 ✗ if(m_valveClosedFlag[rs] == 0) { // else if valve rs is open check if it should close
244 ✗ check_closing_condition(rs, flux_surf); // on the existence of a backflow
245 }
246
247 }
248
249 // ************************************************************************
250
251 ✗ void RISModel::two_valve_case(int rs) {
252
253 ✗ double Rmitral = 0.;
254 ✗ double Raortic = 0.;
255
256 ✗ felInt displ_IterMax_RIS = FelisceParam::instance().timeMaxCaseFile; // to do: it should be taken automatically from the displacement in file.case
257 //Commenting displ_TimeMax as we no longer need it following generalization of RIS model to general opening and closing time
258 //double displ_TimeMax_RIS = displ_IterMax_RIS * m_fstransient-->timeStep;
259
260 ✗ if(m_fstransient->iteration == (m_cycl_RIS+1) * displ_IterMax_RIS) { // to read the file case in cycle
261 ✗ m_cycl_RIS++;
262 }
263
264 ✗ if(rs == 0 && m_verboseRIS > 1 ) { // print just once
265 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n --> m_cycl_RIS = %d ", m_cycl_RIS+1);
266 }
267
268 // 21 June ( Alexandre This ): Modification of the law to be stiffer for the valve closure (40 --> 100)
269 ✗ if(rs == 0 && m_valveClosedFlag[0] == 1 ) {
270
271 ✗ Rmitral = m_R_activeClosed_TwoValves_Intermediate + m_R_activeClosed_TwoValves_Final * exp( -exp( -m_R_FirstValve_ClosingSpeed * ( m_tps - m_closingTime[rs]))); // Exponential evolution of the resistance of the proximal valve
272
273 // Linear evolution value of the resistance of the first valve
274 ✗ if (m_linear_evolution_valves) {
275 ✗ if ( m_closingTime[rs] > 0. && (m_tps - m_closingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[0] * m_fstransient -> timeStep) )
276 ✗ Rmitral = m_rate_linear_transition_firstValve * (m_tps - m_closingTime[rs]) + m_R_inactive_surf; // The proximal valve is closing, linear evolution from m_R_inactive_surf to m_R_active_surf
277 else
278 ✗ Rmitral = m_R_active_surf; // The first valve is completely closed
279 }
280
281 ✗ m_resistances[ m_closedLabels[0] ] = Rmitral;
282 ✗ m_resistances[ m_closedLabels[1] ] = Rmitral;
283 }
284
285 // 21 June ( Alexandre This ): Modification of the law to be stiffer for the valve closure (30 --> 75)
286 ✗ if(rs == 1 && m_valveClosedFlag[1] == 1){// and m_tps > (0.4 + m_cycl_RIS * displ_TimeMax_RIS) ) {
287 // the 3th option is necessary, because in the beginning of a cardiac cycle
288 // the aortic Resitance should be laready greater and not increasing,
289 ✗ Raortic = m_R_activeClosed_TwoValves_Intermediate + m_R_activeClosed_TwoValves_Final * exp( -exp( -m_R_SecondValve_ClosingSpeed * ( m_tps - m_closingTime[rs] ))); // Exponential evolution of the resistance of the second valve
290
291 // Linear evolution in the value of the resistance of the second valve
292 ✗ if (m_linear_evolution_valves) {
293 ✗ if ( m_closingTime[rs] > 0. && (m_tps - m_closingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition_twoValves[1] * m_fstransient -> timeStep) )
294 ✗ Raortic = m_rate_linear_transition_secondValve * (m_tps - m_closingTime[rs]) + m_R_inactive_surf; // The second valve is closing, linear evolution from m_R_inactive_surf to m_R_active_surf
295 else
296 ✗ Raortic = m_R_active_surf; // The second valve is completely closed
297 }
298
299 ✗ m_resistances[ m_closedLabels[0] ] = Raortic;
300 ✗ m_resistances[ m_closedLabels[1] ] = Raortic;
301 }
302
303
304 ✗ if(m_valveClosedFlag.size() > 1) {
305 ✗ if (m_valveClosedFlag[0] == 0 && m_valveClosedFlag[1] == 0)
306 ✗ FEL_ERROR("Valves are open at the same time!");
307 }
308
309 }
310
311 ✗ void RISModel::check_opening_condition(int rs, std::vector<double> pressure_surf) {
312
313 ✗ if( ((pressure_surf[1] - pressure_surf[0]) <= 0) && ((m_tps - m_closingTime[rs]) > m_refractory_time) && (!(preventValvesFromOpeningTooFast(rs)))) { // surf_1 is the distal face, surf_0 the proximal face
314 // negative difference pressure --> the valve must open
315 ✗ double R_value = 0.;
316
317 ✗ if(m_linear_evolution_valves)
318 ✗ R_value = m_R_active_surf;
319 else
320 ✗ R_value = m_R_inactive_surf;
321
322 ✗ m_resistances[ m_closedLabels[0] ] = R_value;
323 ✗ m_resistances[ m_closedLabels[1] ] = R_value;
324
325 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n -------> At time t = %e Valve %d OPENS due to a positive pressure difference (P_0 = %e, P_1 = %e)!\n",
326 ✗ m_tps,rs+1, pressure_surf[0], pressure_surf[1]);
327
328 ✗ m_valveClosedFlag[rs] = 0;
329 ✗ m_openingTime[rs] = m_tps; // update last time of opening
330
331 }
332 else {
333 ✗ if (m_linear_evolution_valves && m_closingTime[rs] > 0. && m_RISModelNum == 1) { // Linear increasing in the value of the resistance of the closed valve/s
334 ✗ double R_value = 0.;
335
336 ✗ if ( (m_tps - m_closingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition * m_fstransient -> timeStep) )
337 ✗ R_value = m_rate_linear_transition * (m_tps - m_closingTime[rs]) + m_R_inactive_surf;
338 else
339 ✗ R_value = m_R_active_surf;
340
341 ✗ m_resistances[ m_closedLabels[0] ] = R_value;
342 ✗ m_resistances[ m_closedLabels[1] ] = R_value;
343 }
344 }
345
346 ✗ if(preventValvesFromOpeningTooFast(rs)) { // print text if a valve has been prevented to open
347
348 ✗ if(m_verboseRIS > 1) {
349 ✗ PetscPrintf(PETSC_COMM_WORLD,"Closed valve of m_ris = %d would have normally opened but has been preventing from!\n", rs);
350 }
351
352 }
353
354 }
355
356 ✗ void RISModel::check_closing_condition(int rs, std::vector<double> flux_surf) {
357
358 ✗ if( flux_surf[0] <= m_flow_ref[rs] && (m_tps - m_openingTime[rs]) > m_refractory_time) { // after "m_refractory_time" seconds the valve is open, it cannot close
359
360 // existence of a backflow --> the valve must close
361 ✗ if (m_RISModelNum == 2 ) {
362 ✗ if(rs == 0 && !m_linear_evolution_valves)
363 ✗ m_R_active_surf = m_R_activeClosed_TwoValves_Intermediate; // In the case of two valves, when the valve closes, its resistance is std::set to 1.e4 and then slightly evolves towards m_R_active_surf in two_valves_cases
364
365 ✗ if(rs == 1 && !m_linear_evolution_valves)
366 ✗ m_R_active_surf = m_R_activeClosed_TwoValves_Intermediate;
367 }
368
369 double R_value;
370
371 ✗ if(m_linear_evolution_valves)
372 ✗ R_value = m_R_inactive_surf; // In the case of linear evolution, when the valve closes, its resistance evolves from m_R_inactive_surf to m_R_active_surf
373 else
374 ✗ R_value = m_R_active_surf; // In the case outside of two valves, when the valves closes, its resistance is directly std::set to m_R_active_surf
375
376 ✗ m_resistances[ m_closedLabels[0] ] = R_value;
377 ✗ m_resistances[ m_closedLabels[1] ] = R_value;
378
379 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n -------> At time t = %e valve %d CLOSES due to a backflow (flow = %e )!\n",
380 ✗ m_tps,rs+1, flux_surf[0] );
381
382 ✗ m_valveClosedFlag[rs] = 1;
383 ✗ m_closingTime[rs] = m_tps; // update last time of closing
384
385 }
386 else {
387 ✗ if (m_linear_evolution_valves && m_openingTime[rs] > 0.) { // Linear decreasing in the value of the resistance of the closed valve/s
388 ✗ double R_value = 0.;
389 ✗ double rate = 0.;
390
391 ✗ if (m_RISModelNum == 2) { // set decreasing rate for the first and the second valve
392 ✗ if (rs == 0)
393 ✗ rate = m_rate_linear_transition_firstValve;
394 ✗ if (rs == 1)
395 ✗ rate = m_rate_linear_transition_secondValve;
396 } else {
397 ✗ rate = m_rate_linear_transition; // set decreasing rate for the valve
398 }
399
400 ✗ if ( (m_tps - m_openingTime[rs]) < std::abs(FelisceParam::instance().nbTimeSteps_linear_transition * m_fstransient -> timeStep) )
401 ✗ R_value = - rate * (m_tps - m_openingTime[rs]) + m_R_active_surf;
402 else
403 ✗ R_value = m_R_inactive_surf;
404
405 ✗ m_resistances[ m_closedLabels[0] ] = R_value;
406 ✗ m_resistances[ m_closedLabels[1] ] = R_value;
407 }
408 }
409 }
410
411 ✗ bool RISModel::preventValvesFromOpeningTooFast(int rs) { // returns 1 if one valve can open rapidly after the other closed, returns 0 if not.
412
413 ✗ if(m_RISModelNum != 2) { // always return 0 if there are not 2 valves
414 ✗ return false;
415 }
416
417 else {
418
419 ✗ if(FelisceParam::instance().preventValvesFromOpeningTooFast == 0) { // prevention not activated
420 ✗ return false;
421 }
422
423 else {
424
425 ✗ if(m_fstransient->iteration > FelisceParam::instance().nbIterationsPreventValvesFromOpeningTooFast) {
426
427 ✗ if(rs == 0) { // if the opening condition is related to the first valve
428
429 ✗ if(std::max(0., m_fstransient->time - m_closingTime[1]) < FelisceParam::instance().nbIterationsPreventValvesFromOpeningTooFast*FelisceParam::instance().timeStep) { // if the other valve (rs = 1) has recently closed, then return true
430 ✗ std::cout << "!/CASE 1: mitral valve prevented from opening!/" << std::endl;
431 ✗ return true;
432 }
433
434 else { // if not returns false
435 ✗ std::cout << "!/CASE 2: mitral valve not prevented from opening it its conditions are met!/" << std::endl;
436 ✗ return false;
437 }
438 }
439
440 else { // if the opening condition is related to the second valve
441
442 ✗ if(std::max(0., m_fstransient->time - m_closingTime[0]) < FelisceParam::instance().nbIterationsPreventValvesFromOpeningTooFast*FelisceParam::instance().timeStep) { // if the other valve (rs = 0) has recently closed, then return true
443 ✗ std::cout << "!/CASE 3: aortic valve prevented from opening!/" << std::endl;
444 ✗ return true;
445 }
446
447 else { // if not returns false
448 ✗ std::cout << "!/CASE 4: aortic valve not prevented from opening if its conditions are met!/" << std::endl;
449 ✗ return false;
450 }
451 }
452 }
453
454 else {
455 ✗ return false;
456 }
457 }
458 }
459
460 return false;
461
462 }
463
464
465 // ******************* Check status admissibility *****************************
466
467 ✗ void RISModel::CheckValveStatus() {
468
469 ✗ for(int rs = 0; rs <m_RISModelNum; rs++) {
470 ✗ m_closedLabels.resize(2,0);
471 ✗ m_closedLabels[0] = m_closed_surf[2*rs];
472 ✗ m_closedLabels[1] = m_closed_surf[2*rs + 1];
473
474 ✗ if( m_valveClosedFlag[rs] == 1) { // the valve is closed
475 ✗ PetscPrintf(PETSC_COMM_WORLD,"******** VALVE MODEL %d is CLOSED\n",rs+1);
476 ✗ check_closed_status(rs);
477 }
478 ✗ if( m_valveClosedFlag[rs] == 0) { // the valve is open
479 ✗ PetscPrintf(PETSC_COMM_WORLD,"******** VALVE MODEL %d is OPEN\n", rs+1);
480 ✗ check_open_status(rs);
481 }
482
483 ✗ m_ModelAdmissibleStatus *= m_ValveAdmissibleStatus[rs]; // admissibility of the model = product of valve admissible status (0 if not correct, 1 if correct)
484
485 }
486
487 }
488
489 ✗ void RISModel::check_closed_status(int rs) {
490
491 ✗ std::vector<double> pressure_surf;
492 ✗ pressure_surf.resize(2, 0.);
493
494 ✗ m_linearProblem[0]->computeMeanQuantity(pressure, m_closedLabels, pressure_surf);
495 ✗ FEL_ASSERT(m_closedLabels.size() > 1);
496
497 ✗ if(m_verboseRIS > 1) {
498 ✗ PetscPrintf(PETSC_COMM_WORLD," Check fissure status: P_0 = %e on label %d, P_1 = %e on label %d\n",
499 ✗ pressure_surf[0], m_closedLabels[0], pressure_surf[1], m_closedLabels[1]);
500 }
501 ✗ if( m_resistances[ m_closedLabels[0] ] > 1 && ( pressure_surf[1] - pressure_surf[0] ) <= 0 and m_tps > 5 * m_refractory_time) {
502
503 ✗ if(m_verboseRIS > 1)
504 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n Status NOT ADMISSIBLE for valve %d\n", rs+1);
505 // positive resistance --> valve is closed
506 // negative pressure difference --> valve should be open
507 // the status is then not admissible
508 ✗ if(FelisceParam::instance().crashIfStatusNotAdmissible) {
509 ✗ m_ValveAdmissibleStatus[rs] = 0; // fissure not ok
510 }
511 else {
512 ✗ m_ValveAdmissibleStatus[rs] = 1; // fissure still ok
513 ✗ PetscPrintf(PETSC_COMM_WORLD," Keep running the simulation even if status not admissible\n");
514 }
515 } else { //the valve is open and fissure OK
516 ✗ if(m_verboseRIS > 1)
517 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n Status ADMISSIBLE for valve %d\n", rs+1);
518 ✗ m_ValveAdmissibleStatus[rs] = 1;
519 }
520
521 }
522
523 ✗ void RISModel::check_open_status(int rs) {
524
525 ✗ std::vector<double> flux_surf;
526 ✗ flux_surf.resize(2, 0.);
527
528 ✗ m_linearProblem[0]->computeMeanQuantity(velocity, m_closedLabels, flux_surf);
529
530 ✗ if(m_verboseRIS > 1) {
531 ✗ PetscPrintf(PETSC_COMM_WORLD," Check fissure status: F_0 = %e on label %d\n",
532 ✗ flux_surf[0], m_closedLabels[0]);
533 }
534
535 ✗ if( (Tools::equal(m_resistances[ m_closedLabels[0] ],0.)) && (flux_surf[0] <= m_flow_ref[0]) ) {
536 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n Status NOT ADMISSIBLE for Valve%d\n", rs+1);
537 // zero resistance --> valve is open
538 // negative flow --> valve should be closed
539 // the status is then not admissible
540 ✗ if(FelisceParam::instance().crashIfStatusNotAdmissible) {
541 ✗ m_ValveAdmissibleStatus[rs] = 0; // fissure not ok
542 }
543 else{
544 ✗ m_ValveAdmissibleStatus[rs] = 1; // fissure still ok
545 ✗ PetscPrintf(PETSC_COMM_WORLD," Keep running the simulation even if status not admissible\n");
546 }
547
548 }
549
550 else { // the valve is open and fissure OK
551
552 ✗ if(m_verboseRIS > 1) {
553 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n Status ADMISSIBLE for valve %d\n", rs+1);
554 }
555
556 ✗ m_ValveAdmissibleStatus[rs] = 1;
557 }
558
559 }
560
561 /* void RISModel::check_open_status(int rs){
562 if(m_verboseRIS>1)
563 PetscPrintf(PETSC_COMM_WORLD," Status Admissible for Valve %d\n", rs+1);
564 m_ValveAdmissibleStatus[rs] = 1;
565 }*/
566
567 ✗ int RISModel::areAllValvesClosed() const {
568
569 ✗ int areAllValvesClosed = 1;
570
571 ✗ for(int i = 0; i < m_RISModelNum; ++i) {
572
573 ✗ if(m_valveClosedFlag[i] == 0) { // if current evaluated valve is open, std::set areAllValvesClosed to 0
574 ✗ areAllValvesClosed = 0;
575 ✗ break;
576 }
577
578 else{ // if current evaluated valve is closed, do nothing
579
580 }
581
582 }
583
584 ✗ return areAllValvesClosed;
585
586 }
587
588 // *********************** Print resistances ***********************************
589
590 ✗ void RISModel::print() const {
591
592 ✗ for(int ic = 0; ic < m_num_closed_surfs/2; ic++) {
593 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n RIS model number = %d --> labels closed surfaces = %d and %d \n",ic+1, m_closed_surf[2*ic], m_closed_surf[2*ic+1]);
594 ✗ PetscPrintf(PETSC_COMM_WORLD," with resistance = %e", m_resistances[ m_closed_surf[2*ic] ]);;
595 }
596 ✗ for(int io = 0; io < m_num_open_surfs/2; io++) {
597 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n RIS model number = %d --> labels open surfaces = %d and %d \n",io+1, m_open_surf[2*io], m_open_surf[2*io+1]);
598 ✗ PetscPrintf(PETSC_COMM_WORLD," with resistance = %e \n",m_resistances[ m_open_surf[2*io] ]);
599 }
600 ✗ for(int is = 0; is < m_num_fake_surfs/2; is++) {
601 ✗ PetscPrintf(PETSC_COMM_WORLD,"\n RIS model number = %d --> labels open surfaces = %d and %d \n",is+1, m_fake_surf[2*is], m_fake_surf[2*is+1]);
602 }
603
604 }
605
606 // ****************************************************************************
607
608 ✗ void RISModel::write(std::string fileName) const{
609
610 static int RISInit = 1;
611 ✗ std::ofstream outputFile;
612 ✗ if(RISInit) {
613 ✗ outputFile.open((fileName.c_str()));
614 } else {
615 ✗ outputFile.open((fileName.c_str()),std::ios::app);
616 }
617
618 ✗ for(int ic = 0; ic < m_num_closed_surfs/2; ic++) {
619 ✗ outputFile << m_resistances[m_closed_surf[2*ic]] << " ;";
620 }
621 ✗ for(int io = 0; io < m_num_open_surfs/2; io++) {
622 ✗ outputFile << m_resistances[ m_open_surf[2*io] ] << " ;" << std::endl;
623 }
624 ✗ RISInit = 0;
625 }
626
627 }
628