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 "Core/felisceTools.hpp" |
21 |
|
|
#include "Solver/cardiacFunction.hpp" |
22 |
|
|
|
23 |
|
|
namespace felisce { |
24 |
|
✗ |
CardiacFunction::CardiacFunction() { |
25 |
|
✗ |
m_fstransient = nullptr; |
26 |
|
|
} |
27 |
|
|
|
28 |
|
✗ |
CardiacFunction::~CardiacFunction() { |
29 |
|
✗ |
m_fstransient = nullptr; |
30 |
|
|
} |
31 |
|
|
|
32 |
|
✗ |
void CardiacFunction::initialize(FelisceTransient::Pointer fstransient) { |
33 |
|
✗ |
m_fstransient = fstransient; |
34 |
|
|
} |
35 |
|
|
|
36 |
|
|
//A class to impose thorso/lungs/bone conductivity. |
37 |
|
✗ |
ThoraxConductivity::ThoraxConductivity() = default; |
38 |
|
|
|
39 |
|
✗ |
ThoraxConductivity::~ThoraxConductivity() = default; |
40 |
|
|
|
41 |
|
✗ |
double ThoraxConductivity::sigmaT(int& label) const { |
42 |
|
|
double conductivity; |
43 |
|
|
|
44 |
|
|
// Thorax conuductivity (mesh region with label == 1) |
45 |
|
|
//if (label == 1){ // Old Zygote |
46 |
|
|
//if (label == 303){ // New systole heart |
47 |
|
✗ |
if (label == 101) { // New end-systole heart && ZygoteOldFinest |
48 |
|
✗ |
conductivity = FelisceParam::instance().sigmaThorax; |
49 |
|
|
} |
50 |
|
|
// Bones conuductivity (mesh region with label == 2) |
51 |
|
|
//else if (label == 2){ // Old Zygote |
52 |
|
|
//else if (label == 302){ // New systole heart |
53 |
|
✗ |
else if (label == 103) { // New end-systole heart && ZygoteOldFinest |
54 |
|
✗ |
conductivity = FelisceParam::instance().sigmaBone; |
55 |
|
|
} |
56 |
|
|
// Lungs conuductivity (mesh region with label == 3) |
57 |
|
|
//else if (label == 3){ // Old Zygote |
58 |
|
|
//else if (label == 301){ // New systole heart |
59 |
|
✗ |
else if (label == 102) { // New end-systole heart && ZygoteOldFinest |
60 |
|
✗ |
conductivity = FelisceParam::instance().sigmaLung; |
61 |
|
|
} else { |
62 |
|
✗ |
if (FelisceParam::verbose() > 10) std::cout << "Warning ! Label " << label << " not found !\n"; |
63 |
|
✗ |
conductivity = FelisceParam::instance().sigmaThorax; |
64 |
|
|
} |
65 |
|
|
|
66 |
|
|
|
67 |
|
✗ |
return conductivity; |
68 |
|
|
} |
69 |
|
|
|
70 |
|
✗ |
double ThoraxConductivity::operator() (int& label) const { |
71 |
|
✗ |
return sigmaT(label); |
72 |
|
|
} |
73 |
|
|
|
74 |
|
|
//A class to impose thorso/lungs/bone conductivity. |
75 |
|
✗ |
HeartConductivity::HeartConductivity() = default; |
76 |
|
|
|
77 |
|
✗ |
HeartConductivity::~HeartConductivity() = default; |
78 |
|
|
|
79 |
|
✗ |
double HeartConductivity::sigmaH(int& label) const { |
80 |
|
|
double conductivity; |
81 |
|
|
|
82 |
|
✗ |
if (label == 2) { |
83 |
|
✗ |
conductivity = 0.0; |
84 |
|
✗ |
} else if (label == 1) { |
85 |
|
✗ |
conductivity = FelisceParam::instance().intraTransvTensor; |
86 |
|
|
} else { |
87 |
|
✗ |
std::cout << "Warning ! Label " << label << " not found !\n"; |
88 |
|
✗ |
conductivity = 0.0; |
89 |
|
|
} |
90 |
|
|
|
91 |
|
✗ |
return conductivity; |
92 |
|
|
} |
93 |
|
|
|
94 |
|
✗ |
double HeartConductivity::operator() (int& label) const { |
95 |
|
✗ |
return sigmaH(label); |
96 |
|
|
} |
97 |
|
|
|
98 |
|
|
//A class to impose a non-pathological applied current to the ellibi-geometry. |
99 |
|
✗ |
AppCurrent::AppCurrent() = default; |
100 |
|
|
|
101 |
|
✗ |
AppCurrent::~AppCurrent() = default; |
102 |
|
|
|
103 |
|
✗ |
double AppCurrent::Iapp(double x, double y, double z, double t, double) const { |
104 |
|
✗ |
double iapp=0.0; |
105 |
|
✗ |
double delay = FelisceParam::instance().delayStim[0]; |
106 |
|
✗ |
double tPeriod=fmod(t,FelisceParam::instance().timePeriodVentricle); |
107 |
|
|
|
108 |
|
|
//current in ventricles: |
109 |
|
✗ |
double pi = acos(-1.0); |
110 |
|
|
//left : left ventricle volume = { (x^2/4^2 + y^2/4^2 + z^2/7.2^2) < 1 } && { (x^2/(4-1.2)^2 + y^2/(4-1.2)^2 + z^2/(7.2-1.2)^2) > 1 } && { z < 2.76 } |
111 |
|
✗ |
double a_L = 4.0; |
112 |
|
✗ |
double b_L = 4.0; |
113 |
|
✗ |
double c_L = 7.2; |
114 |
|
✗ |
double e_L1 = 1.6*0.8; |
115 |
|
✗ |
double e_L2 = 0.0; |
116 |
|
✗ |
double sumL1 = x*x/((a_L-e_L1)*(a_L-e_L1))+y*y/((b_L-e_L1)*(b_L-e_L1))+z*z/((c_L-e_L1)*(c_L-e_L1))-1.0; |
117 |
|
✗ |
double sumL2 = x*x/((a_L-e_L2)*(a_L-e_L2))+y*y/((b_L-e_L2)*(b_L-e_L2))+z*z/((c_L-e_L2)*(c_L-e_L2))-1.0; |
118 |
|
✗ |
double angularVelocityLV = (pi/2.0)/FelisceParam::instance().stimTimeLV; |
119 |
|
✗ |
double alphaLVmax = std::min((tPeriod-delay)*angularVelocityLV,2*pi); |
120 |
|
✗ |
double alphaLV = std::atan((x+1.5)/(2.7586-z)); |
121 |
|
|
|
122 |
|
|
//right : right ventricle volume = { (x^2/7.95^2 + y^2/4^2 + z^2/6.84^2) < 1 } && { (x^2/(7.95-0.7)^2 + y^2/(4-0.7)^2 + z^2/(6.84-0.7)^2) > 1 } && { z < 2.76 } |
123 |
|
✗ |
double a_R = 7.8; |
124 |
|
✗ |
double b_R = 4.0; |
125 |
|
✗ |
double c_R = 7.2*0.95; //=6.84 |
126 |
|
✗ |
double e_R = 1.0*0.8; |
127 |
|
✗ |
double sumR = x*x/((a_R-e_R)*(a_R-e_R))+y*y/((b_R-e_R)*(b_R-e_R))+z*z/((c_R-e_R)*(c_R-e_R))-1.0; |
128 |
|
✗ |
double angularVelocityRV = (pi/2.0)/FelisceParam::instance().stimTimeRV ; |
129 |
|
✗ |
double alphaRVmax = std::min((tPeriod-delay)*angularVelocityRV,2*pi); |
130 |
|
✗ |
double alphaRV = std::atan((-3.7-x)/(2.9586-z)); |
131 |
|
|
|
132 |
|
|
//left |
133 |
|
✗ |
if ((tPeriod>=delay)&&(tPeriod<=(delay+FelisceParam::instance().stimTimeLV))) { //if t \in [delay,delay+t_{actLeft}] |
134 |
|
✗ |
if (sumL1<0.15) { //if (x,y,z) \in S (endocardium of left ventricle) |
135 |
|
✗ |
if (alphaLV < alphaLVmax) { //if std::atan((x-x0)/(z-z0)) < alpha(t) |
136 |
|
✗ |
iapp = (-5.0/2.0*sumL1+0.375)*8.0/200.0; |
137 |
|
|
} |
138 |
|
|
} |
139 |
|
|
} |
140 |
|
|
|
141 |
|
|
//septum (right) |
142 |
|
✗ |
if ((tPeriod>=delay)&&(tPeriod<=(delay+ FelisceParam::instance().stimTimeRV))) { //if t \in [delay,delay+t_{actRight}] |
143 |
|
✗ |
if((sumL2>=-0.2)&&(sumL2<=0.0)&&(x<-2.4)) { //if (x,y,z) \in S (right part of septum) |
144 |
|
✗ |
iapp= (5.0/2.0*sumL2+0.5)*8.0/200.0; |
145 |
|
|
} |
146 |
|
|
} |
147 |
|
|
|
148 |
|
|
//right |
149 |
|
✗ |
if (tPeriod>=delay&&(tPeriod<=(delay+FelisceParam::instance().stimTimeRV))) { //if t \in [delay,delay+t_{actRight}] |
150 |
|
✗ |
if ((x<-2.586)&&(sumL2>-0.2)&&(sumR>=0.0)&&(sumR<0.1)) { //if (x,y,z) \in S (right part of septum or endocardium of right ventricle) |
151 |
|
✗ |
if (alphaRV < alphaRVmax) { //if std::atan((x-x0)/(z-z0)) < alpha(t) |
152 |
|
✗ |
iapp = 0.5*8.0/200.0; |
153 |
|
|
} |
154 |
|
|
} |
155 |
|
|
} |
156 |
|
|
|
157 |
|
✗ |
return iapp; |
158 |
|
|
|
159 |
|
|
} |
160 |
|
|
|
161 |
|
✗ |
double AppCurrent::operator() (double x, double y, double z, double t, double) const { |
162 |
|
✗ |
return Iapp(x,y,z,t); |
163 |
|
|
} |
164 |
|
|
|
165 |
|
|
//A class to impose a non-pathological applied current to the zygote-geometry. |
166 |
|
✗ |
zygoteIapp::zygoteIapp(): |
167 |
|
✗ |
AppCurrent() |
168 |
|
|
{} |
169 |
|
|
|
170 |
|
✗ |
double zygoteIapp::Iapp(double x, double y, double z, double t, double dist) const { |
171 |
|
|
|
172 |
|
✗ |
double iapp = 0.0; |
173 |
|
✗ |
double delay = FelisceParam::instance().delayStim[0]; |
174 |
|
|
|
175 |
|
|
|
176 |
|
|
//APEX, node number: 80185 : 3.316424 31.496351 5.799850 (coordinates of "Nejib heart") |
177 |
|
✗ |
double x_0= 3.316424; |
178 |
|
✗ |
double y_0 = 31.496351; |
179 |
|
✗ |
double z_0 = 5.799850; |
180 |
|
|
|
181 |
|
|
//APEX, node number: 80185 : 2.74681 28.9865 6.20925 (coordinates of "?? heart") |
182 |
|
|
/*double x_0= 2.74681; |
183 |
|
|
double y_0 = 28.9865; |
184 |
|
|
double z_0 = 6.20925; |
185 |
|
|
*/ |
186 |
|
|
//APEX, node number: 80185 : 2.74681 28.9865 6.20925 (coordinates of "?? heart") |
187 |
|
|
/*double x_0= 3.01795; |
188 |
|
|
double y_0 = 28.6617; |
189 |
|
|
double z_0 = 5.27786; |
190 |
|
|
*/ |
191 |
|
|
/* |
192 |
|
|
//APEX, node number: ?? : 29.4715 88.8100 62.8179 (coordinates of heart_024_ventricules_nRef_Vol.mesh) |
193 |
|
|
double x_0= 29.4715; |
194 |
|
|
double y_0 = 88.8100; |
195 |
|
|
double z_0 = 62.8179; |
196 |
|
|
*/ |
197 |
|
|
|
198 |
|
✗ |
double tPeriod=fmod(t,FelisceParam::instance().timePeriodVentricle); |
199 |
|
✗ |
double tPer2 = 4*(tPeriod-delay)*(tPeriod-delay); |
200 |
|
|
|
201 |
|
|
//if simulation starts later than t=0 |
202 |
|
✗ |
if ((tPeriod>=0.0)&&(tPeriod<delay)) { |
203 |
|
✗ |
iapp=0.0; |
204 |
|
|
} |
205 |
|
|
|
206 |
|
✗ |
double sum2 = (x-x_0)*(x-x_0) + (y-y_0)*(y-y_0) + (z-z_0)*(z-z_0); |
207 |
|
|
|
208 |
|
✗ |
if ((tPeriod>=delay)&&(tPeriod<=(delay+ FelisceParam::instance().stimTimeLV))) { //if t \in [delay,delay+t_{actLeft}] |
209 |
|
✗ |
if ( Tools::equal(dist,20.0) || (dist < -0.2) ) { // Left ventricle endocardium |
210 |
|
✗ |
if ( (sum2 >= (tPer2-200.0) ) && (sum2 <= (3.0*tPer2)) ) { |
211 |
|
✗ |
iapp=2.0/500.0; |
212 |
|
|
} |
213 |
|
|
} |
214 |
|
|
} |
215 |
|
|
|
216 |
|
✗ |
if ((tPeriod>=delay)&&(tPeriod<=(delay+ FelisceParam::instance().stimTimeRV))) { //if t \in [delay,delay+t_{actRight}] |
217 |
|
✗ |
if ( Tools::equal(dist,200.0) || Tools::equal(dist,10.0) ) { // Right ventricle endocardium |
218 |
|
✗ |
if ( (sum2 >= (tPer2-200.0) ) && (sum2 <= (3.0*tPer2)) ) { |
219 |
|
✗ |
iapp=2.0/500.0; |
220 |
|
|
} |
221 |
|
|
} |
222 |
|
|
} |
223 |
|
|
|
224 |
|
✗ |
return iapp; |
225 |
|
|
|
226 |
|
|
} |
227 |
|
|
|
228 |
|
✗ |
double zygoteIapp::operator() (double x, double y, double z, double t, double dist) const { |
229 |
|
✗ |
return Iapp(x,y,z,t,dist); |
230 |
|
|
} |
231 |
|
|
|
232 |
|
|
|
233 |
|
|
//A class to impose a non-pathological applied current to the atria. |
234 |
|
✗ |
AppCurrentAtria::AppCurrentAtria() = default; |
235 |
|
|
|
236 |
|
✗ |
AppCurrentAtria::~AppCurrentAtria() = default; |
237 |
|
|
|
238 |
|
✗ |
double AppCurrentAtria::IappAtria(double x, double y, double z, double t) const { |
239 |
|
✗ |
double iapp=0.0; |
240 |
|
✗ |
double delay = FelisceParam::instance().delayStim[0]; |
241 |
|
✗ |
double tPeriod=fmod(t,FelisceParam::instance().timePeriodAtria); |
242 |
|
✗ |
double i0 = 0.01*2.5; |
243 |
|
|
|
244 |
|
|
//current in Atrium |
245 |
|
✗ |
if (FelisceParam::instance().hasHeteroCondAtria) { |
246 |
|
✗ |
i0 = 0.01*2.0; |
247 |
|
|
} |
248 |
|
|
|
249 |
|
|
//double x2879 = -1.926082; |
250 |
|
✗ |
double y2879 = 17.376394; |
251 |
|
✗ |
double z2879 = 1.520698; |
252 |
|
|
|
253 |
|
|
//double x3484 = -1.763522; |
254 |
|
✗ |
double y3484 = 17.849342; |
255 |
|
✗ |
double z3484 = 1.738376; |
256 |
|
|
|
257 |
|
|
|
258 |
|
|
//double xmid = (x2879 + x3484)/2.0; |
259 |
|
✗ |
double ymid = (y2879 + y3484)/2.0; |
260 |
|
✗ |
double zmid = (z2879 + z3484)/2.0; |
261 |
|
|
|
262 |
|
|
//double x3095 = -1.877060; |
263 |
|
✗ |
double y3095 = 17.539717; |
264 |
|
✗ |
double z3095 = 1.705356; |
265 |
|
|
|
266 |
|
|
|
267 |
|
✗ |
double a = std::sqrt((z2879-z3484)*(z2879-z3484)+(y2879-y3484)*(y2879-y3484))/2.0; |
268 |
|
✗ |
double b = std::sqrt((zmid-z3095)*(zmid-z3095)+(ymid-y3095)*(ymid-y3095)); |
269 |
|
|
|
270 |
|
✗ |
double xAxis = zmid - z3095; |
271 |
|
✗ |
double yAxis = ymid - y3095; |
272 |
|
✗ |
xAxis = xAxis/std::sqrt(xAxis*xAxis+yAxis*yAxis); |
273 |
|
✗ |
yAxis = yAxis/std::sqrt(xAxis*xAxis+yAxis*yAxis); |
274 |
|
|
|
275 |
|
✗ |
double theta = acos(xAxis*1.0); |
276 |
|
|
|
277 |
|
✗ |
double X = std::cos(theta)*z-std::sin(theta)*y; |
278 |
|
✗ |
double Y = std::sin(theta)*z+std::cos(theta)*y; |
279 |
|
✗ |
double X0 = std::cos(theta)*zmid-std::sin(theta)*ymid; |
280 |
|
✗ |
double Y0 = std::sin(theta)*zmid+std::cos(theta)*ymid; |
281 |
|
|
// if DA court |
282 |
|
✗ |
if (FelisceParam::instance().dataAssimilation) { |
283 |
|
✗ |
X0 -=0.5; |
284 |
|
✗ |
Y0 +=0.5; |
285 |
|
|
} |
286 |
|
|
|
287 |
|
✗ |
double sumA = (X-X0)*(X-X0)/(a*a)+(Y-Y0)*(Y-Y0)/(b*b); |
288 |
|
|
|
289 |
|
|
// Bachman bundle |
290 |
|
|
// double xB = -0.287020; |
291 |
|
|
// double zB = 1.695450; |
292 |
|
|
// |
293 |
|
|
// double xC = 1.320040; |
294 |
|
|
// double zC = 0.896870; |
295 |
|
|
// |
296 |
|
|
// double xBCmid = (xB+xC)/2.0; |
297 |
|
|
// double zBCmid = (zB+zC)/2.0; |
298 |
|
|
// |
299 |
|
|
// double xxaxis = xBCmid - xB; |
300 |
|
|
// double zzaxis = zBCmid - zB; |
301 |
|
|
// |
302 |
|
|
// double theta2 = acos(xxaxis*1.0); |
303 |
|
|
// |
304 |
|
|
// double XX = std::cos(theta2)*x-std::sin(theta2)*z; |
305 |
|
|
// double ZZ = std::sin(theta2)*x+std::cos(theta2)*z; |
306 |
|
|
// double XX0 = std::cos(theta2)*xBCmid-std::sin(theta2)*zBCmid; |
307 |
|
|
// double ZZ0 = std::sin(theta2)*xBCmid+std::cos(theta2)*zBCmid; |
308 |
|
|
// |
309 |
|
|
// |
310 |
|
|
// double sumB = (XX-XX0)*(XX-XX0)/5.0 + 30.0*(ZZ-ZZ0)*(ZZ-ZZ0); |
311 |
|
|
// |
312 |
|
|
// double sumC = (x-xB)*(x-xB)+(z-zB)*(z-zB); |
313 |
|
|
// |
314 |
|
✗ |
double radius_A = 1.0; |
315 |
|
✗ |
if ((tPeriod>=delay)&&(tPeriod<=(delay+FelisceParam::instance().stimTime))) { //if t \in [delay,delay+t_{act}] |
316 |
|
✗ |
if ((sumA < radius_A)&&(x < -1.40)) { |
317 |
|
✗ |
iapp = i0; |
318 |
|
|
} |
319 |
|
|
} |
320 |
|
|
|
321 |
|
|
// Bachman bundle |
322 |
|
|
// if ((tPeriod>=delay+20.0)&&(tPeriod<=(delay+FelisceParam::instance().stimTime+20.0))){ //if t \in [delay,delay+t_{act}] |
323 |
|
|
// if ((sumB < 0.5)&&(y>=13.5)){ |
324 |
|
|
// iapp = i0; |
325 |
|
|
// } |
326 |
|
|
// } |
327 |
|
|
|
328 |
|
✗ |
return iapp; |
329 |
|
|
} |
330 |
|
|
|
331 |
|
✗ |
double AppCurrentAtria::operator() (double x, double y, double z, double t) const { |
332 |
|
✗ |
return IappAtria(x,y,z,t); |
333 |
|
|
} |
334 |
|
|
|
335 |
|
|
|
336 |
|
|
// Heterogeneous schaf parameters: |
337 |
|
|
|
338 |
|
|
//A class to impose heterogeneous tau_close. |
339 |
|
✗ |
HeteroTauClose::HeteroTauClose() = default; |
340 |
|
|
|
341 |
|
✗ |
HeteroTauClose::~HeteroTauClose() = default; |
342 |
|
|
|
343 |
|
✗ |
double HeteroTauClose::TauClose(double x, double y, double z, double dist, double t) const { |
344 |
|
✗ |
double tauClose = FelisceParam::instance().tauClose; |
345 |
|
✗ |
bool test = true; |
346 |
|
|
|
347 |
|
|
// During a short-long-short QT segment sequence - in order to simulate a Torsades de pointe |
348 |
|
|
// Warning: works only for "zygote" or "heart" ore "ellipseheart" |
349 |
|
✗ |
if ( (FelisceParam::instance().torsade) && (t > FelisceParam::instance().torsadeTimeBegin) && (t < FelisceParam::instance().torsadeTimeEnd) ) { |
350 |
|
✗ |
if ( Tools::equal(dist,200.0) ) { |
351 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRVTorsade; |
352 |
|
✗ |
test = false; |
353 |
|
✗ |
} else if ( ((dist >= -1.0) && (dist <= 0.0)) || Tools::equal(dist,20.0) ) { |
354 |
|
✗ |
tauClose = FelisceParam::instance().tauCloseEndoTorsade; |
355 |
|
✗ |
test = false; |
356 |
|
✗ |
} else if ( (dist > 0.0) && (dist <= 0.33333) ) { |
357 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseCellTorsade; |
358 |
|
✗ |
test = false; |
359 |
|
✗ |
} else if ( (dist > 0.33333) && (dist <= 1.0) ) { |
360 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseEpiTorsade; |
361 |
|
✗ |
test = false; |
362 |
|
✗ |
} else if ( Tools::equal(dist,100.0) || Tools::equal(dist,10.0) ) { |
363 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRVTorsade; |
364 |
|
✗ |
test = false; |
365 |
|
✗ |
} else if ( Tools::equal(dist,-100.0) ) { |
366 |
|
✗ |
tauClose = 1.; |
367 |
|
✗ |
test = false; |
368 |
|
✗ |
} else if ( Tools::equal(dist,-200.0) ) { //atria (complete heart) |
369 |
|
✗ |
tauClose = 1.0; |
370 |
|
✗ |
test = false; |
371 |
|
|
} |
372 |
|
|
|
373 |
|
✗ |
if (test) { |
374 |
|
✗ |
std::cout << "WARNING ! \nPoint (" << x << "," << y << "," << z << ") has distance " << dist << " => tauClose not defined !\n"; |
375 |
|
|
} |
376 |
|
|
} else { |
377 |
|
✗ |
if ( FelisceParam::instance().typeOfAppliedCurrent == "atria" ) { |
378 |
|
✗ |
if ( y > 6.7) { |
379 |
|
✗ |
tauClose = 70.0; |
380 |
|
|
} else { |
381 |
|
✗ |
tauClose = 80.0; |
382 |
|
|
} |
383 |
|
✗ |
} else if ( (FelisceParam::instance().typeOfAppliedCurrent == "zygote") || (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
384 |
|
|
// Right Ventricle |
385 |
|
✗ |
if ( Tools::equal(dist,200.0) ) { |
386 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRV; |
387 |
|
✗ |
test = false; |
388 |
|
|
} |
389 |
|
|
// Left Ventricle -- Endo |
390 |
|
✗ |
else if ( ((dist >= -1.0) && (dist <= 0.0)) || Tools::equal(dist,20.0) ) { |
391 |
|
✗ |
tauClose = FelisceParam::instance().tauCloseEndo; |
392 |
|
✗ |
test = false; |
393 |
|
|
} |
394 |
|
|
// Left Ventricle -- Mcells |
395 |
|
✗ |
else if ( (dist > 0.0) && (dist <= 0.33333) ) { |
396 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseCell; |
397 |
|
✗ |
test = false; |
398 |
|
|
} |
399 |
|
|
// Left Ventricle -- Epi |
400 |
|
✗ |
else if ( (dist > 0.33333) && (dist <= 1.0) ) { |
401 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseEpi; |
402 |
|
✗ |
test = false; |
403 |
|
|
} |
404 |
|
|
// Right Ventricle |
405 |
|
✗ |
else if ( Tools::equal(dist,100.0) || Tools::equal(dist,10.0) ) { |
406 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRV; |
407 |
|
✗ |
test = false; |
408 |
|
|
} |
409 |
|
|
// Right Ventricle |
410 |
|
✗ |
else if ( Tools::equal(dist,-100.0) ) { |
411 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRV; |
412 |
|
✗ |
test = false; |
413 |
|
|
} |
414 |
|
|
// Atria (complete heart) |
415 |
|
✗ |
else if ( Tools::equal(dist,-200.0) ) { |
416 |
|
✗ |
tauClose = 75.0; |
417 |
|
✗ |
test = false; |
418 |
|
|
} |
419 |
|
|
|
420 |
|
✗ |
if (test) { |
421 |
|
✗ |
std::cout << "WARNING ! \nPoint (" << x << "," << y << "," << z << ") has distance " << dist << " => tauClose not defined !\n"; |
422 |
|
|
} |
423 |
|
|
} else { |
424 |
|
✗ |
double a = 4.0; |
425 |
|
✗ |
double b = 4.0; |
426 |
|
✗ |
double c = 7.2; |
427 |
|
✗ |
double e4 = 0.; |
428 |
|
✗ |
double e3 = 0.5*0.8; |
429 |
|
✗ |
double e2 = 1.0*0.8; |
430 |
|
✗ |
double e1 = 1.6*0.8; |
431 |
|
|
|
432 |
|
✗ |
double sum1 = x*x/( (a-e1)*(a-e1) ) + y*y/( (b-e1)*(b-e1) ) + z*z/( (c-e1)*(c-e1) ) -1.0; |
433 |
|
|
|
434 |
|
✗ |
double sum2 = x*x/( (a-e2)*(a-e2) ) + y*y/( (b-e2)*(b-e2) ) + z*z/( (c-e2)*(c-e2) ) -1.0; |
435 |
|
|
|
436 |
|
✗ |
double sum3 = x*x/( (a-e3)*(a-e3) ) + y*y/( (b-e3)*(b-e3) ) + z*z/( (c-e3)*(c-e3) ) -1.0; |
437 |
|
|
|
438 |
|
✗ |
double sum4 = x*x/( (a-e4)*(a-e4) ) + y*y/( (b-e4)*(b-e4) ) + z*z/( (c-e4)*(c-e4) ) -1.0; |
439 |
|
|
|
440 |
|
✗ |
if ( (sum1 >= -0.1 ) && (sum2 < 0) ) { |
441 |
|
✗ |
tauClose = FelisceParam::instance().tauCloseEndo; |
442 |
|
✗ |
test = false; |
443 |
|
✗ |
} else if ( (sum2 >= 0.0) && (sum3 < 0)) { |
444 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseCell; |
445 |
|
✗ |
test = false; |
446 |
|
✗ |
} else if ( (sum3 >= 0.0) && (sum4 <= 0.1) ) { |
447 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseEpi; |
448 |
|
✗ |
test = false; |
449 |
|
|
} |
450 |
|
✗ |
if ( x < -2.344 ) { |
451 |
|
✗ |
tauClose =FelisceParam::instance().tauCloseRV; |
452 |
|
✗ |
test = false; |
453 |
|
|
} |
454 |
|
|
|
455 |
|
✗ |
if(test) { |
456 |
|
✗ |
tauClose = FelisceParam::instance().tauCloseRV; |
457 |
|
|
} |
458 |
|
|
} |
459 |
|
|
} |
460 |
|
✗ |
return tauClose; |
461 |
|
|
} |
462 |
|
|
|
463 |
|
✗ |
double HeteroTauClose::operator() (double x, double y, double z, double dist, double t) const { |
464 |
|
✗ |
return TauClose(x,y,z,dist,t); |
465 |
|
|
} |
466 |
|
|
|
467 |
|
|
|
468 |
|
|
//A class to simulate infarction. |
469 |
|
✗ |
HeteroTauOut::HeteroTauOut() = default; |
470 |
|
|
|
471 |
|
✗ |
HeteroTauOut::~HeteroTauOut() = default; |
472 |
|
|
|
473 |
|
✗ |
double HeteroTauOut::TauOut(double x, double y, double z) const { |
474 |
|
✗ |
double tauOut = FelisceParam::instance().tauOut; |
475 |
|
|
|
476 |
|
✗ |
double x1 = x - FelisceParam::instance().x_infarct; |
477 |
|
✗ |
double y1 = y - FelisceParam::instance().y_infarct; |
478 |
|
✗ |
double z1 = z - FelisceParam::instance().z_infarct; |
479 |
|
|
|
480 |
|
✗ |
if ( FelisceParam::instance().typeOfAppliedCurrent == "circle" ) { |
481 |
|
✗ |
double sum2 = std::sqrt(x1*x1 + y1*y1 + z1*z1); |
482 |
|
✗ |
if (sum2 > FelisceParam::instance().radius_infarct) |
483 |
|
✗ |
tauOut = 1.e+13; |
484 |
|
|
else |
485 |
|
✗ |
tauOut = FelisceParam::instance().tauOut; |
486 |
|
✗ |
} else if ( FelisceParam::instance().hasInfarct) { |
487 |
|
✗ |
double sum2 = std::sqrt(x1*x1 + y1*y1 + z1*z1); |
488 |
|
|
|
489 |
|
✗ |
if (sum2 < FelisceParam::instance().radius_infarct) |
490 |
|
✗ |
tauOut = (FelisceParam::instance().tauOut)/10.0; |
491 |
|
|
else |
492 |
|
✗ |
tauOut = FelisceParam::instance().tauOut; |
493 |
|
|
} |
494 |
|
✗ |
return tauOut; |
495 |
|
|
|
496 |
|
|
} |
497 |
|
|
|
498 |
|
✗ |
double HeteroTauOut::operator() (double x, double y, double z) const { |
499 |
|
✗ |
return TauOut(x,y,z); |
500 |
|
|
} |
501 |
|
|
|
502 |
|
|
//A class to simulate infarction. |
503 |
|
✗ |
HeteroTauIn::HeteroTauIn() = default; |
504 |
|
|
|
505 |
|
✗ |
HeteroTauIn::~HeteroTauIn() = default; |
506 |
|
|
|
507 |
|
✗ |
double HeteroTauIn::TauIn(double x, double y, double z) const { |
508 |
|
✗ |
double tauIn = FelisceParam::instance().tauIn; |
509 |
|
|
|
510 |
|
✗ |
double x1 = x - FelisceParam::instance().x_infarct; |
511 |
|
✗ |
double y1 = y - FelisceParam::instance().y_infarct; |
512 |
|
✗ |
double z1 = z - FelisceParam::instance().z_infarct; |
513 |
|
|
|
514 |
|
✗ |
if ( FelisceParam::instance().typeOfAppliedCurrent == "circle" ) { |
515 |
|
✗ |
double sum2 = std::sqrt(x1*x1 + y1*y1 + z1*z1); |
516 |
|
✗ |
if (sum2 > FelisceParam::instance().radius_infarct) { |
517 |
|
✗ |
tauIn = 1.e+13; |
518 |
|
|
} else { |
519 |
|
✗ |
tauIn = FelisceParam::instance().tauIn; |
520 |
|
|
} |
521 |
|
|
} |
522 |
|
✗ |
return tauIn; |
523 |
|
|
|
524 |
|
|
} |
525 |
|
|
|
526 |
|
✗ |
double HeteroTauIn::operator() (double x, double y, double z) const { |
527 |
|
✗ |
return TauIn(x,y,z); |
528 |
|
|
} |
529 |
|
|
|
530 |
|
|
//A class to impose heterogeneous Ito for the Courtemanche-Ramirez-Nattel model. |
531 |
|
✗ |
HeteroCourtModelIto::HeteroCourtModelIto() = default; |
532 |
|
|
|
533 |
|
✗ |
HeteroCourtModelIto::~HeteroCourtModelIto() = default; |
534 |
|
|
|
535 |
|
✗ |
double HeteroCourtModelIto::CourtCondIto(double x, double y, double z) const { |
536 |
|
✗ |
double courtCondIto = 0.0; |
537 |
|
✗ |
double z1 = 0.5; |
538 |
|
✗ |
double z2 = 4.0; |
539 |
|
✗ |
double y1 = 10.7; |
540 |
|
✗ |
double y2 = 17.4; |
541 |
|
|
|
542 |
|
|
|
543 |
|
✗ |
double z0 = (z1+z2)/2.0; |
544 |
|
✗ |
double y0 = (y1+y2)/2.0; |
545 |
|
|
|
546 |
|
✗ |
double a = (z2-z1)/2.0; |
547 |
|
✗ |
double b = (y2-y1)/2.0; |
548 |
|
|
|
549 |
|
✗ |
if(((z-z0)*(z-z0)/(a*a)+(y-y0)*(y-y0)/(b*b) > 1.1) && (x < -1.6) && ((z-z0)*(z-z0)/(a*a)+(y-y0)*(y-y0)/(b*b) < 1.55) && (z < 1.7) && (y < 17.5)) { |
550 |
|
✗ |
courtCondIto = 0.2215; |
551 |
|
|
} else { |
552 |
|
✗ |
courtCondIto = 0.1652; |
553 |
|
|
} |
554 |
|
|
|
555 |
|
✗ |
return courtCondIto; |
556 |
|
|
|
557 |
|
|
} |
558 |
|
|
|
559 |
|
✗ |
double HeteroCourtModelIto::operator() (double x, double y, double z) const { |
560 |
|
✗ |
return CourtCondIto(x,y,z); |
561 |
|
|
} |
562 |
|
|
|
563 |
|
|
//A class to impose heterogeneous ICaL for the Courtemanche-Ramirez-Nattel model. |
564 |
|
✗ |
HeteroCourtModelICaL::HeteroCourtModelICaL() = default; |
565 |
|
|
|
566 |
|
✗ |
HeteroCourtModelICaL::~HeteroCourtModelICaL() = default; |
567 |
|
|
|
568 |
|
✗ |
double HeteroCourtModelICaL::CourtCondICaL(double x, double y, double z) const { |
569 |
|
✗ |
double courtCondICaL = 0.0; |
570 |
|
✗ |
double z1 = 0.5; |
571 |
|
✗ |
double z2 = 4.0; |
572 |
|
✗ |
double y1 = 10.7; |
573 |
|
✗ |
double y2 = 17.4; |
574 |
|
|
|
575 |
|
|
|
576 |
|
✗ |
double z0 = (z1+z2)/2.0; |
577 |
|
✗ |
double y0 = (y1+y2)/2.0; |
578 |
|
|
|
579 |
|
✗ |
double a = (z2-z1)/2.0; |
580 |
|
✗ |
double b = (y2-y1)/2.0; |
581 |
|
|
|
582 |
|
✗ |
if(((z-z0)*(z-z0)/(a*a)+(y-y0)*(y-y0)/(b*b) > 1.1) && (x < -1.6) && ((z-z0)*(z-z0)/(a*a)+(y-y0)*(y-y0)/(b*b) < 1.55) && (z < 1.7) && (y < 17.5)) { |
583 |
|
✗ |
courtCondICaL = 0.2067; |
584 |
|
|
} else { |
585 |
|
✗ |
courtCondICaL = 0.1238; |
586 |
|
|
} |
587 |
|
|
|
588 |
|
✗ |
return courtCondICaL; |
589 |
|
|
|
590 |
|
|
} |
591 |
|
|
|
592 |
|
✗ |
double HeteroCourtModelICaL::operator() (double x, double y, double z) const { |
593 |
|
✗ |
return CourtCondICaL(x,y,z); |
594 |
|
|
} |
595 |
|
|
|
596 |
|
|
|
597 |
|
|
//A class to impose heterogeneous multiplicative coefficient for the Courtemanche-Ramirez-Nattel model. |
598 |
|
✗ |
HeteroCourtModelMultCoeff::HeteroCourtModelMultCoeff() = default; |
599 |
|
|
|
600 |
|
✗ |
HeteroCourtModelMultCoeff::~HeteroCourtModelMultCoeff() = default; |
601 |
|
|
|
602 |
|
✗ |
double HeteroCourtModelMultCoeff::CourtCondMultCoeff(double x, double y, double z, double ref) const { |
603 |
|
|
IGNORE_UNUSED_ARGUMENT(x); |
604 |
|
|
IGNORE_UNUSED_ARGUMENT(y); |
605 |
|
|
IGNORE_UNUSED_ARGUMENT(z); |
606 |
|
|
|
607 |
|
✗ |
double courtCondMultCoeff = 1.0; |
608 |
|
|
|
609 |
|
✗ |
double g_Na_PM = FelisceParam::instance().g_Na_PM; |
610 |
|
✗ |
double g_Na_CT = FelisceParam::instance().g_Na_CT; |
611 |
|
✗ |
double g_Na_BB = FelisceParam::instance().g_Na_BB; |
612 |
|
|
|
613 |
|
✗ |
if (std::fabs(ref - 31.0) < 1.0e-12) { //SN |
614 |
|
✗ |
courtCondMultCoeff = 1.0; //1.1; |
615 |
|
✗ |
} else if (std::fabs(ref - 32.0) < 1.0e-12) { //CT |
616 |
|
✗ |
if (g_Na_CT < 0.) { |
617 |
|
✗ |
courtCondMultCoeff = 4.0; //3.4; |
618 |
|
|
} |
619 |
|
|
else { |
620 |
|
✗ |
courtCondMultCoeff = g_Na_CT; |
621 |
|
|
} |
622 |
|
✗ |
} else if (std::fabs(ref - 33.0) < 1.0e-12) { //BB |
623 |
|
✗ |
if (FelisceParam::instance().hasPartialBachmannBundleBlock) |
624 |
|
✗ |
courtCondMultCoeff = FelisceParam::instance().valueBBBlock; |
625 |
|
|
else |
626 |
|
✗ |
if (g_Na_BB < 0.) { |
627 |
|
✗ |
courtCondMultCoeff = 6.0; //6.0; |
628 |
|
|
} |
629 |
|
|
else { |
630 |
|
✗ |
courtCondMultCoeff = g_Na_BB; |
631 |
|
|
} |
632 |
|
✗ |
} else if (std::fabs(ref - 34.0) < 1.0e-12) { //RAI |
633 |
|
✗ |
courtCondMultCoeff = 1.0; //0.5; |
634 |
|
✗ |
} else if (std::fabs(ref - 35.0) < 1.0e-12) { //PM |
635 |
|
✗ |
if (g_Na_PM < 0.) { |
636 |
|
✗ |
courtCondMultCoeff = 1.5; //1.3; |
637 |
|
|
} |
638 |
|
|
else { |
639 |
|
✗ |
courtCondMultCoeff = g_Na_PM; |
640 |
|
|
} |
641 |
|
✗ |
} else if (std::fabs(ref - 36.0) < 1.0e-12) { //FO |
642 |
|
✗ |
courtCondMultCoeff = 0.5; |
643 |
|
|
} else { //regular atria |
644 |
|
✗ |
courtCondMultCoeff = 1.0; |
645 |
|
|
} |
646 |
|
|
|
647 |
|
✗ |
return courtCondMultCoeff; |
648 |
|
|
|
649 |
|
|
} |
650 |
|
|
|
651 |
|
✗ |
double HeteroCourtModelMultCoeff::operator() (double x, double y, double z, double ref) const { |
652 |
|
✗ |
return CourtCondMultCoeff(x,y,z,ref); |
653 |
|
|
} |
654 |
|
|
|
655 |
|
|
|
656 |
|
|
// Class to impose variable "s" parameter in FHN |
657 |
|
✗ |
HeteroSparFHN::HeteroSparFHN() = default; |
658 |
|
✗ |
HeteroSparFHN::~HeteroSparFHN() = default; |
659 |
|
|
|
660 |
|
✗ |
double HeteroSparFHN::SparFHN(double x, double y, double z) const { |
661 |
|
|
|
662 |
|
|
// IGNORE_UNUSED_ARGUMENT(z); |
663 |
|
|
// double ratioS = 10.0; |
664 |
|
|
// double y0 = 0.25; |
665 |
|
|
// double y1 = 0.75; |
666 |
|
|
// |
667 |
|
|
// double sPar = FelisceParam::instance().f0; |
668 |
|
|
// if( (x<=0.5)&&(y>=y0)&&(y<=y1)) { |
669 |
|
|
// sPar = sPar*( ((ratioS-1)/ratioS) * (y-y1)/(y1-y0) + 1.0); |
670 |
|
|
// //sPar = 2*sPar*(1.0-ratioS)*y + sPar*(2*ratioS-1); |
671 |
|
|
// } |
672 |
|
|
// return sPar; |
673 |
|
|
|
674 |
|
✗ |
double x1 = x - FelisceParam::instance().x_infarct; |
675 |
|
✗ |
double y1 = y - FelisceParam::instance().y_infarct; |
676 |
|
✗ |
double z1 = z - FelisceParam::instance().z_infarct; |
677 |
|
✗ |
double sum1 = std::sqrt(x1*x1 + y1*y1 + z1*z1); |
678 |
|
|
|
679 |
|
✗ |
double sPar = FelisceParam::instance().f0; |
680 |
|
✗ |
if( sum1 < FelisceParam::instance().radius_infarct) { |
681 |
|
✗ |
sPar = sPar/10.; |
682 |
|
|
} |
683 |
|
✗ |
return sPar; |
684 |
|
|
|
685 |
|
|
} |
686 |
|
|
|
687 |
|
✗ |
double HeteroSparFHN::operator() (double x, double y, double z) const { |
688 |
|
✗ |
return SparFHN(x,y,z); |
689 |
|
|
} |
690 |
|
|
|
691 |
|
|
|
692 |
|
|
// Heterogeneous schaf parameters: |
693 |
|
|
|
694 |
|
|
//A class to impose heterogeneous tau_close. |
695 |
|
✗ |
HeteroMVCoeff::HeteroMVCoeff() = default; |
696 |
|
|
|
697 |
|
✗ |
HeteroMVCoeff::~HeteroMVCoeff() = default; |
698 |
|
|
|
699 |
|
✗ |
int HeteroMVCoeff::MVCoeff(double x, double y, double z, double dist, double) const { |
700 |
|
✗ |
int coeffType = -1; |
701 |
|
✗ |
bool test = true; |
702 |
|
|
|
703 |
|
✗ |
if ( (FelisceParam::instance().typeOfAppliedCurrent == "zygote") || (FelisceParam::instance().typeOfAppliedCurrent == "heart") || (FelisceParam::instance().typeOfAppliedCurrent == "ellipseheart") ) { |
704 |
|
|
// Right Ventricle |
705 |
|
✗ |
if ( Tools::equal(dist,200.0) ) { |
706 |
|
✗ |
coeffType = 0; // RV |
707 |
|
✗ |
test = false; |
708 |
|
|
} |
709 |
|
|
// Left Ventricle -- Endo |
710 |
|
✗ |
else if ( ((dist >= -1.0) && (dist <= /*0.0*/ -0.1)) || Tools::equal(dist,20.0) ) { |
711 |
|
✗ |
coeffType = 1; // Endo |
712 |
|
✗ |
test = false; |
713 |
|
|
} |
714 |
|
|
// Left Ventricle -- Mcells |
715 |
|
✗ |
else if ( (dist > /*0.0*/ -0.1) && (dist <= /*0.33333*/ 0.0 ) ) { |
716 |
|
✗ |
coeffType = 2; // MCell |
717 |
|
✗ |
test = false; |
718 |
|
|
} |
719 |
|
|
// Left Ventricle -- Epi |
720 |
|
✗ |
else if ( (dist > /*0.33333*/ 0.0) && (dist <= 1.0) ) { |
721 |
|
✗ |
coeffType = 3; // Epi |
722 |
|
✗ |
test = false; |
723 |
|
|
} |
724 |
|
|
// Right Ventricle |
725 |
|
✗ |
else if ( Tools::equal(dist,100.0) || Tools::equal(dist,10.0) ) { |
726 |
|
✗ |
coeffType = 0; // RV |
727 |
|
✗ |
test = false; |
728 |
|
|
} |
729 |
|
|
// Right Ventricle |
730 |
|
✗ |
else if ( Tools::equal(dist,-100.0) ) { |
731 |
|
✗ |
coeffType = 0; |
732 |
|
✗ |
test = false; |
733 |
|
|
} |
734 |
|
✗ |
else if ( Tools::equal(dist,-200.0) ) { //atria (complete heart) |
735 |
|
✗ |
coeffType = 4; |
736 |
|
✗ |
test = false; |
737 |
|
|
} |
738 |
|
|
|
739 |
|
✗ |
if (test) { |
740 |
|
✗ |
std::cout << "WARNING ! \nPoint (" << x << "," << y << "," << z << ") has distance " << dist << " => MV coeffType not defined !\n"; |
741 |
|
|
} |
742 |
|
|
} |
743 |
|
✗ |
return coeffType; |
744 |
|
|
} |
745 |
|
|
|
746 |
|
✗ |
int HeteroMVCoeff::operator() (double x, double y, double z, double dist, double t) const { |
747 |
|
✗ |
return MVCoeff(x,y,z,dist,t); |
748 |
|
|
} |
749 |
|
|
} |
750 |
|
|
|