LifeV
ElectroIonicModel.cpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief Base class for ionic models
30 
31  @date 01-2013
32  @author Simone Rossi <simone.rossi@epfl.ch>
33 
34  @contributors
35  @mantainer Simone Rossi <simone.rossi@epfl.ch>
36  @last update 04 - 2014
37  */
38 
39 #include <lifev/electrophysiology/solver/IonicModels/ElectroIonicModel.hpp>
40 
41 
42 namespace LifeV
43 {
44 // ===================================================
45 //! Constructors
46 // ===================================================
47 ElectroIonicModel::ElectroIonicModel() :
48  M_numberOfEquations (0),
49  M_numberOfGatingVariables (0),
50  M_restingConditions (),
51  M_membraneCapacitance (1.),
52  M_appliedCurrent (0.),
53  M_appliedCurrentPtr(),
54  M_pacingProtocol ()
55 {
56 }
57 
58 ElectroIonicModel::ElectroIonicModel ( int n ) :
59  M_numberOfEquations (n),
60  M_numberOfGatingVariables (0),
61  M_restingConditions (std::vector<Real> (M_numberOfEquations, 0.0) ),
62  M_membraneCapacitance (1.),
63  M_appliedCurrent (0.),
64  M_appliedCurrentPtr(),
65  M_pacingProtocol ()
66 {
67 }
68 
69 ElectroIonicModel::ElectroIonicModel ( int n, int g ) :
70  M_numberOfEquations (n),
71  M_numberOfGatingVariables (g),
72  M_restingConditions (std::vector<Real> (M_numberOfEquations, 0.0) ),
73  M_membraneCapacitance (1.),
74  M_appliedCurrent (0.),
75  M_appliedCurrentPtr(),
76  M_pacingProtocol ()
77 {
78 }
79 
80 ElectroIonicModel::ElectroIonicModel ( const ElectroIonicModel& Ionic ) :
81  M_numberOfEquations ( Ionic.Size() ),
82  M_numberOfGatingVariables ( Ionic.numberOfGatingVariables() ),
83  M_restingConditions ( Ionic.restingConditions() ),
84  M_membraneCapacitance ( Ionic.M_membraneCapacitance ),
85  M_appliedCurrent ( Ionic.M_appliedCurrent ),
86  M_pacingProtocol (Ionic.M_pacingProtocol)
87 {
88  if (Ionic.M_appliedCurrentPtr)
89  {
90  M_appliedCurrentPtr.reset (new vector_Type (*Ionic.M_appliedCurrentPtr) );
91  }
92 }
93 
94 // ===================================================
95 //! Methods
96 // ===================================================
97 ElectroIonicModel& ElectroIonicModel::operator = ( const ElectroIonicModel& Ionic )
98 {
99  M_numberOfEquations = Ionic.M_numberOfEquations;
100  M_numberOfGatingVariables = Ionic.M_numberOfGatingVariables;
101  M_restingConditions = Ionic.M_restingConditions;
102  M_membraneCapacitance = Ionic.M_membraneCapacitance;
103  M_appliedCurrent = Ionic.M_appliedCurrent;
104  if (Ionic.M_appliedCurrent)
105  {
106  M_appliedCurrentPtr = Ionic.M_appliedCurrentPtr;
107  }
108  M_pacingProtocol = Ionic.M_pacingProtocol;
109 
110  return *this;
111 }
112 
113 std::vector< std::vector<Real> > ElectroIonicModel::getJac (const std::vector<Real>& v, Real h)
114 {
115  std::vector< std::vector<Real> > J ( M_numberOfEquations, std::vector<Real> (M_numberOfEquations, 0.0) );
116  std::vector<Real> f1 (M_numberOfEquations, 0.0);
117  std::vector<Real> f2 (M_numberOfEquations, 0.0);
118  std::vector<Real> y1 (M_numberOfEquations, 0.0);
119  std::vector<Real> y2 (M_numberOfEquations, 0.0);
120 
121  for (int i = 0; i < M_numberOfEquations; i++)
122  {
123  for (int j = 0; j < M_numberOfEquations; j++)
124  {
125  y1[j] = v[j] + ( (double) (i == j) ) * h;
126  y2[j] = v[j] - ( (double) (i == j) ) * h;
127  }
128  this->computeRhs (y1, f1);
129  f1[0] += M_appliedCurrent;
130  this->computeRhs (y2, f2);
131  f2[0] += M_appliedCurrent;
132 
133  for (int j = 0; j < M_numberOfEquations; j++)
134  {
135  J[j][i] = (f1[j] - f2[j]) / (2.0 * h);
136  }
137  }
138 
139  return J;
140 }
141 
142 MatrixEpetra<Real> ElectroIonicModel::getJac (const vector_Type& v, Real /*h*/)
143 {
144  matrix_Type J (v.map(), M_numberOfEquations, false);
145  // vector<Real> f1(M_numberOfEquations,0.0);
146  // vector<Real> f2(M_numberOfEquations,0.0);
147  // vector< vector<Real> > df ( M_numberOfEquations, vector<Real> (M_numberOfEquations,0.0) );
148  // vector<Real> y1(M_numberOfEquations,0.0);
149  // vector<Real> y2(M_numberOfEquations,0.0);
150  // const Int* k = v.blockMap().MyGlobalElements();
151  //
152  // int* Indices = new int[M_numberOfEquations];
153  // double* Values = new double[M_numberOfEquations];
154  // for(int i=0; i<M_numberOfEquations; i++)
155  // Indices[i] = i;
156  //
157  // for(int i=0; i<M_numberOfEquations; i++)
158  // {
159  // for(int j=0; j<M_numberOfEquations; j++)
160  // {
161  // y1[j] = v[ k[j] ] + ((double)(i==j))*h;
162  // y2[j] = v[ k[j] ] - ((double)(i==j))*h;
163  // }
164  // this->computeRhs(y1, f1);
165  // this->computeRhs(y2, f2);
166  //
167  // for(int j=0; j<M_numberOfEquations; j++)
168  // df[j][i] = (f1[j]-f2[j])/(2.0*h);
169  // }
170  //
171  //
172  // for(int i=0; i<M_numberOfEquations; i++)
173  // {
174  // for(int j=0; j<M_numberOfEquations; j++)
175  // Values[j] = df[i][j];
176  // J.matrixPtr()->InsertGlobalValues (i, M_numberOfEquations, Values, Indices);
177  // }
178  //
179  // J.globalAssemble();
180  //
181  // delete[] Indices;
182  // delete[] Values;
183 
184  return J;
185 }
186 
187 
188 void ElectroIonicModel::computeGatingRhs ( const std::vector<vectorPtr_Type>& v,
189  std::vector<vectorPtr_Type>& rhs )
190 {
191 
192  int nodes = ( * (v.at (1) ) ).epetraVector().MyLength();
193 
194 
195  std::vector<Real> localVec ( M_numberOfEquations, 0.0 );
196  std::vector<Real> localRhs ( M_numberOfEquations - 1, 0.0 );
197 
198  int j (0);
199 
200  for ( int k = 0; k < nodes; k++ )
201  {
202 
203  j = ( * (v.at (1) ) ).blockMap().GID (k);
204 
205  for ( int i = 0; i < M_numberOfEquations; i++ )
206  {
207  localVec.at (i) = ( * ( v.at (i) ) ) [j];
208  }
209 
210  if (M_appliedCurrentPtr)
211  {
212  M_appliedCurrent = (*M_appliedCurrentPtr) [j];
213  }
214  else
215  {
216  M_appliedCurrent = 0.0;
217  }
218 
219  computeGatingRhs ( localVec, localRhs );
220 
221  for ( int i = 1; i < M_numberOfEquations; i++ )
222  {
223  ( * ( rhs.at (i) ) ) [j] = localRhs.at (i - 1);
224  }
225 
226  }
227 
228 }
229 
230 void ElectroIonicModel::computeNonGatingRhs ( const std::vector<vectorPtr_Type>& v,
231  std::vector<vectorPtr_Type>& rhs )
232 {
233 
234  int nodes = ( * (v.at (1) ) ).epetraVector().MyLength();
235 
236 
237  std::vector<Real> localVec ( M_numberOfEquations, 0.0 );
238  int offset = 1 + M_numberOfGatingVariables;
239  std::vector<Real> localRhs ( M_numberOfEquations - offset, 0.0 );
240 
241  int j (0);
242 
243  for ( int k = 0; k < nodes; k++ )
244  {
245 
246  j = ( * (v.at (1) ) ).blockMap().GID (k);
247 
248  for ( int i = 0; i < M_numberOfEquations; i++ )
249  {
250  localVec.at (i) = ( * ( v.at (i) ) ) [j];
251  }
252 
253  if (M_appliedCurrentPtr)
254  {
255  M_appliedCurrent = (*M_appliedCurrentPtr) [j];
256  }
257  else
258  {
259  M_appliedCurrent = 0.0;
260  }
261 
262  computeNonGatingRhs ( localVec, localRhs );
263 
264  for ( int i = offset; i < M_numberOfEquations; i++ )
265  {
266  ( * ( rhs.at (i) ) ) [j] = localRhs.at (i - offset);
267  }
268 
269  }
270 
271 }
272 
273 
274 void ElectroIonicModel::computeRhs ( const std::vector<vectorPtr_Type>& v,
275  std::vector<vectorPtr_Type>& rhs )
276 {
277 
278  int nodes = ( * (v.at (1) ) ).epetraVector().MyLength();
279 
280 
281  std::vector<Real> localVec ( M_numberOfEquations, 0.0 );
282  std::vector<Real> localRhs ( M_numberOfEquations, 0.0 );
283 
284  int j (0);
285 
286  for ( int k = 0; k < nodes; k++ )
287  {
288 
289  j = ( * (v.at (1) ) ).blockMap().GID (k);
290 
291  for ( int i = 0; i < M_numberOfEquations; i++ )
292  {
293  localVec.at (i) = ( * ( v.at (i) ) ) [j];
294  }
295 
296 
297  if (M_appliedCurrentPtr)
298  {
299  M_appliedCurrent = (*M_appliedCurrentPtr) [j];
300  }
301  else
302  {
303  M_appliedCurrent = 0.0;
304  }
305  computeRhs ( localVec, localRhs );
306  addAppliedCurrent (localRhs);
307 
308  for ( int i = 0; i < M_numberOfEquations; i++ )
309  {
310  ( * ( rhs.at (i) ) ) [j] = localRhs.at (i);
311  }
312 
313  }
314 
315 }
316 
317 void ElectroIonicModel::computePotentialRhsICI ( const std::vector<vectorPtr_Type>& v,
318  std::vector<vectorPtr_Type>& rhs,
319  matrix_Type& massMatrix )
320 {
321  int nodes = ( * (v.at (0) ) ).epetraVector().MyLength();
322 
323  ( * ( rhs.at (0) ) ) *= 0.0;
324  std::vector<Real> localVec ( M_numberOfEquations, 0.0 );
325 
326  int j (0);
327 
328  for ( int k = 0; k < nodes; k++ )
329  {
330 
331  j = ( * (v.at (0) ) ).blockMap().GID (k);
332 
333  for ( int i = 0; i < M_numberOfEquations; i++ )
334  {
335  localVec.at (i) = ( * ( v.at (i) ) ) [j];
336  }
337 
338  if (M_appliedCurrentPtr)
339  {
340  M_appliedCurrent = (*M_appliedCurrentPtr) [j];
341  }
342  else
343  {
344  M_appliedCurrent = 0.0;
345  }
346 
347  ( * ( rhs.at (0) ) ) [j] = computeLocalPotentialRhs ( localVec ) + M_appliedCurrent;
348 
349  }
350 
351  ( * ( rhs.at (0) ) ) = massMatrix * ( * ( rhs.at (0) ) );
352 
353 }
354 
355 
356 void ElectroIonicModel::computePotentialRhsSVI ( const std::vector<vectorPtr_Type>& v,
357  std::vector<vectorPtr_Type>& rhs,
358  FESpace<mesh_Type, MapEpetra>& uFESpace )
359 {
360 
361  std::vector<Real> U (M_numberOfEquations, 0.0);
362  Real I (0.0);
363  ( * ( rhs.at (0) ) ) *= 0.0;
364 
365  std::vector<vectorPtr_Type> URepPtr;
366  for ( int k = 0; k < M_numberOfEquations; k++ )
367  {
368  URepPtr.push_back ( vectorPtr_Type ( new VectorEpetra ( * ( v.at (k) ) , Repeated ) ) );
369  }
370 
371  VectorEpetra IappRep ( *M_appliedCurrentPtr, Repeated );
372 
373  std::vector<elvecPtr_Type> elvecPtr;
374  for ( int k = 0; k < M_numberOfEquations; k++ )
375  {
376  elvecPtr.push_back ( elvecPtr_Type ( new VectorElemental ( uFESpace.fe().nbFEDof(), 1 ) ) );
377  }
378 
379  VectorElemental elvec_Iapp ( uFESpace.fe().nbFEDof(), 1 );
380  VectorElemental elvec_Iion ( uFESpace.fe().nbFEDof(), 1 );
381 
382  for (UInt iVol = 0; iVol < uFESpace.mesh()->numVolumes(); ++iVol)
383  {
384 
385  uFESpace.fe().updateJacQuadPt ( uFESpace.mesh()->volumeList ( iVol ) );
386 
387 
388  for ( int k = 0; k < M_numberOfEquations; k++ )
389  {
390  ( * ( elvecPtr.at (k) ) ).zero();
391  }
392  elvec_Iapp.zero();
393  elvec_Iion.zero();
394 
395  UInt eleIDu = uFESpace.fe().currentLocalId();
396  UInt nbNode = ( UInt ) uFESpace.fe().nbFEDof();
397 
398  //! Filling local elvec_u with potential values in the nodes
399  for ( UInt iNode = 0 ; iNode < nbNode ; iNode++ )
400  {
401 
402  Int ig = uFESpace.dof().localToGlobalMap ( eleIDu, iNode );
403 
404  for ( int k = 0; k < M_numberOfEquations; k++ )
405  {
406  ( * ( elvecPtr.at (k) ) ).vec() [iNode] = ( * ( URepPtr.at (k) ) ) [ig];
407  }
408 
409  elvec_Iapp.vec() [ iNode ] = IappRep[ig];
410 
411  }
412 
413  //compute the local vector
414  for ( UInt ig = 0; ig < uFESpace.fe().nbQuadPt(); ig++ )
415  {
416 
417  for ( int k = 0; k < M_numberOfEquations; k++ )
418  {
419  U.at (k) = 0;
420  }
421  I = 0;
422 
423  for ( UInt i = 0; i < uFESpace.fe().nbFEDof(); i++ )
424  {
425 
426  for ( int k = 0; k < M_numberOfEquations; k++ )
427  {
428  U.at (k) += ( * ( elvecPtr.at (k) ) ) (i) * uFESpace.fe().phi ( i, ig );
429  }
430 
431  I += elvec_Iapp (i) * uFESpace.fe().phi ( i, ig );
432 
433  }
434 
435  for ( UInt i = 0; i < uFESpace.fe().nbFEDof(); i++ )
436  {
437 
438  elvec_Iion ( i ) += ( computeLocalPotentialRhs (U) + I ) * uFESpace.fe().phi ( i, ig ) * uFESpace.fe().weightDet ( ig );
439 
440  }
441 
442  }
443 
444  //assembly
445  for ( UInt i = 0 ; i < uFESpace.fe().nbFEDof(); i++ )
446  {
447  Int ig = uFESpace.dof().localToGlobalMap ( eleIDu, i );
448  ( * ( rhs.at (0) ) ).sumIntoGlobalValues (ig, elvec_Iion.vec() [i] );
449  }
450  }
451  rhs.at (0) -> globalAssemble();
452 
453 
454  if (M_appliedCurrentPtr)
455  {
456  M_appliedCurrentPtr -> setMapType (Unique);
457  }
458 
459 }
460 
461 void ElectroIonicModel::computePotentialRhsSVI ( const std::vector<vectorPtr_Type>& v,
462  std::vector<vectorPtr_Type>& rhs,
463  FESpace<mesh_Type, MapEpetra>& uFESpace,
464  const QuadratureRule& qr)
465 {
466  uFESpace.setQuadRule ( qr );
467  computePotentialRhsSVI (v, rhs, uFESpace);
468 }
469 
470 
471 void ElectroIonicModel::computeGatingVariablesWithRushLarsen ( std::vector<vectorPtr_Type>& v, const Real dt )
472 {
473  int nodes = ( * (v.at (0) ) ).epetraVector().MyLength();
474 
475  std::vector<Real> localVec ( M_numberOfEquations, 0.0 );
476 
477  int j (0);
478 
479  for ( int k = 0; k < nodes; k++ )
480  {
481 
482  j = ( * (v.at (0) ) ).blockMap().GID (k);
483 
484  for ( int i = 0; i < M_numberOfEquations; i++ )
485  {
486  localVec.at (i) = ( * ( v.at (i) ) ) [j];
487  }
488 
489  computeGatingVariablesWithRushLarsen (localVec, dt);
490 
491  for ( int i = 0; i < M_numberOfEquations; i++ )
492  {
493  ( * ( v.at (i) ) ) [j] = localVec.at (i);
494  }
495 
496  }
497 
498 }
499 
500 
501 void ElectroIonicModel::initialize ( std::vector<Real>& v )
502 {
503  for (int i (0); i < M_numberOfEquations; i++ )
504  {
505  v.at (i) = M_restingConditions.at (i);
506  }
507 }
508 
509 
510 void ElectroIonicModel::initialize ( std::vector<vectorPtr_Type>& v )
511 {
512  for (int i (0); i < M_numberOfEquations; i++ )
513  {
514  * ( v.at (i) ) = M_restingConditions.at (i);
515  }
516 }
517 
518 
519 }
LifeV::VectorEpetra
VectorEpetra - The Epetra Vector format Wrapper.
Definition: VectorEpetra.hpp:62
LifeV::MatrixEpetra::importFromHDF5
void importFromHDF5(std::string const &fileName, std::string const &matrixName="matrix")
Read a matrix from a HDF5 (.h5) file.
Definition: MatrixEpetra.hpp:1394
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::MapEpetra::importer
Epetra_Import const & importer()
Getter for the Epetra_Import.
Definition: MapEpetra.cpp:394
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::QuadratureRule
QuadratureRule - The basis class for storing and accessing quadrature rules.
Definition: QuadratureRule.hpp:111