EMSolver.hpp 26.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
    }
thomaskummer's avatar
thomaskummer committed
172
173
174
175
176
    
    void setupActivationTimeExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationTimeSolution" )
    {
        EMUtility::setupExporter<Mesh> (*M_activationTimeExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
    }
Simone Rossi's avatar
Simone Rossi committed
177

Simone Rossi's avatar
Simone Rossi committed
178
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Rossi's avatar
Simone Rossi committed
179
    {
Simone Deparis's avatar
Simone Deparis committed
180
181
182
183
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
184
185
    }

Simone Deparis's avatar
Simone Deparis committed
186
187
188
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
thomaskummer's avatar
thomaskummer committed
189
                         std::string activationTimeFileName  = "ActivationTimeSolution",
Simone Deparis's avatar
Simone Deparis committed
190
                         std::string mechanicsFileName  = "MechanicalSolution");
Simone Rossi's avatar
Simone Rossi committed
191

Thomas Kummer's avatar
Thomas Kummer committed
192
193
194
195
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
196

Thomas Kummer's avatar
Thomas Kummer committed
197
198
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
199
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
200

Simone Deparis's avatar
Simone Deparis committed
201
202
203
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
204

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

Thomas Kummer's avatar
Thomas Kummer committed
217
218
219
220
    void setup (GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setup (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
221

Thomas Kummer's avatar
Thomas Kummer committed
222
223
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
224
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
225
    {
Simone Rossi's avatar
Simone Rossi committed
226
227
    	//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
228
229
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
230

Simone Rossi's avatar
Simone Rossi committed
231
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
232
233
234
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
235

Simone Rossi's avatar
Simone Rossi committed
236
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
237
238
239
240
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
241

Simone Rossi's avatar
Simone Rossi committed
242
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
243
244
245
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
246

Simone Rossi's avatar
Simone Rossi committed
247
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
248
249
250
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
251

Simone Rossi's avatar
Simone Rossi committed
252
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
253
254
255
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
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
287
288
289
290
291
292
    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
293
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
294
295
296
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
297

Simone Rossi's avatar
Simone Rossi committed
298
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
299
300
301
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
302
303


Simone Rossi's avatar
Simone Rossi committed
304
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
305
306
307
308
309
310
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
311
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
312
    }
Simone Rossi's avatar
Simone Rossi committed
313

Simone Rossi's avatar
Simone Rossi committed
314
315
316
317
318
319
    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
320
321
322
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
323

Simone Rossi's avatar
Simone Rossi committed
324
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
325
326
327
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
328

Simone Rossi's avatar
Simone Rossi committed
329
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
330
331
332
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
333
334
335
336
337
    
    vectorPtr_Type activationTimePtr()
    {
        return M_activationTimePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
338

thomaskummer's avatar
thomaskummer committed
339
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
340
341
342
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
343

Simone Deparis's avatar
Simone Deparis committed
344
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
345

Simone Deparis's avatar
Simone Deparis committed
346
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
347

Simone Deparis's avatar
Simone Deparis committed
348
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
349

Simone Rossi's avatar
Simone Rossi committed
350
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
351
352
353
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
354
355


Simone Rossi's avatar
Simone Rossi committed
356
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
357
358
359
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
360

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

Simone Rossi's avatar
Simone Rossi committed
363

Simone Rossi's avatar
Simone Rossi committed
364
    void solveActivation (Real dt);
Simone Rossi's avatar
Simone Rossi committed
365

Simone Rossi's avatar
Simone Rossi committed
366
367
368
369
370
371
372
373
374
375
376
377
378
    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
379
380
381
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
    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
397

398
    void showMe() const {}
399
400
401
402
403
404
405
406
407
408
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
409

410
    
Simone Rossi's avatar
Simone Rossi committed
411
protected:
Simone Rossi's avatar
Simone Rossi committed
412
public:
Simone Deparis's avatar
Simone Deparis committed
413
414
415
416
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
417

Simone Deparis's avatar
Simone Deparis committed
418
419
420
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
421
    exporterPtr_Type                     M_activationTimeExporterPtr;
Simone Rossi's avatar
Simone Rossi committed
422

Simone Deparis's avatar
Simone Deparis committed
423
424
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
425
426
    
    vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
427
428
429

    bool                                 M_oneWayCoupling;

Simone Deparis's avatar
Simone Deparis committed
430
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
431

Simone Rossi's avatar
Simone Rossi committed
432
433
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
434

Simone Rossi's avatar
Simone Rossi committed
435
};
Simone Rossi's avatar
Simone Rossi committed
436
437
438

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
439
440
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
441
442
443
444
445
446
447
448
449
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
450
    M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
451
    M_oneWayCoupling     (true),
Simone Rossi's avatar
Simone Rossi committed
452
453
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
454
455
456
457
458
459
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
460
461
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
462
463
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
464
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
465
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
466
467
468
469
470
    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
471
    M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
472
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
Simone Rossi's avatar
Simone Rossi committed
473
474
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
475
476
477
478
479
480
481

{
}


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

Simone Rossi's avatar
Simone Rossi committed
483
484
485
486


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
487
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
488
void
Thomas Kummer's avatar
Thomas Kummer committed
489
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
490
{
Simone Rossi's avatar
Simone Rossi committed
491
492
    M_data.setup (dataFile);
    std::cout << "\nEMSolver - endtime = " << M_data.activationParameter<Real>("endtime");
Thomas Kummer's avatar
Thomas Kummer committed
493
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
494
495
    if (M_commPtr -> MyPID() == 0)
    {
Simone Rossi's avatar
Simone Rossi committed
496
        std::cout << "\nEMS - electro solver setup done! ";
Simone Deparis's avatar
Simone Deparis committed
497
498
499
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
500

Simone Rossi's avatar
Simone Rossi committed
501
502
}

Simone Rossi's avatar
Simone Rossi committed
503
504


Simone Rossi's avatar
Simone Rossi committed
505
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
506
void
Thomas Kummer's avatar
Thomas Kummer committed
507
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
508
{
Simone Deparis's avatar
Simone Deparis committed
509
510
511
512
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
513
514
515
516
517
	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
518
519
520
521
522
523
524
525

    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
526
527
528
529
    M_electroSolverPtr->setParameters();
    M_electroSolverPtr ->showParameters();
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
    M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
530
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
531

Thomas Kummer's avatar
Thomas Kummer committed
532
533
534
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
535
536
537
538
539
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
540
//    }
Simone Deparis's avatar
Simone Deparis committed
541
542
543
544
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
545

Simone Rossi's avatar
Simone Rossi committed
546
547
548
}


Simone Rossi's avatar
Simone Rossi committed
549
550
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
551
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
552
void
Simone Rossi's avatar
Simone Rossi committed
553
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
554
{
Simone Deparis's avatar
Simone Deparis committed
555
556
557
558
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
559
560
561
562
563
564
    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
565
566
567
568
569
570
571
572
573
574
575
576
                                                                  & (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
577
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
Simone Rossi's avatar
Simone Rossi committed
578
579
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);

Simone Rossi's avatar
Simone Rossi committed
580
581
582
583
584

}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
585
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
586
void
Simone Rossi's avatar
Simone Rossi committed
587
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
588
589
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
590
{
Simone Deparis's avatar
Simone Deparis committed
591
592
593
594
595
    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
596
597
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
598
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
599
600
601
602
603
}


/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
604
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
605
void
Simone Rossi's avatar
Simone Rossi committed
606
EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
Simone Deparis's avatar
Simone Deparis committed
607
608
                                                                std::string electroFileName,
                                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
609
                                                                std::string activationTimeFileName,
Simone Deparis's avatar
Simone Deparis committed
610
                                                                std::string mechanicsFileName)
Simone Rossi's avatar
Simone Rossi committed
611
{
Simone Deparis's avatar
Simone Deparis committed
612
613
614
615
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
thomaskummer's avatar
thomaskummer committed
616
617
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
618
619
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
620
621
622
623
624
625
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
626
    // Activation
Simone Deparis's avatar
Simone Deparis committed
627
628
629
630
631
    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
632
                                             M_activationModelPtr -> fiberActivationPtr(),
Simone Deparis's avatar
Simone Deparis committed
633
634
635
                                             UInt (0) );
    M_mechanicsExporterPtr.reset (new exporter_Type() );

thomaskummer's avatar
thomaskummer committed
636
637
638
639
640
641
    // Activation time
    M_activationTimeExporterPtr.reset (new exporter_Type() );
    setupActivationTimeExporter (problemFolder, activationTimeFileName );
    
    M_activationTimePtr.reset (new vector_Type ( M_electroSolverPtr->potentialPtr() -> map() ) );
    *M_activationTimePtr = -1.0;
thomaskummer's avatar
thomaskummer committed
642

thomaskummer's avatar
thomaskummer committed
643
    M_activationTimeExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
644
645
646
647
648
                                            "Activation Time",
                                            M_electroSolverPtr -> feSpacePtr(),
                                            M_activationTimePtr,
                                            UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
649
    // Mechanics
thomaskummer's avatar
thomaskummer committed
650
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
651
652
653
654
655
656
657
658
659
660
661
    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
662
663
664
665
666
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
667
668
}

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

Simone Rossi's avatar
Simone Rossi committed
692
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
693
void
Simone Rossi's avatar
Simone Rossi committed
694
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
695
{
Simone Deparis's avatar
Simone Deparis committed
696
697
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
698
    M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
699
    M_mechanicsExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
700
701
702
703
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
704
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
705
void
Simone Rossi's avatar
Simone Rossi committed
706
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
707
{
Simone Deparis's avatar
Simone Deparis committed
708
709
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
710
711
712
713
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
714
715
}

Simone Rossi's avatar
Simone Rossi committed
716
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
717
void
Simone Rossi's avatar
Simone Rossi committed
718
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
719
{
Simone Deparis's avatar
Simone Deparis committed
720
721
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
722
723
724
725
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
726
727
728
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
729
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
730
void
Simone Rossi's avatar
Simone Rossi committed
731
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
732
{
Simone Deparis's avatar
Simone Deparis committed
733
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
734
735
    auto v = M_electroSolverPtr -> diffusionTensor();
    std::cout << "\nEMS - " << v[0] << ", " << v[1] << ", " << v[2];
Simone Rossi's avatar
Simone Rossi committed
736
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
Thomas Kummer's avatar
Thomas Kummer committed
737
//    M_electroSolverPtr -> solveOneStepGatingVariablesRL();
Simone Rossi's avatar
Simone Rossi committed
738
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
739
    M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
Simone Rossi's avatar
Simone Rossi committed
740
741
}

Simone Rossi's avatar
Simone Rossi committed
742
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
743
void
Simone Rossi's avatar
Simone Rossi committed
744
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
Simone Rossi's avatar
Simone Rossi committed
745
{
Simone Rossi's avatar
Simone Rossi committed
746
    M_activationModelPtr -> solveModel ( dt);
Simone Rossi's avatar
Simone Rossi committed
747
748
}

Simone Rossi's avatar
Simone Rossi committed
749
750
751

} // namespace LifeV

Simone Rossi's avatar
Simone Rossi committed
752
753
754



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