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 "Model/bidomainCurvModel.hpp" |
21 |
|
|
|
22 |
|
|
namespace felisce { |
23 |
|
✗ |
BidomainCurvModel::BidomainCurvModel(): |
24 |
|
|
Model(), |
25 |
|
✗ |
m_ionic(nullptr), |
26 |
|
✗ |
m_schaf(nullptr), |
27 |
|
✗ |
m_fhn(nullptr), |
28 |
|
✗ |
m_court(nullptr), |
29 |
|
✗ |
m_buildIonic(false), |
30 |
|
✗ |
m_buildCourt(false), |
31 |
|
✗ |
m_verbose(0) { |
32 |
|
✗ |
for(std::size_t i=0; i<m_linearProblemBidomainCurv.size(); i++) { |
33 |
|
✗ |
m_linearProblemBidomainCurv[i] = nullptr; |
34 |
|
|
} |
35 |
|
✗ |
m_name = "BidomainCurv"; |
36 |
|
|
} |
37 |
|
|
|
38 |
|
✗ |
BidomainCurvModel::~BidomainCurvModel() { |
39 |
|
✗ |
for(std::size_t i=0; i<m_linearProblemBidomainCurv.size(); i++) { |
40 |
|
✗ |
m_linearProblemBidomainCurv[i] = nullptr; |
41 |
|
|
} |
42 |
|
✗ |
if (m_buildIonic) { |
43 |
|
✗ |
delete m_ionic; |
44 |
|
✗ |
m_ionic = nullptr; |
45 |
|
|
} |
46 |
|
✗ |
if (m_buildCourt) { |
47 |
|
✗ |
delete m_court; |
48 |
|
✗ |
m_court = nullptr; |
49 |
|
|
} |
50 |
|
✗ |
m_solExtraCell.destroy(); |
51 |
|
✗ |
m_extrapolate.destroy(); |
52 |
|
✗ |
if (m_buildCourt) { |
53 |
|
✗ |
m_m.destroy(); |
54 |
|
✗ |
m_h.destroy(); |
55 |
|
✗ |
m_j.destroy(); |
56 |
|
✗ |
m_ao.destroy(); |
57 |
|
✗ |
m_io.destroy(); |
58 |
|
✗ |
m_ua.destroy(); |
59 |
|
✗ |
m_ui.destroy(); |
60 |
|
✗ |
m_xr.destroy(); |
61 |
|
✗ |
m_xs.destroy(); |
62 |
|
✗ |
m_d.destroy(); |
63 |
|
✗ |
m_f.destroy(); |
64 |
|
✗ |
m_fca.destroy(); |
65 |
|
✗ |
m_urel.destroy(); |
66 |
|
✗ |
m_vrel.destroy(); |
67 |
|
✗ |
m_wrel.destroy(); |
68 |
|
✗ |
m_nai.destroy(); |
69 |
|
✗ |
m_nao.destroy(); |
70 |
|
✗ |
m_cao.destroy(); |
71 |
|
✗ |
m_ki.destroy(); |
72 |
|
✗ |
m_ko.destroy(); |
73 |
|
✗ |
m_cai.destroy(); |
74 |
|
✗ |
m_naiont.destroy(); |
75 |
|
✗ |
m_kiont.destroy(); |
76 |
|
✗ |
m_caiont.destroy(); |
77 |
|
✗ |
m_ileak.destroy(); |
78 |
|
✗ |
m_iup.destroy(); |
79 |
|
✗ |
m_itr.destroy(); |
80 |
|
✗ |
m_irel.destroy(); |
81 |
|
✗ |
m_cmdn.destroy(); |
82 |
|
✗ |
m_trpn.destroy(); |
83 |
|
✗ |
m_nsr.destroy(); |
84 |
|
✗ |
m_jsr.destroy(); |
85 |
|
✗ |
m_csqn.destroy(); |
86 |
|
|
} else { |
87 |
|
✗ |
m_W_0.destroy(); |
88 |
|
|
} |
89 |
|
✗ |
delete m_iApp; |
90 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
91 |
|
✗ |
delete m_iAppExt; |
92 |
|
|
} |
93 |
|
|
|
94 |
|
✗ |
void BidomainCurvModel::initializeDerivedLinearProblem() { |
95 |
|
✗ |
for (std::size_t i=0; i<m_linearProblem.size(); i++) { |
96 |
|
✗ |
m_linearProblemBidomainCurv.push_back(static_cast<LinearProblemBidomainCurv*>(m_linearProblem[i])); |
97 |
|
|
} |
98 |
|
|
} |
99 |
|
|
|
100 |
|
|
//! Define m_schaf and m_bdfEdp (use in linearProblem). |
101 |
|
✗ |
void BidomainCurvModel::initializeDerivedModel() { |
102 |
|
✗ |
m_verbose = FelisceParam::verbose(); |
103 |
|
|
|
104 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "schaf") { |
105 |
|
|
//Allocate m_schaf. |
106 |
|
✗ |
m_schaf = new SchafSolver(m_fstransient); |
107 |
|
✗ |
m_ionic = m_schaf; |
108 |
|
|
//Fix order of bdf for schaf solver. |
109 |
|
✗ |
m_schaf->defineOrderBdf(FelisceParam::instance().orderBdfIonic); |
110 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "fhn") { |
111 |
|
|
// Initialize Fitzhugh-Nagumo model |
112 |
|
✗ |
m_fhn = new FhNSolver(m_fstransient); |
113 |
|
✗ |
m_ionic =m_fhn; |
114 |
|
✗ |
m_ionic->defineOrderBdf(FelisceParam::instance().orderBdfIonic,1); |
115 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "court") { |
116 |
|
|
//Allocate m_court. |
117 |
|
✗ |
m_court = new CourtemancheSolver(m_fstransient); |
118 |
|
|
} |
119 |
|
|
|
120 |
|
|
//Define order for bdf used by EDP (linearProblem). |
121 |
|
✗ |
m_bdfEdp.defineOrder(FelisceParam::instance().orderBdfEdp); |
122 |
|
|
|
123 |
|
✗ |
initializeAppCurrent(); |
124 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
125 |
|
✗ |
initializeAppCurrentExt(); |
126 |
|
|
|
127 |
|
|
} |
128 |
|
|
|
129 |
|
✗ |
void BidomainCurvModel::setInitialCondition() { |
130 |
|
✗ |
for (std::size_t iProblem=0; iProblem<m_linearProblemBidomainCurv.size(); iProblem++) { |
131 |
|
✗ |
m_U_0.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
132 |
|
✗ |
m_U_0.set(0.0); |
133 |
|
✗ |
if (FelisceParam::instance().hasInitialCondition) { |
134 |
|
|
|
135 |
|
✗ |
int idVar = -1; |
136 |
|
|
double valueByDofU_0; |
137 |
|
|
|
138 |
|
✗ |
for (std::size_t iListVar = 0; iListVar < m_initialCondition.listVariable().size(); iListVar++) { |
139 |
|
✗ |
for ( std::size_t iUnknown = 0; iUnknown < m_linearProblemBidomainCurv[iProblem]->listUnknown().size(); iUnknown++) { |
140 |
|
✗ |
idVar = m_linearProblemBidomainCurv[iProblem]->listUnknown().idVariable(iUnknown); |
141 |
|
✗ |
if ( m_linearProblemBidomainCurv[iProblem]->listVariable()[idVar].name() == m_initialCondition.listVariable()[iListVar].name() ) { |
142 |
|
✗ |
felInt numDofLpU = m_linearProblemBidomainCurv[iProblem]->numDofLocalPerUnknown(m_linearProblemBidomainCurv[iProblem]->listUnknown()[iUnknown]); |
143 |
|
✗ |
felInt idLocalValue[numDofLpU]; |
144 |
|
✗ |
felInt idGlobalValue[numDofLpU]; |
145 |
|
✗ |
for ( felInt ii = 0; ii < numDofLpU; ii++) { |
146 |
|
✗ |
idLocalValue[ii] = ii; |
147 |
|
|
} |
148 |
|
✗ |
ISLocalToGlobalMappingApply(m_linearProblemBidomainCurv[iProblem]->mappingDofLocalToDofGlobal(m_linearProblemBidomainCurv[iProblem]->listUnknown()[iUnknown]),numDofLpU,&idLocalValue[0],&idGlobalValue[0]); |
149 |
|
✗ |
if (idVar==0) { |
150 |
|
✗ |
for (felInt ii = 0; ii < numDofLpU; ii++) { |
151 |
|
✗ |
valueByDofU_0 = FelisceParam::instance().valueInitCond[iListVar]; |
152 |
|
✗ |
m_U_0.setValue(idGlobalValue[ii],valueByDofU_0,INSERT_VALUES);; |
153 |
|
|
} |
154 |
|
✗ |
} else if (idVar==1) { |
155 |
|
✗ |
for (felInt ii = 0; ii < numDofLpU; ii++) { |
156 |
|
✗ |
if (ii < static_cast<int>(numDofLpU/2)) { |
157 |
|
✗ |
valueByDofU_0 = 0.0; |
158 |
|
|
} else { |
159 |
|
✗ |
valueByDofU_0 = 0.0; |
160 |
|
|
} |
161 |
|
✗ |
m_U_0.setValue(idGlobalValue[ii],valueByDofU_0,INSERT_VALUES); |
162 |
|
|
} |
163 |
|
|
} |
164 |
|
|
} |
165 |
|
|
} |
166 |
|
|
} |
167 |
|
|
} |
168 |
|
✗ |
m_U_0.assembly(); |
169 |
|
|
|
170 |
|
|
//Initialize solution for solver. |
171 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->solution().copyFrom(m_U_0); |
172 |
|
|
|
173 |
|
|
} |
174 |
|
|
|
175 |
|
|
} |
176 |
|
|
|
177 |
|
✗ |
void BidomainCurvModel::initializeAppCurrent() { |
178 |
|
✗ |
m_iApp = new AppCurrentAtria(); |
179 |
|
✗ |
m_iApp->initialize(m_fstransient); |
180 |
|
|
} |
181 |
|
|
|
182 |
|
✗ |
void BidomainCurvModel::initializeAppCurrentExt() { |
183 |
|
✗ |
m_iAppExt = new AppCurrentAtria(); |
184 |
|
✗ |
m_iAppExt->initialize(m_fstransient); |
185 |
|
|
} |
186 |
|
|
|
187 |
|
✗ |
void BidomainCurvModel::evalIapp() { |
188 |
|
✗ |
m_iAppValue.clear(); |
189 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(*m_iApp,m_fstransient->time,m_iAppValue); |
190 |
|
|
} |
191 |
|
|
|
192 |
|
✗ |
void BidomainCurvModel::evalIappExt() { |
193 |
|
✗ |
m_iAppValueExt.clear(); |
194 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(*m_iAppExt,m_fstransient->time,m_iAppValueExt); |
195 |
|
|
} |
196 |
|
|
|
197 |
|
✗ |
void BidomainCurvModel::finalizeRHSDP(std::size_t iProblem) { |
198 |
|
✗ |
int size = m_linearProblemBidomainCurv[iProblem]->numDofLocalPerUnknown(potTransMemb); |
199 |
|
✗ |
double valueForPetsc[size]; |
200 |
|
✗ |
felInt idLocalValue[size]; |
201 |
|
✗ |
felInt idGlobalValue[size]; |
202 |
|
|
|
203 |
|
✗ |
double coefAm = FelisceParam::instance().Am; |
204 |
|
✗ |
double coefCm = FelisceParam::instance().Cm; |
205 |
|
|
|
206 |
|
|
|
207 |
|
|
// RHS = Am*Cm*BDF.RHS (Warning : bdf.RHS containes both problem variables). |
208 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().axpy(1.,m_bdfEdp.vector()); |
209 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().scale(coefAm*coefCm); |
210 |
|
|
|
211 |
|
|
// RHS = RHS + Am*Iapp |
212 |
|
✗ |
for (felInt i = 0; i < size; i++) { |
213 |
|
✗ |
idLocalValue[i] = i; |
214 |
|
|
} |
215 |
|
✗ |
ISLocalToGlobalMappingApply(m_linearProblemBidomainCurv[iProblem]->mappingDofLocalToDofGlobal(potTransMemb),size,&idLocalValue[0],&idGlobalValue[0]); |
216 |
|
✗ |
for (felInt i = 0; i < size; i++) { |
217 |
|
✗ |
felInt indexGlobal=idGlobalValue[i]; |
218 |
|
✗ |
valueForPetsc[i] = coefAm*m_iAppValue[indexGlobal]; |
219 |
|
|
} |
220 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().setValues(size,&idGlobalValue[0],&valueForPetsc[0],ADD_VALUES); |
221 |
|
|
|
222 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().assembly(); |
223 |
|
|
|
224 |
|
|
// For the second equation : RHS = 0 ( or applied exterior current ) |
225 |
|
✗ |
int sizePotExtraCell = m_linearProblemBidomainCurv[iProblem]->numDofLocalPerUnknown(potExtraCell); |
226 |
|
✗ |
double valueForPetscPotExtraCell[sizePotExtraCell]; |
227 |
|
✗ |
felInt idLocalValuePotExtraCell[sizePotExtraCell]; |
228 |
|
✗ |
felInt idGlobalValuePotExtraCell[sizePotExtraCell]; |
229 |
|
|
|
230 |
|
✗ |
for (felInt i = 0; i < sizePotExtraCell; i++) { |
231 |
|
✗ |
idLocalValuePotExtraCell[i] = i; |
232 |
|
|
} |
233 |
|
✗ |
ISLocalToGlobalMappingApply(m_linearProblemBidomainCurv[iProblem]->mappingDofLocalToDofGlobal(potExtraCell),sizePotExtraCell,&idLocalValuePotExtraCell[0],&idGlobalValuePotExtraCell[0]); |
234 |
|
✗ |
for (felInt i = 0; i < sizePotExtraCell; i++) { |
235 |
|
|
// For an applied exterior current |
236 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) { |
237 |
|
✗ |
valueForPetscPotExtraCell[i] = m_iAppValueExt[idGlobalValuePotExtraCell[i]]; |
238 |
|
|
} else { |
239 |
|
✗ |
valueForPetscPotExtraCell[i] = 0.0; |
240 |
|
|
} |
241 |
|
|
} |
242 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().setValues(sizePotExtraCell ,&idGlobalValuePotExtraCell[0],&valueForPetscPotExtraCell[0],INSERT_VALUES); |
243 |
|
|
|
244 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->vector().assembly(); |
245 |
|
|
|
246 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "court") { |
247 |
|
✗ |
m_linearProblemBidomainCurv[0]->vector().axpy(-coefAm*coefCm,m_court->ion()); |
248 |
|
|
} else { |
249 |
|
✗ |
m_linearProblemBidomainCurv[0]->vector().axpy(coefAm,m_ionic->ion()); |
250 |
|
|
} |
251 |
|
|
|
252 |
|
|
// Elisa, Dec. 2014 - Luenberger filter |
253 |
|
✗ |
if (FelisceParam::instance().electrodeControl) { |
254 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->addElectrodeCondtrol(); |
255 |
|
|
} |
256 |
|
|
|
257 |
|
|
|
258 |
|
|
} |
259 |
|
|
|
260 |
|
✗ |
void BidomainCurvModel::preAssemblingMatrixRHS(std::size_t iProblem) { |
261 |
|
|
//Initialize std::vector useful (with m_W_0) with a std::vector with same structure as _RHS in linearProblem. |
262 |
|
|
|
263 |
|
|
// 1) copy structure in m_W_0 (initialize solution at time 0 of EDO in schaf solver). |
264 |
|
✗ |
m_W_0.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
265 |
|
✗ |
m_W_0.zeroEntries(); |
266 |
|
|
|
267 |
|
|
// 2) copy structure in m_extrapolate (contains extrapolate values of linearProblem solution). |
268 |
|
✗ |
m_extrapolate.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
269 |
|
✗ |
m_extrapolate.zeroEntries(); |
270 |
|
|
|
271 |
|
✗ |
m_solExtraCell.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
272 |
|
✗ |
m_solExtraCell.zeroEntries(); |
273 |
|
|
|
274 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "schaf") { |
275 |
|
✗ |
double valueByDofW_0 = 1./((FelisceParam::instance().vMax - FelisceParam::instance().vMin)*(FelisceParam::instance().vMax - FelisceParam::instance().vMin)); |
276 |
|
✗ |
m_W_0.set(valueByDofW_0); |
277 |
|
|
} |
278 |
|
|
|
279 |
|
|
|
280 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "court") { |
281 |
|
✗ |
double valueByDofm = 0.00291; |
282 |
|
✗ |
double valueByDofh = 0.965; |
283 |
|
✗ |
double valueByDofj = 0.978; |
284 |
|
✗ |
double valueByDofao = 0.0304; |
285 |
|
✗ |
double valueByDofio = 0.999; |
286 |
|
✗ |
double valueByDofua = 0.00496; |
287 |
|
✗ |
double valueByDofui = 0.999; |
288 |
|
✗ |
double valueByDofxr = 0.0000329; |
289 |
|
✗ |
double valueByDofxs = 0.0187; |
290 |
|
✗ |
double valueByDofd = 0.000137; |
291 |
|
✗ |
double valueByDoff = 0.999837; |
292 |
|
✗ |
double valueByDoffca = 0.775; |
293 |
|
✗ |
double valueByDofvrel = 1.00; |
294 |
|
✗ |
double valueByDofwrel = 0.999; |
295 |
|
✗ |
double valueByDofnai = 11.2; |
296 |
|
✗ |
double valueByDofnao = 140.0; |
297 |
|
✗ |
double valueByDofcao = 1.8; |
298 |
|
✗ |
double valueByDofki = 139.0; |
299 |
|
✗ |
double valueByDofko = 4.5; |
300 |
|
✗ |
double valueByDofcai = 0.000102; |
301 |
|
✗ |
double valueByDofcmdn = 0.00205; |
302 |
|
✗ |
double valueByDoftrpn = 0.0118; |
303 |
|
✗ |
double valueByDofnsr = 1.49; |
304 |
|
✗ |
double valueByDofjsr = 1.49; |
305 |
|
✗ |
double valueByDofcsqn = 6.51; |
306 |
|
|
|
307 |
|
✗ |
m_m.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
308 |
|
✗ |
m_m.set(valueByDofm); |
309 |
|
✗ |
m_h.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
310 |
|
✗ |
m_h.set(valueByDofh); |
311 |
|
✗ |
m_j.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
312 |
|
✗ |
m_j.set(valueByDofj); |
313 |
|
✗ |
m_ao.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
314 |
|
✗ |
m_ao.set(valueByDofao); |
315 |
|
✗ |
m_io.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
316 |
|
✗ |
m_io.set(valueByDofio); |
317 |
|
✗ |
m_ua.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
318 |
|
✗ |
m_ua.set(valueByDofua); |
319 |
|
✗ |
m_ui.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
320 |
|
✗ |
m_ui.set(valueByDofui); |
321 |
|
✗ |
m_xr.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
322 |
|
✗ |
m_xr.set(valueByDofxr); |
323 |
|
✗ |
m_xs.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
324 |
|
✗ |
m_xs.set(valueByDofxs); |
325 |
|
✗ |
m_d.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
326 |
|
✗ |
m_d.set(valueByDofd); |
327 |
|
✗ |
m_f.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
328 |
|
✗ |
m_f.set(valueByDoff); |
329 |
|
✗ |
m_fca.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
330 |
|
✗ |
m_fca.set(valueByDoffca); |
331 |
|
✗ |
m_urel.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
332 |
|
✗ |
m_urel.set(0.); |
333 |
|
✗ |
m_vrel.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
334 |
|
✗ |
m_vrel.set(valueByDofvrel); |
335 |
|
✗ |
m_wrel.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
336 |
|
✗ |
m_wrel.set(valueByDofwrel); |
337 |
|
✗ |
m_nai.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
338 |
|
✗ |
m_nai.set(valueByDofnai); |
339 |
|
✗ |
m_nao.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
340 |
|
✗ |
m_nao.set(valueByDofnao); |
341 |
|
✗ |
m_cao.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
342 |
|
✗ |
m_cao.set(valueByDofcao); |
343 |
|
✗ |
m_ki.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
344 |
|
✗ |
m_ki.set(valueByDofki); |
345 |
|
✗ |
m_ko.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
346 |
|
✗ |
m_ko.set(valueByDofko); |
347 |
|
✗ |
m_cai.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
348 |
|
✗ |
m_cai.set(valueByDofcai); |
349 |
|
✗ |
m_naiont.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
350 |
|
✗ |
m_naiont.set(0.); |
351 |
|
✗ |
m_kiont.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
352 |
|
✗ |
m_kiont.set(0.); |
353 |
|
✗ |
m_caiont.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
354 |
|
✗ |
m_caiont.set(0.); |
355 |
|
✗ |
m_ileak.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
356 |
|
✗ |
m_ileak.set(0.); |
357 |
|
✗ |
m_iup.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
358 |
|
✗ |
m_iup.set(0.); |
359 |
|
✗ |
m_itr.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
360 |
|
✗ |
m_itr.set(0.); |
361 |
|
✗ |
m_irel.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
362 |
|
✗ |
m_irel.set(0.); |
363 |
|
✗ |
m_cmdn.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
364 |
|
✗ |
m_cmdn.set(valueByDofcmdn); |
365 |
|
✗ |
m_trpn.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
366 |
|
✗ |
m_trpn.set(valueByDoftrpn); |
367 |
|
✗ |
m_nsr.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
368 |
|
✗ |
m_nsr.set(valueByDofnsr); |
369 |
|
✗ |
m_jsr.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
370 |
|
✗ |
m_jsr.set(valueByDofjsr); |
371 |
|
✗ |
m_csqn.duplicateFrom(m_linearProblemBidomainCurv[iProblem]->vector()); |
372 |
|
✗ |
m_csqn.set(valueByDofcsqn); |
373 |
|
|
} |
374 |
|
|
|
375 |
|
|
//Debug print. |
376 |
|
✗ |
if(m_verbose > 10) { |
377 |
|
✗ |
PetscPrintf(MpiInfo::petscComm(),"\n\t\t m_W_0 (==m_solEDO (ionicSolver) at time == 0):"); |
378 |
|
✗ |
m_W_0.view(); |
379 |
|
|
} |
380 |
|
|
|
381 |
|
|
|
382 |
|
|
// Read fibers only for testCase == 1. |
383 |
|
✗ |
if (FelisceParam::instance().testCase == 1) { |
384 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->readData(*io()); |
385 |
|
|
} |
386 |
|
|
|
387 |
|
|
//Initialize heterogeneous parameters for schaf and courtemanche solvers. |
388 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "schaf") { |
389 |
|
✗ |
initializeSchafParameter(); |
390 |
|
|
} |
391 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "fhn") { |
392 |
|
✗ |
initializeFhNParameter(); |
393 |
|
✗ |
} else if (FelisceParam::instance().typeOfIonicModel == "court") { |
394 |
|
✗ |
initializeCourtParameter(); |
395 |
|
|
} |
396 |
|
|
} |
397 |
|
|
|
398 |
|
✗ |
void BidomainCurvModel::postAssemblingMatrixRHS(std::size_t iProblem) { |
399 |
|
✗ |
evalIapp(); |
400 |
|
✗ |
if (FelisceParam::instance().hasAppliedExteriorCurrent) |
401 |
|
✗ |
evalIappExt(); |
402 |
|
✗ |
finalizeRHSDP(iProblem); |
403 |
|
|
// Calculate m_RHS = mass * RHS. |
404 |
|
✗ |
m_linearProblemBidomainCurv[iProblem]->addMatrixRHS(); |
405 |
|
|
} |
406 |
|
|
|
407 |
|
✗ |
void BidomainCurvModel::initializeSchafParameter() { |
408 |
|
|
// Initialize TauClose. |
409 |
|
✗ |
m_schaf->fctTauClose().initialize(m_fstransient); |
410 |
|
✗ |
if (FelisceParam::instance().hasHeteroTauClose) { |
411 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_schaf->fctTauClose(), m_schaf->tauClose()); |
412 |
|
|
} |
413 |
|
|
// Initialize TauOut : heterogeneous in case of infarct. |
414 |
|
✗ |
m_schaf->fctTauOut().initialize(m_fstransient); |
415 |
|
✗ |
if (FelisceParam::instance().hasInfarct) { |
416 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_schaf->fctTauOut(), m_schaf->tauOut()); |
417 |
|
|
} |
418 |
|
|
} |
419 |
|
|
|
420 |
|
✗ |
void BidomainCurvModel::initializeFhNParameter() { |
421 |
|
✗ |
if (FelisceParam::instance().hasInfarct) { |
422 |
|
✗ |
m_fhn->fctSpar().initialize(m_fstransient); |
423 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_fhn->fctSpar(), m_fhn->f0Par()); |
424 |
|
|
} |
425 |
|
|
} |
426 |
|
|
|
427 |
|
✗ |
void BidomainCurvModel::initializeCourtParameter() { |
428 |
|
|
// Initialize gToMax |
429 |
|
✗ |
m_court->fctCourtCondIto().initialize(m_fstransient); |
430 |
|
|
// Initialize gCaLMax |
431 |
|
✗ |
m_court->fctCourtCondICaL().initialize(m_fstransient); |
432 |
|
✗ |
if (FelisceParam::instance().hasHeteroCourtPar) { |
433 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_court->fctCourtCondIto(), m_court->courtCondIto()); |
434 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_court->fctCourtCondICaL(), m_court->courtCondICaL()); |
435 |
|
|
} |
436 |
|
✗ |
if (!FelisceParam::instance().hasHeteroCondAtria) { |
437 |
|
|
// Initialize multiplicative coefficient |
438 |
|
✗ |
m_court->fctCourtCondMultCoeff().initialize(m_fstransient); |
439 |
|
|
// Warning: extraData m_linearProblemBidomainCurv[0]->Reference() size has to be numDof (not number of mesh points). |
440 |
|
✗ |
m_linearProblemBidomainCurv[0]->evalFunctionOnDof(m_court->fctCourtCondMultCoeff(), m_court->courtCondMultCoeff(), m_linearProblemBidomainCurv[0]->Reference()); |
441 |
|
|
} |
442 |
|
|
} |
443 |
|
|
|
444 |
|
|
// Pay attention on the call of this : |
445 |
|
|
// call BidomainModel::writeSolution() function instead of (non-virtual) Model::writeSolution() function |
446 |
|
✗ |
void BidomainCurvModel::writeSolution() { |
447 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
448 |
|
✗ |
if (m_meshIsWritten == false) writeMesh(); |
449 |
|
|
} |
450 |
|
|
|
451 |
|
✗ |
if( (m_fstransient->iteration % FelisceParam::instance().frequencyWriteSolution == 0) or m_hasFinished) { |
452 |
|
✗ |
if(FelisceParam::verbose() > 1) PetscPrintf(MpiInfo::petscComm(),"Write solutions\n"); |
453 |
|
✗ |
for (std::size_t ipb = 0; ipb < m_linearProblem.size(); ipb++) { |
454 |
|
|
|
455 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->writeSolution(MpiInfo::rankProc(), m_ios, m_fstransient->time, m_fstransient->iteration); |
456 |
|
|
|
457 |
|
✗ |
if(FelisceParam::instance().printIonicVar) { |
458 |
|
|
int rankProc; |
459 |
|
✗ |
MPI_Comm_rank(MpiInfo::petscComm(),&rankProc); |
460 |
|
✗ |
if ( m_fstransient->iteration == 0 ) { |
461 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomainCurv[ipb]->numDofPerUnknown(0)]; |
462 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->fromVecToDoubleStar(ionicSolToSave, m_W_0, rankProc, 1, m_linearProblemBidomainCurv[ipb]->numDofPerUnknown(0)); |
463 |
|
✗ |
if (rankProc == 0) { |
464 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), "ionicVar"); |
465 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->writeEnsightCase(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,"ionicVar"); |
466 |
|
|
} |
467 |
|
✗ |
delete [] ionicSolToSave; |
468 |
|
|
|
469 |
|
|
} else { |
470 |
|
✗ |
double * ionicSolToSave = new double[m_linearProblemBidomainCurv[ipb]->numDofPerUnknown(0)]; |
471 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->fromVecToDoubleStar(ionicSolToSave, m_ionic->sol(), rankProc, 1, m_linearProblemBidomainCurv[ipb]->numDofPerUnknown(0)); |
472 |
|
✗ |
if (rankProc == 0) { |
473 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->writeEnsightScalar(ionicSolToSave, static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution), "ionicVar"); |
474 |
|
✗ |
m_linearProblemBidomainCurv[ipb]->writeEnsightCase(static_cast<int>(m_fstransient->iteration/FelisceParam::instance().frequencyWriteSolution)+1,"ionicVar"); |
475 |
|
|
} |
476 |
|
✗ |
delete [] ionicSolToSave; |
477 |
|
|
} |
478 |
|
|
} |
479 |
|
✗ |
if (MpiInfo::rankProc() == 0) { |
480 |
|
✗ |
for(std::size_t iio=0; iio<m_ios.size(); ++iio) { |
481 |
|
✗ |
if ( m_ios[iio]->typeOutput == 1 ) // 1 : ensight{ |
482 |
|
✗ |
m_ios[iio]->postProcess(m_fstransient->time/*, !m_sameGeometricMeshForVariable*/); |
483 |
|
|
} |
484 |
|
|
} |
485 |
|
|
} |
486 |
|
|
} |
487 |
|
|
} |
488 |
|
|
|
489 |
|
✗ |
void BidomainCurvModel::forward() { |
490 |
|
|
//Write solution with ensight. |
491 |
|
✗ |
BidomainCurvModel::writeSolution(); |
492 |
|
|
//Advance time step. |
493 |
|
✗ |
updateTime(FlagMatrixRHS::only_rhs); |
494 |
|
✗ |
printNewTimeIterationBanner(); |
495 |
|
|
|
496 |
|
✗ |
if ( m_fstransient->iteration == 1) { |
497 |
|
✗ |
m_bdfEdp.initialize(m_linearProblemBidomainCurv[0]->solution()); |
498 |
|
✗ |
m_linearProblemBidomainCurv[0]->initializeTimeScheme(&m_bdfEdp); |
499 |
|
|
//m_extrapolate = initial solution. |
500 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
501 |
|
|
|
502 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "court") { |
503 |
|
✗ |
m_court->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomainCurv[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomainCurv[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomainCurv[0]->ao()); |
504 |
|
✗ |
m_court->initializeExtrap(m_extrapolate); |
505 |
|
✗ |
m_court->initializeConcentrationsAndGateConditions(m_m, m_h, m_j, m_ao, m_io, m_ua, m_ui, m_xr, m_xs, m_d, m_f, m_fca, m_urel, m_vrel, m_wrel, m_nai, m_nao, m_cao, m_ki, m_ko, m_cai, m_naiont, m_kiont, m_caiont, m_ileak, m_iup, m_itr, m_irel, m_cmdn, m_trpn, m_nsr, m_jsr, m_csqn); |
506 |
|
✗ |
m_buildCourt = true; |
507 |
|
|
} |
508 |
|
|
else { |
509 |
|
✗ |
m_ionic->defineSizeAndMappingOfIonicProblem(m_linearProblemBidomainCurv[0]->numDofLocalPerUnknown(potTransMemb), m_linearProblemBidomainCurv[0]->mappingDofLocalToDofGlobal(potTransMemb), m_linearProblemBidomainCurv[0]->ao()); |
510 |
|
✗ |
m_ionic->initializeExtrap(m_extrapolate); |
511 |
|
✗ |
m_ionic->initialize(m_W_0); |
512 |
|
✗ |
m_buildIonic = true; |
513 |
|
|
} |
514 |
|
✗ |
if (FelisceParam::instance().dataAssimilation) { |
515 |
|
✗ |
m_linearProblemBidomainCurv[0]->gatherSolution(); |
516 |
|
✗ |
m_linearProblemBidomainCurv[0]->initSolutionTimeBef(); |
517 |
|
✗ |
m_linearProblemBidomainCurv[0]->readDataForDA(*io(),m_fstransient->iteration); |
518 |
|
✗ |
m_linearProblemBidomainCurv[0]->initLevelSet(); |
519 |
|
✗ |
m_linearProblemBidomainCurv[0]->levelSet(FelisceParam::instance().minLSPotTM,FelisceParam::instance().maxLSPotTM); |
520 |
|
|
} |
521 |
|
|
} |
522 |
|
|
else { |
523 |
|
✗ |
m_bdfEdp.update(m_linearProblemBidomainCurv[0]->solution()); |
524 |
|
✗ |
m_bdfEdp.extrapolate(m_extrapolate); |
525 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "court") { |
526 |
|
✗ |
m_court->updateExtrap(m_extrapolate); |
527 |
|
|
} |
528 |
|
|
else { |
529 |
|
✗ |
m_ionic->update(m_extrapolate); |
530 |
|
|
} |
531 |
|
|
} |
532 |
|
|
|
533 |
|
✗ |
if (FelisceParam::instance().dataAssimilation) { |
534 |
|
✗ |
m_linearProblemBidomainCurv[0]->gatherSolution(); |
535 |
|
✗ |
m_linearProblemBidomainCurv[0]->solutionTimeBef(); |
536 |
|
✗ |
m_linearProblemBidomainCurv[0]->readDataForDA(*io(),m_fstransient->iteration); |
537 |
|
✗ |
m_linearProblemBidomainCurv[0]->levelSet(FelisceParam::instance().minLSPotTM,FelisceParam::instance().maxLSPotTM); |
538 |
|
|
} |
539 |
|
|
|
540 |
|
✗ |
if (FelisceParam::instance().typeOfIonicModel == "court") { |
541 |
|
✗ |
m_court->computeIon(); |
542 |
|
|
} |
543 |
|
|
else { |
544 |
|
✗ |
m_ionic->computeRHS(); |
545 |
|
✗ |
m_ionic->solveEDO(); |
546 |
|
✗ |
m_ionic->computeIon(); |
547 |
|
✗ |
m_ionic->updateBdf(); |
548 |
|
|
} |
549 |
|
|
|
550 |
|
✗ |
m_bdfEdp.computeRHSTime(m_fstransient->timeStep); |
551 |
|
|
|
552 |
|
|
//double timeDA =fmod(m_fstransient->iteration,1); |
553 |
|
|
|
554 |
|
|
// Assemble Matrix of linearProblem. |
555 |
|
✗ |
if ( m_fstransient->iteration == 1 ) { |
556 |
|
✗ |
m_linearProblemBidomainCurv[0]->assembleMatrixRHSBD(MpiInfo::rankProc()); |
557 |
|
✗ |
} else if (FelisceParam::instance().dataAssimilation) { //&&(timeDA == 0)) |
558 |
|
✗ |
m_linearProblemBidomainCurv[0]->assembleMatrixRHSBD(MpiInfo::rankProc(), FlagMatrixRHS::only_rhs); |
559 |
|
|
} |
560 |
|
|
|
561 |
|
|
//Specific operations before solve the system : add static matrix _Matrix[1] to the dynamic matrix of the system _Matrix. |
562 |
|
✗ |
postAssemblingMatrixRHS(); |
563 |
|
✗ |
m_linearProblemBidomainCurv[0]->solve(MpiInfo::rankProc(), MpiInfo::numProc()); |
564 |
|
|
|
565 |
|
|
} |
566 |
|
|
} |
567 |
|
|
|