EMSolver.hpp 22.2 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
#include <lifev/em/solver/electrophysiology/IonicModelsList.hpp>
Simone Rossi's avatar
Simone Rossi committed
51
52
53
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>


Simone Rossi's avatar
Simone Rossi committed
54
55
//#include <lifev/em/solver/activation/activeStressModels/ActiveStressActivation.hpp>
#include <lifev/em/solver/activation/ActivationModelsList.hpp>
Simone Rossi's avatar
Simone Rossi committed
56

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>
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 Rossi's avatar
Simone Rossi committed
110
    typedef boost::shared_ptr<Activation>                     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
                                    const Real& z,
                                    const ID&   /*i*/ ) >       function_Type;



Simone Rossi's avatar
Simone Rossi committed
140
    EMSolver(commPtr_Type comm);
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 Rossi's avatar
Simone Rossi committed
146
    void loadMesh (std::string meshName, std::string meshPath)
Simone Deparis's avatar
Simone Deparis committed
147
    {
Simone Rossi's avatar
Simone Rossi committed
148
        std::cout << "EMS - Loading mesh\n";
Simone Rossi's avatar
Simone Rossi committed
149
        M_fullMeshPtr.reset( new Mesh() );
Simone Deparis's avatar
Simone Deparis committed
150
        MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
Simone Rossi's avatar
Simone Rossi committed
151
152
153
154
155
156
157
158
159
160
        if(M_commPtr)
        {
			M_localMeshPtr->setComm(M_commPtr);
			M_fullMeshPtr->setComm(M_commPtr);
        }
        else
        {
        	M_commPtr = M_localMeshPtr -> comm();
        }

Simone Deparis's avatar
Simone Deparis committed
161
    }
Simone Rossi's avatar
Simone Rossi committed
162

Simone Rossi's avatar
Simone Rossi committed
163
    void setupElectroExporter ( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
Simone Deparis's avatar
Simone Deparis committed
164
165
166
    {
        M_electroSolverPtr -> setupExporter (*M_electroExporterPtr, outputFileName, problemFolder);
    }
Simone Rossi's avatar
Simone Rossi committed
167

Simone Rossi's avatar
Simone Rossi committed
168
    void setupActivationExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
Simone Rossi's avatar
Simone Rossi committed
169
    {
Simone Deparis's avatar
Simone Deparis committed
170
        EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
Simone Rossi's avatar
Simone Rossi committed
171
172
    }

Simone Rossi's avatar
Simone Rossi committed
173
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Rossi's avatar
Simone Rossi committed
174
    {
Simone Deparis's avatar
Simone Deparis committed
175
176
177
178
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
179
180
    }

Simone Deparis's avatar
Simone Deparis committed
181
182
183
184
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
                         std::string mechanicsFileName  = "MechanicalSolution");
Simone Rossi's avatar
Simone Rossi committed
185

Simone Rossi's avatar
Simone Rossi committed
186
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list);
Simone Rossi's avatar
Simone Rossi committed
187

Simone Deparis's avatar
Simone Deparis committed
188
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
189

Simone Deparis's avatar
Simone Deparis committed
190
191
192
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
193

Simone Rossi's avatar
Simone Rossi committed
194
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
195
196
197
198
199
    {
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "EMS - setting up activation solver\n";
        }
Simone Rossi's avatar
Simone Rossi committed
200
201
202
203
        M_activationModelPtr.reset ( Activation::EMActivationFactory::instance().createObject ( M_data.activationParameter<std::string>( "ActivationModel" ) ) );
        M_activationModelPtr->setup(M_data, map);
        M_activationModelPtr->setVariablesPtr(*M_electroSolverPtr);
        M_activationModelPtr->setI4fPtr( M_EMStructuralOperatorPtr -> I4fPtr() );
Simone Deparis's avatar
Simone Deparis committed
204
    }
Simone Rossi's avatar
Simone Rossi committed
205

Simone Rossi's avatar
Simone Rossi committed
206
    void setup (GetPot& dataFile, Teuchos::ParameterList& list);
Simone Rossi's avatar
Simone Rossi committed
207

Simone Rossi's avatar
Simone Rossi committed
208
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
209
    {
Simone Rossi's avatar
Simone Rossi committed
210
211
    	//Here we call the buildSystem Of the Structural operator
    	// the coefficient is the density in front of the mass matrix
Simone Deparis's avatar
Simone Deparis committed
212
213
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
214

Simone Rossi's avatar
Simone Rossi committed
215
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
216
217
218
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
219

Simone Rossi's avatar
Simone Rossi committed
220
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
221
222
223
224
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
225

Simone Rossi's avatar
Simone Rossi committed
226
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
227
228
229
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
230

Simone Rossi's avatar
Simone Rossi committed
231
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
232
233
234
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
235

Simone Rossi's avatar
Simone Rossi committed
236
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
237
238
239
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
240
241


Simone Rossi's avatar
Simone Rossi committed
242
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
243
244
245
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
246

Simone Rossi's avatar
Simone Rossi committed
247
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
248
249
250
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
251
252


Simone Rossi's avatar
Simone Rossi committed
253
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
254
255
256
257
258
259
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
260
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
261
    }
Simone Rossi's avatar
Simone Rossi committed
262

Simone Rossi's avatar
Simone Rossi committed
263
264
265
266
267
268
    void setupSheetVector ( Real sx, Real sy, Real sz )
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupSheetVector(sx, sy, sz);
    }

    electroSolverPtr_Type electroSolverPtr()
Simone Deparis's avatar
Simone Deparis committed
269
270
271
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
272

Simone Rossi's avatar
Simone Rossi committed
273
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
274
275
276
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
277

Simone Rossi's avatar
Simone Rossi committed
278
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
279
280
281
    {
        return M_activationModelPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
282

Simone Deparis's avatar
Simone Deparis committed
283
    void saveSolution (Real time);
Simone Rossi's avatar
Simone Rossi committed
284

Simone Deparis's avatar
Simone Deparis committed
285
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
286

Simone Deparis's avatar
Simone Deparis committed
287
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
288

Simone Deparis's avatar
Simone Deparis committed
289
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
290

Simone Rossi's avatar
Simone Rossi committed
291
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
292
293
294
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
295
296


Simone Rossi's avatar
Simone Rossi committed
297
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
298
299
300
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
301

Simone Rossi's avatar
Simone Rossi committed
302
303
    void solveElectrophysiology (function_Type& stimulus, Real time = 0.0);

Simone Rossi's avatar
Simone Rossi committed
304

Simone Rossi's avatar
Simone Rossi committed
305
    void solveActivation (Real dt);
Simone Rossi's avatar
Simone Rossi committed
306

Simone Rossi's avatar
Simone Rossi committed
307
308
309
310
311
312
313
314
315
316
317
318
319
    vectorPtr_Type getElectroFibers()
    {
    	return M_electroSolverPtr -> fiberPtr();
    }
    vectorPtr_Type getMechanicsFibers()
    {
    	return M_EMStructuralOperatorPtr -> EMMaterial() -> fiberVectorPtr();
    }
    vectorPtr_Type getMechanicsSheets()
    {
    	return M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr();
    }
    //  bcInterface_Type bcInterface()
Simone Deparis's avatar
Simone Deparis committed
320
321
322
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
    EMData& data()
    {
    	return M_data;
    }


    commPtr_Type comm() const
    {
    	return M_commPtr;
    }

    void setComm(commPtr_Type comm)
    {
    	M_commPtr = comm;
    }
Simone Rossi's avatar
Simone Rossi committed
338
339

protected:
Simone Rossi's avatar
Simone Rossi committed
340
public:
Simone Deparis's avatar
Simone Deparis committed
341
342
343
344
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
345

Simone Deparis's avatar
Simone Deparis committed
346
347
348
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
Simone Rossi's avatar
Simone Rossi committed
349

Simone Deparis's avatar
Simone Deparis committed
350
351
352
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
    vectorPtr_Type                       M_activationTime;
Simone Rossi's avatar
Simone Rossi committed
353
354
355

    bool                                 M_oneWayCoupling;

Simone Deparis's avatar
Simone Deparis committed
356
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
357

Simone Rossi's avatar
Simone Rossi committed
358
359
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
360

Simone Rossi's avatar
Simone Rossi committed
361
};
Simone Rossi's avatar
Simone Rossi committed
362
363
364

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
365
366
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
367
368
369
370
371
372
373
374
375
376
377
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
    M_activationTime     ( ),
    M_oneWayCoupling     (true),
Simone Rossi's avatar
Simone Rossi committed
378
379
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
380
381
382
383
384
385
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
386
387
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
388
389
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
390
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
391
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
392
393
394
395
396
397
398
    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),
Simone Rossi's avatar
Simone Rossi committed
399
400
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
401
402
403
404
405
406
407

{
}


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
408

Simone Rossi's avatar
Simone Rossi committed
409
410
411
412


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
413
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
414
void
Simone Rossi's avatar
Simone Rossi committed
415
416
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile,
										Teuchos::ParameterList& list)
Simone Rossi's avatar
Simone Rossi committed
417
{
Simone Rossi's avatar
Simone Rossi committed
418
419
420
    M_data.setup (dataFile);
    std::cout << "\nEMSolver - endtime = " << M_data.activationParameter<Real>("endtime");
    setupElectroSolver ( dataFile, list );
Simone Deparis's avatar
Simone Deparis committed
421
422
    if (M_commPtr -> MyPID() == 0)
    {
Simone Rossi's avatar
Simone Rossi committed
423
        std::cout << "\nEMS - electro solver setup done! ";
Simone Deparis's avatar
Simone Deparis committed
424
425
426
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
427

Simone Rossi's avatar
Simone Rossi committed
428
429
}

Simone Rossi's avatar
Simone Rossi committed
430
431


Simone Rossi's avatar
Simone Rossi committed
432
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
433
void
Simone Rossi's avatar
Simone Rossi committed
434
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
Simone Rossi's avatar
Simone Rossi committed
435
{
Simone Deparis's avatar
Simone Deparis committed
436
437
438
439
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
440
441
442
443
444
	ionicModelPtr_Type ionicModelPtr;
	std::string ionicModelName = M_data.electroParameter<std::string>("IonicModel");
	ionicModelPtr.reset (ionicModel_Type::IonicModelFactory::instance().createObject ( ionicModelName ) );

//    M_electroSolverPtr.reset( new ElectroSolver ( meshName, meshPath, dataFile , ionicModelPtr ) );
Simone Deparis's avatar
Simone Deparis committed
445
446
447
448
449
450
451
452

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

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
Simone Rossi's avatar
Simone Rossi committed
453
454
455
456
457
458
    M_electroSolverPtr->setParameters();
    M_electroSolverPtr ->showParameters();
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
    M_electroSolverPtr ->showParameters();
    M_electroSolverPtr->init (M_commPtr);

Simone Deparis's avatar
Simone Deparis committed
459
460
461
462
463
464
465
466
467
468
469
470
471
    if (M_localMeshPtr)
    {
        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
    }
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
472

Simone Rossi's avatar
Simone Rossi committed
473
474
475
}


Simone Rossi's avatar
Simone Rossi committed
476
477
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
478
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
479
void
Simone Rossi's avatar
Simone Rossi committed
480
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
481
{
Simone Deparis's avatar
Simone Deparis committed
482
483
484
485
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
486
487
488
489
490
491
    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
492
493
494
495
496
497
498
499
500
501
502
503
                                                                  & (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
504
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
Simone Rossi's avatar
Simone Rossi committed
505
506
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);

Simone Rossi's avatar
Simone Rossi committed
507
508
509
510
511

}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
512
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
513
void
Simone Rossi's avatar
Simone Rossi committed
514
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
515
516
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
517
{
Simone Deparis's avatar
Simone Deparis committed
518
519
520
521
522
    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
523
524
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
Simone Deparis's avatar
Simone Deparis committed
525
    M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
526
527
528
529
530
}


/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
531
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
532
void
Simone Rossi's avatar
Simone Rossi committed
533
EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
Simone Deparis's avatar
Simone Deparis committed
534
535
536
                                                                std::string electroFileName,
                                                                std::string activationFileName,
                                                                std::string mechanicsFileName)
Simone Rossi's avatar
Simone Rossi committed
537
{
Simone Deparis's avatar
Simone Deparis committed
538
539
540
541
542
543
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
544
545
546
547
548
549
550
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );


Simone Deparis's avatar
Simone Deparis committed
551
552
553
554
555
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
                                             M_electroSolverPtr -> feSpacePtr(),
Simone Rossi's avatar
Simone Rossi committed
556
                                             M_activationModelPtr -> fiberActivationPtr(),
Simone Deparis's avatar
Simone Deparis committed
557
558
559
560
561
562
563
564
565
566
567
568
569
570
                                             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
571
572
573
574
575
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
576
577
}

Simone Rossi's avatar
Simone Rossi committed
578
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
579
void
Simone Rossi's avatar
Simone Rossi committed
580
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time)
Simone Rossi's avatar
Simone Rossi committed
581
{
Simone Deparis's avatar
Simone Deparis committed
582
    M_electroExporterPtr -> postProcess (time);
Simone Rossi's avatar
Simone Rossi committed
583
584
    if(M_activationExporterPtr) std::cout << "\nActivation exporter available.";
    if(M_activationModelPtr -> fiberActivationPtr()) std::cout << "\nFiber activation exporter available.";
Simone Deparis's avatar
Simone Deparis committed
585
    M_activationExporterPtr -> postProcess (time);
Simone Rossi's avatar
Simone Rossi committed
586
    if(M_mechanicsExporterPtr) std::cout << "\nMechanics exporter available.";
Simone Deparis's avatar
Simone Deparis committed
587
    M_mechanicsExporterPtr -> postProcess (time);
Simone Rossi's avatar
Simone Rossi committed
588
589
}

Simone Rossi's avatar
Simone Rossi committed
590
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
591
void
Simone Rossi's avatar
Simone Rossi committed
592
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
593
{
Simone Deparis's avatar
Simone Deparis committed
594
595
596
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
    M_mechanicsExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
597
598
599
600
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
601
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
602
void
Simone Rossi's avatar
Simone Rossi committed
603
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
604
{
Simone Deparis's avatar
Simone Deparis committed
605
606
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
607
608
609
610
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
611
612
}

Simone Rossi's avatar
Simone Rossi committed
613
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
614
void
Simone Rossi's avatar
Simone Rossi committed
615
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
616
{
Simone Deparis's avatar
Simone Deparis committed
617
618
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
619
620
621
622
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
623
624
625
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
626
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
627
void
Simone Rossi's avatar
Simone Rossi committed
628
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
629
{
Simone Deparis's avatar
Simone Deparis committed
630
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
631
632
    auto v = M_electroSolverPtr -> diffusionTensor();
    std::cout << "\nEMS - " << v[0] << ", " << v[1] << ", " << v[2];
Simone Rossi's avatar
Simone Rossi committed
633
634
635
636
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
    M_electroSolverPtr -> solveOneICIStep();
}

Simone Rossi's avatar
Simone Rossi committed
637
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
638
void
Simone Rossi's avatar
Simone Rossi committed
639
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
Simone Rossi's avatar
Simone Rossi committed
640
{
Simone Rossi's avatar
Simone Rossi committed
641
    M_activationModelPtr -> solveModel ( dt);
Simone Rossi's avatar
Simone Rossi committed
642
643
}

Simone Rossi's avatar
Simone Rossi committed
644
645
646

} // namespace LifeV

Simone Rossi's avatar
Simone Rossi committed
647
648
649



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