EMSolver.hpp 25.8 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

Thomas Kummer's avatar
Thomas Kummer committed
186
187
188
189
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
190

Thomas Kummer's avatar
Thomas Kummer committed
191
192
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
193
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
194

Simone Deparis's avatar
Simone Deparis committed
195
196
197
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
198

Simone Rossi's avatar
Simone Rossi committed
199
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
200
201
202
203
204
    {
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "EMS - setting up activation solver\n";
        }
Simone Rossi's avatar
Simone Rossi committed
205
206
207
208
        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
209
    }
Simone Rossi's avatar
Simone Rossi committed
210

Thomas Kummer's avatar
Thomas Kummer committed
211
212
213
214
    void setup (GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setup (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
215

Thomas Kummer's avatar
Thomas Kummer committed
216
217
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
218
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
219
    {
Simone Rossi's avatar
Simone Rossi committed
220
221
    	//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
222
223
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
224

Simone Rossi's avatar
Simone Rossi committed
225
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
226
227
228
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
229

Simone Rossi's avatar
Simone Rossi committed
230
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
231
232
233
234
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
235

Simone Rossi's avatar
Simone Rossi committed
236
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
237
238
239
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
240

Simone Rossi's avatar
Simone Rossi committed
241
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
242
243
244
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
245

Simone Rossi's avatar
Simone Rossi committed
246
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
247
248
249
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
250
251


252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
    void setupFiberVector( const std::string& fileName,
						   const std::string& fieldName,
						   const std::string& postDir = "./",
						   const std::string& polynomialDegree = "P1"  )
    {
    	setupMechanicalFiberVector(fileName, fieldName, postDir, polynomialDegree);
    	M_electroSolverPtr->setFiberPtr(getMechanicsFibers());
    }

    void setupMechanicalFiberVector( const std::string& fileName,
    		                         const std::string& fieldName,
									 const std::string& postDir = "./",
                                     const std::string& polynomialDegree = "P1"  )
    {
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
    }

    void setupMechanicalSheetVector( const std::string& fileName,
    		                         const std::string& fieldName,
									 const std::string& postDir = "./",
                                     const std::string& polynomialDegree = "P1"  )
    {
        ElectrophysiologyUtility::importVectorField (getMechanicsSheets(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
    }

    void setupElectroFiberVector( const std::string& fileName,
								  const std::string& fieldName,
								  const std::string& postDir = "./",
								  const std::string& polynomialDegree = "P1"  )
    {
        ElectrophysiologyUtility::importVectorField (getElectroFibers(), fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
    }



Simone Rossi's avatar
Simone Rossi committed
287
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
288
289
290
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
291

Simone Rossi's avatar
Simone Rossi committed
292
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
293
294
295
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
296
297


Simone Rossi's avatar
Simone Rossi committed
298
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
299
300
301
302
303
304
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
305
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
306
    }
Simone Rossi's avatar
Simone Rossi committed
307

Simone Rossi's avatar
Simone Rossi committed
308
309
310
311
312
313
    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
314
315
316
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
317

Simone Rossi's avatar
Simone Rossi committed
318
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
319
320
321
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
322

Simone Rossi's avatar
Simone Rossi committed
323
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
324
325
326
    {
        return M_activationModelPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
327

thomaskummer's avatar
thomaskummer committed
328
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
329
330
331
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
332

Simone Deparis's avatar
Simone Deparis committed
333
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
334

Simone Deparis's avatar
Simone Deparis committed
335
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
336

Simone Deparis's avatar
Simone Deparis committed
337
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
338

Simone Rossi's avatar
Simone Rossi committed
339
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
340
341
342
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
343
344


Simone Rossi's avatar
Simone Rossi committed
345
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
346
347
348
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
349

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

Simone Rossi's avatar
Simone Rossi committed
352

Simone Rossi's avatar
Simone Rossi committed
353
    void solveActivation (Real dt);
Simone Rossi's avatar
Simone Rossi committed
354

Simone Rossi's avatar
Simone Rossi committed
355
356
357
358
359
360
361
362
363
364
365
366
367
    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
368
369
370
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
    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
386

387
    void showMe() const {}
388
389
390
391
392
393
394
395
396
397
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
398

399
    
Simone Rossi's avatar
Simone Rossi committed
400
protected:
Simone Rossi's avatar
Simone Rossi committed
401
public:
Simone Deparis's avatar
Simone Deparis committed
402
403
404
405
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
406

Simone Deparis's avatar
Simone Deparis committed
407
408
409
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
410
    exporterPtr_Type                     M_activationTimeExporterPtr;
Simone Rossi's avatar
Simone Rossi committed
411

Simone Deparis's avatar
Simone Deparis committed
412
413
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
414
415
    
    vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
416
417
418

    bool                                 M_oneWayCoupling;

Simone Deparis's avatar
Simone Deparis committed
419
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
420

Simone Rossi's avatar
Simone Rossi committed
421
422
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
423

Simone Rossi's avatar
Simone Rossi committed
424
};
Simone Rossi's avatar
Simone Rossi committed
425
426
427

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
428
429
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
430
431
432
433
434
435
436
437
438
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
439
    M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
440
    M_oneWayCoupling     (true),
Simone Rossi's avatar
Simone Rossi committed
441
442
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
443
444
445
446
447
448
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
449
450
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
451
452
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
453
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
454
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
455
456
457
458
459
    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),
thomaskummer's avatar
thomaskummer committed
460
    M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
461
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
Simone Rossi's avatar
Simone Rossi committed
462
463
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
464
465
466
467
468
469
470

{
}


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

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


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
476
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
477
void
Thomas Kummer's avatar
Thomas Kummer committed
478
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
479
{
Simone Rossi's avatar
Simone Rossi committed
480
481
    M_data.setup (dataFile);
    std::cout << "\nEMSolver - endtime = " << M_data.activationParameter<Real>("endtime");
Thomas Kummer's avatar
Thomas Kummer committed
482
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
483
484
    if (M_commPtr -> MyPID() == 0)
    {
Simone Rossi's avatar
Simone Rossi committed
485
        std::cout << "\nEMS - electro solver setup done! ";
Simone Deparis's avatar
Simone Deparis committed
486
487
488
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
489

Simone Rossi's avatar
Simone Rossi committed
490
491
}

Simone Rossi's avatar
Simone Rossi committed
492
493


Simone Rossi's avatar
Simone Rossi committed
494
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
495
void
Thomas Kummer's avatar
Thomas Kummer committed
496
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
497
{
Simone Deparis's avatar
Simone Deparis committed
498
499
500
501
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
502
503
504
505
506
	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
507
508
509
510
511
512
513
514

    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
515
516
517
518
    M_electroSolverPtr->setParameters();
    M_electroSolverPtr ->showParameters();
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
    M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
519
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
520

Thomas Kummer's avatar
Thomas Kummer committed
521
522
523
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
524
525
526
527
528
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
529
//    }
Simone Deparis's avatar
Simone Deparis committed
530
531
532
533
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
534

Simone Rossi's avatar
Simone Rossi committed
535
536
537
}


Simone Rossi's avatar
Simone Rossi committed
538
539
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
540
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
541
void
Simone Rossi's avatar
Simone Rossi committed
542
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
543
{
Simone Deparis's avatar
Simone Deparis committed
544
545
546
547
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
548
549
550
551
552
553
    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
554
555
556
557
558
559
560
561
562
563
564
565
                                                                  & (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
566
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
Simone Rossi's avatar
Simone Rossi committed
567
568
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);

Simone Rossi's avatar
Simone Rossi committed
569
570
571
572
573

}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
574
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
575
void
Simone Rossi's avatar
Simone Rossi committed
576
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
577
578
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
579
{
Simone Deparis's avatar
Simone Deparis committed
580
581
582
583
584
    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
585
586
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
587
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
588
589
590
591
592
}


/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
593
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
594
void
Simone Rossi's avatar
Simone Rossi committed
595
EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
Simone Deparis's avatar
Simone Deparis committed
596
597
598
                                                                std::string electroFileName,
                                                                std::string activationFileName,
                                                                std::string mechanicsFileName)
Simone Rossi's avatar
Simone Rossi committed
599
{
Simone Deparis's avatar
Simone Deparis committed
600
601
602
603
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
thomaskummer's avatar
thomaskummer committed
604
605
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
606
607
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
608
609
610
611
612
613
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
614
    // Activation
Simone Deparis's avatar
Simone Deparis committed
615
616
617
618
619
    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
620
                                             M_activationModelPtr -> fiberActivationPtr(),
Simone Deparis's avatar
Simone Deparis committed
621
622
623
                                             UInt (0) );
    M_mechanicsExporterPtr.reset (new exporter_Type() );

thomaskummer's avatar
thomaskummer committed
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
//    // Activation time
//    M_activationTimeExporterPtr.reset (new exporter_Type() );
//    setupActivationExporter (problemFolder, activationFileName );
//    
//    vectorPtr_Type M_activationTime ( new vector_Type ( M_electroSolverPtr()->potentialPtr() -> map() ) );
//    *activationTimeVector = -1.0;
//
//    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
//                                            "Activation Time",
//                                            M_electroSolverPtr -> feSpacePtr(),
//                                            M_activationModelPtr -> fiberActivationPtr(),
//                                            UInt (0) );
//    M_mechanicsExporterPtr.reset (new exporter_Type() );
//    
//    
//    activationTimeExporter.addVariable (ExporterData<mesh_Type>::ScalarField, "Activation Time", solver.electroSolverPtr()->feSpacePtr(), activationTimeVector, UInt (0) );
//    
//    activationTimeExporter.setMeshProcId (solver.localMeshPtr(), solver.comm()->MyPID() );
//    activationTimeExporter.setPrefix ("ActivationTime");
//    activationTimeExporter.setPostDir (problemFolder);

    // Mechanics
Simone Deparis's avatar
Simone Deparis committed
646
647
648
649
650
651
652
653
654
655
656
    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
657
658
659
660
661
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
662
663
}

thomaskummer's avatar
restart    
thomaskummer committed
664
665
666
667
668
669
670
671
672
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::setTimeIndex (const UInt& time)
{
    M_electroExporterPtr -> setTimeIndex (time);
    M_activationExporterPtr -> setTimeIndex (time);
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
673
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
674
void
thomaskummer's avatar
thomaskummer committed
675
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
676
{
thomaskummer's avatar
thomaskummer committed
677
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
678
679
    //if(M_activationExporterPtr) std::cout << "\nActivation exporter available.";
    //if(M_activationModelPtr -> fiberActivationPtr()) std::cout << "\nFiber activation exporter available.";
thomaskummer's avatar
thomaskummer committed
680
    M_activationExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
681
    //if(M_mechanicsExporterPtr) std::cout << "\nMechanics exporter available.";
thomaskummer's avatar
thomaskummer committed
682
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
683
684
}

Simone Rossi's avatar
Simone Rossi committed
685
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
686
void
Simone Rossi's avatar
Simone Rossi committed
687
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
688
{
Simone Deparis's avatar
Simone Deparis committed
689
690
691
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
    M_mechanicsExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
692
693
694
695
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
696
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
697
void
Simone Rossi's avatar
Simone Rossi committed
698
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
699
{
Simone Deparis's avatar
Simone Deparis committed
700
701
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
702
703
704
705
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
706
707
}

Simone Rossi's avatar
Simone Rossi committed
708
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
709
void
Simone Rossi's avatar
Simone Rossi committed
710
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
711
{
Simone Deparis's avatar
Simone Deparis committed
712
713
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
714
715
716
717
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
718
719
720
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
721
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
722
void
Simone Rossi's avatar
Simone Rossi committed
723
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
724
{
Simone Deparis's avatar
Simone Deparis committed
725
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
726
727
    auto v = M_electroSolverPtr -> diffusionTensor();
    std::cout << "\nEMS - " << v[0] << ", " << v[1] << ", " << v[2];
Simone Rossi's avatar
Simone Rossi committed
728
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
Thomas Kummer's avatar
Thomas Kummer committed
729
//    M_electroSolverPtr -> solveOneStepGatingVariablesRL();
Simone Rossi's avatar
Simone Rossi committed
730
731
732
    M_electroSolverPtr -> solveOneICIStep();
}

Simone Rossi's avatar
Simone Rossi committed
733
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
734
void
Simone Rossi's avatar
Simone Rossi committed
735
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
Simone Rossi's avatar
Simone Rossi committed
736
{
Simone Rossi's avatar
Simone Rossi committed
737
    M_activationModelPtr -> solveModel ( dt);
Simone Rossi's avatar
Simone Rossi committed
738
739
}

Simone Rossi's avatar
Simone Rossi committed
740
741
742

} // namespace LifeV

Simone Rossi's avatar
Simone Rossi committed
743
744
745



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