GCC Code Coverage Report

Directory: ./
File: Solver/cardiacCycle.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 61 0.0%
Branches: 0 80 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
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/cardiacCycle.hpp"
21
22 namespace felisce
23 {
24 ✗ CardiacCycle::CardiacCycle(FelisceTransient::Pointer fstransient,
25 ✗ RISModel* risModel, std::vector<LinearProblem*> linearProblem):
26 ✗ m_fstransient(fstransient),
27 ✗ m_ris(risModel),
28 ✗ m_linearProblem(linearProblem) {
29
30 ✗ if (! FelisceParam::instance().RISModels) {
31 ✗ FEL_ERROR("Valve model not initialized");
32 }
33
34 ✗ m_inflow.resize(3,0.);
35 ✗ _CardiacCycleLabel.resize(1,0);
36 ✗ m_closed_surf.resize(FelisceParam::instance().closed_surf.size(),0);
37 ✗ m_p0_isovol_relax.resize(1, 0.);
38 ✗ m_flag_valve_was_open = false;
39
40 ✗ m_dt = m_fstransient->timeStep ;
41 ✗ _CardiacCycleLabel[0] = FelisceParam::instance().CardiacCycleLabel;
42 ✗ m_duration_ejection = FelisceParam::instance(). duration_ejection;
43 ✗ m_duration_isovol_relax = FelisceParam::instance().duration_isovol_relax;
44 ✗ m_duration_filling = FelisceParam::instance().duration_filling;
45 ✗ m_duration_isovol_contract = FelisceParam::instance().duration_isovol_contract;
46 ✗ m_t0_ejection = 0.;
47
48 ✗ FEL_ASSERT(FelisceParam::instance().closed_surf.size() > 1);
49 ✗ for(std::size_t ic = 0; ic < FelisceParam::instance().closed_surf.size()/2; ic++) {
50 ✗ m_closed_surf[2*ic] = FelisceParam::instance().closed_surf[2*ic];
51 ✗ m_closed_surf[2*ic+1] = FelisceParam::instance().closed_surf[2*ic+1];
52 }
53
54 }
55
56
57
58
59 ✗ void CardiacCycle::ApplyCardiacCycleInput() {
60 ✗ InflowCondition();
61 ✗ PressureInput();
62 }
63
64
65 ✗ void CardiacCycle::InflowCondition() {
66
67 ✗ double tps = 0.;
68 ✗ tps = m_fstransient->time;
69
70 ✗ if(m_ris->valve_is_closed() == 1) { // Valve is closed: Neumann BC.
71 // the last time instant when the valve was closed:
72 ✗ m_t0_ejection = tps;
73 }
74
75 ✗ if(m_ris->valve_is_closed() == 0 ) { // Valve is open
76
77 ✗ PetscPrintf(PETSC_COMM_WORLD,"####### SYSTOLE\n");
78 ✗ PetscPrintf(PETSC_COMM_WORLD, "******* Cardiac phase: Ejection\n\n" );
79
80 // the inlet is supposed to be oriented along y-axis : to change ! when some one will have more time than me. Sv
81 ✗ m_inflow[1] = 110* std::sin( (tps - m_t0_ejection) * M_PI / m_duration_ejection ) ;
82 ✗ m_penalisationParam = 10e30;
83 }
84 }
85
86
87 ✗ void CardiacCycle::PressureInput() {
88
89 ✗ std::vector<double> pres_distal_valve(m_closed_surf.size(), 0.);
90 ✗ double tphase = 0.;
91 ✗ m_pressure_in.resize(1, 0.);
92
93 ✗ if(m_ris->valve_is_closed() == 0) { // Valve is open: Dirichlet BC.
94 ✗ m_flag_valve_was_open = true; // to remember that the valve was open
95 ✗ m_t0_isovol_relax = m_fstransient->time; //t0_isovol_relax will remember the last time instant when the valve was open
96 }
97
98 ✗ if (m_ris->valve_is_closed() == 1 ) { // Valve is closed
99 ✗ m_penalisationParam = 0.;
100 ✗ tphase = m_fstransient->time - m_t0_isovol_relax;
101
102 ✗ if(m_flag_valve_was_open) {
103 // we enter here only the first time step after the valve closes
104 // the inlet mean pressure at the first time instant the valve close
105 ✗ m_linearProblem[0]->computeMeanQuantity(pressure, _CardiacCycleLabel, m_p0_isovol_relax);
106 ✗ m_flag_valve_was_open = false;
107 }
108
109 ✗ if(tphase < m_duration_isovol_relax +m_dt ) {
110 ✗ PetscPrintf(PETSC_COMM_WORLD, "####### DIASTOLE \n");
111 ✗ PetscPrintf(PETSC_COMM_WORLD, "******* Cardiac phase: Isovolumetric relaxation\n\n");
112 ✗ m_pressure_in[0] = m_p0_isovol_relax[0] * exp( -tphase / (m_duration_isovol_relax/4.) );
113 }
114
115 ✗ if (tphase >= (m_duration_isovol_relax + m_dt) && tphase < (m_duration_isovol_relax + m_duration_filling +m_dt)) {
116 ✗ PetscPrintf(PETSC_COMM_WORLD, "####### DIASTOLE \n");
117 ✗ PetscPrintf(PETSC_COMM_WORLD, "******* Cardiac phase: Filling\n\n");
118 }
119
120 ✗ if(tphase > (m_duration_isovol_relax + m_duration_filling)) {
121 ✗ PetscPrintf(PETSC_COMM_WORLD, "####### SYSTOLE \n");
122 ✗ PetscPrintf(PETSC_COMM_WORLD, "******* Cardiac phase: Isovolumetric contraction\n\n");
123
124 ✗ m_linearProblem[0]->computeMeanQuantity(pressure, m_closed_surf, pres_distal_valve);
125 ✗ m_pressure_in[0] = 1.2 * pres_distal_valve[1] * std::pow( tphase - (m_duration_isovol_relax+m_duration_filling), 6) / std::pow(m_duration_isovol_contract, 6) + pres_distal_valve[0] ;
126 }
127
128 // PetscPrintf(PETSC_COMM_WORLD,"******* Inlet BC : pressure\n");
129 // PetscPrintf(PETSC_COMM_WORLD, " p_in = %e\n", m_pressure_in[0]);
130
131 }
132
133 }
134
135
136 }
137
138