LifeV
HeartFunctors.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 Heart Functors for the Luo-Rudy Kinetics
30 
31  @date 04−2010
32  @author
33 
34  @contributor Simone Rossi <simone.rossi@epfl.ch>, Ricardo Ruiz-Baier <ricardo.ruiz@epfl.ch>
35  @mantainer Simone Rossi <simone.rossi@epfl.ch>
36  */
37 
38 
39 
40 
41 #include <lifev/heart/solver/HeartFunctors.hpp>
42 
43 using namespace LifeV;
44 
45 // ===================================================
46 // Constructors & Destructor
47 // ===================================================
48 HeartFunctors::HeartFunctors() :
49  M_dataFile ( ),
50  M_stimulusSource ( ),
51  M_stimulusPeriod1 ( ),
52  M_stimulusPeriod2 ( ),
53  M_stimulusPeriod3 ( ),
54  M_stimulusPeriod4 ( ),
55  M_stimulusPeriod5 ( ),
56  M_stimulusPeriod6 ( ),
57  M_stimulusStart1 ( ),
58  M_stimulusStop1 ( ),
59  M_stimulusValue1 ( ),
60  M_stimulusRadius1 ( ),
61  M_stimulusCenter1 (3),
62  M_stimulusStart2 ( ),
63  M_stimulusStop2 ( ),
64  M_stimulusValue2 ( ),
65  M_stimulusRadius2 ( ),
66  M_stimulusCenter2 (3),
67  M_stimulusStart3 ( ),
68  M_stimulusStop3 ( ),
69  M_stimulusValue3 ( ),
70  M_stimulusRadius3 ( ),
71  M_stimulusCenter3 (3),
72  M_stimulusStart4 ( ),
73  M_stimulusStop4 ( ),
74  M_stimulusValue4 ( ),
75  M_stimulusRadius4 ( ),
76  M_stimulusCenter4 (3),
77  M_stimulusStart5 ( ),
78  M_stimulusStop5 ( ),
79  M_stimulusValue5 ( ),
80  M_stimulusRadius5 ( ),
81  M_stimulusCenter5 (3),
82  M_stimulusStart6 ( ),
83  M_stimulusStop6 ( ),
84  M_stimulusValue6 ( ),
85  M_stimulusRadius6 ( ),
86  M_stimulusCenter6 (3),
87  M_sphereX ( ),
88  M_sphereY ( ),
89  M_sphereZ ( ),
90  M_sphereR ( ),
91  M_conductivityReduction (2),
92  M_cylinderX ( ),
93  M_cylinderY ( ),
94  M_cylinderZ ( ),
95  M_cylinderA ( ),
96  M_cylinderB ( ),
97  M_cylinderC ( ),
98  M_cylinderR ( ),
99  M_minimumCylinderX ( ),
100  M_maximumCylinderX ( ),
101  M_minimumBoxX ( ),
102  M_minimumBoxY ( ),
103  M_minimumBoxZ ( ),
104  M_maximumBoxX ( ),
105  M_maximumBoxY ( ),
106  M_maximumBoxZ ( ),
107  //parametre Iapp IappZygote
108  M_timePeriod ( ),
109  M_appliedCurrentRightVentriculeAngle ( ),
110  M_appliedCurrentLeftVentriculeAngle ( ),
111  M_appliedCurrentStimulusTimeRightVentricule ( ),
112  M_appliedCurrentStimulusTimeLeftVentrivcule ( ),
113  M_ventricularFibrillation ( ),
114  M_nb_fibrillationSources ( ),
115  M_fibrillationSources ( ),
116  M_restPotential ( )
117 {
118 }
119 
120 
121 HeartFunctors::HeartFunctors ( GetPot& dataFile ) :
122  M_dataFile (dataFile),
123  M_stimulusSource (dataFile ("electric/physics/stim_source", 1) ),
124  M_stimulusPeriod1 (dataFile ("electric/physics/stim_period_1", 200.) ),
125  M_stimulusPeriod2 (dataFile ("electric/physics/stim_period_2", 200.) ),
126  M_stimulusPeriod3 (dataFile ("electric/physics/stim_period_3", 200.) ),
127  M_stimulusPeriod4 (dataFile ("electric/physics/stim_period_4", 200.) ),
128  M_stimulusPeriod5 (dataFile ("electric/physics/stim_period_5", 200.) ),
129  M_stimulusPeriod6 (dataFile ("electric/physics/stim_period_6", 200.) ),
130  M_stimulusStart1 (dataFile ("electric/physics/stim_start_1", 0.) ),
131  M_stimulusStop1 (dataFile ("electric/physics/stim_stop_1", 0.) ),
132  M_stimulusValue1 (dataFile ("electric/physics/stim_value_1", 0.) ),
133  M_stimulusRadius1 (dataFile ("electric/physics/stim_radius_1", 0.) ),
134  M_stimulusCenter1 (3),
135  M_stimulusStart2 (dataFile ("electric/physics/stim_start_2", 0.) ),
136  M_stimulusStop2 (dataFile ("electric/physics/stim_stop_2", 0.) ),
137  M_stimulusValue2 (dataFile ("electric/physics/stim_value_2", 0.) ),
138  M_stimulusRadius2 (dataFile ("electric/physics/stim_radius_2", 0.) ),
139  M_stimulusCenter2 (3),
140  M_stimulusStart3 (dataFile ("electric/physics/stim_start_3", 0.) ),
141  M_stimulusStop3 (dataFile ("electric/physics/stim_stop_3", 0.) ),
142  M_stimulusValue3 (dataFile ("electric/physics/stim_value_3", 0.) ),
143  M_stimulusRadius3 (dataFile ("electric/physics/stim_radius_3", 0.) ),
144  M_stimulusCenter3 (3),
145  M_stimulusStart4 (dataFile ("electric/physics/stim_start_4", 0.) ),
146  M_stimulusStop4 (dataFile ("electric/physics/stim_stop_4", 0.) ),
147  M_stimulusValue4 (dataFile ("electric/physics/stim_value_4", 0.) ),
148  M_stimulusRadius4 (dataFile ("electric/physics/stim_radius_4", 0.) ),
149  M_stimulusCenter4 (3),
150  M_stimulusStart5 (dataFile ("electric/physics/stim_start_5", 0.) ),
151  M_stimulusStop5 (dataFile ("electric/physics/stim_stop_5", 0.) ),
152  M_stimulusValue5 (dataFile ("electric/physics/stim_value_5", 0.) ),
153  M_stimulusRadius5 (dataFile ("electric/physics/stim_radius_5", 0.) ),
154  M_stimulusCenter5 (3),
155  M_stimulusStart6 (dataFile ("electric/physics/stim_start_6", 0.) ),
156  M_stimulusStop6 (dataFile ("electric/physics/stim_stop_6", 0.) ),
157  M_stimulusValue6 (dataFile ("electric/physics/stim_value_6", 0.) ),
158  M_stimulusRadius6 (dataFile ("electric/physics/stim_radius_6", 0.) ),
159  M_stimulusCenter6 (3),
160  M_sphereX (dataFile ("electric/physics/sphere_center", 0., 0) ),
161  M_sphereY (dataFile ("electric/physics/sphere_center", 0., 1) ),
162  M_sphereZ (dataFile ("electric/physics/sphere_center", 0., 2) ),
163  M_sphereR (dataFile ("electric/physics/sphere_radius", 0.) ),
164  M_conductivityReduction (2),
165  M_cylinderX (dataFile ("electric/physics/x_cylinder", 0.) ),
166  M_cylinderY (dataFile ("electric/physics/y_cylinder", 0.) ),
167  M_cylinderZ (dataFile ("electric/physics/z_cylinder", 0.) ),
168  M_cylinderA (dataFile ("electric/physics/a_cylinder", 0.) ),
169  M_cylinderB (dataFile ("electric/physics/b_cylinder", 0.) ),
170  M_cylinderC (dataFile ("electric/physics/c_cylinder", 0.) ),
171  M_cylinderR (dataFile ("electric/physics/r_cylinder", 0.) ),
172  M_minimumCylinderX (dataFile ("electric/physics/xmin_cylinder", 0.) ),
173  M_maximumCylinderX (dataFile ("electric/physics/xmax_cylinder", 0.) ),
174  M_minimumBoxX (dataFile ("electric/physics/box_vertex_min", 0., 0) ),
175  M_minimumBoxY (dataFile ("electric/physics/box_vertex_min", 0., 1) ),
176  M_minimumBoxZ (dataFile ("electric/physics/box_vertex_min", 0., 2) ),
177  M_maximumBoxX (dataFile ("electric/physics/box_vertex_max", 0., 0) ),
178  M_maximumBoxY (dataFile ("electric/physics/box_vertex_max", 0., 1) ),
179  M_maximumBoxZ (dataFile ("electric/physics/box_vertex_max", 0., 2) ),
180  //parametre Iapp IappZygote
181  M_timePeriod (dataFile ("electric/physics/Time_period", 700.0) ),
182  M_appliedCurrentRightVentriculeAngle (dataFile ("electric/physics/Iapp_RV_angle", 360.) ),
183  M_appliedCurrentLeftVentriculeAngle (dataFile ("electric/physics/Iapp_LV_angle", 360.) ),
184  M_appliedCurrentStimulusTimeRightVentricule (dataFile ("electric/physics/Iapp_stim_time_RV", 6.) ),
185  M_appliedCurrentStimulusTimeLeftVentrivcule (dataFile ("electric/physics/Iapp_stim_time_LV", 10.) ),
186  M_ventricularFibrillation (dataFile ("electric/physics/Ventricular_Fibrillation", 0) ),
187  M_nb_fibrillationSources (dataFile ("electric/physics/nb_fibrillation_sources", 20) ),
188  M_fibrillationSources (dataFile ("electric/physics/fibrillation_sources", 0) ),
189  M_restPotential (dataFile ("electric/physics/u0", 0.) )
190 {
191  M_conductivityReduction (0) = dataFile ("electric/physics/sigma_reduction", 1., 0);
192  M_conductivityReduction (1) = dataFile ("electric/physics/sigma_reduction", 1., 1);
193  M_stimulusCenter1 (0) = dataFile ("electric/physics/stim_center_1", 0., 0);
194  M_stimulusCenter1 (1) = dataFile ("electric/physics/stim_center_1", 0., 1);
195  M_stimulusCenter1 (2) = dataFile ("electric/physics/stim_center_1", 0., 2);
196  M_stimulusCenter2 (0) = dataFile ("electric/physics/stim_center_2", 0., 0);
197  M_stimulusCenter2 (1) = dataFile ("electric/physics/stim_center_2", 0., 1);
198  M_stimulusCenter2 (2) = dataFile ("electric/physics/stim_center_2", 0., 2);
199  M_stimulusCenter3 (0) = dataFile ("electric/physics/stim_center_3", 0., 0);
200  M_stimulusCenter3 (1) = dataFile ("electric/physics/stim_center_3", 0., 1);
201  M_stimulusCenter3 (2) = dataFile ("electric/physics/stim_center_3", 0., 2);
202  M_stimulusCenter4 (0) = dataFile ("electric/physics/stim_center_4", 0., 0);
203  M_stimulusCenter4 (1) = dataFile ("electric/physics/stim_center_4", 0., 1);
204  M_stimulusCenter4 (2) = dataFile ("electric/physics/stim_center_4", 0., 2);
205  M_stimulusCenter5 (0) = dataFile ("electric/physics/stim_center_5", 0., 0);
206  M_stimulusCenter5 (1) = dataFile ("electric/physics/stim_center_5", 0., 1);
207  M_stimulusCenter5 (2) = dataFile ("electric/physics/stim_center_5", 0., 2);
208  M_stimulusCenter6 (0) = dataFile ("electric/physics/stim_center_6", 0., 0);
209  M_stimulusCenter6 (1) = dataFile ("electric/physics/stim_center_6", 0., 1);
210  M_stimulusCenter6 (2) = dataFile ("electric/physics/stim_center_6", 0., 2);
211 }
212 
213 // ===================================================
214 // Set Methods
215 // ===================================================
216 Real
217 HeartFunctors::setAppliedCurrent ( const Real& x, const Real& y, const Real& z, const Real& t ) const
218 {
219  Real appliedCurrent = 0.0;
220  Real pi = std::acos ( -1.0 );
221  Real AngularVelocityLeftVentricule = ( pi / 2 ) / M_appliedCurrentStimulusTimeLeftVentrivcule ; // rd/s
222  Real AngularVelocityRightVentricule = ( pi / 2 ) / M_appliedCurrentStimulusTimeRightVentricule ;
223  Real sumL1, sumL4, aL, bL, cL, e1, e4;
224  //exitation istantane' pi=0;
225  aL = 40;
226  bL = 40;
227  cL = 72;
228  e4 = 0;
229  e1 = 16 * 0.8;
230 
231  sumL1 = x * x / ( ( aL - e1 ) * ( aL - e1 ) ) * 100 +
232  y * y / ( ( bL - e1 ) * ( bL - e1 ) ) * 100 +
233  z * z / ( ( cL - e1 ) * ( cL - e1 ) ) * 100 - 1;
234 
235  sumL4 = x * x / ( ( aL - e4 ) * ( aL - e4 ) ) * 100 +
236  y * y / ( ( bL - e4 ) * ( bL - e4 ) ) * 100 +
237  z * z / ( ( cL - e4 ) * ( cL - e4 ) ) * 100 - 1;
238 
239  Real sumR2, aR, bR, cR, e2;
240  aR = 78;
241  bR = 40;
242  cR = 72 * 0.95;
243  e2 = 10 * 0.8;
244 
245  sumR2 = x * x / ( ( aR - e2 ) * ( aR - e2 ) ) * 100 +
246  y * y / ( ( bR - e2 ) * ( bR - e2 ) ) * 100 +
247  z * z / ( ( cR - e2 ) * ( cR - e2 ) ) * 100 - 1;
248 
249  //============ Coment if BBG
250  if ( fmod (t, M_timePeriod) >= 0 &&
251  fmod (t, M_timePeriod) <= M_appliedCurrentStimulusTimeLeftVentrivcule )
252  {
253  if ( sumL1 < 0.15 &&
254  ( std::atan ( ( x + 1.5e+00) / ( 2.7586e+00 - z) ) < fmod (t, M_timePeriod) * AngularVelocityLeftVentricule ) ) // -0.05<sumL1< 0.15
255  {
256  // if ( ( atan(( x + 1.5e+00) /( 2.7586e+00 - z) ) < pi/4 ) ) // BBG Partiel angle pi/4
257  if ( ( std::atan ( ( x + 1.5e+00) / ( 2.7586e+00 - z) ) < M_appliedCurrentLeftVentriculeAngle * pi / 180 ) ) //
258  {
259  appliedCurrent = ( -5 / 2 * sumL1 + 0.375) * 8 ;
260  }
261  }
262  }
263 
264  //=====================septum droit
265  if ( fmod (t, M_timePeriod) >= 0 &&
266  fmod (t, M_timePeriod) <= M_appliedCurrentStimulusTimeLeftVentrivcule &&
267  sumL4 >= -0.2 && sumL4 <= 0.0 && x < -2.4 )
268  {
269  appliedCurrent = (5 / 2 * sumL4 + 0.5) * 8 ;
270  }
271 
272  //============ Coment if BBD
273  if ( fmod (t, M_timePeriod) >= 0 &&
274  fmod (t, M_timePeriod) <= M_appliedCurrentStimulusTimeLeftVentrivcule )
275  {
276  if (x < -2.586e+00)
277  if (sumL4 > -0.2 )
278  {
279  if ( sumR2 >= 0.0 && sumR2 < 0.1 &&
280  std::atan ( (- 3.7e+00 - x ) / ( 2.9586e+00 - z) ) < fmod (t, M_timePeriod) * AngularVelocityRightVentricule )
281  // if ( ( atan((- 3.7e+00 - x) /( 2.9586e+00 - z) ) < pi/6 )) //BBD Partiel angle pi/6
282  if ( std::atan ( (- 3.7e+00 - x) / ( 2.9586e+00 - z) ) < M_appliedCurrentRightVentriculeAngle * pi / 180 ) //
283  {
284  appliedCurrent = 0.5 * 8 ; //( -5 * sumR2 + 0.5) * 8 ;//0.3*4.;
285  }
286  }
287  }
288 
289  return appliedCurrent;
290 
291 }
292 
293 
294 Real
295 HeartFunctors::setAppliedCurrentZygote (const double& t,
296  const double& x,
297  const double& y,
298  const double& z,
299  const ID& /*i*/,
300  const markerID_Type& ref)
301 {
302  // double pi = acos(-1.0);
303  Real appliedCurrent = 0.0;
304  Real x0 = 3.316424;
305  Real y0 = 31.496351;
306  Real z0 = 5.799850; //APEX, node number: 80185 : 3.316424 31.496351 5.799850
307 
308  if ( fmod (t, M_timePeriod) <= 25 )
309  {
310  if ( ref == 2 ) // ventricule gauche if (atan( d1/d2) < pi* t/20 )
311  {
312  if ( ( ( x - x0 ) * ( x - x0 ) +
313  ( y - y0 ) * ( y - y0 ) +
314  ( z - z0 ) * ( z - z0 ) >= fmod (t, M_timePeriod) * fmod (t, M_timePeriod) - 100 )
315  &&
316  ( ( x - x0 ) * ( x - x0 ) +
317  ( y - y0 ) * ( y - y0 ) +
318  ( z - z0 ) * ( z - z0 ) <= 3 * fmod (t, M_timePeriod) * fmod (t, M_timePeriod) ) )
319  {
320  appliedCurrent = 2;
321  }
322  }
323  if ( ( ref == 1 ) || ( ref == 20 ) ) // if (ref==20) //BBD
324  {
325  //venticule droit if (atan( d1/d2) < pi* t/20 )
326  if ( ( ( x - x0 ) * ( x - x0 ) +
327  ( y - y0 ) * ( y - y0 ) +
328  ( z - z0 ) * ( z - z0 ) >= fmod (t, M_timePeriod) * fmod (t, M_timePeriod) - 100 )
329  &&
330  ( ( x - x0 ) * ( x - x0 ) +
331  ( y - y0 ) * ( y - y0 ) +
332  ( z - z0 ) * ( z - z0 ) <= 3 * fmod (t, M_timePeriod) * fmod (t, M_timePeriod) ) )
333  {
334  appliedCurrent = 2;
335  }
336  }
337  }
338 
339  return appliedCurrent;
340 }
341 
342 
343 Real
344 HeartFunctors::setStimulus ( const Real& t, const Real& x, const Real& y, const Real& z, const ID& id) const
345 {
346  Real returnValue1;
347  Real returnValue2;
348  Real returnValue3;
349  Real returnValue4;
350  Real returnValue5;
351  Real returnValue6;
352  Real timeReset1 (t - static_cast<int> (t / M_stimulusPeriod1) * M_stimulusPeriod1);
353  Real timeReset2 (t - static_cast<int> (t / M_stimulusPeriod2) * M_stimulusPeriod2);
354  Real timeReset3 (t - static_cast<int> (t / M_stimulusPeriod3) * M_stimulusPeriod3);
355  Real timeReset4 (t - static_cast<int> (t / M_stimulusPeriod4) * M_stimulusPeriod4);
356  Real timeReset5 (t - static_cast<int> (t / M_stimulusPeriod5) * M_stimulusPeriod5);
357  Real timeReset6 (t - static_cast<int> (t / M_stimulusPeriod6) * M_stimulusPeriod6);
358 
359  if ( (timeReset1 >= M_stimulusStart1 && timeReset1 <= M_stimulusStop1) &&
360  ( ( ( x - M_stimulusCenter1 (0) ) * ( x - M_stimulusCenter1 (0) ) +
361  ( y - M_stimulusCenter1 (1) ) * ( y - M_stimulusCenter1 (1) ) +
362  ( z - M_stimulusCenter1 (2) ) * ( z - M_stimulusCenter1 (2) ) ) <= ( M_stimulusRadius1 * M_stimulusRadius1 ) ) )
363  {
364  returnValue1 = M_stimulusValue1;
365  }
366  else
367  {
368  returnValue1 = 0.;
369  }
370 
371  if ( (timeReset2 >= M_stimulusStart2 && timeReset2 <= M_stimulusStop2) &&
372  ( ( ( x - M_stimulusCenter2 (0) ) * ( x - M_stimulusCenter2 (0) ) +
373  ( y - M_stimulusCenter2 (1) ) * ( y - M_stimulusCenter2 (1) ) +
374  ( z - M_stimulusCenter2 (2) ) * ( z - M_stimulusCenter2 (2) ) ) <= ( M_stimulusRadius2 * M_stimulusRadius2 ) ) )
375  {
376  returnValue2 = M_stimulusValue2;
377  }
378  else
379  {
380  returnValue2 = 0.;
381  }
382 
383  if ( (timeReset3 >= M_stimulusStart3 && timeReset3 <= M_stimulusStop3) &&
384  ( ( ( x - M_stimulusCenter3 (0) ) * ( x - M_stimulusCenter3 (0) ) +
385  ( y - M_stimulusCenter3 (1) ) * ( y - M_stimulusCenter3 (1) ) +
386  ( z - M_stimulusCenter3 (2) ) * ( z - M_stimulusCenter3 (2) ) ) <= ( M_stimulusRadius3 * M_stimulusRadius3 ) ) )
387  {
388  returnValue3 = M_stimulusValue3;
389  }
390  else
391  {
392  returnValue3 = 0.;
393  }
394 
395  if ( (timeReset4 >= M_stimulusStart4 && timeReset4 <= M_stimulusStop4) &&
396  ( ( ( x - M_stimulusCenter4 (0) ) * ( x - M_stimulusCenter4 (0) ) +
397  ( y - M_stimulusCenter4 (1) ) * ( y - M_stimulusCenter4 (1) ) +
398  ( z - M_stimulusCenter4 (2) ) * ( z - M_stimulusCenter4 (2) ) ) <= ( M_stimulusRadius4 * M_stimulusRadius4 ) ) )
399  {
400  returnValue4 = M_stimulusValue4;
401  }
402  else
403  {
404  returnValue4 = 0.;
405  }
406 
407  if ( (timeReset5 >= M_stimulusStart5 && timeReset5 <= M_stimulusStop5) &&
408  ( ( ( x - M_stimulusCenter5 (0) ) * ( x - M_stimulusCenter5 (0) ) +
409  ( y - M_stimulusCenter5 (1) ) * ( y - M_stimulusCenter5 (1) ) +
410  ( z - M_stimulusCenter5 (2) ) * ( z - M_stimulusCenter5 (2) ) ) <= ( M_stimulusRadius5 * M_stimulusRadius5 ) ) )
411  {
412  returnValue5 = M_stimulusValue5;
413  }
414  else
415  {
416  returnValue5 = 0.;
417  }
418 
419  if ( (timeReset6 >= M_stimulusStart6 && timeReset6 <= M_stimulusStop6) &&
420  ( ( ( x - M_stimulusCenter6 (0) ) * ( x - M_stimulusCenter6 (0) ) +
421  ( y - M_stimulusCenter6 (1) ) * ( y - M_stimulusCenter6 (1) ) +
422  ( z - M_stimulusCenter6 (2) ) * ( z - M_stimulusCenter6 (2) ) ) <= ( M_stimulusRadius6 * M_stimulusRadius6 ) ) )
423  {
424  returnValue6 = M_stimulusValue6;
425  }
426  else
427  {
428  returnValue6 = 0.;
429  }
430 
431  Real returnValue = returnValue1 + returnValue2 + returnValue3 + returnValue4 + returnValue5 + returnValue6;
432 
433  if (id == 0)
434  {
435  return returnValue;
436  }
437  else
438  {
439  return returnValue;
440  }
441 
442 }
443 
444 
445 
446 
447 
448 
449 
450 
451 Real
452 HeartFunctors::setReducedConductivitySphere ( const Real& x,
453  const Real& y,
454  const Real& z,
455  const Real& /*t*/,
456  const ID& id,
457  const Real& sigma) const
458 {
459  if ( ( ( x - M_sphereX ) * ( x - M_sphereX ) +
460  ( y - M_sphereY ) * ( y - M_sphereY ) +
461  ( z - M_sphereZ ) * ( z - M_sphereZ ) ) < M_sphereR * M_sphereR )
462  {
463  return sigma * M_conductivityReduction (id);
464  }
465  else
466  {
467  return sigma;
468  }
469 }
470 
471 
472 
473 Real
474 HeartFunctors::setReducedConductivityCylinder ( const Real& x,
475  const Real& y,
476  const Real& z,
477  const Real& /*t*/,
478  const ID& id,
479  const Real& sigma ) const
480 {
481  Real distance2, distanceX, distanceY, distanceZ;
482  distanceX = ( ( M_cylinderB * M_cylinderB + M_cylinderC * M_cylinderC ) * ( M_cylinderX - x ) -
483  ( M_cylinderA * M_cylinderC ) * ( M_cylinderZ - z ) -
484  ( M_cylinderA * M_cylinderB ) * ( M_cylinderY - y ) ) /
485  ( M_cylinderA * M_cylinderA + M_cylinderB * M_cylinderB + M_cylinderC * M_cylinderC );
486 
487  distanceY = ( ( M_cylinderA * M_cylinderA + M_cylinderC * M_cylinderC ) * ( M_cylinderY - y ) -
488  ( M_cylinderB * M_cylinderC ) * ( M_cylinderZ - z ) -
489  ( M_cylinderA * M_cylinderB ) * ( M_cylinderX - x ) ) /
490  ( M_cylinderA * M_cylinderA + M_cylinderB * M_cylinderB + M_cylinderC * M_cylinderC );
491 
492  distanceZ = ( ( M_cylinderA * M_cylinderA + M_cylinderB * M_cylinderB ) * ( M_cylinderZ - z ) -
493  ( M_cylinderA * M_cylinderC ) * ( M_cylinderX - x ) -
494  ( M_cylinderC * M_cylinderB ) * ( M_cylinderY - y ) ) /
495  ( M_cylinderA * M_cylinderA + M_cylinderB * M_cylinderB + M_cylinderC * M_cylinderC );
496 
497  distance2 = std::pow (distanceX, 2) + std::pow (distanceY, 2) + std::pow (distanceZ, 2);
498 
499  if ( ( distance2 < M_cylinderR * M_cylinderR ) && ( x < M_maximumCylinderX ) && ( x > M_minimumCylinderX ) )
500  {
501  return sigma * M_conductivityReduction (id);
502  }
503  else
504  {
505  return sigma;
506  }
507 }
508 
509 
510 
511 
512 Real
513 HeartFunctors::setReducedConductivityBox ( const Real& x,
514  const Real& y,
515  const Real& z,
516  const Real& /*t*/,
517  const ID& id,
518  const Real& sigma ) const
519 {
520  if ( ( x > M_minimumBoxX ) && ( x < M_maximumBoxX ) &&
521  ( y > M_minimumBoxY ) && ( y < M_maximumBoxY ) &&
522  ( z > M_minimumBoxZ ) && ( z < M_maximumBoxZ ) )
523  {
524  return sigma * M_conductivityReduction (id);
525  }
526  else
527  {
528  return sigma;
529  }
530 }
531 
532 Real
533 HeartFunctors::setInitialScalar ( const Real& /*t*/ ,
534  const Real& /*x*/,
535  const Real& /*y*/,
536  const Real& /*z*/,
537  const ID& /*id*/ )
538 {
539  return M_restPotential;
540 }
541 
542 
543 
544 Real
545 HeartFunctors::setZeroScalar ( const Real& /*t*/,
546  const Real& /*x*/,
547  const Real& /*y*/,
548  const Real& /*z*/,
549  const ID& /*id*/ )
550 {
551  return 0.;
552 }
553 // ===================================================
554 // Get Methods
555 // ===================================================
556 
557 
558 
559 HeartFunctors::region1_Type
560 HeartFunctors::appliedCurrent()
561 {
562  region1_Type f;
563  f = std::bind (&HeartFunctors::setAppliedCurrent, this, 1, 2, 3, 4 );
564  return f;
565 }
566 
567 HeartFunctors::region1_Type
568 HeartFunctors::stimulus()
569 {
570  region1_Type f;
571  f = std::bind (&HeartFunctors::setStimulus, this, 1, 2, 3, 4, 5);
572  return f;
573 }
574 
575 const HeartFunctors::region_Type
576 HeartFunctors::reducedConductivitySphere()
577 {
578  region_Type f;
579  f = std::bind (&HeartFunctors::setReducedConductivitySphere, this, 1, 2, 3, 4, 5, 6);
580  return f;
581 }
582 
583 const HeartFunctors::region_Type
584 HeartFunctors::reducedConductivityCylinder()
585 {
586  region_Type f;
587  f = std::bind (&HeartFunctors::setReducedConductivityCylinder, this, 1, 2, 3, 4, 5, 6);
588  return f;
589 }
590 
591 const HeartFunctors::region_Type
592 HeartFunctors::reducedConductivityBox()
593 {
594  region_Type f;
595  f = std::bind (&HeartFunctors::setReducedConductivityBox, this, 1, 2, 3, 4, 5, 6);
596  return f;
597 }
598 
599 const HeartFunctors::region1_Type
600 HeartFunctors::initialScalar()
601 {
602  region1_Type f;
603  f = std::bind (&HeartFunctors::setInitialScalar, this, 1, 2, 3, 4, 5);
604  return f;
605 }
606 
607 const HeartFunctors::region1_Type
608 HeartFunctors::zeroScalar()
609 {
610  region1_Type f;
611  f = std::bind (&HeartFunctors::setZeroScalar, this, 1, 2, 3, 4, 5);
612  return f;
613 }
GetPot::operator()
double operator()(const char *VarName, const double &Default, unsigned Idx) const
Definition: GetPot.hpp:2154
LifeV::HeartFunctors::setStimulus
Real setStimulus(const Real &t, const Real &x, const Real &y, const Real &z, const ID &id) const
Definition: HeartFunctors.cpp:344
GetPot::operator()
double operator()(const char *VarName, const double &Default) const
Definition: GetPot.hpp:2033
LifeV::HeartFunctors::setAppliedCurrentZygote
Real setAppliedCurrentZygote(const double &t, const double &x, const double &y, const double &z, const ID &, const markerID_Type &ref)
Definition: HeartFunctors.cpp:295
LifeV::HeartFunctors::setInitialScalar
Real setInitialScalar(const Real &t, const Real &x, const Real &y, const Real &z, const ID &id)
Definition: HeartFunctors.cpp:533
LifeV::HeartFunctors::region_Type
std::function< Real(Real const &x, Real const &y, Real const &z, Real const &, ID const &id, Real const &) > region_Type
Definition: HeartFunctors.hpp:64
LifeV::markerID_Type
ID markerID_Type
markerID_Type is the type used to store the geometric entity marker IDs
Definition: Marker.hpp:81
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::HeartFunctors::reducedConductivitySphere
const region_Type reducedConductivitySphere()
Definition: HeartFunctors.cpp:576
LifeV::HeartFunctors::appliedCurrent
region1_Type appliedCurrent()
Definition: HeartFunctors.cpp:560
LifeV::HeartFunctors::setReducedConductivityCylinder
Real setReducedConductivityCylinder(const Real &x, const Real &y, const Real &z, const Real &, const ID &id, const Real &sigma) const
Reduces the conductivity in a cylinder.
Definition: HeartFunctors.cpp:474
LifeV::ID
uint32_type ID
IDs.
Definition: LifeV.hpp:194
LifeV::HeartFunctors::region1_Type
std::function< Real(Real const &x, Real const &y, Real const &z, Real const &, ID const &id) > region1_Type
Definition: HeartFunctors.hpp:65
LifeV::HeartFunctors::M_appliedCurrentRightVentriculeAngle
Real M_appliedCurrentRightVentriculeAngle
Definition: HeartFunctors.hpp:218
LifeV::HeartFunctors::initialScalar
const region1_Type initialScalar()
Definition: HeartFunctors.cpp:600
LifeV::HeartFunctors::setAppliedCurrent
Real setAppliedCurrent(const Real &x, const Real &y, const Real &z, const Real &t) const
current volume source
Definition: HeartFunctors.cpp:217
LifeV::HeartFunctors::HeartFunctors
HeartFunctors(GetPot &dataFile)
Definition: HeartFunctors.cpp:121
LifeV::HeartFunctors::setReducedConductivitySphere
Real setReducedConductivitySphere(const Real &x, const Real &y, const Real &z, const Real &, const ID &id, const Real &sigma) const
Reduces the conductivity in a sphere.
Definition: HeartFunctors.cpp:452
LifeV::HeartFunctors::setZeroScalar
Real setZeroScalar(const Real &t, const Real &x, const Real &y, const Real &z, const ID &id)
Definition: HeartFunctors.cpp:545
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
GetPot::operator()
int operator()(const char *VarName, int Default) const
Definition: GetPot.hpp:2021
LifeV::HeartFunctors::M_appliedCurrentLeftVentriculeAngle
Real M_appliedCurrentLeftVentriculeAngle
Definition: HeartFunctors.hpp:219
LifeV::HeartFunctors::M_appliedCurrentStimulusTimeRightVentricule
Real M_appliedCurrentStimulusTimeRightVentricule
Definition: HeartFunctors.hpp:220
LifeV::HeartFunctors::zeroScalar
const region1_Type zeroScalar()
Definition: HeartFunctors.cpp:608
LifeV::HeartFunctors::M_appliedCurrentStimulusTimeLeftVentrivcule
Real M_appliedCurrentStimulusTimeLeftVentrivcule
Definition: HeartFunctors.hpp:221
LifeV::HeartFunctors::reducedConductivityBox
const region_Type reducedConductivityBox()
Definition: HeartFunctors.cpp:592
LifeV::HeartFunctors::stimulus
region1_Type stimulus()
Definition: HeartFunctors.cpp:568
LifeV::HeartFunctors::reducedConductivityCylinder
const region_Type reducedConductivityCylinder()
Definition: HeartFunctors.cpp:584
LifeV::HeartFunctors::setReducedConductivityBox
Real setReducedConductivityBox(const Real &x, const Real &y, const Real &z, const Real &, const ID &id, const Real &sigma) const
Definition: HeartFunctors.cpp:513