Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Bidomain.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File Bidomain.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Bidomain cardiac electrophysiology homogenised model.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 
38 #include <CardiacEPSolver/EquationSystems/Bidomain.h>
39 #include <CardiacEPSolver/Filters/FilterCheckpointCellModel.h>
40 
41 namespace Nektar
42 {
43  /**
44  * @class Bidomain
45  *
46  * Base model of cardiac electrophysiology of the form
47  * \f{align*}{
48  * \frac{\partial u}{\partial t} = \nabla^2 u + J_{ion},
49  * \f}
50  * where the reaction term, \f$J_{ion}\f$ is defined by a specific cell
51  * model.
52  *
53  * This implementation, at present, treats the reaction terms explicitly
54  * and the diffusive element implicitly.
55  */
56 
57  /**
58  * Registers the class with the Factory.
59  */
60  string Bidomain::className
61  = GetEquationSystemFactory().RegisterCreatorFunction(
62  "Bidomain",
63  Bidomain::create,
64  "Bidomain model of cardiac electrophysiology with 3D diffusion.");
65 
66 
67  /**
68  *
69  */
70  Bidomain::Bidomain(
71  const LibUtilities::SessionReaderSharedPtr& pSession)
72  : UnsteadySystem(pSession)
73  {
74  }
75 
76  void Bidomain::v_InitObject()
77  {
78  UnsteadySystem::v_InitObject();
79  m_session->LoadParameter("Chi", m_chi);
80  m_session->LoadParameter("Cm", m_capMembrane);
81 
82  std::string vCellModel;
83  m_session->LoadSolverInfo("CELLMODEL", vCellModel, "");
84 
85  ASSERTL0(vCellModel != "", "Cell Model not specified.");
86 
87  m_cell = GetCellModelFactory().CreateInstance(vCellModel, m_session, m_fields[0]);
88  m_intVariables.push_back(0);
89  m_intVariables.push_back(1);
90 
91  // Load variable coefficients
92  StdRegions::VarCoeffType varCoeffEnum[3] = {
93  StdRegions::eVarCoeffD00,
94  StdRegions::eVarCoeffD11,
95  StdRegions::eVarCoeffD22
96  };
97  std::string varName[3] = {
98  "AnisotropicConductivityX",
99  "AnisotropicConductivityY",
100  "AnisotropicConductivityZ"
101  };
102 
103 
104  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
105  {
106  for (int i = 0; i < m_spacedim; ++i)
107  {
108  int nq = m_fields[0]->GetNpoints();
109  Array<OneD,NekDouble> x0(nq);
110  Array<OneD,NekDouble> x1(nq);
111  Array<OneD,NekDouble> x2(nq);
112 
113  // get the coordinates
114  m_fields[0]->GetCoords(x0,x1,x2);
115  tmp1 = Array<OneD, const Array<OneD, NekDouble> >(nq);
116  tmp2 = Array<OneD, const Array<OneD, NekDouble> >(nq);
117  tmp3 = Array<OneD, const Array<OneD, NekDouble> >(nq);
118  tmp1[i] = Array<OneD, NekDouble>(nq);
119  tmp2[i] = Array<OneD, NekDouble>(nq);
120  tmp3[i] = Array<OneD, NekDouble>(nq);
121 
122  LibUtilities::EquationSharedPtr ifunc1
123  = m_session->GetFunction("IntracellularConductivity", varName[i]);
124  LibUtilities::EquationSharedPtr ifunc2
125  = m_session->GetFunction("ExtracellularConductivity", varName[i]);
126  for(int j = 0; j < nq; j++)
127  {
128  tmp1[i][j] = ifunc1->Evaluate(x0[j],x1[j],x2[j],0.0);
129  tmp2[i][j] = ifunc2->Evaluate(x0[j],x1[j],x2[j],0.0);
130  }
131  Vmath::Vadd(nq, tmp1[i], 1, tmp2[i], 1, tmp3[i], 1);
132  m_vardiffi[varCoeffEnum[i]] = tmp1[i];
133  m_vardiffie[varCoeffEnum[i]] = tmp3[i];
134  }
135  }
136 
137 
138  if (m_session->DefinesParameter("StimulusDuration"))
139  {
140  ASSERTL0(m_session->DefinesFunction("Stimulus", "u"),
141  "Stimulus function not defined.");
142  m_session->LoadParameter("StimulusDuration", m_stimDuration);
143  }
144  else
145  {
146  m_stimDuration = 0;
147  }
148 
149 
150  // Search through the loaded filters and pass the cell model to any
151  // CheckpointCellModel filters loaded.
152  int k = 0;
153  const LibUtilities::FilterMap& f = m_session->GetFilters();
154  LibUtilities::FilterMap::const_iterator x;
155  for (x = f.begin(); x != f.end(); ++x, ++k)
156  {
157  if (x->first == "CheckpointCellModel")
158  {
159  boost::shared_ptr<FilterCheckpointCellModel> c
160  = boost::dynamic_pointer_cast<FilterCheckpointCellModel>(
161  m_filters[k]);
162  c->SetCellModel(m_cell);
163  }
164  }
165 
166  if (!m_explicitDiffusion)
167  {
168  m_ode.DefineImplicitSolve (&Bidomain::DoImplicitSolve, this);
169  }
170  m_ode.DefineOdeRhs(&Bidomain::DoOdeRhs, this);
171  }
172 
173 
174  /**
175  *
176  */
177  Bidomain::~Bidomain()
178  {
179 
180  }
181 
182 
183  /**
184  * @param inarray Input array.
185  * @param outarray Output array.
186  * @param time Current simulation time.
187  * @param lambda Timestep.
188  */
189  void Bidomain::DoImplicitSolve(
190  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
191  Array<OneD, Array<OneD, NekDouble> >&outarray,
192  const NekDouble time,
193  const NekDouble lambda)
194  {
195  int nvariables = inarray.num_elements();
196  int nq = m_fields[0]->GetNpoints();
197 
198  Array<OneD, NekDouble> grad0(nq), grad1(nq), grad2(nq), grad(nq);
199  Array<OneD, NekDouble> ggrad0(nq), ggrad1(nq), ggrad2(nq), ggrad(nq), temp(nq);
200 
201  // We solve ( \sigma\nabla^2 - HHlambda ) Y[i] = rhs [i]
202  // inarray = input: \hat{rhs} -> output: \hat{Y}
203  // outarray = output: nabla^2 \hat{Y}
204  // where \hat = modal coeffs
205  for (int i = 0; i < nvariables; ++i)
206  {
207  // Only apply diffusion to first variable.
208  if (i > 1) {
209  Vmath::Vcopy(nq, &inarray[i][0], 1, &outarray[i][0], 1);
210  continue;
211  }
212  if (i == 0) {
213  StdRegions::ConstFactorMap factors;
214  factors[StdRegions::eFactorLambda] = (1.0/lambda)*(m_capMembrane*m_chi);
215  if (m_spacedim==1) {
216  // Take first partial derivative
217  m_fields[i]->PhysDeriv(inarray[1],ggrad0);
218  // Take second partial derivative
219  m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
220  // Multiply by Intracellular-Conductivity
221  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
222  {
223  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), ggrad0, 1, ggrad0, 1);
224  }
225  // Add partial derivatives together
226  Vmath::Vcopy(nq, ggrad0, 1, ggrad, 1);
227  Vmath::Smul(nq, -1.0, ggrad, 1, ggrad, 1);
228  // Multiply 1.0/timestep/lambda
229  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1, temp, 1);
230  Vmath::Vadd(nq, ggrad, 1, temp, 1, m_fields[i]->UpdatePhys(), 1);
231  // Solve a system of equations with Helmholtz solver and transform
232  // back into physical space.
233  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(),NullFlagList,factors);
234  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
235  m_fields[i]->SetPhysState(true);
236  // Copy the solution vector (required as m_fields must be set).
237  outarray[i] = m_fields[i]->GetPhys();
238  }
239 
240  if (m_spacedim==2) {
241  // Take first partial derivative
242  m_fields[i]->PhysDeriv(inarray[1],ggrad0,ggrad1);
243  // Take second partial derivative
244  m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
245  m_fields[i]->PhysDeriv(1,ggrad1,ggrad1);
246  // Multiply by Intracellular-Conductivity
247  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
248  {
249  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), ggrad0, 1, ggrad0, 1);
250  Vmath::Smul(nq, m_session->GetParameter("sigmaiy"), ggrad1, 1, ggrad1, 1);
251  }
252  // Add partial derivatives together
253  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
254  Vmath::Smul(nq, -1.0, ggrad, 1, ggrad, 1);
255  // Multiply 1.0/timestep/lambda
256  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1, temp, 1);
257  Vmath::Vadd(nq, ggrad, 1, temp, 1, m_fields[i]->UpdatePhys(), 1);
258  // Solve a system of equations with Helmholtz solver and transform
259  // back into physical space.
260  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(),NullFlagList,factors,m_vardiffi);
261  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
262  m_fields[i]->SetPhysState(true);
263  // Copy the solution vector (required as m_fields must be set).
264  outarray[i] = m_fields[i]->GetPhys();
265  }
266 
267  if (m_spacedim==3) {
268  // Take first partial derivative
269  m_fields[i]->PhysDeriv(inarray[1],ggrad0,ggrad1,ggrad2);
270  // Take second partial derivative
271  m_fields[i]->PhysDeriv(0,ggrad0,ggrad0);
272  m_fields[i]->PhysDeriv(1,ggrad1,ggrad1);
273  m_fields[i]->PhysDeriv(2,ggrad2,ggrad2);
274  // Multiply by Intracellular-Conductivity
275  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
276  {
277  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), ggrad0, 1, ggrad0, 1);
278  Vmath::Smul(nq, m_session->GetParameter("sigmaiy"), ggrad1, 1, ggrad1, 1);
279  Vmath::Smul(nq, m_session->GetParameter("sigmaiz"), ggrad2, 1, ggrad2, 1);
280  }
281  // Add partial derivatives together
282  Vmath::Vadd(nq, ggrad0, 1, ggrad1, 1, ggrad, 1);
283  Vmath::Vadd(nq, ggrad2, 1, ggrad, 1, ggrad, 1);
284  Vmath::Smul(nq, -1.0, ggrad, 1, ggrad, 1);
285  // Multiply 1.0/timestep/lambda
286  Vmath::Smul(nq, -factors[StdRegions::eFactorLambda], inarray[i], 1, temp, 1);
287  Vmath::Vadd(nq, ggrad, 1, temp, 1, m_fields[i]->UpdatePhys(), 1);
288  // Solve a system of equations with Helmholtz solver and transform
289  // back into physical space.
290  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(),NullFlagList,factors,m_vardiffi);
291  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
292  m_fields[i]->SetPhysState(true);
293  // Copy the solution vector (required as m_fields must be set).
294  outarray[i] = m_fields[i]->GetPhys();
295  }
296 
297  }
298  if (i == 1) {
299  StdRegions::ConstFactorMap factors;
300  factors[StdRegions::eFactorLambda] = 0.0;
301  if (m_spacedim==1) {
302  // Take first partial derivative
303  m_fields[i]->PhysDeriv(m_fields[0]->UpdatePhys(),grad0);
304  // Take second derivative
305  m_fields[i]->PhysDeriv(0,grad0,grad0);
306  // Multiply by Intracellular-Conductivity
307  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
308  {
309  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), grad0, 1, grad0, 1);
310  }
311  // and sum terms
312  Vmath::Vcopy(nq, grad0, 1, grad, 1);
313  Vmath::Smul(nq, (-1.0*m_session->GetParameter("sigmaix"))/(m_session->GetParameter("sigmaix")+m_session->GetParameter("sigmaix")), grad, 1, grad, 1);
314  // Now solve Poisson problem for \phi_e
315  m_fields[i]->SetPhys(grad);
316  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(), NullFlagList, factors);
317  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
318  m_fields[i]->SetPhysState(true);
319  // Copy the solution vector (required as m_fields must be set).
320  outarray[i] = m_fields[i]->GetPhys();
321  }
322 
323  if (m_spacedim==2) {
324  // Take first partial derivative
325  m_fields[i]->PhysDeriv(m_fields[0]->UpdatePhys(),grad0,grad1);
326  // Take second derivative
327  m_fields[i]->PhysDeriv(0,grad0,grad0);
328  m_fields[i]->PhysDeriv(1,grad1,grad1);
329  // Multiply by Intracellular-Conductivity
330  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
331  {
332  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), grad0, 1, grad0, 1);
333  Vmath::Smul(nq, m_session->GetParameter("sigmaiy"), grad1, 1, grad1, 1);
334  }
335  // and sum terms
336  Vmath::Vadd(nq, grad0, 1, grad1, 1, grad, 1);
337  Vmath::Smul(nq, -1.0, grad, 1, grad, 1);
338  // Now solve Poisson problem for \phi_e
339  m_fields[i]->SetPhys(grad);
340  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(), NullFlagList, factors, m_vardiffie);
341  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
342  m_fields[i]->SetPhysState(true);
343  // Copy the solution vector (required as m_fields must be set).
344  outarray[i] = m_fields[i]->GetPhys();
345  }
346 
347  if (m_spacedim==3) {
348  // Take first partial derivative
349  m_fields[i]->PhysDeriv(m_fields[0]->UpdatePhys(),grad0,grad1,grad2);
350  // Take second derivative
351  m_fields[i]->PhysDeriv(0,grad0,grad0);
352  m_fields[i]->PhysDeriv(1,grad1,grad1);
353  m_fields[i]->PhysDeriv(2,grad2,grad2);
354  // Multiply by Intracellular-Conductivity
355  if (m_session->DefinesFunction("IntracellularConductivity") && m_session->DefinesFunction("ExtracellularConductivity"))
356  {
357  Vmath::Smul(nq, m_session->GetParameter("sigmaix"), grad0, 1, grad0, 1);
358  Vmath::Smul(nq, m_session->GetParameter("sigmaiy"), grad1, 1, grad1, 1);
359  Vmath::Smul(nq, m_session->GetParameter("sigmaiz"), grad2, 1, grad2, 1);
360  }
361  // and sum terms
362  Vmath::Vadd(nq, grad0, 1, grad1, 1, grad, 1);
363  Vmath::Vadd(nq, grad2, 1, grad, 1, grad, 1);
364  Vmath::Smul(nq, -1.0, grad, 1, grad, 1);
365  // Now solve Poisson problem for \phi_e
366  m_fields[i]->SetPhys(grad);
367  m_fields[i]->HelmSolve(m_fields[i]->GetPhys(), m_fields[i]->UpdateCoeffs(), NullFlagList, factors, m_vardiffie);
368  m_fields[i]->BwdTrans( m_fields[i]->GetCoeffs(), m_fields[i]->UpdatePhys());
369  m_fields[i]->SetPhysState(true);
370  // Copy the solution vector (required as m_fields must be set).
371  outarray[i] = m_fields[i]->GetPhys();
372  }
373 
374  }
375  }
376  }
377 
378 
379  void Bidomain::DoOdeRhs(
380  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
381  Array<OneD, Array<OneD, NekDouble> >&outarray,
382  const NekDouble time)
383  {
384  int nq = m_fields[0]->GetNpoints();
385  m_cell->TimeIntegrate(inarray, outarray, time);
386  if (m_stimDuration > 0 && time < m_stimDuration)
387  {
388  Array<OneD,NekDouble> x0(nq);
389  Array<OneD,NekDouble> x1(nq);
390  Array<OneD,NekDouble> x2(nq);
391  Array<OneD,NekDouble> result(nq);
392 
393  // get the coordinates
394  m_fields[0]->GetCoords(x0,x1,x2);
395 
396  LibUtilities::EquationSharedPtr ifunc
397  = m_session->GetFunction("Stimulus", "u");
398  ifunc->Evaluate(x0,x1,x2,time, result);
399 
400  Vmath::Vadd(nq, outarray[0], 1, result, 1, outarray[0], 1);
401  }
402  Vmath::Smul(nq, 1.0/m_capMembrane, outarray[0], 1, outarray[0], 1);
403  }
404 
405 
406  void Bidomain::v_SetInitialConditions(NekDouble initialtime,
407  bool dumpInitialConditions,
408  const int domain)
409  {
410  EquationSystem::v_SetInitialConditions(initialtime, dumpInitialConditions, domain);
411  m_cell->Initialise();
412  }
413 
414  /**
415  *
416  */
417  void Bidomain::v_GenerateSummary(SummaryList& s)
418  {
419  UnsteadySystem::v_GenerateSummary(s);
420 
421  /// @TODO Update summary
422  ASSERTL0(false, "Update the generate summary");
423 //
424 // out << "\tChi : " << m_chi << endl;
425 // out << "\tCm : " << m_capMembrane << endl;
426 // if (m_session->DefinesFunction("IntracellularConductivity", "AnisotropicConductivityX") &&
427 // m_session->GetFunctionType("IntracellularConductivity", "AnisotropicConductivityX") == LibUtilities::eFunctionTypeExpression)
428 // {
429 // out << "\tIntra-Diffusivity-x : "
430 // << m_session->GetFunction("IntracellularConductivity", "AnisotropicConductivityX")->GetExpression()
431 // << endl;
432 // }
433 // if (m_session->DefinesFunction("IntracellularConductivity", "AnisotropicConductivityY") &&
434 // m_session->GetFunctionType("IntracellularConductivity", "AnisotropicConductivityY") == LibUtilities::eFunctionTypeExpression)
435 // {
436 // out << "\tIntra-Diffusivity-y : "
437 // << m_session->GetFunction("IntracellularConductivity", "AnisotropicConductivityY")->GetExpression()
438 // << endl;
439 // }
440 // if (m_session->DefinesFunction("IntracellularConductivity", "AnisotropicConductivityZ") &&
441 // m_session->GetFunctionType("IntracellularConductivity", "AnisotropicConductivityZ") == LibUtilities::eFunctionTypeExpression)
442 // {
443 // out << "\tIntra-Diffusivity-z : "
444 // << m_session->GetFunction("IntracellularConductivity", "AnisotropicConductivityZ")->GetExpression()
445 // << endl;
446 // }
447 // if (m_session->DefinesFunction("ExtracellularConductivity", "AnisotropicConductivityX") &&
448 // m_session->GetFunctionType("ExtracellularConductivity", "AnisotropicConductivityX") == LibUtilities::eFunctionTypeExpression)
449 // {
450 // out << "\tExtra-Diffusivity-x : "
451 // << m_session->GetFunction("ExtracellularConductivity", "AnisotropicConductivityX")->GetExpression()
452 // << endl;
453 // }
454 // if (m_session->DefinesFunction("ExtracellularConductivity", "AnisotropicConductivityY") &&
455 // m_session->GetFunctionType("ExtracellularConductivity", "AnisotropicConductivityY") == LibUtilities::eFunctionTypeExpression)
456 // {
457 // out << "\tExtra-Diffusivity-y : "
458 // << m_session->GetFunction("ExtracellularConductivity", "AnisotropicConductivityY")->GetExpression()
459 // << endl;
460 // }
461 // if (m_session->DefinesFunction("ExtracellularConductivity", "AnisotropicConductivityZ") &&
462 // m_session->GetFunctionType("ExtracellularConductivity", "AnisotropicConductivityZ") == LibUtilities::eFunctionTypeExpression)
463 // {
464 // out << "\tExtra-Diffusivity-z : "
465 // << m_session->GetFunction("ExtracellularConductivity", "AnisotropicConductivityZ")->GetExpression()
466 // << endl;
467 // }
468  m_cell->GenerateSummary(s);
469  }
470 
471 }
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2