Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AlievPanfilov.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File AlievPanfilov.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: Aliev-Panfilov phenomological cell model.
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
36 #include <iostream>
37 #include <string>
38 
39 #include <LibUtilities/BasicUtils/Vmath.hpp>
40 #include <CardiacEPSolver/CellModels/AlievPanfilov.h>
41 
42 namespace Nektar
43 {
44  /**
45  * Registers the class with the Factory.
46  */
47  std::string CellModelAlievPanfilov::className
48  = GetCellModelFactory().RegisterCreatorFunction(
49  "AlievPanfilov",
50  CellModelAlievPanfilov::create,
51  "Phenomological model of canine cardiac electrophysiology.");
52 
53  CellModelAlievPanfilov::CellModelAlievPanfilov(
54  const LibUtilities::SessionReaderSharedPtr& pSession,
55  const MultiRegions::ExpListSharedPtr& pField)
56  : CellModel(pSession, pField)
57  {
58  pSession->LoadParameter("k", m_k, 0.0);
59  pSession->LoadParameter("a", m_a, 0.0);
60  pSession->LoadParameter("mu1", m_mu1, 0.0);
61  pSession->LoadParameter("mu2", m_mu2, 0.0);
62  pSession->LoadParameter("eps", m_eps, 0.0);
63 
64  m_uu = Array<OneD, NekDouble>(m_nq, 0.0);
65  m_uuu = Array<OneD, NekDouble>(m_nq, 0.0);
66  m_tmp1 = Array<OneD, NekDouble>(m_nq, 0.0);
67  m_tmp2 = Array<OneD, NekDouble>(m_nq, 0.0);
68 
69  m_nvar = 2;
70  m_concentrations.push_back(1);
71  }
72 
73 
74  void CellModelAlievPanfilov::v_Update(
75  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
76  Array<OneD, Array<OneD, NekDouble> >&outarray,
77  const NekDouble time)
78  {
79  // inarray[0] holds initial physical u values throughout
80  // inarray[1] holds initial physical v values throughout
81 
82  // compute u^2: m_u = u*u
83  Vmath::Vmul(m_nq, &inarray[0][0], 1, &inarray[0][0], 1, &m_uu[0], 1);
84 
85  // compute u^3: m_u = u*u*u
86  Vmath::Vmul(m_nq, &inarray[0][0], 1, &m_uu[0], 1, &m_uuu[0], 1);
87 
88  // --------------------------------------
89  // Compute reaction term f(u,v)
90  // --------------------------------------
91 // if (m_spatialParameters->Exists("a"))
92 // {
93 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
94 // &inarray[0][0], 1, &m_tmp1[0], 1);
95 //
96 // Vmath::Vvtvm(m_nq, &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
97 // &m_uu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
98 //
99 // Vmath::Svtvm(m_nq, -1.0, &m_uu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
100 // }
101 // else
102 // {
103  // Ru = au
104  Vmath::Smul(m_nq, m_a, &inarray[0][0], 1, &m_tmp1[0], 1);
105  // Ru = (-1-a)u*u + au
106  Vmath::Svtvp(m_nq, (-1.0-m_a), &m_uu[0], 1, &m_tmp1[0], 1,
107  &m_tmp1[0], 1);
108 // }
109  // Ru = u*u*u - (1+a)u*u + au
110  Vmath::Vadd(m_nq, &m_uuu[0], 1, &m_tmp1[0], 1, &m_tmp1[0], 1);
111  // Ru = k(u*u*u - (1+a)u*u + au)
112 // if (m_spatialParameters->Exists("k"))
113 // {
114 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("k")->GetPhys()[0], 1,
115 // &m_tmp1[0], 1, &m_tmp1[0], 1);
116 // }
117 // else
118 // {
119  Vmath::Smul(m_nq, m_k, &m_tmp1[0], 1, &m_tmp1[0], 1);
120 // }
121 
122  // Ru = k(u*u*u - (1+a)u*u + au) + I_stim
123  Vmath::Vadd(m_nq, &outarray[0][0], 1, &m_tmp1[0], 1, &outarray[0][0], 1);
124 
125  // Ru = k(u*u*u - (1+a)u*u + au) + uv + I_stim
126  Vmath::Vvtvp(m_nq, &inarray[0][0], 1, &inarray[1][0], 1, &m_tmp1[0], 1,
127  &outarray[0][0], 1);
128  // Ru = -k(u*u*u - (1+a)u*u + au) - uv - I_stim
129  Vmath::Neg(m_nq, &outarray[0][0], 1);
130 
131 
132  // --------------------------------------
133  // Compute reaction term g(u,v)
134  // --------------------------------------
135  // tmp2 = mu2 + u
136  Vmath::Sadd(m_nq, m_mu2, &inarray[0][0], 1, &m_tmp2[0], 1);
137 
138  // tmp2 = v/(mu2 + u)
139  Vmath::Vdiv(m_nq, &inarray[1][0], 1, &m_tmp2[0], 1, &m_tmp2[0], 1);
140 
141  // tmp2 = mu1*v/(mu2 + u)
142  Vmath::Smul(m_nq, m_mu1, &m_tmp2[0], 1, &m_tmp2[0], 1);
143 
144  // tmp1 = Eps + mu1*v/(mu2+u)
145  Vmath::Sadd(m_nq, m_eps, &m_tmp2[0], 1, &m_tmp2[0], 1);
146 
147  // tmp1 = (-a-1) + u
148 // if (m_spatialParameters->Exists("a"))
149 // {
150 // Vmath::Vsub(m_nq, &inarray[0][0], 1,
151 // &m_spatialParameters->GetData("a")->GetPhys()[0], 1,
152 // &m_tmp1[0], 1);
153 //
154 // Vmath::Sadd(m_nq, -1.0, &inarray[0][0], 1, &m_tmp1[0], 1);
155 // }
156 // else
157 // {
158  Vmath::Sadd(m_nq, (-m_a-1), &inarray[0][0], 1, &m_tmp1[0], 1);
159 // }
160 
161  // tmp1 = k(u-a-1)
162 // if (m_spatialParameters->Exists("k"))
163 // {
164 // Vmath::Vmul(m_nq, &m_spatialParameters->GetData("k")->GetPhys()[0], 1,
165 // &m_tmp1[0], 1, &m_tmp1[0], 1);
166 // }
167 // else
168 // {
169  Vmath::Smul(m_nq, m_k, &m_tmp1[0], 1, &m_tmp1[0], 1);
170 // }
171 
172  // tmp1 = ku(u-a-1) + v
173  Vmath::Vvtvp(m_nq, &inarray[0][0], 1, &m_tmp1[0], 1, &inarray[1][0], 1,
174  &m_tmp1[0], 1);
175 
176  // tmp1 = -ku(u-a-1)-v
177  Vmath::Neg(m_nq, &m_tmp1[0], 1);
178 
179  // outarray = [Eps + mu1*v/(mu2+u)] * [-ku(u-a-1)-v]
180  Vmath::Vmul(m_nq, &m_tmp1[0], 1, &m_tmp2[0], 1, &outarray[1][0], 1);
181  }
182 
183  /**
184  *
185  */
186  void CellModelAlievPanfilov::v_GenerateSummary(SummaryList& s)
187  {
188  SolverUtils::AddSummaryItem(s, "Cell model","Aliev-Panfilov");
189  SolverUtils::AddSummaryItem(s, "k", m_k);
190  SolverUtils::AddSummaryItem(s, "a", m_a);
191  SolverUtils::AddSummaryItem(s, "eps", m_eps);
192  SolverUtils::AddSummaryItem(s, "mu1", m_mu1);
193  SolverUtils::AddSummaryItem(s, "mu2", m_mu2);
194  }
195 
196 
197  /**
198  *
199  */
200  void CellModelAlievPanfilov::v_SetInitialConditions()
201  {
202  Vmath::Fill(m_nq, 0.0, m_cellSol[0], 1);
203  Vmath::Fill(m_nq, 0.0, m_cellSol[1], 1);
204  }
205 }
Generated on Sun Mar 15 2015 00:11:29 for Nektar++ by   doxygen 1.8.1.2