GCC Code Coverage Report

Directory: ./
File: Solver/courtemancheSolver.cpp
Date: 2024-04-14 07:32:34
Exec Total Coverage
Lines: 0 466 0.0%
Branches: 0 354 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: A. Collin
13 //
14
15 // System includes
16
17 // External includes
18
19 // Project includes
20 #include "Solver/courtemancheSolver.hpp"
21
22 namespace felisce {
23 ✗ CourtemancheSolver::CourtemancheSolver(FelisceTransient::Pointer fstransient):
24 ✗ m_fstransient(fstransient),
25 ✗ m_size(0),
26 ✗ m_aoPetsc(FELISCE_PETSC_NULLPTR){
27 }
28
29
30 ✗ CourtemancheSolver::~CourtemancheSolver() {
31 ✗ m_courtCondIto.clear();
32 ✗ m_uExtrap.destroy();
33 ✗ m_ion.destroy();
34 ✗ m_m.destroy();
35 ✗ m_h.destroy();
36 ✗ m_j.destroy();
37 ✗ m_ao.destroy();
38 ✗ m_io.destroy();
39 ✗ m_ua.destroy();
40 ✗ m_ui.destroy();
41 ✗ m_xr.destroy();
42 ✗ m_xs.destroy();
43 ✗ m_d.destroy();
44 ✗ m_f.destroy();
45 ✗ m_fca.destroy();
46 ✗ m_urel.destroy();
47 ✗ m_vrel.destroy();
48 ✗ m_wrel.destroy();
49 ✗ m_nai.destroy();
50 ✗ m_nao.destroy();
51 ✗ m_cao.destroy();
52 ✗ m_ki.destroy();
53 ✗ m_ko.destroy();
54 ✗ m_cai.destroy();
55 ✗ m_naiont.destroy();
56 ✗ m_kiont.destroy();
57 ✗ m_caiont.destroy();
58 ✗ m_ileak.destroy();
59 ✗ m_iup.destroy();
60 ✗ m_itr.destroy();
61 ✗ m_irel.destroy();
62 ✗ m_cmdn.destroy();
63 ✗ m_trpn.destroy();
64 ✗ m_nsr.destroy();
65 ✗ m_jsr.destroy();
66 ✗ m_csqn.destroy();
67 }
68
69
70 //Size of equation not necessarily V_0 size.
71 ✗ void CourtemancheSolver::defineSizeAndMappingOfIonicProblem(felInt size, ISLocalToGlobalMapping& mapping, AO ao) {
72 ✗ m_size = size;
73 ✗ m_localDofToGlobalDof = mapping;
74
75 ✗ m_aoPetsc = ao;
76 //Print the mapping.
77 //ISLocalToGlobalMappingView(m_localDofToGlobalDof, PETSC_VIEWER_STDOUT_WORLD);
78 }
79
80
81 ✗ void CourtemancheSolver::initializeExtrap(PetscVector& V_0) {
82 //Initialize u_extrapolate.
83 ✗ m_uExtrap.duplicateFrom(V_0);
84 ✗ m_uExtrap.copyFrom(V_0);
85
86 //Initialize ion.
87 ✗ m_ion.duplicateFrom(V_0);
88 ✗ m_ion.zeroEntries();
89
90 }
91
92
93 ✗ void CourtemancheSolver::updateExtrap(PetscVector& V_1) {
94 //Update u_extrapolate.
95 ✗ swap(V_1,m_uExtrap);
96 }
97
98 ✗ void CourtemancheSolver::initializeConcentrationsAndGateConditions(PetscVector& m,PetscVector& h, PetscVector& j, PetscVector& ao, PetscVector& io, PetscVector& ua, PetscVector& ui, PetscVector& xr, PetscVector& xs, PetscVector& d, PetscVector& f, PetscVector& fca, PetscVector& urel, PetscVector& vrel, PetscVector& wrel, PetscVector& nai, PetscVector& nao, PetscVector& cao, PetscVector& ki, PetscVector& ko, PetscVector& cai, PetscVector& naiont, PetscVector& kiont, PetscVector& caiont, PetscVector& ileak, PetscVector& iup, PetscVector& itr, PetscVector& irel, PetscVector& cmdn, PetscVector& trpn, PetscVector& nsr, PetscVector& jsr, PetscVector& csqn) {
99 ✗ m_m.duplicateFrom(m);
100 ✗ m_m.copyFrom(m);
101 ✗ m_h.duplicateFrom(h);
102 ✗ m_h.copyFrom(h);
103 ✗ m_j.duplicateFrom(j);
104 ✗ m_j.copyFrom(j);
105 ✗ m_ao.duplicateFrom(ao);
106 ✗ m_ao.copyFrom(ao);
107 ✗ m_io.duplicateFrom(io);
108 ✗ m_io.copyFrom(io);
109 ✗ m_ua.duplicateFrom(ua);
110 ✗ m_ua.copyFrom(ua);
111 ✗ m_ui.duplicateFrom(ui);
112 ✗ m_ui.copyFrom(ui);
113 ✗ m_xr.duplicateFrom(xr);
114 ✗ m_xr.copyFrom(xr);
115 ✗ m_xs.duplicateFrom(xs);
116 ✗ m_xs.copyFrom(xs);
117 ✗ m_d.duplicateFrom(d);
118 ✗ m_d.copyFrom(d);
119 ✗ m_f.duplicateFrom(f);
120 ✗ m_f.copyFrom(f);
121 ✗ m_fca.duplicateFrom(fca);
122 ✗ m_fca.copyFrom(fca);
123 ✗ m_urel.duplicateFrom(urel);
124 ✗ m_urel.copyFrom(urel);
125 ✗ m_vrel.duplicateFrom(vrel);
126 ✗ m_vrel.copyFrom(vrel);
127 ✗ m_wrel.duplicateFrom(wrel);
128 ✗ m_wrel.copyFrom(wrel);
129 ✗ m_nai.duplicateFrom(nai);
130 ✗ m_nai.copyFrom(nai);
131 ✗ m_nao.duplicateFrom(nao);
132 ✗ m_nao.copyFrom(nao);
133 ✗ m_cao.duplicateFrom(cao);
134 ✗ m_cao.copyFrom(cao);
135 ✗ m_ki.duplicateFrom(ki);
136 ✗ m_ki.copyFrom(ki);
137 ✗ m_ko.duplicateFrom(ko);
138 ✗ m_ko.copyFrom(ko);
139 ✗ m_cai.duplicateFrom(cai);
140 ✗ m_cai.copyFrom(cai);
141 ✗ m_naiont.duplicateFrom(naiont);
142 ✗ m_naiont.copyFrom(naiont);
143 ✗ m_kiont.duplicateFrom(kiont);
144 ✗ m_kiont.copyFrom(kiont);
145 ✗ m_caiont.duplicateFrom(caiont);
146 ✗ m_caiont.copyFrom(caiont);
147 ✗ m_ileak.duplicateFrom(ileak);
148 ✗ m_ileak.copyFrom(ileak);
149 ✗ m_iup.duplicateFrom(iup);
150 ✗ m_iup.copyFrom(iup);
151 ✗ m_itr.duplicateFrom(itr);
152 ✗ m_itr.copyFrom(itr);
153 ✗ m_irel.duplicateFrom(irel);
154 ✗ m_irel.copyFrom(irel);
155 ✗ m_cmdn.duplicateFrom(cmdn);
156 ✗ m_cmdn.copyFrom(cmdn);
157 ✗ m_trpn.duplicateFrom(trpn);
158 ✗ m_trpn.copyFrom(trpn);
159 ✗ m_nsr.duplicateFrom(nsr);
160 ✗ m_nsr.copyFrom(nsr);
161 ✗ m_jsr.duplicateFrom(jsr);
162 ✗ m_jsr.copyFrom(jsr);
163 ✗ m_csqn.duplicateFrom(csqn);
164 ✗ m_csqn.copyFrom(csqn);
165 }
166
167
168
169 ✗ double CourtemancheSolver::iIonTotal(double v, felInt pos) {
170
171 ✗ double iIonTot = 0.0;
172 ✗ double R = 8.3143;
173 ✗ double temp = 310.0;
174 ✗ double frdy = 96.4867;
175 ✗ double dt = m_fstransient->timeStep;
176
177
178 // Compute INa: ina
179 ✗ double ina = 0.0;
180 ✗ double gna = 7.8;
181 ✗ double ena = ((R*temp)/frdy)*std::log(m_value_nao/m_value_nai);
182 ✗ double am = 3.2;
183
184 ✗ if (std::abs(v - 47.13) < 1.0e-7) {
185 ✗ am = 0.32*(v+47.13)/(1.0-exp(-0.1*(v+47.13)));
186 }
187
188 ✗ double bm = 0.08*exp(-v/11.0);
189 double ah, bh, aj, bj;
190
191 ✗ if (v < -40.0) {
192 ✗ ah = 0.135*exp((80.0+v)/-6.8);
193 ✗ bh = 3.56*exp(0.079*v)+310000.0*exp(0.35*v);
194 ✗ aj = (-127140.0*exp(0.2444*v)-0.00003474*exp(-0.04391*v))*((v+37.78)/(1.0+exp(0.311*(v+79.23))));
195 ✗ bj = (0.1212*exp(-0.01052*v))/(1.0+exp(-0.1378*(v+40.14)));
196 } else {
197 ✗ ah = 0.0;
198 ✗ bh = 1.0/(0.13*(1.0+exp((v+10.66)/-11.1)));
199 ✗ aj = 0.0;
200 ✗ bj = (0.3*exp(-0.0000002535*v))/(1.0+exp(-0.1*(v+32.0)));
201 }
202
203 ✗ m_value_h = ah/(ah+bh)-((ah/(ah+bh))-m_value_h)*exp(-dt/(1.0/(ah+bh)));
204 ✗ m_value_j = aj/(aj+bj)-((aj/(aj+bj))-m_value_j)*exp(-dt/(1.0/(aj+bj)));
205 ✗ m_value_m = am/(am+bm)-((am/(am+bm))-m_value_m)*exp(-dt/(1.0/(am+bm)));
206
207 ✗ if (!FelisceParam::instance().hasHeteroCondAtria) {
208 ✗ ina = gna*m_value_m*m_value_m*m_value_m*m_value_h*m_value_j*(v-ena)*m_courtCondMultCoeff[pos];
209 } else {
210 ✗ ina = gna*m_value_m*m_value_m*m_value_m*m_value_h*m_value_j*(v-ena);
211 }
212
213 // Compute ICaL: ical
214 ✗ double ical = 0.0;
215 ✗ double dss = 1.0/(1.0+exp(-(v+10.0)/8.0));
216 ✗ double taud = (1.0-exp((v+10.0)/-6.24))/(0.035*(v+10.0)*(1.0+exp((v+10.0)/-6.24)));
217 ✗ double fss = 1.0/(1.0+exp((v+28.0)/6.9));
218 ✗ double tauf = 9.0/(0.0197*exp(-std::pow((0.0337*(v+10.0)),2))+0.02);
219 ✗ double fcass = 1.0/(1.0+m_value_cai/0.00035);
220 ✗ double taufca = 2.0;
221
222 ✗ m_value_d = dss-(dss-m_value_d)*exp(-dt/taud);
223 ✗ m_value_f = fss-(fss-m_value_f)*exp(-dt/tauf);
224 ✗ m_value_fca = fcass-(fcass-m_value_fca)*exp(-dt/taufca);
225
226 ✗ double gcalbar = 0.1238;
227 ✗ if (FelisceParam::instance().hasHeteroCourtPar) {
228 ✗ gcalbar = m_courtCondICaL[pos];
229 }
230
231 ✗ ical = gcalbar*m_value_d*m_value_f*m_value_fca*(v-65.0);
232
233
234 // Compute IKr: ikr
235 ✗ double ikr = 0.0;
236
237 ✗ double gkr = 0.0294;//*std::sqrt(ko/5.4);
238 ✗ double ekr = ((R*temp)/frdy)*std::log(m_value_ko/m_value_ki);
239
240 ✗ double xrss = 1.0/(1.0+exp(-(v+14.1)/6.5));
241 ✗ double tauxr = 1.0/(0.0003*(v+14.1)/(1.0-exp(-(v+14.1)/5.0))+0.000073898*(v-3.3328)/(exp((v-3.3328)/5.1237)-1.0));
242
243 ✗ m_value_xr = xrss-(xrss-m_value_xr)*exp(-dt/tauxr);
244
245 ✗ double r = 1.0/(1.0+exp((v+15.0)/22.4));
246
247 ✗ ikr = gkr*m_value_xr*r*(v-ekr);
248
249
250 // Compute IKs: iks
251 ✗ double iks = 0.0;
252
253 ✗ double gks = 0.129;
254 ✗ double eks = ((R*temp)/frdy)*std::log(m_value_ko/m_value_ki);
255 ✗ double axs = 4.0e-5*(v-19.9)/(1.0-exp(-(v-19.9)/17.0));
256 ✗ double bxs = 3.5e-5*(v-19.9)/(exp((v-19.9)/9.0)-1.0);
257 ✗ double tauxs = 1.0/(2.0*(axs+bxs));
258 ✗ double xsss = 1.0/std::pow((1.0+exp(-(v-19.9)/12.7)),0.5);
259 ✗ m_value_xs = xsss-(xsss-m_value_xs)*exp(-dt/tauxs);
260
261 //Warning ! The parameter of the CRN model gks becomes 5.0*gks in order to accelerate the atrial repolarization.
262 ✗ if (FelisceParam::instance().torsade) {
263 ✗ if ( (m_fstransient->time > FelisceParam::instance().torsadeTimeBegin) && (m_fstransient->time < FelisceParam::instance().torsadeTimeEnd) ) {
264 ✗ iks = gks*m_value_xs*m_value_xs*(v-eks);
265 }
266 } else {
267 ✗ iks = 5./*3.0*/*gks*m_value_xs*m_value_xs*(v-eks);
268 }
269
270 // Compute IKi: iki
271 ✗ double iki = 0.0;
272
273 ✗ double gki = 0.09;//*std::pow(ko/5.4,0.4);
274 ✗ double eki = ((R*temp)/frdy)*std::log(m_value_ko/m_value_ki);
275
276 ✗ double kin = 1.0/(1.0+exp(0.07*(v+80.0)));
277
278 ✗ iki = gki*kin*(v-eki);
279
280
281 // Compute Ikur: ikur
282 ✗ double ikur = 0.0;
283
284 ✗ double gkur = 0.005+0.05/(1+exp(-(v-15.0)/13.0));
285 ✗ double ekur = ((R*temp)/frdy)*std::log(m_value_ko/m_value_ki);
286 ✗ double alphauakur = 0.65/(exp(-(v+10.0)/8.5)+exp(-(v-30.0)/59.0));
287 ✗ double betauakur = 0.65/(2.5+exp((v+82.0)/17.0));
288 ✗ double tauuakur = 1.0/(3.0*(alphauakur+betauakur));
289 ✗ double uakurss = 1.0/(1.0+exp(-(v+30.3)/9.6));
290 ✗ double alphauikur = 1.0/(21.0+exp(-(v-185.0)/28.0));
291 ✗ double betauikur = exp((v-158.0)/16.0);
292 ✗ double tauuikur = 1.0/(3.0*(alphauikur+betauikur));
293 ✗ double uikurss = 1.0/(1.0+exp((v-99.45)/27.48));
294
295 ✗ m_value_ua = uakurss-(uakurss-m_value_ua)*exp(-dt/tauuakur);
296 ✗ m_value_ui = uikurss-(uikurss-m_value_ui)*exp(-dt/tauuikur);
297
298 ✗ ikur = gkur*m_value_ua*m_value_ua*m_value_ua*m_value_ui*(v-ekur);
299
300
301 // Compute ITo: ito
302 ✗ double ito = 0.0;
303 ✗ double gito = 0.1652;
304 ✗ if (FelisceParam::instance().hasHeteroCourtPar) {
305 ✗ gito = m_courtCondIto[pos];
306 }
307
308 ✗ double erevto = ((R*temp)/frdy)*std::log(m_value_ko/m_value_ki);
309
310 ✗ double alphaato = 0.65/(exp(-(v+10.0)/8.5)+exp(-(v-30.0)/59.0));
311 ✗ double betaato = 0.65/(2.5+exp((v+82.0)/17.0));
312 ✗ double tauato = 1.0/(3.0*(alphaato+betaato));
313 ✗ double atoss = 1.0/(1.0+exp(-(v+20.47)/17.54));
314 ✗ m_valuem_ao = atoss-(atoss-m_valuem_ao)*exp(-dt/tauato);
315
316 ✗ double alphaiito = 1.0/(18.53+exp((v+113.7)/10.95));
317 ✗ double betaiito = 1.0/(35.56+exp(-(v+1.26)/7.44));
318 ✗ double tauiito = 1.0/(3.0*(alphaiito+betaiito));
319 ✗ double iitoss = 1.0/(1.0+exp((v+43.1)/5.3));
320 ✗ m_value_io = iitoss-(iitoss-m_value_io)*exp(-dt/tauiito);
321
322 ✗ ito = gito*m_valuem_ao*m_valuem_ao*m_valuem_ao*m_value_io*(v-erevto);
323
324
325 // Compute INaCa: inaca
326 ✗ double inaca = 0.0;
327
328 ✗ double gammas = 0.35;
329 ✗ double kmnancx = 87.5;
330 ✗ double kmcancx = 1.38;
331 ✗ double ksatncx = 0.1;
332
333 ✗ inaca = 1600.0*(exp(gammas*frdy*v/(R*temp))*m_value_nai*m_value_nai*m_value_nai*m_value_cao-exp((gammas-1.0)*frdy*v/(R*temp))*m_value_nao*m_value_nao*m_value_nao*m_value_cai)/((std::pow(kmnancx,3)+std::pow(m_value_nao,3))*(kmcancx+m_value_cao)*(1.0+ksatncx*exp((gammas-1.0)*frdy*v/(R*temp))));
334
335
336 // Compute INaK: inak
337 ✗ double inak = 0.0;
338 ✗ double sigma = (exp(m_value_nao/67.3)-1.0)/7.0;
339 ✗ double ibarnak = 0.6;
340 ✗ double kmko = 1.5;
341 ✗ double kmnai = 10.0;
342
343 ✗ double fnak = 1.0/(1.0+0.1245*exp((-0.1*v*frdy)/(R*temp))+0.0365*sigma*exp((-v*frdy)/(R*temp)));
344
345 ✗ inak = ibarnak*fnak*(1.0/(1.0+std::pow((kmnai/m_value_nai),1.5)))*(m_value_ko/(m_value_ko+kmko));
346
347
348 // Compute IPCa: ipca
349 ✗ double ipca = 0.0;
350 ✗ double ibarpca = 0.275;
351 ✗ double kmpca = 0.0005;
352
353 ✗ ipca = (ibarpca*m_value_cai)/(kmpca+m_value_cai);
354
355
356 // Compute IbCa: ibca
357 ✗ double ibca = 0.0;
358 ✗ double gcab = 0.00113;
359 ✗ double ecan = ((R*temp)/frdy/2.0)*std::log(m_value_cao/m_value_cai);
360
361 ✗ ibca = gcab*(v-ecan);
362
363
364 // Compute IbNa: ibna
365 ✗ double ibna = 0.0;
366 ✗ double gnab = 0.000674;
367 ✗ double enan = ((R*temp)/frdy)*std::log(m_value_nao/m_value_nai);
368
369 ✗ ibna = gnab*(v-enan);
370
371
372 // Compute ILeak: ileak
373 ✗ m_value_ileak = 0.0;
374 ✗ double caupmax = 15.0;
375 ✗ double iupmax = 0.005;
376
377 ✗ m_value_ileak = m_value_nsr*iupmax/caupmax;
378
379
380 // Compute Iup: iup
381 ✗ m_value_iup = 0.0;
382 ✗ double kup = 0.00092;
383 ✗ m_value_iup = iupmax/(1.0+(kup/m_value_cai)) ;
384
385
386 // Compute Itr: itr
387 ✗ m_value_itr = 0.0;
388 ✗ double tautr = 180.0;
389 ✗ m_value_itr = (m_value_nsr-m_value_jsr)/tautr;
390
391
392 // Compute Irel: irel
393 ✗ m_value_irel = 0.0;
394 ✗ double krel = 30.0;
395 ✗ double constvrel = 96.48;
396 ✗ double fn = constvrel*(1.0e-12)*m_value_irel-(5.0e-13)*(ical/2.0-inaca/5.0)/frdy;
397
398 ✗ double tauurel = 8.0;
399 ✗ double urelss = 1.0/(1.0+exp(-(fn-3.4175e-13)/13.67e-16));
400
401 ✗ double tauvrel = 1.91+2.09/(1.0+exp(-(fn-3.4175e-13)/13.67e-16));
402 ✗ double vrelss = 1.0-1.0/(1.0+exp(-(fn-6.835e-14)/13.67e-16));
403
404 ✗ double tauwrel = 6.0*(1.0-exp(-(v-7.9)/5.0))/((1.0+0.3*exp(-(v-7.9)/5.0))*(v-7.9));
405 ✗ double wrelss = 1.0-1.0/(1.0+exp(-(v-40.0)/17.0));
406
407 ✗ m_value_urel = urelss-(urelss-m_value_urel)*exp(-dt/tauurel);
408 ✗ m_value_vrel = vrelss-(vrelss-m_value_vrel)*exp(-dt/tauvrel);
409 ✗ m_value_wrel = wrelss-(wrelss-m_value_wrel)*exp(-dt/tauwrel);
410
411 ✗ m_value_irel = krel*m_value_urel*m_value_urel*m_value_vrel*m_value_wrel*(m_value_jsr-m_value_cai);
412
413 // Compute total current for Na:
414 ✗ m_value_naiont = ina+ibna+3.0*inak+3.0*inaca;
415 // Compute total current for K:
416 ✗ m_value_kiont = ikr+iks+iki-2.0*inak+ito+ikur;
417 // Compute total current for Ca:
418 ✗ m_value_caiont = ical+ibca+ipca-2.0*inaca;
419
420 ✗ iIonTot = ical + ibca + ibna + inaca + ipca + inak + ikr + iks + ikur + ito + iki + ina;
421
422 ✗ return iIonTot;
423 }
424
425 ✗ void CourtemancheSolver::updateConcentrations() {
426 ✗ double dt = m_fstransient->timeStep;
427 ✗ double frdy = 96.4867;
428 ✗ double vi = 13668.0;
429
430
431 // Update Na concentration
432 ✗ double dnai = -dt*m_value_naiont/(vi*frdy);
433 ✗ m_value_nai = dnai + m_value_nai;
434
435 // Update K concentration
436 ✗ double dki = -dt*m_value_kiont/(vi*frdy);
437 ✗ m_value_ki = dki + m_value_ki;
438
439 // Update Ca concentration
440 ✗ double vup = 1109.52;
441 ✗ double constvrel = 96.48;
442
443 ✗ double b1cai = -m_value_caiont/(2.0*frdy*vi)+(vup*(m_value_ileak-m_value_iup)+constvrel*m_value_irel)/vi;
444
445 ✗ double cmdnmax = 0.05;
446 ✗ double trpnmax = 0.07;
447 ✗ double kmcmdn = 0.00238;
448 ✗ double kmtrpn = 0.0005;
449 ✗ double b2cai = 1.0+trpnmax*kmtrpn/std::pow((m_value_cai+kmtrpn),2)+cmdnmax*kmcmdn/std::pow((m_value_cai+kmcmdn),2);
450
451 ✗ double dcai = dt*b1cai/b2cai;
452
453 ✗ m_value_cai = dcai+m_value_cai;
454 }
455
456 ✗ PetscVector& CourtemancheSolver::ionicVariable(int iVar) {
457 ✗ switch (iVar+1) {
458 ✗ case 1:
459 ✗ return m_m;
460 ✗ case 2:
461 ✗ return m_h;
462 ✗ case 3:
463 ✗ return m_j;
464 ✗ case 4:
465 ✗ return m_ao;
466 ✗ case 5:
467 ✗ return m_io;
468 ✗ case 6:
469 ✗ return m_ua;
470 ✗ case 7:
471 ✗ return m_ui;
472 ✗ case 8:
473 ✗ return m_xr;
474 ✗ case 9:
475 ✗ return m_xs;
476 ✗ case 10:
477 ✗ return m_d;
478 ✗ case 11:
479 ✗ return m_f;
480 ✗ case 12:
481 ✗ return m_fca;
482 ✗ case 13:
483 ✗ return m_urel;
484 ✗ case 14:
485 ✗ return m_vrel;
486 ✗ case 15:
487 ✗ return m_wrel;
488 ✗ case 16:
489 ✗ return m_nai;
490 ✗ case 17:
491 ✗ return m_nao;
492 ✗ case 18:
493 ✗ return m_cao;
494 ✗ case 19:
495 ✗ return m_ki;
496 ✗ case 20:
497 ✗ return m_ko;
498 ✗ case 21:
499 ✗ return m_cai;
500 ✗ case 22:
501 ✗ return m_naiont;
502 ✗ case 23:
503 ✗ return m_kiont;
504 ✗ case 24:
505 ✗ return m_caiont;
506 ✗ case 25:
507 ✗ return m_cmdn;
508 ✗ case 26:
509 ✗ return m_trpn;
510 ✗ case 27:
511 ✗ return m_nsr;
512 ✗ case 28:
513 ✗ return m_jsr;
514 ✗ case 29:
515 ✗ return m_csqn;
516 ✗ default:
517 ✗ FEL_ERROR("Courtemanche ionic variable not defined");
518 ✗ return m_csqn;
519 }
520 }
521
522
523
524 ✗ void CourtemancheSolver::computeIon() {
525
526 felInt pos;
527 double value_uExtrap;
528 double value_ion;
529
530 ✗ for (felInt i = 0; i < m_size; i++) {
531 ✗ ISLocalToGlobalMappingApply(m_localDofToGlobalDof,1,&i,&pos);
532 // value_uExtrap = m_uExtrap(i);
533 ✗ m_uExtrap.getValues(1,&pos,&value_uExtrap);
534
535 ✗ m_m.getValues(1,&pos,&m_value_m);
536 ✗ m_h.getValues(1,&pos,&m_value_h);
537 ✗ m_j.getValues(1,&pos,&m_value_j);
538 ✗ m_ao.getValues(1,&pos,&m_valuem_ao);
539 ✗ m_io.getValues(1,&pos,&m_value_io);
540 ✗ m_ua.getValues(1,&pos,&m_value_ua);
541 ✗ m_ui.getValues(1,&pos,&m_value_ui);
542 ✗ m_xr.getValues(1,&pos,&m_value_xr);
543 ✗ m_xs.getValues(1,&pos,&m_value_xs);
544 ✗ m_d.getValues(1,&pos,&m_value_d);
545 ✗ m_f.getValues(1,&pos,&m_value_f);
546 ✗ m_fca.getValues(1,&pos,&m_value_fca);
547 ✗ m_urel.getValues(1,&pos,&m_value_urel);
548 ✗ m_vrel.getValues(1,&pos,&m_value_vrel);
549 ✗ m_wrel.getValues(1,&pos,&m_value_wrel);
550 ✗ m_nai.getValues(1,&pos,&m_value_nai);
551 ✗ m_nao.getValues(1,&pos,&m_value_nao);
552 ✗ m_cao.getValues(1,&pos,&m_value_cao);
553 ✗ m_ki.getValues(1,&pos,&m_value_ki);
554 ✗ m_ko.getValues(1,&pos,&m_value_ko);
555 ✗ m_cai.getValues(1,&pos,&m_value_cai);
556 ✗ m_naiont.getValues(1,&pos,&m_value_naiont);
557 ✗ m_kiont.getValues(1,&pos,&m_value_kiont);
558 ✗ m_caiont.getValues(1,&pos,&m_value_caiont);
559 ✗ m_ileak.getValues(1,&pos,&m_value_ileak);
560 ✗ m_iup.getValues(1,&pos,&m_value_iup);
561 ✗ m_itr.getValues(1,&pos,&m_value_itr);
562 ✗ m_irel.getValues(1,&pos,&m_value_irel);
563 ✗ m_cmdn.getValues(1,&pos,&m_value_cmdn);
564 ✗ m_trpn.getValues(1,&pos,&m_value_trpn);
565 ✗ m_nsr.getValues(1,&pos,&m_value_nsr);
566 ✗ m_jsr.getValues(1,&pos,&m_value_jsr);
567 ✗ m_csqn.getValues(1,&pos,&m_value_csqn);
568
569 ✗ value_ion = iIonTotal(value_uExtrap,pos);
570
571 ✗ felInt meshId = pos;
572 ✗ AOPetscToApplication(m_aoPetsc,1,&meshId);
573 ✗ felInt sizeVent = 45580;
574 ✗ if (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart")
575 ✗ sizeVent = 32320;
576 ✗ if ( (FelisceParam::instance().typeOfIonicModel == "courtAtriaSchafVent") || (FelisceParam::instance().typeOfIonicModel == "courtAtriaMvVent") ) {
577 ✗ if (meshId < sizeVent)
578 ✗ value_ion = 0.0;
579 }
580
581 ✗ updateConcentrations();
582
583
584 ✗ m_ion.setValue(pos,value_ion, INSERT_VALUES);
585 ✗ m_m.setValue(pos,m_value_m, INSERT_VALUES);
586 ✗ m_h.setValue(pos,m_value_h, INSERT_VALUES);
587 ✗ m_j.setValue(pos,m_value_j, INSERT_VALUES);
588 ✗ m_ao.setValue(pos,m_valuem_ao, INSERT_VALUES);
589 ✗ m_io.setValue(pos,m_value_io, INSERT_VALUES);
590 ✗ m_ua.setValue(pos,m_value_ua, INSERT_VALUES);
591 ✗ m_ui.setValue(pos,m_value_ui, INSERT_VALUES);
592 ✗ m_xr.setValue(pos,m_value_xr, INSERT_VALUES);
593 ✗ m_xs.setValue(pos,m_value_xs, INSERT_VALUES);
594 ✗ m_d.setValue(pos,m_value_d, INSERT_VALUES);
595 ✗ m_f.setValue(pos,m_value_f, INSERT_VALUES);
596 ✗ m_fca.setValue(pos,m_value_fca, INSERT_VALUES);
597 ✗ m_urel.setValue(pos,m_value_urel, INSERT_VALUES);
598 ✗ m_vrel.setValue(pos,m_value_vrel, INSERT_VALUES);
599 ✗ m_wrel.setValue(pos,m_value_wrel, INSERT_VALUES);
600 ✗ m_nai.setValue(pos,m_value_nai, INSERT_VALUES);
601 ✗ m_nao.setValue(pos,m_value_nao, INSERT_VALUES);
602 ✗ m_cao.setValue(pos,m_value_cao, INSERT_VALUES);
603 ✗ m_ki.setValue(pos,m_value_ki, INSERT_VALUES);
604 ✗ m_ko.setValue(pos,m_value_ko, INSERT_VALUES);
605 ✗ m_cai.setValue(pos,m_value_cai, INSERT_VALUES);
606 ✗ m_naiont.setValue(pos,m_value_naiont, INSERT_VALUES);
607 ✗ m_kiont.setValue(pos,m_value_kiont, INSERT_VALUES);
608 ✗ m_caiont.setValue(pos,m_value_caiont, INSERT_VALUES);
609 ✗ m_ileak.setValue(pos,m_value_ileak, INSERT_VALUES);
610 ✗ m_iup.setValue(pos,m_value_iup, INSERT_VALUES);
611 ✗ m_itr.setValue(pos,m_value_itr, INSERT_VALUES);
612 ✗ m_irel.setValue(pos,m_value_irel, INSERT_VALUES);
613 ✗ m_cmdn.setValue(pos,m_value_cmdn, INSERT_VALUES);
614 ✗ m_trpn.setValue(pos,m_value_trpn, INSERT_VALUES);
615 ✗ m_nsr.setValue(pos,m_value_nsr, INSERT_VALUES);
616 ✗ m_jsr.setValue(pos,m_value_jsr, INSERT_VALUES);
617 ✗ m_csqn.setValue(pos,m_value_csqn, INSERT_VALUES);
618
619
620 }
621
622 ✗ m_ion.assembly();
623 ✗ m_m.assembly();
624 ✗ m_h.assembly();
625 ✗ m_j.assembly();
626 ✗ m_ao.assembly();
627 ✗ m_io.assembly();
628 ✗ m_ua.assembly();
629 ✗ m_ui.assembly();
630 ✗ m_xr.assembly();
631 ✗ m_xs.assembly();
632 ✗ m_d.assembly();
633 ✗ m_f.assembly();
634 ✗ m_fca.assembly();
635 ✗ m_urel.assembly();
636 ✗ m_vrel.assembly();
637 ✗ m_wrel.assembly();
638 ✗ m_nai.assembly();
639 ✗ m_nao.assembly();
640 ✗ m_cao.assembly();
641 ✗ m_ki.assembly();
642 ✗ m_ko.assembly();
643 ✗ m_cai.assembly();
644 ✗ m_naiont.assembly();
645 ✗ m_kiont.assembly();
646 ✗ m_caiont.assembly();
647 ✗ m_ileak.assembly();
648 ✗ m_iup.assembly();
649 ✗ m_itr.assembly();
650 ✗ m_irel.assembly();
651 ✗ m_cmdn.assembly();
652 ✗ m_trpn.assembly();
653 ✗ m_nsr.assembly();
654 ✗ m_jsr.assembly();
655 ✗ m_csqn.assembly();
656
657 }
658
659
660 }
661