EMSolver.hpp 21 KB
Newer Older
Simone Rossi's avatar
Simone Rossi committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
//@HEADER
/*
 *******************************************************************************

 Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
 Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University

 This file is part of LifeV.

 LifeV is free software; you can redistribute it and/or modify
 it under the terms of the GNU Lesser General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 LifeV is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public License
 along with LifeV.  If not, see <http://www.gnu.org/licenses/>.

 *******************************************************************************
 */
//@HEADER
/*!
 @file
 @brief Class for solving the Monodomain equations in electrophysiology.

 @date 02-2013
 @author Simone Rossi <simone.rossi@epfl.ch>

 @last update 02-2013

 This class provides interfaces to solve the monodomain equation
 ( reaction diffusion equation ) using the ETA framework.
 The solution can be performed using three different methods:
 -operator splitting method (at this point available only with forward Euler
 for the reaction step and backward Euler for the diffusion step. );
 -Ionic Currents Interpolation (at this point only forward Euler);
 -State Variable interpolation (at this point only forward Euler).
 */

#ifndef _EMSOLVER_H_
#define _EMSOLVER_H_

Simone Rossi's avatar
Simone Rossi committed
47
#include <lifev/core/mesh/MeshLoadingUtility.hpp>
Simone Rossi's avatar
Simone Rossi committed
48

Simone Rossi's avatar
Simone Rossi committed
49
#include <lifev/em/solver/electrophysiology/EMMonodomainSolver.hpp>
Simone Rossi's avatar
Simone Rossi committed
50
51
52
53
54
55
56
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>
#include <lifev/em/solver/mechanics/EMStructuralConstitutiveLaw.hpp>
#include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>


#include <lifev/em/solver/activation/activeStressModels/ActiveStressActivation.hpp>

Simone Rossi's avatar
Simone Rossi committed
57
58
59
60
61
62
63
64
65
66
67
//
//#include <lifev/em/solver/EMStructuralOperator.hpp>
//#include <lifev/em/solver/EMGeneralizedActiveHolzapfelOgdenMaterial.hpp>
//#include <lifev/em/solver/EMActiveStrainSolver.hpp>
//#include <lifev/core/interpolation/RBFlocallyRescaledVectorial.hpp>
//#include <lifev/core/interpolation/RBFlocallyRescaledScalar.hpp>
//#include <lifev/core/interpolation/RBFrescaledVectorial.hpp>
//#include <lifev/core/interpolation/RBFrescaledScalar.hpp>
////#include <lifev/core/interpolation/RBFscalar.hpp>
//#include <lifev/core/interpolation/RBFvectorial.hpp>
//
Simone Rossi's avatar
Simone Rossi committed
68
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
Simone Rossi's avatar
Simone Rossi committed
69
70
71
//
//
//#include <lifev/em/solver/EMEvaluate.hpp>
Simone Rossi's avatar
Simone Rossi committed
72
73
74
75
76
77

namespace LifeV
{

//! EMSolver - Class featuring the solution of the electromechanical problem with monodomain equation

Simone Rossi's avatar
Simone Rossi committed
78
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
Simone Rossi's avatar
Simone Rossi committed
79
80
class EMSolver
{
Simone Rossi's avatar
Simone Rossi committed
81
public:
Simone Deparis's avatar
Simone Deparis committed
82
    typedef Mesh                                              mesh_Type;
Simone Rossi's avatar
Simone Rossi committed
83

Simone Deparis's avatar
Simone Deparis committed
84
    typedef boost::shared_ptr<mesh_Type>                      meshPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
85

Simone Deparis's avatar
Simone Deparis committed
86
    typedef Epetra_Comm                                       comm_Type;
Simone Rossi's avatar
Simone Rossi committed
87

Simone Deparis's avatar
Simone Deparis committed
88
    typedef boost::shared_ptr<Epetra_Comm>                    commPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
89

Simone Deparis's avatar
Simone Deparis committed
90
    typedef VectorEpetra                                      vector_Type;
Simone Rossi's avatar
Simone Rossi committed
91

Simone Deparis's avatar
Simone Deparis committed
92
    typedef StructuralConstitutiveLawData                     structureData_Type;
Simone Rossi's avatar
Simone Rossi committed
93

Simone Deparis's avatar
Simone Deparis committed
94
    typedef boost::shared_ptr<structureData_Type>             structureDataPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
95

Simone Deparis's avatar
Simone Deparis committed
96
    typedef EMStructuralOperator< mesh_Type >                 structuralOperator_Type;
Simone Rossi's avatar
Simone Rossi committed
97

Simone Deparis's avatar
Simone Deparis committed
98
    typedef boost::shared_ptr< structuralOperator_Type >      structuralOperatorPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
99

Simone Deparis's avatar
Simone Deparis committed
100
    typedef BCHandler                                          bc_Type;
Simone Rossi's avatar
Simone Rossi committed
101

Simone Deparis's avatar
Simone Deparis committed
102
    typedef boost::shared_ptr< bc_Type >                       bcPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
103

Simone Deparis's avatar
Simone Deparis committed
104
    typedef StructuralOperator< mesh_Type >                   physicalSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
105

Simone Deparis's avatar
Simone Deparis committed
106
    typedef BCInterface3D< bc_Type, physicalSolver_Type >  bcInterface_Type;
Simone Rossi's avatar
Simone Rossi committed
107

Simone Deparis's avatar
Simone Deparis committed
108
    typedef boost::shared_ptr< bcInterface_Type >              bcInterfacePtr_Type;
Simone Rossi's avatar
Simone Rossi committed
109

Simone Deparis's avatar
Simone Deparis committed
110
    typedef boost::shared_ptr<ActivationModel>                 activationModelPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
111

Simone Deparis's avatar
Simone Deparis committed
112
    typedef ElectroSolver                                      electroSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
113

Simone Deparis's avatar
Simone Deparis committed
114
    typedef boost::shared_ptr<electroSolver_Type>              electroSolverPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
115

Simone Deparis's avatar
Simone Deparis committed
116
    typedef ElectroIonicModel                                  ionicModel_Type;
Simone Rossi's avatar
Simone Rossi committed
117

Simone Deparis's avatar
Simone Deparis committed
118
    typedef boost::shared_ptr<ionicModel_Type>                 ionicModelPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
119

Simone Deparis's avatar
Simone Deparis committed
120
    typedef ExporterHDF5< Mesh >                               exporter_Type;
Simone Rossi's avatar
Simone Rossi committed
121

Simone Deparis's avatar
Simone Deparis committed
122
    typedef boost::shared_ptr<ExporterHDF5< Mesh > >           exporterPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
123
124
125
126
127
128
129
130
131
132

    typedef FESpace< RegionMesh<LinearTetra>, MapEpetra >      solidFESpace_Type;

    typedef boost::shared_ptr<solidFESpace_Type>                solidFESpacePtr_Type;

    typedef ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > solidETFESpace_Type;

    typedef boost::shared_ptr<solidETFESpace_Type>              solidETFESpacePtr_Type;

    typedef boost::function < Real (const Real& t,
Simone Deparis's avatar
Simone Deparis committed
133
134
                                    const Real&    x,
                                    const Real&    y,
Simone Rossi's avatar
Simone Rossi committed
135
136
137
138
139
140
                                    const Real& z,
                                    const ID&   /*i*/ ) >       function_Type;



    EMSolver();
Simone Rossi's avatar
Simone Rossi committed
141

Simone Deparis's avatar
Simone Deparis committed
142
    EMSolver (const EMSolver& solver);
Simone Rossi's avatar
Simone Rossi committed
143

Simone Rossi's avatar
Simone Rossi committed
144

Simone Rossi's avatar
Simone Rossi committed
145

Simone Deparis's avatar
Simone Deparis committed
146
147
    inline void loadMesh (std::string meshName, std::string meshPath)
    {
Simone Rossi's avatar
Simone Rossi committed
148
        std::cout << "EMS - Loading mesh\n";
Simone Deparis's avatar
Simone Deparis committed
149
150
        MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
    }
Simone Rossi's avatar
Simone Rossi committed
151

Simone Deparis's avatar
Simone Deparis committed
152
153
154
155
    inline void setupElectroExporter ( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
    {
        M_electroSolverPtr -> setupExporter (*M_electroExporterPtr, outputFileName, problemFolder);
    }
Simone Rossi's avatar
Simone Rossi committed
156

Simone Deparis's avatar
Simone Deparis committed
157
    inline void setupActivationExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
Simone Rossi's avatar
Simone Rossi committed
158
    {
Simone Deparis's avatar
Simone Deparis committed
159
        EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
Simone Rossi's avatar
Simone Rossi committed
160
161
    }

Simone Deparis's avatar
Simone Deparis committed
162
    inline void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Rossi's avatar
Simone Rossi committed
163
    {
Simone Deparis's avatar
Simone Deparis committed
164
165
166
167
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
168
169
    }

Simone Deparis's avatar
Simone Deparis committed
170
171
172
173
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
                         std::string mechanicsFileName  = "MechanicalSolution");
Simone Rossi's avatar
Simone Rossi committed
174

Simone Deparis's avatar
Simone Deparis committed
175
176
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize);
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list, std::string ionicModelName);
Simone Rossi's avatar
Simone Rossi committed
177

Simone Deparis's avatar
Simone Deparis committed
178
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
179

Simone Deparis's avatar
Simone Deparis committed
180
181
182
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
183

Simone Deparis's avatar
Simone Deparis committed
184
185
186
187
188
189
190
191
    inline void setupActivation (const MapEpetra& map)
    {
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "EMS - setting up activation solver\n";
        }
        M_activationModelPtr.reset ( new ActivationModel (map) );
    }
Simone Rossi's avatar
Simone Rossi committed
192

Simone Deparis's avatar
Simone Deparis committed
193
194
    void setup (GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize, commPtr_Type commPtr);
    void setup (GetPot& dataFile, Teuchos::ParameterList& list, std::string ionicModelName, commPtr_Type commPtr);
Simone Rossi's avatar
Simone Rossi committed
195

Simone Deparis's avatar
Simone Deparis committed
196
197
198
199
    inline void buildMechanicalSystem()
    {
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
200

Simone Deparis's avatar
Simone Deparis committed
201
202
203
204
    inline void buildElectroSystem()
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
205

Simone Deparis's avatar
Simone Deparis committed
206
207
208
209
210
    inline  void buildSystem()
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
211

Simone Deparis's avatar
Simone Deparis committed
212
213
214
215
    inline void initializeElectroVariables()
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
216

Simone Deparis's avatar
Simone Deparis committed
217
218
219
220
    inline void initialize()
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
221

Simone Deparis's avatar
Simone Deparis committed
222
223
224
225
    inline bcInterfacePtr_Type bcInterfacePtr()
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
226
227


Simone Deparis's avatar
Simone Deparis committed
228
229
230
231
    inline void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
232

Simone Deparis's avatar
Simone Deparis committed
233
234
235
236
    inline void setupElectroFiberVector ( VectorSmall<3>& fibers)
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
237
238


Simone Deparis's avatar
Simone Deparis committed
239
240
241
242
243
244
245
246
247
    inline void setupFiberVector ( Real fx, Real fy, Real fz )
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
        M_EMStructuralOperatorPtr -> EMMaterial() -> setFiberVectorPtr ( M_electroSolverPtr -> fiberPtr() );
    }
Simone Rossi's avatar
Simone Rossi committed
248

Simone Deparis's avatar
Simone Deparis committed
249
250
251
252
    inline electroSolverPtr_Type electroSolverPtr()
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
253

Simone Deparis's avatar
Simone Deparis committed
254
255
256
257
    inline structuralOperator_Type structuralOperatorPtr()
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
258

Simone Deparis's avatar
Simone Deparis committed
259
260
261
262
    inline activationModelPtr_Type  activationModelPtr()
    {
        return M_activationModelPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
263

Simone Deparis's avatar
Simone Deparis committed
264
    void saveSolution (Real time);
Simone Rossi's avatar
Simone Rossi committed
265

Simone Deparis's avatar
Simone Deparis committed
266
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
267

Simone Deparis's avatar
Simone Deparis committed
268
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
269

Simone Deparis's avatar
Simone Deparis committed
270
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
271

Simone Deparis's avatar
Simone Deparis committed
272
273
274
275
    inline void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
276
277


Simone Deparis's avatar
Simone Deparis committed
278
279
280
281
    inline void solveMechanics()
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
282

Simone Deparis's avatar
Simone Deparis committed
283
    inline void solveElectrophysiology (function_Type& stimulus, Real time = 0.0);
Simone Rossi's avatar
Simone Rossi committed
284
285


Simone Deparis's avatar
Simone Deparis committed
286
287
288
289
290
    inline void solveActivation (UInt index, Real dt);
    //  inline bcInterface_Type bcInterface()
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
291
292

protected:
Simone Deparis's avatar
Simone Deparis committed
293
294
295
296
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
297

Simone Deparis's avatar
Simone Deparis committed
298
299
300
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
Simone Rossi's avatar
Simone Rossi committed
301

Simone Deparis's avatar
Simone Deparis committed
302
303
304
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
    vectorPtr_Type                       M_activationTime;
Simone Rossi's avatar
Simone Rossi committed
305
306
307

    bool                                 M_oneWayCoupling;

Simone Deparis's avatar
Simone Deparis committed
308
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
309
310


Simone Rossi's avatar
Simone Rossi committed
311
};
Simone Rossi's avatar
Simone Rossi committed
312
313
314

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
315
316
317
318
319
320
321
322
323
324
325
326
327
328
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
EMSolver<Mesh, ElectroSolver, ActivationModel>::EMSolver() :
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
    M_activationTime     ( ),
    M_oneWayCoupling     (true),
    M_commPtr               ( )
Simone Rossi's avatar
Simone Rossi committed
329
330
331
332
333
334
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
335
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
Simone Deparis's avatar
Simone Deparis committed
336
337
338
EMSolver<Mesh, ElectroSolver, ActivationModel>::EMSolver (const EMSolver& solver) :
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
339
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
340
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
    M_electroExporterPtr ( solver.M_electroExporterPtr),
    M_activationExporterPtr ( solver.M_activationExporterPtr),
    M_mechanicsExporterPtr  ( solver.M_mechanicsExporterPtr),
    M_localMeshPtr      ( solver.M_localMeshPtr),
    M_fullMeshPtr      ( solver.M_fullMeshPtr),
    M_activationTime     ( solver.M_activationTime),
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
    M_commPtr               ( solver.M_commPtr)

{
}


/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
358
359
360
361
EMSolver<Mesh, ElectroSolver, ActivationModel>::setup ( GetPot& dataFile,
                                                        Teuchos::ParameterList& list,
                                                        short int ionicModelSize,
                                                        commPtr_Type commPtr)
Simone Rossi's avatar
Simone Rossi committed
362
{
Simone Deparis's avatar
Simone Deparis committed
363
364
365
366
    M_commPtr = commPtr;
    setupElectroSolver ( dataFile, list, ionicModelSize );
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
367
368
369
370
371
372
373
}


/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
374
375
376
377
EMSolver<Mesh, ElectroSolver, ActivationModel>::setup ( GetPot& dataFile,
                                                        Teuchos::ParameterList& list,
                                                        std::string ionicModelName,
                                                        commPtr_Type commPtr)
Simone Rossi's avatar
Simone Rossi committed
378
{
Simone Deparis's avatar
Simone Deparis committed
379
380
381
382
383
384
385
386
    M_commPtr = commPtr;
    setupElectroSolver ( dataFile, list, ionicModelName );
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - electro solver setup done! ";
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
387
388
}

Simone Rossi's avatar
Simone Rossi committed
389

Simone Rossi's avatar
Simone Rossi committed
390
/////////////////////
Simone Rossi's avatar
Simone Rossi committed
391
392
393
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
394
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list, short int ionicModelSize)
Simone Rossi's avatar
Simone Rossi committed
395
396
397
398
399
{
}

template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
400
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list,  std::string ionicModelName)
Simone Rossi's avatar
Simone Rossi committed
401
{
Simone Deparis's avatar
Simone Deparis committed
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
    ionicModelPtr_Type ionicModelPtr;
    ionicModelPtr.reset (ionicModel_Type::IonicModelFactory::instance().createObject ( ionicModelName ) );

    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up electrophysiology solver ";
    }

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
    if (M_localMeshPtr)
    {
        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
    }
    M_electroSolverPtr -> setParameters ( list );
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
430

Simone Rossi's avatar
Simone Rossi committed
431
432
433
}


Simone Rossi's avatar
Simone Rossi committed
434
435
436
437
/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
438
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
439
{
Simone Deparis's avatar
Simone Deparis committed
440
441
442
443
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
444
445
446
447
448
449
    boost::shared_ptr<StructuralConstitutiveLawData> dataStructure (new StructuralConstitutiveLawData( ) );
    dataStructure->setup (dataFile);

    std::string dOrder =  dataFile ( "solid/space_discretization/order", "P1");
    solidFESpacePtr_Type dFESpace ( new solidFESpace_Type (M_localMeshPtr, dOrder, 3, M_commPtr) );
    solidETFESpacePtr_Type dETFESpace ( new solidETFESpace_Type ( M_localMeshPtr,
Simone Deparis's avatar
Simone Deparis committed
450
451
452
453
454
455
456
457
458
459
460
461
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

    setupMechanicalBC (data_file_name, "solid",  dFESpace);
    M_EMStructuralOperatorPtr.reset (new structuralOperator_Type() );
    M_EMStructuralOperatorPtr->setup ( dataStructure,
                                       dFESpace,
                                       dETFESpace,
                                       M_bcInterfacePtr->handler(),
                                       M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
462
463
464
465
466
467
468
469
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);

}

/////////////////////
// Setting up the electrophysiology solver
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
470
471
472
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupMechanicalBC (std::string data_file_name,
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
473
{
Simone Deparis's avatar
Simone Deparis committed
474
475
476
477
478
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up bc interface\n";
    }
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
479
480
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
Simone Deparis's avatar
Simone Deparis committed
481
    M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
482
483
484
485
486
487
488
}


/////////////////////
//Setup exporters
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
489
490
491
492
EMSolver<Mesh, ElectroSolver, ActivationModel>::setupExporters (std::string problemFolder,
                                                                std::string electroFileName,
                                                                std::string activationFileName,
                                                                std::string mechanicsFileName)
Simone Rossi's avatar
Simone Rossi committed
493
{
Simone Deparis's avatar
Simone Deparis committed
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
                                             M_electroSolverPtr -> feSpacePtr(),
                                             M_activationModelPtr -> activationPtr(),
                                             UInt (0) );
    M_mechanicsExporterPtr.reset (new exporter_Type() );

    setupMechanicsExporter (problemFolder, mechanicsFileName);
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "displacement",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> displacementPtr(),
                                            UInt (0) );
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> fiberVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
520
521
522
523
}

template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
524
EMSolver<Mesh, ElectroSolver, ActivationModel>::saveSolution (Real time)
Simone Rossi's avatar
Simone Rossi committed
525
{
Simone Deparis's avatar
Simone Deparis committed
526
527
528
    M_electroExporterPtr -> postProcess (time);
    M_activationExporterPtr -> postProcess (time);
    M_mechanicsExporterPtr -> postProcess (time);
Simone Rossi's avatar
Simone Rossi committed
529
530
531
532
533
534
}

template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::closeExporters()
{
Simone Deparis's avatar
Simone Deparis committed
535
536
537
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
    M_mechanicsExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
538
539
540
541
542
543
544
545
}


///////////////////////////////
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::oneWayCoupling()
{
Simone Deparis's avatar
Simone Deparis committed
546
547
548
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
    M_EMStructuralOperatorPtr -> EMMaterial() -> setActivationPtr ( M_activationModelPtr -> activationPtr() );
Simone Rossi's avatar
Simone Rossi committed
549
550
551
552
553
554
}

template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
EMSolver<Mesh, ElectroSolver, ActivationModel>::twoWayCoupling()
{
Simone Deparis's avatar
Simone Deparis committed
555
556
557
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
    M_EMStructuralOperatorPtr -> EMMaterial() -> setActivationPtr ( M_activationModelPtr -> activationPtr() );
Simone Rossi's avatar
Simone Rossi committed
558
559
560
561
562
}

////////////////////////////
template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
563
EMSolver<Mesh, ElectroSolver, ActivationModel>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
564
{
Simone Deparis's avatar
Simone Deparis committed
565
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
566
567
568
569
570
571
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
    M_electroSolverPtr -> solveOneICIStep();
}

template<typename Mesh , typename ElectroSolver, typename ActivationModel>
void
Simone Deparis's avatar
Simone Deparis committed
572
EMSolver<Mesh, ElectroSolver, ActivationModel>::solveActivation (UInt index, Real dt)
Simone Rossi's avatar
Simone Rossi committed
573
{
Simone Deparis's avatar
Simone Deparis committed
574
    M_activationModelPtr -> solveModel ( * ( (M_electroSolverPtr -> globalSolution() ) [index] ), dt);
Simone Rossi's avatar
Simone Rossi committed
575
576
}

Simone Rossi's avatar
Simone Rossi committed
577
578
579

} // namespace LifeV

Simone Rossi's avatar
Simone Rossi committed
580
581
582



Simone Rossi's avatar
Simone Rossi committed
583
#endif //_MONODOMAINSOLVER_H_