GCC Code Coverage Report

Directory: ./
File: Solver/MVSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 907 0.0%
Branches: 0 672 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: M. Boulakia
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/MVSolver.hpp"
21
22 namespace felisce {
23 ✗ MVSolver::MVSolver(FelisceTransient::Pointer fstransient):
24 ✗ IonicSolver(fstransient) {
25 ✗ m_nComp = 3;
26 ✗ m_assimilation = false;
27 ✗ m_cellTypeCut = 4; // E.Tixier: Test MV model for atria. This is ugly bu necessary (for now) to be able to pass the Buildbot tests.
28
29 ✗ if (FelisceParam::instance().CellsType != "hetero") {
30 ✗ m_uo.resize(1);
31 ✗ m_uu.resize(1);
32 ✗ m_thetav.resize(1);
33 ✗ m_thetaw.resize(1);
34 ✗ m_thetavm.resize(1);
35 ✗ m_thetao.resize(1);
36 ✗ m_tauvm1.resize(1);
37 ✗ m_tauvm2.resize(1);
38 ✗ m_tauvp.resize(1);
39 ✗ m_tauwm1.resize(1);
40 ✗ m_tauwm2.resize(1);
41 ✗ m_kwm.resize(1);
42 ✗ m_uwm.resize(1);
43 ✗ m_tauwm.resize(1);
44 ✗ m_taufi.resize(1);
45 ✗ m_tauo1.resize(1);
46 ✗ m_tauo2.resize(1);
47 ✗ m_tauso1.resize(1);
48 ✗ m_tauso2.resize(1);
49 ✗ m_kso.resize(1);
50 ✗ m_uso.resize(1);
51 ✗ m_taus1.resize(1);
52 ✗ m_taus2.resize(1);
53 ✗ m_ks.resize(1);
54 ✗ m_us.resize(1);
55 ✗ m_tausi.resize(1);
56 ✗ m_tauwinf.resize(1);
57 ✗ m_winfstar.resize(1);
58 ✗ m_gfi.resize(1);
59 ✗ m_gso.resize(1);
60 ✗ m_gsi.resize(1);
61 ✗ m_gfi[0] = 1.;m_gso[0] = 1.;m_gsi[0] = 1.;
62 ✗ if (FelisceParam::instance().CellsType == "endo") {
63 ✗ m_uo[0] = 0.; // u_o
64 ✗ m_uu[0] = 1.56; // u_u
65 ✗ m_thetav[0] = 0.3; // theta_v
66 ✗ m_thetaw[0] = 0.13; // theta_w
67 ✗ m_thetavm[0] = 0.2; // theta_v^-
68 ✗ m_thetao[0] = 0.006; // theta_o
69 ✗ m_tauvm1[0] = 75.; // tau_v1^-
70 ✗ m_tauvm2[0] = 10.; // tau_v2^-
71 ✗ m_tauvp[0] = 1.4506; // tau_v^+
72 ✗ m_tauwm1[0] = 6.; // tau_w1^-
73 ✗ m_tauwm2[0] = 140.; // tau_w2^+
74 ✗ m_kwm[0] = 200.; // k_w^-
75 ✗ m_uwm[0] = 0.016; // u_w^-
76 ✗ m_tauwm[0] = 280.; // tau_w^+
77 ✗ m_taufi[0] = 0.1; // tau_fi
78 ✗ m_tauo1[0] = 470.; // tau_o1
79 ✗ m_tauo2[0] = 6.; // tau_o2
80 ✗ m_tauso1[0] = 40.; // tau_so1
81 ✗ m_tauso2[0] = 1.2; // tau_so2
82 ✗ m_kso[0] = 2.; // k_so
83 ✗ m_uso[0] = 0.65; // u_so
84 ✗ m_taus1[0] = 2.7342; // tau_s1
85 ✗ m_taus2[0] = 2.; // tau_s2
86 ✗ m_ks[0] = 2.0994; // k_s
87 ✗ m_us[0] = 0.9087; // u_s
88 ✗ m_tausi[0] = 2.9013; // tau_si
89 ✗ m_tauwinf[0] = 0.0273;// tau_{w inf}
90 ✗ m_winfstar[0] = 0.78; // w^*_{inf}
91 }
92 ✗ else if (FelisceParam::instance().CellsType == "M") {
93 ✗ m_uo[0] = 0.; // u_o
94 ✗ m_uu[0] = 1.61; // u_u
95 ✗ m_thetav[0] = 0.3; // theta_v
96 ✗ m_thetaw[0] = 0.13; // theta_w
97 ✗ m_thetavm[0] = 0.1; // theta_v^-
98 ✗ m_thetao[0] = 0.005; // theta_o
99 ✗ m_tauvm1[0] = 80.; // tau_v1^-
100 ✗ m_tauvm2[0] = 1.4506; // tau_v2^-
101 ✗ m_tauvp[0] = 1.4506; // tau_v^+
102 ✗ m_tauwm1[0] = 70.; // tau_w1^-
103 ✗ m_tauwm2[0] = 8.; // tau_w2^+
104 ✗ m_kwm[0] = 200.; // k_w^-
105 ✗ m_uwm[0] = 0.016; // u_w^-
106 ✗ m_tauwm[0] = 280.; // tau_w^+
107 ✗ m_taufi[0] = 0.078; // tau_fi
108 ✗ m_tauo1[0] = 410.; // tau_o1
109 ✗ m_tauo2[0] = 7.; // tau_o2
110 ✗ m_tauso1[0] = 91.; // tau_so1
111 ✗ m_tauso2[0] = 0.8; // tau_so2
112 ✗ m_kso[0] = 2.1; // k_so
113 ✗ m_uso[0] = 0.6; // u_so
114 ✗ m_taus1[0] = 2.7342; // tau_s1
115 ✗ m_taus2[0] = 4.; // tau_s2
116 ✗ m_ks[0] = 2.0994; // k_s
117 ✗ m_us[0] = 0.9087; // u_s
118 ✗ m_tausi[0] = 3.3849; // tau_si
119 ✗ m_tauwinf[0] = 0.01; // tau_{w inf}
120 ✗ m_winfstar[0] = 0.5; // w^*_{inf}
121 }
122 ✗ else if (FelisceParam::instance().CellsType == "epi") {
123 ✗ m_uo[0] = 0.; // u_o
124 ✗ m_uu[0] = 1.55; // u_u
125 ✗ m_thetav[0] = 0.3; // theta_v
126 ✗ m_thetaw[0] = 0.13; // theta_w
127 ✗ m_thetavm[0] = 0.006; // theta_v^-
128 ✗ m_thetao[0] = 0.006; // theta_o
129 ✗ m_tauvm1[0] = 60.; // tau_v1^-
130 ✗ m_tauvm2[0] = 1150.; // tau_v2^-
131 ✗ m_tauvp[0] = 1.4506; // tau_v^+
132 ✗ m_tauwm1[0] = 60.; // tau_w1^-
133 ✗ m_tauwm2[0] = 15.; // tau_w2^+
134 ✗ m_kwm[0] = 65.; // k_w^-
135 ✗ m_uwm[0] = 0.03; // u_w^-
136 ✗ m_tauwm[0] = 200.; // tau_w^+
137 ✗ m_taufi[0] = 0.11; // tau_fi
138 ✗ m_tauo1[0] = 400.; // tau_o1
139 ✗ m_tauo2[0] = 6.; // tau_o2
140 ✗ m_tauso1[0] = 30.0181;// tau_so1
141 ✗ m_tauso2[0] = 0.9957; // tau_so2
142 ✗ m_kso[0] = 2.0458; // k_so
143 ✗ m_uso[0] = 0.65; // u_so
144 ✗ m_taus1[0] = 2.7342; // tau_s1
145 ✗ m_taus2[0] = 16.; // tau_s2
146 ✗ m_ks[0] = 2.0994; // k_s
147 ✗ m_us[0] = 0.9087; // u_s
148 ✗ m_tausi[0] = 1.8875; // tau_si
149 ✗ m_tauwinf[0] = 0.07; // tau_{w inf}
150 ✗ m_winfstar[0] = 0.94; // w^*_{inf}
151 }
152 ✗ else if (FelisceParam::instance().CellsType == "PB") {
153 ✗ m_uo[0] = 0.; // u_o
154 ✗ m_uu[0] = 1.45; // u_u
155 ✗ m_thetav[0] = 0.35; // theta_v
156 ✗ m_thetaw[0] = 0.13; // theta_w
157 ✗ m_thetavm[0] = 0.175; // theta_v^-
158 ✗ m_thetao[0] = 0.006; // theta_o
159 ✗ m_tauvm1[0] = 10.; // tau_v1^-
160 ✗ m_tauvm2[0] = 1150.; // tau_v2^-
161 ✗ m_tauvp[0] = 1.4506; // tau_v^+
162 ✗ m_tauwm1[0] = 140.; // tau_w1^-
163 ✗ m_tauwm2[0] = 6.25; // tau_w2^+
164 ✗ m_kwm[0] = 65.; // k_w^-
165 ✗ m_uwm[0] = 0.015; // u_w^-
166 ✗ m_tauwm[0] = 326.; // tau_w^+
167 ✗ m_taufi[0] = 0.105; // tau_fi
168 ✗ m_tauo1[0] = 400.; // tau_o1
169 ✗ m_tauo2[0] = 6.; // tau_o2
170 ✗ m_tauso1[0] = 30.0181;// tau_so1
171 ✗ m_tauso2[0] = 0.9957; // tau_so2
172 ✗ m_kso[0] = 2.0458; // k_so
173 ✗ m_uso[0] = 0.65; // u_so
174 ✗ m_taus1[0] = 2.7342; // tau_s1
175 ✗ m_taus2[0] = 16.; // tau_s2
176 ✗ m_ks[0] = 2.0994; // k_s
177 ✗ m_us[0] = 0.9087; // u_s
178 ✗ m_tausi[0] = 1.8875; // tau_si
179 ✗ m_tauwinf[0] = 0.175; // tau_{w inf}
180 ✗ m_winfstar[0] = 0.9; // w^*_{inf}
181 }
182 ✗ else if (FelisceParam::instance().CellsType == "TNNP") {
183 ✗ m_uo[0] = 0.; // u_o
184 ✗ m_uu[0] = 1.58; // u_u
185 ✗ m_thetaw[0] = 0.015; // theta_w
186 ✗ m_thetavm[0] = 0.015; // theta_v^-
187 ✗ m_thetao[0] = 0.006; // theta_o
188 ✗ m_tauvm1[0] = 60.; // tau_v1^-
189 ✗ m_tauvm2[0] = 1150.; // tau_v2^-
190 ✗ m_tauvp[0] = 1.4506; // tau_v^+
191 ✗ m_tauwm1[0] = 70.; // tau_w1^-
192 ✗ m_tauwm2[0] = 20.; // tau_w2^+
193 ✗ m_kwm[0] = 65.; // k_w^-
194 ✗ m_uwm[0] = 0.03; // u_w^-
195 ✗ m_tauwm[0] = 280.; // tau_w^+
196 ✗ m_taufi[0] = 0.11; // tau_fi
197 ✗ m_tauo1[0] = 6.; // tau_o1
198 ✗ m_tauo2[0] = 6.; // tau_o2
199 ✗ m_tauso1[0] = 43.; // tau_so1
200 ✗ m_tauso2[0] = 0.2; // tau_so2
201 ✗ m_kso[0] = 2.; // k_so
202 ✗ m_uso[0] = 0.65; // u_so
203 ✗ m_taus1[0] = 2.7342; // tau_s1
204 ✗ m_taus2[0] = 3.; // tau_s2
205 ✗ m_ks[0] = 2.0994; // k_s
206 ✗ m_us[0] = 0.9087; // u_s
207 ✗ m_tausi[0] = 2.8723; // tau_si
208 ✗ m_tauwinf[0] = 0.07; // tau_{w inf}
209 ✗ m_winfstar[0] = 0.94; // w^*_{inf}
210 }
211
212 }
213 ✗ else if (FelisceParam::instance().CellsType == "hetero") {
214 ✗ m_uo.resize(m_cellTypeCut);
215 ✗ m_uu.resize(m_cellTypeCut);
216 ✗ m_thetav.resize(m_cellTypeCut);
217 ✗ m_thetaw.resize(m_cellTypeCut);
218 ✗ m_thetavm.resize(m_cellTypeCut);
219 ✗ m_thetao.resize(m_cellTypeCut);
220 ✗ m_tauvm1.resize(m_cellTypeCut);
221 ✗ m_tauvm2.resize(m_cellTypeCut);
222 ✗ m_tauvp.resize(m_cellTypeCut);
223 ✗ m_tauwm1.resize(m_cellTypeCut);
224 ✗ m_tauwm2.resize(m_cellTypeCut);
225 ✗ m_kwm.resize(m_cellTypeCut);
226 ✗ m_uwm.resize(m_cellTypeCut);
227 ✗ m_tauwm.resize(m_cellTypeCut);
228 ✗ m_taufi.resize(m_cellTypeCut);
229 ✗ m_tauo1.resize(m_cellTypeCut);
230 ✗ m_tauo2.resize(m_cellTypeCut);
231 ✗ m_tauso1.resize(m_cellTypeCut);
232 ✗ m_tauso2.resize(m_cellTypeCut);
233 ✗ m_kso.resize(m_cellTypeCut);
234 ✗ m_uso.resize(m_cellTypeCut);
235 ✗ m_taus1.resize(m_cellTypeCut);
236 ✗ m_taus2.resize(m_cellTypeCut);
237 ✗ m_ks.resize(m_cellTypeCut);
238 ✗ m_us.resize(m_cellTypeCut);
239 ✗ m_tausi.resize(m_cellTypeCut);
240 ✗ m_tauwinf.resize(m_cellTypeCut);
241 ✗ m_winfstar.resize(m_cellTypeCut);
242 ✗ m_gfi.resize(m_cellTypeCut);
243 ✗ m_gso.resize(m_cellTypeCut);
244 ✗ m_gsi.resize(m_cellTypeCut);
245 // Right Ventricle (RV)
246
247 ✗ m_uo[0] = 0.; // u_o
248 ✗ m_uu[0] = 1.55; // u_u
249 ✗ m_thetav[0] = 0.3; // theta_v
250 ✗ m_thetaw[0] = 0.13; // theta_w
251 ✗ m_thetavm[0] = 0.006; // theta_v^-
252 ✗ m_thetao[0] = 0.006; // theta_o
253 ✗ m_tauvm1[0] = 60.; // tau_v1^-
254 ✗ m_tauvm2[0] = 1150.; // tau_v2^-
255 ✗ m_tauvp[0] = 1.4506; // tau_v^+
256 ✗ m_tauwm1[0] = 60.; // tau_w1^-
257 ✗ m_tauwm2[0] = 15.; // tau_w2^+
258 ✗ m_kwm[0] = 65.; // k_w^-
259 ✗ m_uwm[0] = 0.03; // u_w^-
260 ✗ m_tauwm[0] = 200.; // tau_w^+
261 ✗ m_taufi[0] = 0.11; // tau_fi
262 ✗ m_tauo1[0] = 400.; // tau_o1
263 ✗ m_tauo2[0] = 6.; // tau_o2
264 ✗ if (FelisceParam::instance().tau_so_1_rv < 0.) {
265 ✗ m_tauso1[0] = 20.; //26.; // tau_so1 (30.0181)
266 }
267 else {
268 ✗ m_tauso1[0] = FelisceParam::instance().tau_so_1_rv;
269 }
270 ✗ m_tauso2[0] = 0.9957; // tau_so2
271 ✗ m_kso[0] = 2.0458; // k_so
272 ✗ m_uso[0] = 0.65; // u_so
273 ✗ m_taus1[0] = 2.7342; // tau_s1
274 ✗ m_taus2[0] = 16.; // tau_s2
275 ✗ m_ks[0] = 2.0994; // k_s
276 ✗ m_us[0] = 0.9087; // u_s
277 ✗ m_tausi[0] = 1.8875; // tau_si
278 ✗ m_tauwinf[0] = 0.07; // tau_{w inf}
279 ✗ m_winfstar[0] = 0.94; // w^*_{inf}
280 ✗ m_gfi[0] = FelisceParam::instance().gfi_rv;
281 ✗ m_gso[0] = FelisceParam::instance().gso_rv;
282 ✗ m_gsi[0] = FelisceParam::instance().gsi_rv;
283 // Endo
284
285 ✗ m_uo[1] = 0.; // u_o
286 ✗ m_uu[1] = 1.56; // u_u
287 ✗ m_thetav[1] = 0.3; // theta_v
288 ✗ m_thetaw[1] = 0.13; // theta_w
289 ✗ m_thetavm[1] = 0.2; // theta_v^-
290 ✗ m_thetao[1] = 0.006; // theta_o
291 ✗ m_tauvm1[1] = 75.; // tau_v1^-
292 ✗ m_tauvm2[1] = 10.; // tau_v2^-
293 ✗ m_tauvp[1] = 1.4506; // tau_v^+
294 ✗ m_tauwm1[1] = 6.; // tau_w1^-
295 ✗ m_tauwm2[1] = 140.; // tau_w2^+
296 ✗ m_kwm[1] = 200.; // k_w^-
297 ✗ m_uwm[1] = 0.016; // u_w^-
298 ✗ m_tauwm[1] = 280.; // tau_w^+
299 ✗ m_taufi[1] = 0.1; // tau_fi
300 ✗ m_tauo1[1] = 470.; // tau_o1
301 ✗ m_tauo2[1] = 6.; // tau_o2
302 ✗ if (FelisceParam::instance().tau_so_1_endo < 0.) {
303 ✗ m_tauso1[1] = 30.; //40.; // tau_so1 (40.)
304 }
305 else {
306 ✗ m_tauso1[1] = FelisceParam::instance().tau_so_1_endo;
307 }
308 ✗ m_tauso2[1] = 1.2; // tau_so2
309 ✗ m_kso[1] = 2.; // k_so
310 ✗ m_uso[1] = 0.65; // u_so
311 ✗ m_taus1[1] = 2.7342; // tau_s1
312 ✗ m_taus2[1] = 2.; // tau_s2
313 ✗ m_ks[1] = 2.0994; // k_s
314 ✗ m_us[1] = 0.9087; // u_s
315 ✗ m_tausi[1] = 2.9013; // tau_si
316 ✗ m_tauwinf[1] = 0.0273;// tau_{w inf}
317 ✗ m_winfstar[1] = 0.78; // w^*_{inf}
318 ✗ m_gfi[1] = FelisceParam::instance().gfi_endo;
319 ✗ m_gso[1] = FelisceParam::instance().gso_endo;
320 ✗ m_gsi[1] = FelisceParam::instance().gsi_endo;
321
322 // M-Cells
323
324 ✗ m_uo[2] = 0.; // u_o
325 ✗ m_uu[2] = 1.61; // u_u
326 ✗ m_thetav[2] = 0.3; // theta_v
327 ✗ m_thetaw[2] = 0.13; // theta_w
328 ✗ m_thetavm[2] = 0.1; // theta_v^-
329 ✗ m_thetao[2] = 0.005; // theta_o
330 ✗ m_tauvm1[2] = 80.; // tau_v1^-
331 ✗ m_tauvm2[2] = 1.4506; // tau_v2^-
332 ✗ m_tauvp[2] = 1.4506; // tau_v^+
333 ✗ m_tauwm1[2] = 70.; // tau_w1^-
334 ✗ m_tauwm2[2] = 8.; // tau_w2^+
335 ✗ m_kwm[2] = 200.; // k_w^-
336 ✗ m_uwm[2] = 0.016; // u_w^-
337 ✗ m_tauwm[2] = 280.; // tau_w^+
338 ✗ m_taufi[2] = 0.078; // tau_fi
339 ✗ m_tauo1[2] = 410.; // tau_o1
340 ✗ m_tauo2[2] = 7.; // tau_o2
341 ✗ if (FelisceParam::instance().tau_so_1_mcel < 0.) {
342 ✗ m_tauso1[2] = 45.; //61.; // tau_so1 (91.)
343 }
344 else {
345 ✗ m_tauso1[2] = FelisceParam::instance().tau_so_1_mcel;
346 }
347 ✗ m_tauso2[2] = 0.8; // tau_so2
348 ✗ m_kso[2] = 2.1; // km_so
349 ✗ m_uso[2] = 0.6; // u_so
350 ✗ m_taus1[2] = 2.7342; // tau_s1
351 ✗ m_taus2[2] = 4.; // tau_s2
352 ✗ m_ks[2] = 2.0994; // k_s
353 ✗ m_us[2] = 0.9087; // u_s
354 ✗ m_tausi[2] = 3.3849; // tau_si
355 ✗ m_tauwinf[2] = 0.01; // tau_{w inf}
356 ✗ m_winfstar[2] = 0.5; // w^*_{inf}
357 ✗ m_gfi[2] = FelisceParam::instance().gfi_mid;
358 ✗ m_gso[2] = FelisceParam::instance().gso_mid;
359 ✗ m_gsi[2] = FelisceParam::instance().gsi_mid;
360 // Epi
361
362 ✗ m_uo[3] = 0.; // u_o
363 ✗ m_uu[3] = 1.55; // u_u
364 ✗ m_thetav[3] = 0.3; // theta_v
365 ✗ m_thetaw[3] = 0.13; // theta_w
366 ✗ m_thetavm[3] = 0.006; // theta_v^-
367 ✗ m_thetao[3] = 0.006; // theta_o
368 ✗ m_tauvm1[3] = 60.; // tau_v1^-
369 ✗ m_tauvm2[3] = 1150.; // tau_v2^-
370 ✗ m_tauvp[3] = 1.4506; // tau_v^+
371 ✗ m_tauwm1[3] = 60.; // tau_w1^-
372 ✗ m_tauwm2[3] = 15.; // tau_w2^+
373 ✗ m_kwm[3] = 65.; // k_w^-
374 ✗ m_uwm[3] = 0.03; // u_w^-
375 ✗ m_tauwm[3] = 200.; // tau_w^+
376 ✗ m_taufi[3] = 0.11; // tau_fi
377 ✗ m_tauo1[3] = 400.; // tau_o1
378 ✗ m_tauo2[3] = 6.; // tau_o2
379 ✗ if (FelisceParam::instance().tau_so_1_epi < 0.) {
380 ✗ m_tauso1[3] = 19.; //25.; // tau_so1 (30.0181)
381 }
382 else {
383 ✗ m_tauso1[3] = FelisceParam::instance().tau_so_1_epi;
384 }
385 ✗ m_tauso2[3] = 0.9957; // tau_so2
386 ✗ m_kso[3] = 2.0458; // k_so
387 ✗ m_uso[3] = 0.65; // u_so
388 ✗ m_taus1[3] = 2.7342; // tau_s1
389 ✗ m_taus2[3] = 16.; // tau_s2
390 ✗ m_ks[3] = 2.0994; // k_s
391 ✗ m_us[3] = 0.9087; // u_s
392 ✗ m_tausi[3] = 1.8875; // tau_si
393 ✗ m_tauwinf[3] = 0.07; // tau_{w inf}
394 ✗ m_winfstar[3] = 0.94; // w^*_{inf}
395 ✗ m_gfi[3] = FelisceParam::instance().gfi_epi;
396 ✗ m_gso[3] = FelisceParam::instance().gso_epi;
397 ✗ m_gsi[3] = FelisceParam::instance().gsi_epi;
398
399 // Atria (Atria-ventricular MV model)
400 ✗ if (m_cellTypeCut > 4){
401 ✗ m_uo[4] = 0.; // u_o
402 ✗ m_uu[4] = 1.55; // u_u
403 ✗ m_thetav[4] = 0.3; // theta_v
404 ✗ m_thetaw[4] = 0.13; // theta_w
405 ✗ m_thetavm[4] = 0.006; // theta_v^-
406 ✗ m_thetao[4] = 0.006; // theta_o
407 ✗ m_tauvm1[4] = 60.; // tau_v1^-
408 ✗ m_tauvm2[4] = 1150.; // tau_v2^-
409 ✗ m_tauvp[4] = 1.4506; // tau_v^+
410 ✗ m_tauwm1[4] = 60.; // tau_w1^-
411 ✗ m_tauwm2[4] = 15.; // tau_w2^+
412 ✗ m_kwm[4] = 65.; // k_w^-
413 ✗ m_uwm[4] = 0.03; // u_w^-
414 ✗ m_tauwm[4] = 200.; // tau_w^+
415 ✗ m_taufi[4] = 0.11; // tau_fi
416 ✗ m_tauo1[4] = 400.; // tau_o1
417 ✗ m_tauo2[4] = 6.; // tau_o2
418 ✗ m_tauso1[4] = 30.; // tau_so1
419 ✗ m_tauso2[4] = 0.9957; // tau_so2
420 ✗ m_kso[4] = 2.0458; // k_so
421 ✗ m_uso[4] = 0.65; // u_so
422 ✗ m_taus1[4] = 2.7342; // tau_s1
423 ✗ m_taus2[4] = 16.; // tau_s2
424 ✗ m_ks[4] = 2.0994; // k_s
425 ✗ m_us[4] = 0.9087; // u_s
426 ✗ m_tausi[4] = 1.8875; // tau_si
427 ✗ m_tauwinf[4] = 0.07; // tau_{w inf}
428 ✗ m_winfstar[4] = 0.94; // w^*_{inf}
429 ✗ m_gfi[4] = 1.;
430 ✗ m_gso[4] = 1.;
431 ✗ m_gsi[4] = 1.;
432 }
433 }
434 else {
435 ✗ FEL_ERROR("Cells type not defined.");
436 }
437
438 }
439
440 ✗ MVSolver::~MVSolver()
441 = default;
442
443 ✗ void MVSolver::computeRHS() {
444
445 ✗ double& dt = m_fstransient->timeStep;
446
447 ✗ m_bdf.computeRHSTime(dt, m_vecRHS);
448
449 double value_uExtrap;
450 double value_RHSv;
451 double value_RHSw;
452 double value_RHSs;
453 double hv;
454 double hw;
455 double hm;
456 double ho;
457 double vinf;
458 double tauvm;
459 double winf;
460 double tauwm;
461 double taus;
462
463 ✗ double Vmin = FelisceParam::instance().vMin;
464 ✗ double Vmax = FelisceParam::instance().vMax;
465
466 ✗ int iType = 0;
467
468 felInt pos;
469 ✗ for (felInt i = 0; i < m_size; i++) {
470 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
471
472 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
473
474 //________________________________________________________________________________________________________________________________________________________________
475
476 ✗ if(m_assimilation ){
477
478 //For assimilation
479 ✗ if(m_fstransient->iteration == 1) {
480 ✗ m_vEDO.resize(m_size);
481 ✗ m_wEDO.resize(m_size);
482 ✗ m_sEDO.resize(m_size);
483
484 ✗ p_vEDO.resize(m_size);
485 ✗ p_wEDO.resize(m_size);
486 ✗ p_sEDO.resize(m_size);
487 }
488
489 ✗ if(m_fstransient->iteration == 1) {
490 ✗ m_Vm.resize(m_size);
491 ✗ p_Vm.resize(m_size);
492 ✗ m_Vm[pos] = (m_Vm[pos] - Vmin)/(Vmax-Vmin);
493 ✗ p_Vm[pos] = (p_Vm[pos] - Vmin)/(Vmax-Vmin);
494 }
495
496 ✗ if(m_fstransient->iteration > 1) {
497 ✗ if(m_assimilation) {
498 ✗ value_uExtrap = m_Vm[pos];
499 }
500 }
501
502 }
503 //________________________________________________________________________________________________________________________________________________________________
504
505 ✗ value_uExtrap = (value_uExtrap - Vmin)/(Vmax-Vmin); // \in [0,1]
506
507 ✗ if (FelisceParam::instance().CellsType == "hetero")
508 ✗ iType = m_cellType[pos];
509
510 ✗ if (iType < m_cellTypeCut) {
511
512 // hm = H(u-theta_v^-)
513 ✗ if (value_uExtrap > m_thetavm[iType])
514 ✗ hm=1.;
515 else
516 ✗ hm=0.;
517
518 // v_inf = 1 - H(u-theta_v^-)
519 ✗ vinf = 1-hm;
520
521 // hv = H(u-theta_v)
522 ✗ if (value_uExtrap > m_thetav[iType])
523 ✗ hv=1.;
524 else
525 ✗ hv=0.;
526
527 // hw = H(u-theta_w)
528 ✗ if (value_uExtrap > m_thetaw[iType])
529 ✗ hw=1.;
530 else
531 ✗ hw=0.;
532
533 // ho = H(u-theta_o)
534 ✗ if(value_uExtrap > m_thetao[iType])
535 ✗ ho = 1.;
536 else
537 ✗ ho = 0.;
538
539 // w_inf = (1 - H(u-theta_o)) * (1 / u/tau_w_inf) + H(u-theat_o)*w_inf^*
540 ✗ winf = (1-ho)*(1-value_uExtrap/m_tauwinf[iType]) + ho*m_winfstar[iType];
541
542 // tau_v^- = (1 - H(u-theta_v^+))*tau_v1^- + H(u-theta_v^-)*tau_v2^-
543 ✗ tauvm = (1-hm)*m_tauvm1[iType] + hm*m_tauvm2[iType];
544 // dv/dt = (1 - H(u-theta_v))*(v_inf - v)/tau_v^- - H(u-theta_v)*v/tau_v^+
545 ✗ value_RHSv = (1-hv)*vinf/tauvm;
546 ✗ m_vecRHS[0].setValue(pos,value_RHSv, ADD_VALUES);
547
548 // tau_w^- = tau_w1^- + (tau_w2^- - tau_W1^-) * (1 + tanh(k_w^- * (u-u-u_w^-)) )/2
549 ✗ tauwm = m_tauwm1[iType] + (m_tauwm2[iType]-m_tauwm1[iType]) * (1 + tanh(m_kwm[iType]*(value_uExtrap-m_uwm[iType])) )/2;
550 // dw/dt = (1 - H(u-theta_w))*(w_inf - w)/tau_w^- - H(u-theta_w)*w/tau_w^+
551 ✗ value_RHSw = (1-hw)*winf/tauwm;
552 ✗ m_vecRHS[1].setValue(pos,value_RHSw, ADD_VALUES);
553
554 // tau_s = (1 - H(u-theta_w))*tau_s1 + H(u-theta_w)*tau_s2
555 ✗ taus = (1-hw)*m_taus1[iType] + hw*m_taus2[iType];
556 // ds/dt = ((1 + tanh(k_s * (u-u_s)) )/2 - s)/tau_s
557 ✗ value_RHSs = (1 + tanh(m_ks[iType]*(value_uExtrap-m_us[iType])))/(2*taus);
558 ✗ m_vecRHS[2].setValue(pos,value_RHSs, ADD_VALUES);
559 }
560 ✗ else if (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent" ){
561 ✗ felInt meshId = pos;
562 ✗ AOPetscToApplication(m_ao,1,&meshId);
563 ✗ felInt sizeVent = 45580;
564 ✗ if (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart")
565 ✗ sizeVent = 32320;
566 ✗ if (meshId > (sizeVent-1) ) {
567 ✗ m_vecRHS[0].setValue(pos,0., ADD_VALUES);
568 ✗ m_vecRHS[1].setValue(pos,0., ADD_VALUES);
569 ✗ m_vecRHS[2].setValue(pos,0., ADD_VALUES);
570 }
571 else {
572 ✗ std::ostringstream msg;
573 ✗ msg << "Problem in MVsolver atria/ventricles nodes : node " << pos << " not defined, " << iType << "." << std::endl;
574 ✗ FEL_ERROR(msg.str());
575 }
576 }
577 else {
578 ✗ m_vecRHS[0].setValue(pos,0., ADD_VALUES);
579 ✗ m_vecRHS[1].setValue(pos,0., ADD_VALUES);
580 ✗ m_vecRHS[2].setValue(pos,0., ADD_VALUES);
581 }
582 }
583
584 ✗ m_vecRHS[0].assembly();
585 ✗ m_vecRHS[1].assembly();
586 ✗ m_vecRHS[2].assembly();
587 }
588
589 ✗ void MVSolver::solveEDO() {
590 ✗ double& coeffDeriv = m_bdf.coeffDeriv0();
591 ✗ double& dt = m_fstransient->timeStep;
592 ✗ int& numIt = m_fstransient->iteration;
593 double value_uExtrap;
594
595 double value_RHSv;
596 double value_RHSw;
597 double value_RHSs;
598 double valuem_solEDOv;
599 double valuem_solEDOw;
600 double valuem_solEDOs;
601
602 double beta1;
603 double beta2;
604 double beta3;
605
606 double tauvm;
607 double tauwm;
608 double taus;
609 double hm;
610 double hv;
611 double hw;
612
613 ✗ double Vmin = FelisceParam::instance().vMin;
614 ✗ double Vmax = FelisceParam::instance().vMax;
615
616 ✗ int iType = 0;
617
618 felInt pos;
619 ✗ for (felInt i = 0; i < m_size; i++) {
620 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
621
622 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
623 ✗ m_vecRHS[0].getValues(1,&pos,&value_RHSv);//value_RHS = m_RHS(i)
624 ✗ m_vecRHS[1].getValues(1,&pos,&value_RHSw);//value_RHS = m_RHS(i)
625 ✗ m_vecRHS[2].getValues(1,&pos,&value_RHSs);//value_RHS = m_RHS(i)
626
627 //________________________________________________________________________________________________________________________________________________________________
628
629 ✗ if(m_assimilation) {
630 ✗ if(numIt > 1) {
631 ✗ value_uExtrap = m_Vm[pos];
632 }
633 }
634
635 //________________________________________________________________________________________________________________________________________________________________
636
637 ✗ value_uExtrap = (value_uExtrap - Vmin)/(Vmax-Vmin); // \in [0,1]
638
639 ✗ if (FelisceParam::instance().CellsType == "hetero")
640 ✗ iType = m_cellType[pos];
641
642 ✗ if (iType < m_cellTypeCut) {
643
644 // hm = H(u-theta_v^-)
645 ✗ if (value_uExtrap > m_thetavm[iType])
646 ✗ hm=1.;
647 else
648 ✗ hm=0.;
649
650 // hv = H(u-theta_v)
651 ✗ if (value_uExtrap > m_thetav[iType])
652 ✗ hv=1.;
653 else
654 ✗ hv=0.;
655
656 // hw = H(u-theta_w)
657 ✗ if (value_uExtrap > m_thetaw[iType])
658 ✗ hw=1.;
659 else
660 ✗ hw=0.;
661
662 // tau_v^- = (1 - H(u-theta_v^-))*tau_v1^- + H(u-theta_v^-)*tau_v2^-
663 ✗ tauvm = (1-hm)*m_tauvm1[iType] + hm*m_tauvm2[iType];
664 ✗ beta1 = coeffDeriv/dt + (1-hv)/tauvm + hv/m_tauvp[iType];
665
666 // tau_w^- = tau_w1^- + (tau_w2^- - tau_W1^-) * (1 + tanh(k_w^- * (u-u-u_w^-)) )/2
667 ✗ tauwm = m_tauwm1[iType] + (m_tauwm2[iType]-m_tauwm1[iType]) * (1 + tanh(m_kwm[iType]*(value_uExtrap-m_uwm[iType])) )/2;
668 ✗ beta2=coeffDeriv/dt + (1-hw)/tauwm + hw/m_tauwm[iType];
669
670 // tau_s = (1 - H(u-theta_w))*tau_s1 + H(u-theta_w)*tau_s2
671 ✗ taus = (1-hw)*m_taus1[iType] + hw*m_taus2[iType];
672 ✗ beta3 = coeffDeriv/dt + 1./taus;
673
674 // v^n+1 = RHSv /beta1
675 ✗ valuem_solEDOv = value_RHSv * 1./beta1;
676 // w^n+1 = RHSw / beta1
677 ✗ valuem_solEDOw = value_RHSw * 1./beta2;
678 // s^n+1 = RHSs / beta3
679 ✗ valuem_solEDOs = value_RHSs * 1./beta3;
680
681 ✗ if(m_assimilation)
682 {
683
684 ✗ if(numIt == 1) {
685 ✗ valuem_solEDOv = 1.;
686 ✗ valuem_solEDOw = 1.;
687 ✗ valuem_solEDOs = 0.;
688
689 ✗ m_vEDO[pos] = valuem_solEDOv;
690 ✗ m_wEDO[pos] = valuem_solEDOw;
691 ✗ m_sEDO[pos] = valuem_solEDOs;
692
693 }
694
695 //________________________________________________________________________________________________________________________________________________________________
696
697
698 ✗ p_vEDO[pos] = valuem_solEDOv;
699 ✗ p_wEDO[pos] = valuem_solEDOw;
700 ✗ p_sEDO[pos] = valuem_solEDOs;
701
702 ✗ valuem_solEDOv = m_vEDO[pos];
703 ✗ valuem_solEDOw = m_wEDO[pos];
704 ✗ valuem_solEDOs = m_sEDO[pos];
705 }
706
707 //________________________________________________________________________________________________________________________________________________________________
708
709 ✗ m_vecSolEDO[0].setValue(pos,valuem_solEDOv, INSERT_VALUES);
710 ✗ m_vecSolEDO[1].setValue(pos,valuem_solEDOw, INSERT_VALUES);
711 ✗ m_vecSolEDO[2].setValue(pos,valuem_solEDOs, INSERT_VALUES);
712 }
713 ✗ else if (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent" ){
714 ✗ felInt meshId = pos;
715 ✗ AOPetscToApplication(m_ao,1,&meshId);
716 ✗ felInt sizeVent = 45580;
717 ✗ if (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart")
718 ✗ sizeVent = 32320;
719 ✗ if (meshId > (sizeVent-1) ) {
720 ✗ m_vecSolEDO[0].setValue(pos,0., INSERT_VALUES);
721 ✗ m_vecSolEDO[1].setValue(pos,0., INSERT_VALUES);
722 ✗ m_vecSolEDO[2].setValue(pos,0., INSERT_VALUES);
723 }
724 else {
725 ✗ std::ostringstream msg;
726 ✗ msg << "Problem in MVsolver atria/ventricles nodes : node " << pos << " not defined, " << iType << "." << std::endl;
727 ✗ FEL_ERROR(msg.str());
728 }
729 }
730 else {
731 ✗ m_vecSolEDO[0].setValue(pos,0., INSERT_VALUES);
732 ✗ m_vecSolEDO[1].setValue(pos,0., INSERT_VALUES);
733 ✗ m_vecSolEDO[2].setValue(pos,0., INSERT_VALUES);
734 }
735
736 }
737
738 ✗ m_vecSolEDO[0].assembly();
739 ✗ m_vecSolEDO[1].assembly();
740 ✗ m_vecSolEDO[2].assembly();
741 }
742
743 ✗ void MVSolver::computeIon() {
744 double value_uExtrap;
745 double valuem_solEDOv;
746 double valuem_solEDOw;
747 double valuem_solEDOs;
748 double value_ion;
749
750 double Jfi;
751 double Jso;
752 double Jsi;
753 double hv;
754 double hw;
755 double ho;
756 double tauo;
757 double tauso;
758
759 ✗ double Vmin = FelisceParam::instance().vMin;
760 ✗ double Vmax = FelisceParam::instance().vMax;
761 double correctionCoef;
762 ✗ int iType = 0;
763
764 felInt pos;
765 ✗ for (felInt i = 0; i < m_size; i++) {
766 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
767
768 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);//value_uExtrap = m_uExtrap(i)
769 ✗ m_vecSolEDO[0].getValues(1,&pos,&valuem_solEDOv);//valuem_solEDOv = m_vecSolEDO[0](i)
770 ✗ m_vecSolEDO[1].getValues(1,&pos,&valuem_solEDOw);//valuem_solEDOw = m_vecSolEDO[1](i)
771 ✗ m_vecSolEDO[2].getValues(1,&pos,&valuem_solEDOs);//valuem_solEDOs = m_vecSolEDO[2](i)
772
773
774 //________________________________________________________________________________________________________________________________________________________________
775
776 ✗ if(m_assimilation) {
777 ✗ if(m_fstransient->iteration > 1) {
778 ✗ value_uExtrap = m_Vm[pos];
779 }
780 }
781
782 //________________________________________________________________________________________________________________________________________________________________
783
784
785 ✗ value_uExtrap = (value_uExtrap - Vmin)/(Vmax-Vmin); // \in [0,1]
786
787 ✗ if (FelisceParam::instance().CellsType == "hetero")
788 ✗ iType = m_cellType[pos];
789
790 ✗ if (iType < m_cellTypeCut) {
791
792 ✗ if (value_uExtrap > m_thetav[iType])
793 ✗ hv=1.;
794 else
795 ✗ hv=0.;
796
797 ✗ if (value_uExtrap > m_thetaw[iType])
798 ✗ hw=1.;
799 else
800 ✗ hw=0.;
801
802 ✗ if (value_uExtrap > m_thetao[iType])
803 ✗ ho=1.;
804 else
805 ✗ ho=0.;
806
807 // tau_o = (1 - H(u-thea_o))
808 ✗ tauo = (1-ho)*m_tauo1[iType] + ho*m_tauo2[iType];
809 // tau_so = tau_so1 + (tau_so2 - tau_so1)*(1 + tanh(k_so * (u-u_so)) )/2
810 ✗ tauso = m_tauso1[iType] + (m_tauso2[iType] - m_tauso1[iType])*(1 + tanh(m_kso[iType]*(value_uExtrap-m_uso[iType])) )/2;
811
812 // J_fi = -v*H(u-theta_v)*(u-theta_v)*(u_u-u)/tau_fi
813 ✗ Jfi = - valuem_solEDOv * hv * (value_uExtrap-m_thetav[iType]) * (m_uu[iType]-value_uExtrap) / m_taufi[iType];
814 ✗ Jfi *= m_gfi[iType];
815 // J_so = (u-u_o)*(1 - H(u-theta_w))/tau_o + H(u-theta_w)/tau_so
816 ✗ Jso = (value_uExtrap - m_uo[iType]) * (1 - hw)/tauo + hw/tauso;
817 ✗ Jso *= m_gso[iType];
818 // J_si = -H(u-theta_w)*w*s/tau_si
819 ✗ Jsi = - hw * valuem_solEDOw * valuem_solEDOs / m_tausi[iType];
820 ✗ Jsi *= m_gsi[iType];
821 ✗ correctionCoef = FelisceParam::instance().Cm * (FelisceParam::instance().vMax - FelisceParam::instance().vMin);
822 // Infarct
823 ✗ if (FelisceParam::instance().hasInfarct) {
824 ✗ Jfi *= m_tauOut[pos];
825 }
826 ✗ value_ion = - ( Jfi + Jso + Jsi )*correctionCoef;
827 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
828 }
829 ✗ else if (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent" ){
830 ✗ felInt meshId = pos;
831 ✗ AOPetscToApplication(m_ao,1,&meshId);
832 ✗ felInt sizeVent = 45580;
833 ✗ if (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart")
834 ✗ sizeVent = 32320;
835 ✗ if (meshId > (sizeVent-1) ) {
836 ✗ m_ion.setValue(pos,0., INSERT_VALUES);
837 }
838 else {
839 ✗ std::ostringstream msg;
840 ✗ msg << "Problem in MVsolver atria/ventricles nodes : node " << pos << " not defined, " << iType << "." << std::endl;
841 ✗ FEL_ERROR(msg.str());
842 }
843 }
844 else {
845 ✗ m_ion.setValue(pos,0., INSERT_VALUES);
846 }
847 }
848 ✗ m_ion.assembly();
849 }
850
851
852 /////////////////////////
853 // Setters and Getters //
854 /////////////////////////
855
856
857 ✗ void MVSolver::setParameters(std::vector<double> parameters) {
858
859
860 ✗ if(parameters.size() == 28) {
861
862 ✗ m_uo.resize(1);
863 ✗ m_uu.resize(1);
864 ✗ m_thetav.resize(1);
865 ✗ m_thetaw.resize(1);
866 ✗ m_thetavm.resize(1);
867 ✗ m_thetao.resize(1);
868 ✗ m_tauvm1.resize(1);
869 ✗ m_tauvm2.resize(1);
870 ✗ m_tauvp.resize(1);
871 ✗ m_tauwm1.resize(1);
872 ✗ m_tauwm2.resize(1);
873 ✗ m_kwm.resize(1);
874 ✗ m_uwm.resize(1);
875 ✗ m_tauwm.resize(1);
876 ✗ m_taufi.resize(1);
877 ✗ m_tauo1.resize(1);
878 ✗ m_tauo2.resize(1);
879 ✗ m_tauso1.resize(1);
880 ✗ m_tauso2.resize(1);
881 ✗ m_kso.resize(1);
882 ✗ m_uso.resize(1);
883 ✗ m_taus1.resize(1);
884 ✗ m_taus2.resize(1);
885 ✗ m_ks.resize(1);
886 ✗ m_us.resize(1);
887 ✗ m_tausi.resize(1);
888 ✗ m_tauwinf.resize(1);
889 ✗ m_winfstar.resize(1);
890
891
892
893 ✗ m_uo[0] = parameters[0];
894 ✗ m_uu[0] = parameters[1];
895 ✗ m_thetav[0] = parameters[2];
896 ✗ m_thetaw[0] = parameters[3];
897 ✗ m_thetavm[0] = parameters[4];
898 ✗ m_thetao[0] = parameters[5];
899 ✗ m_tauvm1[0] = parameters[6];
900 ✗ m_tauvm2[0] = parameters[7];
901 ✗ m_tauvp[0] = parameters[8];
902 ✗ m_tauwm1[0] = parameters[9];
903 ✗ m_tauwm2[0] = parameters[10];
904 ✗ m_kwm[0] = parameters[11];
905 ✗ m_uwm[0] = parameters[12];
906 ✗ m_tauwm[0] = parameters[13];
907 ✗ m_taufi[0] = parameters[14];
908 ✗ m_tauo1[0] = parameters[15];
909 ✗ m_tauo2[0] = parameters[16];
910 ✗ m_tauso1[0] = parameters[17];
911 ✗ m_tauso2[0] = parameters[18];
912 ✗ m_kso[0] = parameters[19];
913 ✗ m_uso[0] = parameters[20];
914 ✗ m_taus1[0] = parameters[21];
915 ✗ m_taus2[0] = parameters[22];
916 ✗ m_ks[0] = parameters[23];
917 ✗ m_us[0] = parameters[24];
918 ✗ m_tausi[0] = parameters[25];
919 ✗ m_tauwinf[0] = parameters[26];
920 ✗ m_winfstar[0] = parameters[27];
921
922 } else {
923
924 ✗ for(unsigned int i=0 ; i<m_posparam.size() ; i++) {
925
926 ✗ const int posp = m_posparam[i];
927
928 ✗ switch(posp) {
929
930 ✗ case 0:
931 ✗ m_uo.resize(1);
932 ✗ m_uo[0] = parameters[i];
933 ✗ break;
934
935 ✗ case 1:
936 ✗ m_uu.resize(1);
937 ✗ m_uu[0] = parameters[i];
938 ✗ break;
939
940 ✗ case 2:
941 ✗ m_thetav.resize(1);
942 ✗ m_thetav[0] = parameters[i];
943 ✗ break;
944
945 ✗ case 3:
946 ✗ m_thetaw.resize(1);
947 ✗ m_thetaw[0] = parameters[i];
948 ✗ break;
949
950 ✗ case 4:
951 ✗ m_thetavm.resize(1);
952 ✗ m_thetavm[0] = parameters[i];
953 ✗ break;
954
955 ✗ case 5:
956 ✗ m_thetao.resize(1);
957 ✗ m_thetao[0] = parameters[i];
958 ✗ break;
959
960 ✗ case 6:
961 ✗ m_tauvm1.resize(1);
962 ✗ m_tauvm1[0] = parameters[i];
963 ✗ break;
964
965 ✗ case 7:
966 ✗ m_tauvm2.resize(1);
967 ✗ m_tauvm2[0] = parameters[i];
968 ✗ break;
969
970 ✗ case 8:
971 ✗ m_tauvp.resize(1);
972 ✗ m_tauvp[0] = parameters[i];
973 ✗ break;
974
975 ✗ case 9:
976 ✗ m_tauwm1.resize(1);
977 ✗ m_tauwm1[0] = parameters[i];
978 ✗ break;
979
980 ✗ case 10:
981 ✗ m_tauwm2.resize(1);
982 ✗ m_tauwm2[0] = parameters[i];
983 ✗ break;
984
985 ✗ case 11:
986 ✗ m_kwm.resize(1);
987 ✗ m_kwm[0] = parameters[i];
988 ✗ break;
989
990 ✗ case 12:
991 ✗ m_uwm.resize(1);
992 ✗ m_uwm[0] = parameters[i];
993 ✗ break;
994
995 ✗ case 13:
996 ✗ m_tauwm.resize(1);
997 ✗ m_tauwm[0] = parameters[i];
998 ✗ break;
999
1000 ✗ case 14:
1001 ✗ m_taufi.resize(1);
1002 ✗ m_taufi[0] = parameters[i];
1003 ✗ break;
1004
1005 ✗ case 15:
1006 ✗ m_tauo1.resize(1);
1007 ✗ m_tauo1[0] = parameters[i];
1008 ✗ break;
1009
1010 ✗ case 16:
1011 ✗ m_tauo2.resize(1);
1012 ✗ m_tauo2[0] = parameters[i];
1013 ✗ break;
1014
1015 ✗ case 17:
1016 ✗ m_tauso1.resize(1);
1017 ✗ m_tauso1[0] = parameters[i];
1018 ✗ break;
1019
1020 ✗ case 18:
1021 ✗ m_tauso2.resize(1);
1022 ✗ m_tauso2[0] = parameters[i];
1023 ✗ break;
1024
1025 ✗ case 19:
1026 ✗ m_kso.resize(1);
1027 ✗ m_kso[0] = parameters[i];
1028 ✗ break;
1029
1030 ✗ case 20:
1031 ✗ m_uso.resize(1);
1032 ✗ m_uso[0] = parameters[i];
1033 ✗ break;
1034
1035 ✗ case 21:
1036 ✗ m_taus1.resize(1);
1037 ✗ m_taus1[0] = parameters[i];
1038 ✗ break;
1039
1040 ✗ case 22:
1041 ✗ m_taus2.resize(1);
1042 ✗ m_taus2[0] = parameters[i];
1043 ✗ break;
1044
1045 ✗ case 23:
1046 ✗ m_ks.resize(1);
1047 ✗ m_ks[0] = parameters[i];
1048 ✗ break;
1049
1050 ✗ case 24:
1051 ✗ m_us.resize(1);
1052 ✗ m_us[0] = parameters[i];
1053 ✗ break;
1054
1055 ✗ case 25:
1056 ✗ m_tausi.resize(1);
1057 ✗ m_tausi[0] = parameters[i];
1058 ✗ break;
1059
1060 ✗ case 26:
1061 ✗ m_tauwinf.resize(1);
1062 ✗ m_tauwinf[0] = parameters[i];
1063 ✗ break;
1064
1065 ✗ case 27:
1066 ✗ m_winfstar.resize(1);
1067 ✗ m_winfstar[0] = parameters[i];
1068 ✗ break;
1069
1070 ✗ default:
1071 ✗ std::cout << "Numerotation error in setParameters of MVSolver" << std::endl;
1072 ✗ break;
1073 }
1074
1075 }//End for
1076
1077 }//End else
1078
1079 }
1080
1081
1082
1083 ✗ void MVSolver::getParameters(std::vector<double>& p) {
1084
1085 ✗ std::vector<double> param;
1086
1087 ✗ if(m_posparam.size() == 28) {
1088
1089 ✗ param.push_back(m_uo[0]);
1090 ✗ param.push_back(m_uu[0]);
1091 ✗ param.push_back(m_thetav[0]);
1092 ✗ param.push_back(m_thetaw[0]);
1093 ✗ param.push_back(m_thetavm[0]);
1094 ✗ param.push_back(m_thetao[0]);
1095 ✗ param.push_back(m_tauvm1[0]);
1096 ✗ param.push_back(m_tauvm2[0]);
1097 ✗ param.push_back(m_tauvp[0]);
1098 ✗ param.push_back(m_tauwm1[0]);
1099 ✗ param.push_back(m_tauwm2[0]);
1100 ✗ param.push_back(m_kwm[0]);
1101 ✗ param.push_back(m_uwm[0]);
1102 ✗ param.push_back(m_tauwm[0]);
1103 ✗ param.push_back(m_taufi[0]);
1104 ✗ param.push_back(m_tauo1[0]);
1105 ✗ param.push_back(m_tauo2[0]);
1106 ✗ param.push_back(m_tauso1[0]);
1107 ✗ param.push_back(m_tauso2[0]);
1108 ✗ param.push_back(m_kso[0]);
1109 ✗ param.push_back(m_uso[0]);
1110 ✗ param.push_back(m_taus1[0]);
1111 ✗ param.push_back(m_taus2[0]);
1112 ✗ param.push_back(m_ks[0]);
1113 ✗ param.push_back(m_us[0]);
1114 ✗ param.push_back(m_tausi[0]);
1115 ✗ param.push_back(m_tauwinf[0]);
1116 ✗ param.push_back(m_winfstar[0]);
1117
1118 } else {
1119
1120 ✗ for(unsigned int i=0 ; i<m_posparam.size() ; i++) {
1121
1122 ✗ const int posp = m_posparam[i];
1123
1124 ✗ switch(posp) {
1125
1126 ✗ case 0:
1127 ✗ param.push_back(m_uo[0]);
1128 ✗ break;
1129
1130 ✗ case 1:
1131 ✗ param.push_back(m_uu[0]);
1132 ✗ break;
1133
1134 ✗ case 2:
1135 ✗ param.push_back(m_thetav[0]);
1136 ✗ break;
1137
1138 ✗ case 3:
1139 ✗ param.push_back(m_thetaw[0]);
1140 ✗ break;
1141
1142 ✗ case 4:
1143 ✗ param.push_back(m_thetavm[0]);
1144 ✗ break;
1145
1146 ✗ case 5:
1147 ✗ param.push_back(m_thetao[0]);
1148 ✗ break;
1149
1150 ✗ case 6:
1151 ✗ param.push_back(m_tauvm1[0]);
1152 ✗ break;
1153
1154 ✗ case 7:
1155 ✗ param.push_back(m_tauvm2[0]);
1156 ✗ break;
1157
1158 ✗ case 8:
1159 ✗ param.push_back(m_tauvp[0]);
1160 ✗ break;
1161
1162 ✗ case 9:
1163 ✗ param.push_back(m_tauwm1[0]);
1164 ✗ break;
1165
1166 ✗ case 10:
1167 ✗ param.push_back(m_tauwm2[0]);
1168 ✗ break;
1169
1170 ✗ case 11:
1171 ✗ param.push_back(m_kwm[0]);
1172 ✗ break;
1173
1174 ✗ case 12:
1175 ✗ param.push_back(m_uwm[0]);
1176 ✗ break;
1177
1178 ✗ case 13:
1179 ✗ param.push_back(m_tauwm[0]);
1180 ✗ break;
1181
1182 ✗ case 14:
1183 ✗ param.push_back(m_taufi[0]);
1184 ✗ break;
1185
1186 ✗ case 15:
1187 ✗ param.push_back(m_tauo1[0]);
1188 ✗ break;
1189
1190 ✗ case 16:
1191 ✗ param.push_back(m_tauo2[0]);
1192 ✗ break;
1193
1194 ✗ case 17:
1195 ✗ param.push_back(m_tauso1[0]);
1196 ✗ break;
1197
1198 ✗ case 18:
1199 ✗ param.push_back(m_tauso2[0]);
1200 ✗ break;
1201
1202 ✗ case 19:
1203 ✗ param.push_back(m_kso[0]);
1204 ✗ break;
1205
1206 ✗ case 20:
1207 ✗ param.push_back(m_uso[0]);
1208 ✗ break;
1209
1210 ✗ case 21:
1211 ✗ param.push_back(m_taus1[0]);
1212 ✗ break;
1213
1214 ✗ case 22:
1215 ✗ param.push_back(m_taus2[0]);
1216 ✗ break;
1217
1218 ✗ case 23:
1219 ✗ param.push_back(m_ks[0]);
1220 ✗ break;
1221
1222 ✗ case 24:
1223 ✗ param.push_back(m_us[0]);
1224 ✗ break;
1225
1226 ✗ case 25:
1227 ✗ param.push_back(m_tausi[0]);
1228 ✗ break;
1229
1230 ✗ case 26:
1231 ✗ param.push_back(m_tauwinf[0]);
1232 ✗ break;
1233
1234 ✗ case 27:
1235 ✗ param.push_back(m_winfstar[0]);
1236 ✗ break;
1237
1238 ✗ default:
1239 ✗ std::cout << "Numerotation error in setParameters of MVSolver" << std::endl;
1240 ✗ break;
1241
1242 }
1243
1244 }//End for
1245
1246 }//End else
1247
1248 ✗ p = param;
1249
1250 }
1251
1252
1253
1254 ✗ void MVSolver::getEDOvalues(std::vector<double>& v,std::vector<double>& w,std::vector<double>& s) {
1255
1256 ✗ v.resize(m_vEDO.size());
1257 ✗ w.resize(m_wEDO.size());
1258 ✗ s.resize(m_sEDO.size());
1259
1260 ✗ v = m_vEDO;
1261 ✗ w = m_wEDO;
1262 ✗ s = m_sEDO;
1263
1264 }
1265
1266 ✗ void MVSolver::getEDOpvalues(std::vector<double>& v,std::vector<double>& w,std::vector<double>& s) {
1267
1268 ✗ v.resize(p_vEDO.size());
1269 ✗ w.resize(p_wEDO.size());
1270 ✗ s.resize(p_sEDO.size());
1271
1272 ✗ v = p_vEDO;
1273 ✗ w = p_wEDO;
1274 ✗ s = p_sEDO;
1275
1276 }
1277
1278
1279 ✗ void MVSolver::setEDOvalues(std::vector<double> v,std::vector<double> w,std::vector<double> s) {
1280
1281 ✗ m_vEDO = v;
1282 ✗ m_wEDO = w;
1283 ✗ m_sEDO = s;
1284
1285 }
1286
1287 ✗ void MVSolver::getVmpvalues(std::vector<double>& Vm) {
1288
1289 ✗ Vm.resize(p_vEDO.size());
1290
1291 ✗ Vm = p_Vm;
1292
1293 }
1294
1295 ✗ void MVSolver::getVmvalues(std::vector<double>& Vm) {
1296
1297 ✗ Vm.resize(m_vEDO.size());
1298
1299 ✗ Vm = m_Vm;
1300
1301 }
1302
1303
1304 ✗ void MVSolver::setVmvalues(std::vector<double> Vm) {
1305
1306 ✗ m_Vm = Vm;
1307
1308 }
1309
1310 ✗ void MVSolver::setAssimilation(bool assim) {
1311
1312 ✗ m_assimilation = assim;
1313
1314 }
1315
1316
1317 ✗ void MVSolver::setPosParam(std::vector<int> pp) {
1318
1319 ✗ m_posparam.clear();
1320
1321 ✗ for(unsigned ipp=0 ; ipp<pp.size() ; ipp++) {
1322 ✗ m_posparam.push_back(pp[ipp]);
1323 }
1324
1325 }
1326
1327
1328
1329
1330 }
1331
1332
1333
1334