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


thomaskummer's avatar
thomaskummer committed
54
55
#include <lifev/structure/solver/WallTensionEstimator.hpp>

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

Simone Rossi's avatar
Simone Rossi committed
59
60
61
62
63
64
65
66
67
68
69
//
//#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
70
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
Simone Rossi's avatar
Simone Rossi committed
71
72
73
//
//
//#include <lifev/em/solver/EMEvaluate.hpp>
Simone Rossi's avatar
Simone Rossi committed
74
75
76
77
78
79

namespace LifeV
{

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

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

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

Simone Deparis's avatar
Simone Deparis committed
88
    typedef Epetra_Comm                                       comm_Type;
Simone Rossi's avatar
Simone Rossi committed
89

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

Simone Deparis's avatar
Simone Deparis committed
92
    typedef VectorEpetra                                      vector_Type;
thomaskummer's avatar
thomaskummer committed
93
    
Simone Deparis's avatar
Simone Deparis committed
94
    typedef StructuralConstitutiveLawData                     structureData_Type;
Simone Rossi's avatar
Simone Rossi committed
95

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

Simone Deparis's avatar
Simone Deparis committed
98
    typedef EMStructuralOperator< mesh_Type >                 structuralOperator_Type;
Simone Rossi's avatar
Simone Rossi committed
99

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

Simone Deparis's avatar
Simone Deparis committed
102
    typedef BCHandler                                          bc_Type;
Simone Rossi's avatar
Simone Rossi committed
103

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

Simone Deparis's avatar
Simone Deparis committed
106
    typedef StructuralOperator< mesh_Type >                   physicalSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
107

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

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

Simone Rossi's avatar
Simone Rossi committed
112
    typedef boost::shared_ptr<Activation>                     activationModelPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
113

Simone Deparis's avatar
Simone Deparis committed
114
    typedef ElectroSolver                                      electroSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
115

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

Simone Deparis's avatar
Simone Deparis committed
118
    typedef ElectroIonicModel                                  ionicModel_Type;
Simone Rossi's avatar
Simone Rossi committed
119

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

Simone Deparis's avatar
Simone Deparis committed
122
    typedef ExporterHDF5< Mesh >                               exporter_Type;
Simone Rossi's avatar
Simone Rossi committed
123

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

    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
135
136
                                    const Real&    x,
                                    const Real&    y,
Simone Rossi's avatar
Simone Rossi committed
137
138
139
140
141
                                    const Real& z,
                                    const ID&   /*i*/ ) >       function_Type;



Simone Rossi's avatar
Simone Rossi committed
142
    EMSolver(commPtr_Type comm);
Simone Rossi's avatar
Simone Rossi committed
143

Simone Deparis's avatar
Simone Deparis committed
144
    EMSolver (const EMSolver& solver);
Simone Rossi's avatar
Simone Rossi committed
145

Simone Rossi's avatar
Simone Rossi committed
146

Simone Rossi's avatar
Simone Rossi committed
147

Simone Rossi's avatar
Simone Rossi committed
148
    void loadMesh (std::string meshName, std::string meshPath)
Simone Deparis's avatar
Simone Deparis committed
149
    {
Simone Rossi's avatar
Simone Rossi committed
150
        std::cout << "EMS - Loading mesh\n";
Simone Rossi's avatar
Simone Rossi committed
151
        M_fullMeshPtr.reset( new Mesh() );
Simone Deparis's avatar
Simone Deparis committed
152
        MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
Simone Rossi's avatar
Simone Rossi committed
153
154
155
156
157
158
159
160
161
162
        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
163
    }
Simone Rossi's avatar
Simone Rossi committed
164

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

Simone Rossi's avatar
Simone Rossi committed
170
    void setupActivationExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
Simone Rossi's avatar
Simone Rossi committed
171
    {
Simone Deparis's avatar
Simone Deparis committed
172
        EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
Simone Rossi's avatar
Simone Rossi committed
173
    }
thomaskummer's avatar
thomaskummer committed
174
175
176
177
178
    
    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
179

thomaskummer's avatar
thomaskummer committed
180
181
182
183
    void setupVonMisesStressExporter ( std::string problemFolder = "./", std::string outputFileName = "VonMisesStress" )
    {
        EMUtility::setupExporter<Mesh> (*M_vonMisesStressExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
    }
thomaskummer's avatar
thomaskummer committed
184
185
186
187
188
189
190
191
192
193
    
    void setupVonMisesStressExporterP ( std::string problemFolder = "./", std::string outputFileName = "VonMisesStressP" )
    {
        EMUtility::setupExporter<Mesh> (*M_vonMisesStressExporterPtrP, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
    }
    
    void setupVonMisesStressExporterA ( std::string problemFolder = "./", std::string outputFileName = "VonMisesStressA" )
    {
        EMUtility::setupExporter<Mesh> (*M_vonMisesStressExporterPtrA, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
    }
thomaskummer's avatar
thomaskummer committed
194

Simone Rossi's avatar
Simone Rossi committed
195
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Rossi's avatar
Simone Rossi committed
196
    {
Simone Deparis's avatar
Simone Deparis committed
197
198
199
200
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
201
202
    }

203
204
    void importHdf5 ();

Simone Deparis's avatar
Simone Deparis committed
205
206
207
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
thomaskummer's avatar
thomaskummer committed
208
                         std::string activationTimeFileName  = "ActivationTimeSolution",
thomaskummer's avatar
thomaskummer committed
209
                         std::string mechanicsFileName  = "MechanicalSolution",
thomaskummer's avatar
thomaskummer committed
210
211
212
                         std::string vonMisesStressFileName  = "VonMisesStress",
                         std::string vonMisesStressFileNameP  = "VonMisesStressP",
                         std::string vonMisesStressFileNameA  = "VonMisesStressA");
Simone Rossi's avatar
Simone Rossi committed
213

Thomas Kummer's avatar
Thomas Kummer committed
214
215
216
217
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
218

Thomas Kummer's avatar
Thomas Kummer committed
219
220
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
221
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
222

Simone Deparis's avatar
Simone Deparis committed
223
224
225
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
226

Simone Rossi's avatar
Simone Rossi committed
227
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
228
229
230
231
232
    {
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "EMS - setting up activation solver\n";
        }
Simone Rossi's avatar
Simone Rossi committed
233
234
235
236
        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
237
    }
Simone Rossi's avatar
Simone Rossi committed
238

Thomas Kummer's avatar
Thomas Kummer committed
239
240
241
242
    void setup (GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setup (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
243

Thomas Kummer's avatar
Thomas Kummer committed
244
245
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
246
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
247
    {
Simone Rossi's avatar
Simone Rossi committed
248
249
    	//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
250
251
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
252

Simone Rossi's avatar
Simone Rossi committed
253
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
254
255
256
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
257

Simone Rossi's avatar
Simone Rossi committed
258
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
259
260
261
262
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
263

Simone Rossi's avatar
Simone Rossi committed
264
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
265
266
267
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
268

Simone Rossi's avatar
Simone Rossi committed
269
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
270
271
272
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
273

Simone Rossi's avatar
Simone Rossi committed
274
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
275
276
277
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
278
279


280
281
282
283
284
285
286
287
288
289
290
291
292
293
    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"  )
    {
thomaskummer's avatar
thomaskummer committed
294
        
295
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
296
297
298
299
300
301
302
303
304
        
//        if ( polynomialDegree == "P1" )
//        {
//            FESpace<RegionMesh<LinearTetra> , MapEpetra > p2FESpace ( fullMeshPtr(), "P2", 1, fullMeshPtr() -> comm() );
//            vector_Type p2FibersRep (*getMechanicsFibers(), Repeated);
//            vector_Type p1FibersRep = M_EMStructuralOperatorPtr -> dispFESpacePtr() -> feToFEInterpolate(p2FESpace, p2FibersRep);
//            getMechanicsFibers().reset ( new vector_Type (p1FibersRep, Unique) );
//        }

305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
    }

    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
325
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
326
327
328
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
329

Simone Rossi's avatar
Simone Rossi committed
330
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
331
332
333
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
334
335


Simone Rossi's avatar
Simone Rossi committed
336
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
337
338
339
340
341
342
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
343
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
344
    }
Simone Rossi's avatar
Simone Rossi committed
345

Simone Rossi's avatar
Simone Rossi committed
346
347
348
349
350
351
    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
352
353
354
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
355

Simone Rossi's avatar
Simone Rossi committed
356
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
357
358
359
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
360

Simone Rossi's avatar
Simone Rossi committed
361
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
362
363
364
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
365
366
367
368
369
    
    vectorPtr_Type activationTimePtr()
    {
        return M_activationTimePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
370

thomaskummer's avatar
thomaskummer committed
371
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
372
373
374
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
375

Simone Deparis's avatar
Simone Deparis committed
376
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
377

Simone Deparis's avatar
Simone Deparis committed
378
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
379

Simone Deparis's avatar
Simone Deparis committed
380
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
381

Simone Rossi's avatar
Simone Rossi committed
382
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
383
384
385
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
386
387


Simone Rossi's avatar
Simone Rossi committed
388
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
389
390
391
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
392

393
394
395
396
397
    void solveMechanicsLin()
    {
        M_EMStructuralOperatorPtr -> solveLin ();
    }

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

Simone Rossi's avatar
Simone Rossi committed
400

Simone Rossi's avatar
Simone Rossi committed
401
    void solveActivation (Real dt);
Simone Rossi's avatar
Simone Rossi committed
402

Simone Rossi's avatar
Simone Rossi committed
403
404
405
406
407
408
409
410
411
412
413
414
415
    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
416
417
418
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
    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
434

435
    void showMe() const {}
436
437
438
439
440
441
442
443
444
445
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
446

447
    
Simone Rossi's avatar
Simone Rossi committed
448
protected:
Simone Rossi's avatar
Simone Rossi committed
449
public:
Simone Deparis's avatar
Simone Deparis committed
450
451
452
453
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
454

Simone Deparis's avatar
Simone Deparis committed
455
456
457
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
458
    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
459
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
thomaskummer's avatar
thomaskummer committed
460
461
    exporterPtr_Type                     M_vonMisesStressExporterPtrP;
    exporterPtr_Type                     M_vonMisesStressExporterPtrA;
thomaskummer's avatar
thomaskummer committed
462
    
Simone Deparis's avatar
Simone Deparis committed
463
464
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
465
466
    
    vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
467
468

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
469
    
thomaskummer's avatar
thomaskummer committed
470
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
thomaskummer's avatar
thomaskummer committed
471
472
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
thomaskummer's avatar
thomaskummer committed
473

Simone Rossi's avatar
Simone Rossi committed
474

Simone Deparis's avatar
Simone Deparis committed
475
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
476

Simone Rossi's avatar
Simone Rossi committed
477
478
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
479

Simone Rossi's avatar
Simone Rossi committed
480
};
Simone Rossi's avatar
Simone Rossi committed
481
482
483

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
484
485
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
486
487
488
489
490
491
492
493
494
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
495
    M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
496
    M_oneWayCoupling     (true),
thomaskummer's avatar
thomaskummer committed
497
    M_wteTotal ( ),
thomaskummer's avatar
thomaskummer committed
498
499
//    M_wtePassive ( ),
//    M_wteActive ( ),
Simone Rossi's avatar
Simone Rossi committed
500
501
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
502
503
504
505
506
507
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
508
509
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
510
511
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
512
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
513
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
514
515
516
517
518
    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
519
    M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
520
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
thomaskummer's avatar
thomaskummer committed
521
    M_wteTotal                   (solver.M_wteTotal),
thomaskummer's avatar
thomaskummer committed
522
523
//    M_wtePassive                   (solver.M_wtePassive),
//    M_wteActive                   (solver.M_wteActive),
Simone Rossi's avatar
Simone Rossi committed
524
525
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
526
527
528
529
530
531
532

{
}


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

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


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
538
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
539
void
Thomas Kummer's avatar
Thomas Kummer committed
540
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
541
{
Simone Rossi's avatar
Simone Rossi committed
542
543
    M_data.setup (dataFile);
    std::cout << "\nEMSolver - endtime = " << M_data.activationParameter<Real>("endtime");
Thomas Kummer's avatar
Thomas Kummer committed
544
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
545
546
    if (M_commPtr -> MyPID() == 0)
    {
Simone Rossi's avatar
Simone Rossi committed
547
        std::cout << "\nEMS - electro solver setup done! ";
Simone Deparis's avatar
Simone Deparis committed
548
549
550
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
551

Simone Rossi's avatar
Simone Rossi committed
552
553
}

Simone Rossi's avatar
Simone Rossi committed
554
555


Simone Rossi's avatar
Simone Rossi committed
556
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
557
void
Thomas Kummer's avatar
Thomas Kummer committed
558
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
559
{
Simone Deparis's avatar
Simone Deparis committed
560
561
562
563
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
564
565
566
567
568
	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
569
570
571
572
573
574
575
576

    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
577
578
579
580
    M_electroSolverPtr->setParameters();
    M_electroSolverPtr ->showParameters();
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
    M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
581
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
582

Thomas Kummer's avatar
Thomas Kummer committed
583
584
585
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
586
587
588
589
590
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
591
//    }
Simone Deparis's avatar
Simone Deparis committed
592
593
594
595
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
596

Simone Rossi's avatar
Simone Rossi committed
597
598
599
}


Simone Rossi's avatar
Simone Rossi committed
600
601
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
602
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
603
void
Simone Rossi's avatar
Simone Rossi committed
604
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
605
{
Simone Deparis's avatar
Simone Deparis committed
606
607
608
609
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
610
611
612
613
614
615
    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
616
617
618
619
620
621
622
623
624
625
626
627
                                                                  & (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
628
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
Simone Rossi's avatar
Simone Rossi committed
629
630
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);

thomaskummer's avatar
thomaskummer committed
631
    M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
thomaskummer's avatar
thomaskummer committed
632
633
//    M_wtePassive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "passive");
//    M_wteActive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "active");
Simone Rossi's avatar
Simone Rossi committed
634
635
636
637
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
638
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
639
void
Simone Rossi's avatar
Simone Rossi committed
640
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
641
642
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
643
{
Simone Deparis's avatar
Simone Deparis committed
644
645
646
647
648
    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
649
650
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
651
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
652
653
}

654
655
656
657
658
659
660
    
/////////////////////
//restart hdf5
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::importHdf5 ()
{
thomaskummer's avatar
thomaskummer committed
661
    M_electroExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
662
    M_activationExporterPtr -> importHdf5 (0.);
thomaskummer's avatar
thomaskummer committed
663
664
665
666
667
    M_activationTimeExporterPtr -> importHdf5 (30);
    M_vonMisesStressExporterPtr -> importHdf5 (30);
    M_vonMisesStressExporterPtrP -> importHdf5 (30);
    M_vonMisesStressExporterPtrA -> importHdf5 (30);
    M_mechanicsExporterPtr -> importHdf5 (30);
668
669
}
                                                   
Simone Rossi's avatar
Simone Rossi committed
670
671
672

/////////////////////
//Setup exporters
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
676
677
678
679
EMSolver<Mesh, ElectroSolver>::setupExporters ( std::string problemFolder,
                                                std::string electroFileName,
                                                std::string activationFileName,
                                                std::string activationTimeFileName,
                                                std::string mechanicsFileName,
thomaskummer's avatar
thomaskummer committed
680
681
682
                                                std::string vonMisesStressFileName,
                                               std::string vonMisesStressFileNameP,
                                               std::string vonMisesStressFileNameA)
Simone Rossi's avatar
Simone Rossi committed
683
{
Simone Deparis's avatar
Simone Deparis committed
684
685
686
687
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
thomaskummer's avatar
thomaskummer committed
688
689
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
690
691
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
692
693
694
695
696
697
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
698
    // Activation
Simone Deparis's avatar
Simone Deparis committed
699
700
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
thomaskummer's avatar
thomaskummer committed
701
702
703
704
705
//    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
//                                             "Activation",
//                                              M_electroSolverPtr -> feSpacePtr(),
//                                             M_activationModelPtr -> fiberActivationPtr(),
//                                             UInt (0) );
Simone Deparis's avatar
Simone Deparis committed
706

thomaskummer's avatar
thomaskummer committed
707
708
709
710
711
712
    // 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
713

thomaskummer's avatar
thomaskummer committed
714
    M_activationTimeExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
715
716
717
718
719
                                            "Activation Time",
                                            M_electroSolverPtr -> feSpacePtr(),
                                            M_activationTimePtr,
                                            UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
720
721
722
723
724
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
725
726
727
728
729
                                                "Von Mises Stress Total",
                                                M_electroSolverPtr -> feSpacePtr(),
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
    M_vonMisesStressExporterPtrP.reset (new exporter_Type() );
    setupVonMisesStressExporterP (problemFolder, vonMisesStressFileNameP );
    
    M_vonMisesStressExporterPtrP -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                                "Von Mises Stress Total P",
                                                M_electroSolverPtr -> feSpacePtr(),
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );

    M_vonMisesStressExporterPtrA.reset (new exporter_Type() );
    setupVonMisesStressExporterA (problemFolder, vonMisesStressFileNameA );
    
    M_vonMisesStressExporterPtrA -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                                "Von Mises Stress Total A",
                                                M_electroSolverPtr -> feSpacePtr(),
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );

thomaskummer's avatar
thomaskummer committed
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "X Stress Total",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteTotal.sigmaXPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Y Stress Total",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteTotal.sigmaYPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Z Stress Total",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteTotal.sigmaZPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
//                                                "Von Mises Stress Passive",
//                                                M_electroSolverPtr -> feSpacePtr(),
//                                                M_wtePassive.vonMisesStressPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "X Stress Passive",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wtePassive.sigmaXPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Y Stress Passive",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wtePassive.sigmaYPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Z Stress Passive",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wtePassive.sigmaZPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
//                                                "Von Mises Stress Active",
//                                                M_electroSolverPtr -> feSpacePtr(),
//                                                M_wteActive.vonMisesStressPtr(),
//                                                UInt (0) );
//    
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "X Stress Active",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteActive.sigmaXPtr(),
//                                                UInt (0) );
//
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Y Stress Active",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteActive.sigmaYPtr(),
//                                                UInt (0) );
//    
//    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
//                                                "Z Stress Active",
//                                                M_EMStructuralOperatorPtr -> dispFESpacePtr(),
//                                                M_wteActive.sigmaZPtr(),
//                                                UInt (0) );
thomaskummer's avatar
thomaskummer committed
813
    
thomaskummer's avatar
thomaskummer committed
814
    // Mechanics
thomaskummer's avatar
thomaskummer committed
815
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
816
817
818
819
820
821
822
823
824
825
826
    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
827
828
829
830
831
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
832
833
}

thomaskummer's avatar
restart    
thomaskummer committed
834
835
836
837
838
839
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
840
    M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
841
    M_vonMisesStressExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
842
843
    M_vonMisesStressExporterPtrP -> setTimeIndex (time);
    M_vonMisesStressExporterPtrA -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
844
845
846
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
847
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
848
void
thomaskummer's avatar
thomaskummer committed
849
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
850
{
thomaskummer's avatar
thomaskummer committed
851
852
853
854
855
856
857
858
859
860
861
862
863
//    M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
//    
//    M_wteTotal.setStressType ( "total" );
//    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
//    M_vonMisesStressExporterPtr -> postProcess (time);
//    
//    M_wteTotal.setStressType ( "passiv" );
//    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
//    M_vonMisesStressExporterPtrP -> postProcess (time);
//
//    M_wteTotal.setStressType ( "active" );
//    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
//    M_vonMisesStressExporterPtrA -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
864

thomaskummer's avatar
thomaskummer committed
865
866
867
868
//    M_wtePassive.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
//    M_wtePassive.analyzeTensionsRecoveryVonMisesStress();
//    M_wteActive.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
//    M_wteActive.analyzeTensionsRecoveryVonMisesStress();
thomaskummer's avatar
thomaskummer committed
869

thomaskummer's avatar
thomaskummer committed
870
    
thomaskummer's avatar
thomaskummer committed
871
    //M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
872
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
873
874
    //M_activationTimeExporterPtr -> postProcess (time);
    //M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
875
876
}

Simone Rossi's avatar
Simone Rossi committed
877
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
878
void
Simone Rossi's avatar
Simone Rossi committed
879
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
880
{
Simone Deparis's avatar
Simone Deparis committed
881
882
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
883
    M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
884
    M_mechanicsExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
885
    M_vonMisesStressExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
886
887
888
889
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
890
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
891
void
Simone Rossi's avatar
Simone Rossi committed
892
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
893
{
Simone Deparis's avatar
Simone Deparis committed
894
895
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
896
897
898
899
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
900
901
}

Simone Rossi's avatar
Simone Rossi committed
902
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
903
void
Simone Rossi's avatar
Simone Rossi committed
904
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
905
{
Simone Deparis's avatar
Simone Deparis committed
906
907
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
908
909
910
911
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
912
913
914
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
915
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
916
void
Simone Rossi's avatar
Simone Rossi committed
917
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
918
{
Simone Deparis's avatar
Simone Deparis committed
919
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
920
    auto v = M_electroSolverPtr -> diffusionTensor();
thomaskummer's avatar
thomaskummer committed
921
    //std::cout << "\nEMS - " << v[0] << ", " << v[1] << ", " << v[2];
Simone Rossi's avatar
Simone Rossi committed
922
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
Thomas Kummer's avatar
Thomas Kummer committed
923
//    M_electroSolverPtr -> solveOneStepGatingVariablesRL();
Simone Rossi's avatar
Simone Rossi committed
924
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
925
    M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
Simone Rossi's avatar
Simone Rossi committed
926
927
}

Simone Rossi's avatar
Simone Rossi committed
928
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
929
void
Simone Rossi's avatar
Simone Rossi committed
930
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
Simone Rossi's avatar
Simone Rossi committed
931
{
Simone Rossi's avatar
Simone Rossi committed
932
    M_activationModelPtr -> solveModel ( dt);
Simone Rossi's avatar
Simone Rossi committed
933
934
}

Simone Rossi's avatar
Simone Rossi committed
935
936
937

} // namespace LifeV

Simone Rossi's avatar
Simone Rossi committed
938

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