EMSolver.hpp 40.4 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
    {
thomaskummer's avatar
thomaskummer committed
150
151
        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
152
            std::cout << "\n\nEMSolver: loadMesh ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
153
        }
thomaskummer's avatar
thomaskummer committed
154

Simone Rossi's avatar
Simone Rossi committed
155
        M_fullMeshPtr.reset( new Mesh() );
Simone Deparis's avatar
Simone Deparis committed
156
        MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
Simone Rossi's avatar
Simone Rossi committed
157
158
159
160
161
162
163
164
165
        if(M_commPtr)
        {
			M_localMeshPtr->setComm(M_commPtr);
			M_fullMeshPtr->setComm(M_commPtr);
        }
        else
        {
        	M_commPtr = M_localMeshPtr -> comm();
        }
thomaskummer's avatar
thomaskummer committed
166
167
168
        
        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
169
            std::cout << "EMSolver: loadMesh - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
170
        }
thomaskummer's avatar
thomaskummer committed
171

Simone Deparis's avatar
Simone Deparis committed
172
    }
Simone Rossi's avatar
Simone Rossi committed
173

thomaskummer's avatar
thomaskummer committed
174
    void setupElectroExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Deparis's avatar
Simone Deparis committed
175
176
177
    {
        M_electroSolverPtr -> setupExporter (*M_electroExporterPtr, outputFileName, problemFolder);
    }
Simone Rossi's avatar
Simone Rossi committed
178

Simone Rossi's avatar
Simone Rossi committed
179
    void setupActivationExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
Simone Rossi's avatar
Simone Rossi committed
180
    {
Simone Deparis's avatar
Simone Deparis committed
181
        EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
Simone Rossi's avatar
Simone Rossi committed
182
    }
thomaskummer's avatar
thomaskummer committed
183
    
thomaskummer's avatar
thomaskummer committed
184
185
186
187
//    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
188

thomaskummer's avatar
thomaskummer committed
189
190
191
192
    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
193
    
thomaskummer's avatar
thomaskummer committed
194
195
196
197
198
199
200
201
202
//    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
203

thomaskummer's avatar
thomaskummer committed
204
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
Simone Rossi's avatar
Simone Rossi committed
205
    {
Simone Deparis's avatar
Simone Deparis committed
206
207
208
209
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
210
211
    }

212
213
    void importHdf5 ();

Simone Deparis's avatar
Simone Deparis committed
214
215
216
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
thomaskummer's avatar
thomaskummer committed
217
                         //std::string activationTimeFileName  = "ActivationTimeSolution",
thomaskummer's avatar
thomaskummer committed
218
                         std::string mechanicsFileName  = "MechanicalSolution",
thomaskummer's avatar
thomaskummer committed
219
                         std::string vonMisesStressFileName  = "VonMisesStress");
Simone Rossi's avatar
Simone Rossi committed
220

Thomas Kummer's avatar
Thomas Kummer committed
221
222
223
224
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
225

Thomas Kummer's avatar
Thomas Kummer committed
226
227
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
228
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
229

Simone Deparis's avatar
Simone Deparis committed
230
231
232
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
233

Simone Rossi's avatar
Simone Rossi committed
234
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
235
236
237
    {
        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
238
            std::cout << "\nEMSolver: setupActivation ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
239
        }
thomaskummer's avatar
thomaskummer committed
240
        
Simone Rossi's avatar
Simone Rossi committed
241
242
243
244
        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() );
thomaskummer's avatar
thomaskummer committed
245
246
247

        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
248
            std::cout << "EMSolver: setupActivation - done " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
249
250
        }

Simone Deparis's avatar
Simone Deparis committed
251
    }
Simone Rossi's avatar
Simone Rossi committed
252

Thomas Kummer's avatar
Thomas Kummer committed
253
254
255
256
    void setup (GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setup (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
257

Thomas Kummer's avatar
Thomas Kummer committed
258
259
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
260
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
261
    {
thomaskummer's avatar
thomaskummer committed
262
263
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: buildMechanicalSystem ... " << '\r' << std::flush;
        
Simone Rossi's avatar
Simone Rossi committed
264
265
    	//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
266
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
thomaskummer's avatar
thomaskummer committed
267
268
        
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: buildMechanicalSystem - done" << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
269
    }
Simone Rossi's avatar
Simone Rossi committed
270

Simone Rossi's avatar
Simone Rossi committed
271
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
272
    {
thomaskummer's avatar
thomaskummer committed
273
274
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: buildElectroSystem ... " << '\r' << std::flush;

Simone Deparis's avatar
Simone Deparis committed
275
        M_electroSolverPtr -> setupMatrices();
thomaskummer's avatar
thomaskummer committed
276
277
        
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: buildElectroSystem - done" << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
278
    }
Simone Rossi's avatar
Simone Rossi committed
279

Simone Rossi's avatar
Simone Rossi committed
280
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
281
282
283
284
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
285

Simone Rossi's avatar
Simone Rossi committed
286
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
287
    {
thomaskummer's avatar
thomaskummer committed
288
289
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: initializeElectroVariables ... " << '\r' << std::flush;

Simone Deparis's avatar
Simone Deparis committed
290
        M_electroSolverPtr -> setInitialConditions();
thomaskummer's avatar
thomaskummer committed
291
292
293
        
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: initializeElectroVariables - done" << '\r' << std::flush;

Simone Deparis's avatar
Simone Deparis committed
294
    }
Simone Rossi's avatar
Simone Rossi committed
295

Simone Rossi's avatar
Simone Rossi committed
296
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
297
298
299
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
300

Simone Rossi's avatar
Simone Rossi committed
301
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
302
303
304
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
305
306


307
308
309
310
311
    void setupFiberVector( const std::string& fileName,
						   const std::string& fieldName,
						   const std::string& postDir = "./",
						   const std::string& polynomialDegree = "P1"  )
    {
thomaskummer's avatar
thomaskummer committed
312
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: setupFiberVector ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
313

314
315
    	setupMechanicalFiberVector(fileName, fieldName, postDir, polynomialDegree);
    	M_electroSolverPtr->setFiberPtr(getMechanicsFibers());
thomaskummer's avatar
thomaskummer committed
316
        
thomaskummer's avatar
thomaskummer committed
317
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: setupFiberVector - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
318

319
320
321
322
323
324
325
    }

    void setupMechanicalFiberVector( const std::string& fileName,
    		                         const std::string& fieldName,
									 const std::string& postDir = "./",
                                     const std::string& polynomialDegree = "P1"  )
    {
thomaskummer's avatar
thomaskummer committed
326
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: setupMechanicalFiberVector ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
327
        
328
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
329
330
331
332
333
334
335
336
        
//        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) );
//        }
thomaskummer's avatar
thomaskummer committed
337
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: setupMechanicalFiberVector - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
338

339
340
341
342
343
344
345
    }

    void setupMechanicalSheetVector( const std::string& fileName,
    		                         const std::string& fieldName,
									 const std::string& postDir = "./",
                                     const std::string& polynomialDegree = "P1"  )
    {
thomaskummer's avatar
thomaskummer committed
346
        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: setupMechanicalSheetVector ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
347

348
        ElectrophysiologyUtility::importVectorField (getMechanicsSheets(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
349
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: setupMechanicalSheetVector - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
350

351
352
353
354
355
356
357
358
359
360
361
362
    }

    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
363
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
364
365
366
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
367

Simone Rossi's avatar
Simone Rossi committed
368
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
369
370
371
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
372
373


Simone Rossi's avatar
Simone Rossi committed
374
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
375
376
377
378
379
380
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
381
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
382
    }
Simone Rossi's avatar
Simone Rossi committed
383

Simone Rossi's avatar
Simone Rossi committed
384
385
386
387
388
389
    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
390
391
392
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
393

Simone Rossi's avatar
Simone Rossi committed
394
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
395
396
397
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
398

Simone Rossi's avatar
Simone Rossi committed
399
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
400
401
402
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
403
    
thomaskummer's avatar
thomaskummer committed
404
405
406
407
//    vectorPtr_Type activationTimePtr()
//    {
//        return M_activationTimePtr;
//    }
Simone Rossi's avatar
Simone Rossi committed
408

thomaskummer's avatar
thomaskummer committed
409
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
410
411
412
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
413

Simone Deparis's avatar
Simone Deparis committed
414
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
415

Simone Deparis's avatar
Simone Deparis committed
416
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
417

Simone Deparis's avatar
Simone Deparis committed
418
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
419

Simone Rossi's avatar
Simone Rossi committed
420
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
421
422
423
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
424
425


Simone Rossi's avatar
Simone Rossi committed
426
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
427
428
429
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
430

431
432
433
434
435
    void solveMechanicsLin()
    {
        M_EMStructuralOperatorPtr -> solveLin ();
    }

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

Simone Rossi's avatar
Simone Rossi committed
438

Simone Rossi's avatar
Simone Rossi committed
439
    void solveActivation (Real dt);
thomaskummer's avatar
thomaskummer committed
440
    
thomaskummer's avatar
thomaskummer committed
441
    void computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
Simone Rossi's avatar
Simone Rossi committed
442

Simone Rossi's avatar
Simone Rossi committed
443
444
445
446
447
448
449
450
451
452
453
454
455
    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
456
457
458
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
    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
474

475
    void showMe() const {}
476
477
478
479
480
481
482
483
484
485
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
486

thomaskummer's avatar
thomaskummer committed
487
    std::vector<meshPtr_Type> mesh()
thomaskummer's avatar
thomaskummer committed
488
    {
thomaskummer's avatar
thomaskummer committed
489
        std::vector<meshPtr_Type> meshVector;
thomaskummer's avatar
thomaskummer committed
490
491
        meshVector.push_back(M_fullMeshPtr);
        meshVector.push_back(M_localMeshPtr);
thomaskummer's avatar
thomaskummer committed
492
        return meshVector;
thomaskummer's avatar
thomaskummer committed
493
    }
494
    
Simone Rossi's avatar
Simone Rossi committed
495
protected:
Simone Rossi's avatar
Simone Rossi committed
496
public:
Simone Deparis's avatar
Simone Deparis committed
497
498
499
500
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
501

Simone Deparis's avatar
Simone Deparis committed
502
503
504
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
505
//    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
506
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
thomaskummer's avatar
thomaskummer committed
507
508
//    exporterPtr_Type                     M_vonMisesStressExporterPtrP;
//    exporterPtr_Type                     M_vonMisesStressExporterPtrA;
thomaskummer's avatar
thomaskummer committed
509
    
Simone Deparis's avatar
Simone Deparis committed
510
511
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
512
    
thomaskummer's avatar
thomaskummer committed
513
    //vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
514
515

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
516
    
thomaskummer's avatar
thomaskummer committed
517
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
thomaskummer's avatar
thomaskummer committed
518
519
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
thomaskummer's avatar
thomaskummer committed
520

Simone Rossi's avatar
Simone Rossi committed
521

Simone Deparis's avatar
Simone Deparis committed
522
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
523

Simone Rossi's avatar
Simone Rossi committed
524
525
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
526

Simone Rossi's avatar
Simone Rossi committed
527
};
Simone Rossi's avatar
Simone Rossi committed
528
529
530

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
531
532
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
533
534
535
536
537
538
539
540
541
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
542
    //M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
543
    M_oneWayCoupling     (true),
thomaskummer's avatar
thomaskummer committed
544
    M_wteTotal ( ),
thomaskummer's avatar
thomaskummer committed
545
546
//    M_wtePassive ( ),
//    M_wteActive ( ),
Simone Rossi's avatar
Simone Rossi committed
547
548
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
549
550
551
552
553
554
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
555
556
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
557
558
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
559
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
560
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
561
562
563
564
565
    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
566
    //M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
567
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
thomaskummer's avatar
thomaskummer committed
568
    M_wteTotal                   (solver.M_wteTotal),
thomaskummer's avatar
thomaskummer committed
569
570
//    M_wtePassive                   (solver.M_wtePassive),
//    M_wteActive                   (solver.M_wteActive),
Simone Rossi's avatar
Simone Rossi committed
571
572
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
573
574
575
576
577
578
579

{
}


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

Simone Rossi's avatar
Simone Rossi committed
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
Thomas Kummer's avatar
Thomas Kummer committed
587
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
588
{
Simone Rossi's avatar
Simone Rossi committed
589
    M_data.setup (dataFile);
thomaskummer's avatar
thomaskummer committed
590
    
Simone Deparis's avatar
Simone Deparis committed
591
592
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
593
        std::cout << "\n\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " ...";
Simone Deparis's avatar
Simone Deparis committed
594
    }
thomaskummer's avatar
thomaskummer committed
595
596
    
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
597
598
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
thomaskummer's avatar
thomaskummer committed
599
600
601
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
602
        std::cout << "\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
603
604
    }

Simone Rossi's avatar
Simone Rossi committed
605
606
}

Simone Rossi's avatar
Simone Rossi committed
607
608


Simone Rossi's avatar
Simone Rossi committed
609
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
610
void
Thomas Kummer's avatar
Thomas Kummer committed
611
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
612
{
Simone Deparis's avatar
Simone Deparis committed
613
614
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
615
        std::cout << "\nEMSolver: setupElectroSolver ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
616
    }
Simone Rossi's avatar
Simone Rossi committed
617
618
619
620
621
	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
622
623
624

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
Simone Rossi's avatar
Simone Rossi committed
625
    M_electroSolverPtr->setParameters();
thomaskummer's avatar
thomaskummer committed
626
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Simone Rossi's avatar
Simone Rossi committed
627
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
thomaskummer's avatar
thomaskummer committed
628
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
629
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
630

Thomas Kummer's avatar
Thomas Kummer committed
631
632
633
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
634
635
636
637
638
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
639
//    }
thomaskummer's avatar
thomaskummer committed
640
    
Simone Deparis's avatar
Simone Deparis committed
641
642
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
643
        std::cout << "\nEMSolver: setupElectroSolver - done" << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
644
    }
Simone Rossi's avatar
Simone Rossi committed
645

Simone Rossi's avatar
Simone Rossi committed
646
647
648
}


Simone Rossi's avatar
Simone Rossi committed
649
650
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
651
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
652
void
Simone Rossi's avatar
Simone Rossi committed
653
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
654
{
Simone Deparis's avatar
Simone Deparis committed
655
656
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
657
        std::cout << "\nEMSolver: setupMechanicalSolver ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
658
    }
thomaskummer's avatar
thomaskummer committed
659
    
Simone Rossi's avatar
Simone Rossi committed
660
661
662
663
664
665
    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
666
667
668
669
670
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

thomaskummer's avatar
thomaskummer committed
671
    dFESpace->setQuadRule(quadRuleTetra4pt);
thomaskummer's avatar
thomaskummer committed
672
673
    //if (M_commPtr -> MyPID() == 0) std::cout << "\nPolynomial degree: " << dFESpace->polynomialDegree();
    //if (M_commPtr -> MyPID() == 0) std::cout << "\nDof: " << dFESpace->dof() << std::endl;
thomaskummer's avatar
thomaskummer committed
674
    
Simone Deparis's avatar
Simone Deparis committed
675
676
677
678
679
680
681
    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
682
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
thomaskummer's avatar
thomaskummer committed
683
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data, dataFile);
Simone Rossi's avatar
Simone Rossi committed
684

thomaskummer's avatar
thomaskummer committed
685
    M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
thomaskummer's avatar
thomaskummer committed
686
687
//    M_wtePassive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "passive");
//    M_wteActive.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, "active");
thomaskummer's avatar
thomaskummer committed
688
689
690
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
691
        std::cout << "\nEMSolver: setupMechanicalSolver - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
692
693
    }

Simone Rossi's avatar
Simone Rossi committed
694
695
696
697
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
698
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
699
void
Simone Rossi's avatar
Simone Rossi committed
700
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
701
702
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
703
{
Simone Deparis's avatar
Simone Deparis committed
704
705
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
706
        std::cout << "\nEMSolver: setupMechanicalBC ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
707
    }
thomaskummer's avatar
thomaskummer committed
708

Simone Deparis's avatar
Simone Deparis committed
709
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
710
711
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
712
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
thomaskummer's avatar
thomaskummer committed
713
714
715
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
716
        std::cout << "EMSolver: setupMechanicalBC - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
717
    }
Simone Rossi's avatar
Simone Rossi committed
718
719
}

720
721
722
723
724
725
726
    
/////////////////////
//restart hdf5
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::importHdf5 ()
{
thomaskummer's avatar
thomaskummer committed
727
728
    M_electroExporterPtr -> import (30);
    M_activationExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
729
    //M_activationTimeExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
730
    M_vonMisesStressExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
731
732
    //M_vonMisesStressExporterPtrP -> importHdf5 (30);
    //M_vonMisesStressExporterPtrA -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
733
    M_mechanicsExporterPtr -> importHdf5 (30);
734
735
}
                                                   
Simone Rossi's avatar
Simone Rossi committed
736
737
738

/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
739
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
740
void
thomaskummer's avatar
thomaskummer committed
741
742
743
EMSolver<Mesh, ElectroSolver>::setupExporters ( std::string problemFolder,
                                                std::string electroFileName,
                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
744
                                                //std::string activationTimeFileName,
thomaskummer's avatar
thomaskummer committed
745
                                                std::string mechanicsFileName,
thomaskummer's avatar
thomaskummer committed
746
                                                std::string vonMisesStressFileName)
Simone Rossi's avatar
Simone Rossi committed
747
{
Simone Deparis's avatar
Simone Deparis committed
748
749
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
750
        std::cout << "\nEMSolver: setupExporters ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
751
    }
thomaskummer's avatar
thomaskummer committed
752
753
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
754
755
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
756
757
758
759
760
761
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
762
    // Activation
Simone Deparis's avatar
Simone Deparis committed
763
764
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
thomaskummer's avatar
thomaskummer committed
765
766
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
thomaskummer's avatar
thomaskummer committed
767
                                             M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
768
769
                                             M_activationModelPtr -> fiberActivationPtr(),
                                             UInt (0) );
Simone Deparis's avatar
Simone Deparis committed
770

thomaskummer's avatar
thomaskummer committed
771
772
773
774
775
776
777
778
779
780
781
782
//    // 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;
//
//    M_activationTimeExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
//                                            "Activation Time",
//                                            M_electroSolverPtr -> feSpacePtr(),
//                                            M_activationTimePtr,
//                                            UInt (0) );
thomaskummer's avatar
thomaskummer committed
783
    
thomaskummer's avatar
thomaskummer committed
784
785
786
787
788
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
789
                                                "Von Mises Stress",
thomaskummer's avatar
thomaskummer committed
790
                                                M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
791
792
793
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
//    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
811

thomaskummer's avatar
thomaskummer committed
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
//    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
877
    
thomaskummer's avatar
thomaskummer committed
878
    // Mechanics
thomaskummer's avatar
thomaskummer committed
879
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
880
881
882
883
884
885
886
887
888
889
890
    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
891
892
893
894
895
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
thomaskummer's avatar
thomaskummer committed
896
897
898
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
899
        std::cout << "EMSolver: setupExporters - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
900
    }
Simone Rossi's avatar
Simone Rossi committed
901
902
}

thomaskummer's avatar
restart    
thomaskummer committed
903
904
905
906
907
908
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
909
    //M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
910
    M_vonMisesStressExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
911
912
    //M_vonMisesStressExporterPtrP -> setTimeIndex (time);
    //M_vonMisesStressExporterPtrA -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
913
914
915
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
916
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
917
void
thomaskummer's avatar
thomaskummer committed
918
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
919
{
thomaskummer's avatar
thomaskummer committed
920
921
922
923
    M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
    
    M_vonMisesStressExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
924
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
925
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
926
    //M_activationTimeExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
927
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
928
929
}

Simone Rossi's avatar
Simone Rossi committed
930
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
931
void
Simone Rossi's avatar
Simone Rossi committed
932
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
933
{
Simone Deparis's avatar
Simone Deparis committed
934
    M_electroExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
935
936
    M_activationExporterPtr -> closeFile();
    //M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
937
    M_mechanicsExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
938
    M_vonMisesStressExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
939
940
941
942
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
943
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
944
void
Simone Rossi's avatar
Simone Rossi committed
945
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
946
{
Simone Deparis's avatar
Simone Deparis committed
947
948
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
949
950
951
952
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
953
954
}

Simone Rossi's avatar
Simone Rossi committed
955
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
956
void
Simone Rossi's avatar
Simone Rossi committed
957
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
958
{
Simone Deparis's avatar
Simone Deparis committed
959
960
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
961
962
963
964
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
965
966
967
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
968
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
969
void
Simone Rossi's avatar
Simone Rossi committed
970
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
971
{
thomaskummer's avatar
thomaskummer committed
972
973
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
974
        std::cout << "\nEMSolver: solveElectrophysiology ... " << '\r' << std::flush;