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

thomaskummer's avatar
thomaskummer committed
165
    void setupElectroExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
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
    
thomaskummer's avatar
thomaskummer committed
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
    
thomaskummer's avatar
thomaskummer committed
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

thomaskummer's avatar
thomaskummer committed
195
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
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
                         std::string vonMisesStressFileName  = "VonMisesStress");
Simone Rossi's avatar
Simone Rossi committed
211

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

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

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

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

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

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

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

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

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

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

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


278
279
280
281
282
283
284
285
286
287
288
289
290
291
    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
292
        
293
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
294
295
296
297
298
299
300
301
302
        
//        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) );
//        }

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

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

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


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

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

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

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

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

Simone Rossi's avatar
Simone Rossi committed
373

Simone Deparis's avatar
Simone Deparis committed
374
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
375

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

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

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


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

391
392
393
394
395
    void solveMechanicsLin()
    {
        M_EMStructuralOperatorPtr -> solveLin ();
    }

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

Simone Rossi's avatar
Simone Rossi committed
398

Simone Rossi's avatar
Simone Rossi committed
399
    void solveActivation (Real dt);
thomaskummer's avatar
thomaskummer committed
400
    
thomaskummer's avatar
thomaskummer committed
401
    void computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
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

thomaskummer's avatar
thomaskummer committed
447
    std::vector<meshPtr_Type> mesh()
thomaskummer's avatar
thomaskummer committed
448
    {
thomaskummer's avatar
thomaskummer committed
449
        std::vector<meshPtr_Type> meshVector;
thomaskummer's avatar
thomaskummer committed
450
451
        meshVector.push_back(M_fullMeshPtr);
        meshVector.push_back(M_localMeshPtr);
thomaskummer's avatar
thomaskummer committed
452
        return meshVector;
thomaskummer's avatar
thomaskummer committed
453
    }
454
    
Simone Rossi's avatar
Simone Rossi committed
455
protected:
Simone Rossi's avatar
Simone Rossi committed
456
public:
Simone Deparis's avatar
Simone Deparis committed
457
458
459
460
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
461

Simone Deparis's avatar
Simone Deparis committed
462
463
464
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
465
//    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
466
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
thomaskummer's avatar
thomaskummer committed
467
468
//    exporterPtr_Type                     M_vonMisesStressExporterPtrP;
//    exporterPtr_Type                     M_vonMisesStressExporterPtrA;
thomaskummer's avatar
thomaskummer committed
469
    
Simone Deparis's avatar
Simone Deparis committed
470
471
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
472
    
thomaskummer's avatar
thomaskummer committed
473
    //vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
474
475

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
476
    
thomaskummer's avatar
thomaskummer committed
477
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
thomaskummer's avatar
thomaskummer committed
478
479
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
thomaskummer's avatar
thomaskummer committed
480

Simone Rossi's avatar
Simone Rossi committed
481

Simone Deparis's avatar
Simone Deparis committed
482
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
483

Simone Rossi's avatar
Simone Rossi committed
484
485
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
486

Simone Rossi's avatar
Simone Rossi committed
487
};
Simone Rossi's avatar
Simone Rossi committed
488
489
490

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


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

{
}


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

Simone Rossi's avatar
Simone Rossi committed
541
542
543
544


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
545
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
546
void
Thomas Kummer's avatar
Thomas Kummer committed
547
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
548
{
Simone Rossi's avatar
Simone Rossi committed
549
    M_data.setup (dataFile);
thomaskummer's avatar
thomaskummer committed
550
    
Simone Deparis's avatar
Simone Deparis committed
551
552
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
553
554
        std::cout << "\nEMSolver: setup ... ";
        std::cout << "endtime = " << M_data.activationParameter<Real>("endtime");
Simone Deparis's avatar
Simone Deparis committed
555
    }
thomaskummer's avatar
thomaskummer committed
556
557
    
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
558
559
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
560
561
}

Simone Rossi's avatar
Simone Rossi committed
562
563


Simone Rossi's avatar
Simone Rossi committed
564
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
565
void
Thomas Kummer's avatar
Thomas Kummer committed
566
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
567
{
Simone Deparis's avatar
Simone Deparis committed
568
569
570
571
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
572
573
574
575
576
	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
577
578
579
580
581
582
583
584

    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
585
    M_electroSolverPtr->setParameters();
thomaskummer's avatar
thomaskummer committed
586
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Simone Rossi's avatar
Simone Rossi committed
587
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
thomaskummer's avatar
thomaskummer committed
588
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
589
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
590

Thomas Kummer's avatar
Thomas Kummer committed
591
592
593
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
594
595
596
597
598
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
599
//    }
Simone Deparis's avatar
Simone Deparis committed
600
601
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
602
        std::cout << "... done\n";
Simone Deparis's avatar
Simone Deparis committed
603
    }
Simone Rossi's avatar
Simone Rossi committed
604

Simone Rossi's avatar
Simone Rossi committed
605
606
607
}


Simone Rossi's avatar
Simone Rossi committed
608
609
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
610
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
611
void
Simone Rossi's avatar
Simone Rossi committed
612
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
613
{
Simone Deparis's avatar
Simone Deparis committed
614
615
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
616
        std::cout << "\nEMSolver: setupMechanicalSolver ... ";
Simone Deparis's avatar
Simone Deparis committed
617
    }
thomaskummer's avatar
thomaskummer committed
618
    
Simone Rossi's avatar
Simone Rossi committed
619
620
621
622
623
624
    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
625
626
627
628
629
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

thomaskummer's avatar
thomaskummer committed
630
    dFESpace->setQuadRule(quadRuleTetra4pt);
thomaskummer's avatar
thomaskummer committed
631
632
    //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
633
    
Simone Deparis's avatar
Simone Deparis committed
634
635
636
637
638
639
640
    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
641
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
thomaskummer's avatar
thomaskummer committed
642
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data, dataFile);
Simone Rossi's avatar
Simone Rossi committed
643

thomaskummer's avatar
thomaskummer committed
644
    M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
thomaskummer's avatar
thomaskummer committed
645
646
//    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
647
648
649
650
651
652
    
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "done";
    }

Simone Rossi's avatar
Simone Rossi committed
653
654
655
656
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
657
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
658
void
Simone Rossi's avatar
Simone Rossi committed
659
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
660
661
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
662
{
Simone Deparis's avatar
Simone Deparis committed
663
664
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
665
        std::cout << "\nEMSolver: setupMechanicalBC ... ";
Simone Deparis's avatar
Simone Deparis committed
666
    }
thomaskummer's avatar
thomaskummer committed
667

Simone Deparis's avatar
Simone Deparis committed
668
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
669
670
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
671
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
thomaskummer's avatar
thomaskummer committed
672
673
674
675
676
    
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "done";
    }
Simone Rossi's avatar
Simone Rossi committed
677
678
}

679
680
681
682
683
684
685
    
/////////////////////
//restart hdf5
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::importHdf5 ()
{
thomaskummer's avatar
thomaskummer committed
686
687
    M_electroExporterPtr -> import (30);
    M_activationExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
688
    //M_activationTimeExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
689
    M_vonMisesStressExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
690
691
    //M_vonMisesStressExporterPtrP -> importHdf5 (30);
    //M_vonMisesStressExporterPtrA -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
692
    M_mechanicsExporterPtr -> importHdf5 (30);
693
694
}
                                                   
Simone Rossi's avatar
Simone Rossi committed
695
696
697

/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
698
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
699
void
thomaskummer's avatar
thomaskummer committed
700
701
702
EMSolver<Mesh, ElectroSolver>::setupExporters ( std::string problemFolder,
                                                std::string electroFileName,
                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
703
                                                //std::string activationTimeFileName,
thomaskummer's avatar
thomaskummer committed
704
                                                std::string mechanicsFileName,
thomaskummer's avatar
thomaskummer committed
705
                                                std::string vonMisesStressFileName)
Simone Rossi's avatar
Simone Rossi committed
706
{
Simone Deparis's avatar
Simone Deparis committed
707
708
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
709
        std::cout << "\nEMSolver: setupExporters ... ";
Simone Deparis's avatar
Simone Deparis committed
710
    }
thomaskummer's avatar
thomaskummer committed
711
712
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
713
714
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
715
716
717
718
719
720
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
721
    // Activation
Simone Deparis's avatar
Simone Deparis committed
722
723
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
thomaskummer's avatar
thomaskummer committed
724
725
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
thomaskummer's avatar
thomaskummer committed
726
                                             M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
727
728
                                             M_activationModelPtr -> fiberActivationPtr(),
                                             UInt (0) );
Simone Deparis's avatar
Simone Deparis committed
729

thomaskummer's avatar
thomaskummer committed
730
731
732
733
734
735
736
737
738
739
740
741
//    // 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
742
    
thomaskummer's avatar
thomaskummer committed
743
744
745
746
747
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
748
                                                "Von Mises Stress",
thomaskummer's avatar
thomaskummer committed
749
                                                M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
750
751
752
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
//    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
770

thomaskummer's avatar
thomaskummer committed
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
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
//    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
836
    
thomaskummer's avatar
thomaskummer committed
837
    // Mechanics
thomaskummer's avatar
thomaskummer committed
838
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
839
840
841
842
843
844
845
846
847
848
849
    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
850
851
852
853
854
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
thomaskummer's avatar
thomaskummer committed
855
856
857
858
859
    
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "done";
    }
Simone Rossi's avatar
Simone Rossi committed
860
861
}

thomaskummer's avatar
restart    
thomaskummer committed
862
863
864
865
866
867
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
868
    //M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
869
    M_vonMisesStressExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
870
871
    //M_vonMisesStressExporterPtrP -> setTimeIndex (time);
    //M_vonMisesStressExporterPtrA -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
872
873
874
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
875
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
876
void
thomaskummer's avatar
thomaskummer committed
877
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
878
{
thomaskummer's avatar
thomaskummer committed
879
880
881
882
    M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
    
    M_vonMisesStressExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
883
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
884
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
885
    //M_activationTimeExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
886
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
887
888
}

Simone Rossi's avatar
Simone Rossi committed
889
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
890
void
Simone Rossi's avatar
Simone Rossi committed
891
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
892
{
Simone Deparis's avatar
Simone Deparis committed
893
    M_electroExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
894
895
    M_activationExporterPtr -> closeFile();
    //M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
896
    M_mechanicsExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
897
    M_vonMisesStressExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
898
899
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>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
905
{
Simone Deparis's avatar
Simone Deparis committed
906
907
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
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
}

Simone Rossi's avatar
Simone Rossi committed
914
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
915
void
Simone Rossi's avatar
Simone Rossi committed
916
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
917
{
Simone Deparis's avatar
Simone Deparis committed
918
919
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
920
921
922
923
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
924
925
926
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
927
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
928
void
Simone Rossi's avatar
Simone Rossi committed
929
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
930
{
thomaskummer's avatar
thomaskummer committed
931
932
933
934
935
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "\nEMSolver: solveElectrophysiology ... ";
    }
    
Simone Deparis's avatar
Simone Deparis committed
936
    setAppliedCurrent ( stimulus, time );
thomaskummer's avatar
thomaskummer committed
937

Simone Rossi's avatar
Simone Rossi committed
938
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
thomaskummer's avatar
thomaskummer committed
939
    //M_electroSolverPtr -> solveOneStepGatingVariablesRL();
thomaskummer's avatar
thomaskummer committed
940
    
Simone Rossi's avatar
Simone Rossi committed
941
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
942
    
thomaskummer's avatar
thomaskummer committed
943
    //M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
thomaskummer's avatar
thomaskummer committed
944
945
946
947
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "done";
    }
Simone Rossi's avatar
Simone Rossi committed
948
949
}

Simone Rossi's avatar
Simone Rossi committed
950
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
951
void
Simone Rossi's avatar
Simone Rossi committed
952
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
thomaskummer's avatar
thomaskummer committed
953
    
thomaskummer's avatar
thomaskummer committed
954
{
thomaskummer's avatar
thomaskummer committed
955
956
957
958
959
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "\nEMSolver: solveActivation ... ";
    }
    
thomaskummer's avatar
thomaskummer committed
960
    computeI4f (M_activationModelPtr->I4f(), *M_EMStructuralOperatorPtr->EMMaterial()->fiberVectorPtr(), *M_EMStructuralOperatorPtr->displacementPtr(), M_EMStructuralOperatorPtr->dispFESpacePtr());
thomaskummer's avatar
thomaskummer committed
961

thomaskummer's avatar
thomaskummer committed
962
    M_activationModelPtr -> solveModelPathology ( dt, M_fullMeshPtr, M_EMStructuralOperatorPtr -> dispFESpacePtr() );
thomaskummer's avatar
thomaskummer committed
963
964
965
966
967
    
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "done";
    }
Simone Rossi's avatar
Simone Rossi committed
968
969
}

thomaskummer's avatar
thomaskummer committed
970
971
template<typename Mesh , typename ElectroSolver>
void
thomaskummer's avatar
thomaskummer committed
972
EMSolver<Mesh, ElectroSolver>::computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr)
thomaskummer's avatar
thomaskummer committed
973
974
975
976
{
    VectorEpetra dUdx (disp);
    VectorEpetra dUdy (disp);
    VectorEpetra dUdz (disp);
thomaskummer's avatar
thomaskummer committed
977
    
thomaskummer's avatar
thomaskummer committed
978
979
980
    dUdx = GradientRecovery::ZZGradient (feSpacePtr, disp, 0);
    dUdy = GradientRecovery::ZZGradient (feSpacePtr, disp, 1);
    dUdz = GradientRecovery::ZZGradient (feSpacePtr, disp, 2);
thomaskummer's avatar
thomaskummer committed
981
    
thomaskummer's avatar
thomaskummer committed
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
    int n = i4f.epetraVector().MyLength();
    int i (0); int j (0); int k (0);
    MatrixSmall<3,3> F; VectorSmall<3> f0;
    
    for (int p (0); p < n; p++)
    {
        i = dUdx.blockMap().GID (p);
        j = dUdx.blockMap().GID (p + n);
        k = dUdx.blockMap().GID (p + 2 * n);
        
        F(0,0) = 1.0 + dUdx[i];
        F(0,1) =       dUdy[i];
        F(0,2) =       dUdz[i];
        F(1,0) =       dUdx[j];
        F(1,1) = 1.0 + dUdy[j];
        F(1,2) =       dUdz[j];
        F(2,0) =       dUdx[k];
        F(2,1) =       dUdy[k];
        F(2,2) = 1.0 + dUdz[k];