EMSolver.hpp 38.9 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
152
153
154
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "\nEMSolver: loadMesh ...";
        }
        
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
169
170
        
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << " done";
        }
Simone Deparis's avatar
Simone Deparis committed
171
    }
Simone Rossi's avatar
Simone Rossi committed
172

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

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

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

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

211
212
    void importHdf5 ();

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

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

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

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

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

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

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

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

Thomas Kummer's avatar
Thomas Kummer committed
257
258
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
259
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
260
    {
Simone Rossi's avatar
Simone Rossi committed
261
262
    	//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
263
264
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
265

Simone Rossi's avatar
Simone Rossi committed
266
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
267
268
269
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
270

Simone Rossi's avatar
Simone Rossi committed
271
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
272
273
274
275
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
276

Simone Rossi's avatar
Simone Rossi committed
277
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
278
279
280
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
281

Simone Rossi's avatar
Simone Rossi committed
282
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
283
284
285
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
286

Simone Rossi's avatar
Simone Rossi committed
287
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
288
289
290
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
291
292


293
294
295
296
297
298
299
300
301
302
303
304
305
306
    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
307
        
308
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
309
310
311
312
313
314
315
316
317
        
//        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) );
//        }

318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
    }

    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
338
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
339
340
341
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
342

Simone Rossi's avatar
Simone Rossi committed
343
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
344
345
346
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
347
348


Simone Rossi's avatar
Simone Rossi committed
349
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
350
351
352
353
354
355
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
356
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
357
    }
Simone Rossi's avatar
Simone Rossi committed
358

Simone Rossi's avatar
Simone Rossi committed
359
360
361
362
363
364
    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
365
366
367
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
368

Simone Rossi's avatar
Simone Rossi committed
369
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
370
371
372
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
373

Simone Rossi's avatar
Simone Rossi committed
374
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
375
376
377
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
378
    
thomaskummer's avatar
thomaskummer committed
379
380
381
382
//    vectorPtr_Type activationTimePtr()
//    {
//        return M_activationTimePtr;
//    }
Simone Rossi's avatar
Simone Rossi committed
383

thomaskummer's avatar
thomaskummer committed
384
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
385
386
387
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
388

Simone Deparis's avatar
Simone Deparis committed
389
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
390

Simone Deparis's avatar
Simone Deparis committed
391
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
392

Simone Deparis's avatar
Simone Deparis committed
393
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
394

Simone Rossi's avatar
Simone Rossi committed
395
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
396
397
398
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
399
400


Simone Rossi's avatar
Simone Rossi committed
401
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
402
403
404
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
405

406
407
408
409
410
    void solveMechanicsLin()
    {
        M_EMStructuralOperatorPtr -> solveLin ();
    }

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

Simone Rossi's avatar
Simone Rossi committed
413

Simone Rossi's avatar
Simone Rossi committed
414
    void solveActivation (Real dt);
thomaskummer's avatar
thomaskummer committed
415
    
thomaskummer's avatar
thomaskummer committed
416
    void computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
Simone Rossi's avatar
Simone Rossi committed
417

Simone Rossi's avatar
Simone Rossi committed
418
419
420
421
422
423
424
425
426
427
428
429
430
    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
431
432
433
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
    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
449

450
    void showMe() const {}
451
452
453
454
455
456
457
458
459
460
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
461

thomaskummer's avatar
thomaskummer committed
462
    std::vector<meshPtr_Type> mesh()
thomaskummer's avatar
thomaskummer committed
463
    {
thomaskummer's avatar
thomaskummer committed
464
        std::vector<meshPtr_Type> meshVector;
thomaskummer's avatar
thomaskummer committed
465
466
        meshVector.push_back(M_fullMeshPtr);
        meshVector.push_back(M_localMeshPtr);
thomaskummer's avatar
thomaskummer committed
467
        return meshVector;
thomaskummer's avatar
thomaskummer committed
468
    }
469
    
Simone Rossi's avatar
Simone Rossi committed
470
protected:
Simone Rossi's avatar
Simone Rossi committed
471
public:
Simone Deparis's avatar
Simone Deparis committed
472
473
474
475
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
476

Simone Deparis's avatar
Simone Deparis committed
477
478
479
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
480
//    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
481
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
thomaskummer's avatar
thomaskummer committed
482
483
//    exporterPtr_Type                     M_vonMisesStressExporterPtrP;
//    exporterPtr_Type                     M_vonMisesStressExporterPtrA;
thomaskummer's avatar
thomaskummer committed
484
    
Simone Deparis's avatar
Simone Deparis committed
485
486
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
487
    
thomaskummer's avatar
thomaskummer committed
488
    //vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
489
490

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
491
    
thomaskummer's avatar
thomaskummer committed
492
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
thomaskummer's avatar
thomaskummer committed
493
494
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
thomaskummer's avatar
thomaskummer committed
495

Simone Rossi's avatar
Simone Rossi committed
496

Simone Deparis's avatar
Simone Deparis committed
497
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
498

Simone Rossi's avatar
Simone Rossi committed
499
500
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
501

Simone Rossi's avatar
Simone Rossi committed
502
};
Simone Rossi's avatar
Simone Rossi committed
503
504
505

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
506
507
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
508
509
510
511
512
513
514
515
516
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
517
    //M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
518
    M_oneWayCoupling     (true),
thomaskummer's avatar
thomaskummer committed
519
    M_wteTotal ( ),
thomaskummer's avatar
thomaskummer committed
520
521
//    M_wtePassive ( ),
//    M_wteActive ( ),
Simone Rossi's avatar
Simone Rossi committed
522
523
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
524
525
526
527
528
529
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
530
531
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
532
533
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
534
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
535
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
536
537
538
539
540
    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
541
    //M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
542
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
thomaskummer's avatar
thomaskummer committed
543
    M_wteTotal                   (solver.M_wteTotal),
thomaskummer's avatar
thomaskummer committed
544
545
//    M_wtePassive                   (solver.M_wtePassive),
//    M_wteActive                   (solver.M_wteActive),
Simone Rossi's avatar
Simone Rossi committed
546
547
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
548
549
550
551
552
553
554

{
}


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

Simone Rossi's avatar
Simone Rossi committed
556
557
558
559


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
560
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
561
void
Thomas Kummer's avatar
Thomas Kummer committed
562
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
563
{
Simone Rossi's avatar
Simone Rossi committed
564
    M_data.setup (dataFile);
thomaskummer's avatar
thomaskummer committed
565
    
Simone Deparis's avatar
Simone Deparis committed
566
567
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
568
        std::cout << "\n\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " ...";
Simone Deparis's avatar
Simone Deparis committed
569
    }
thomaskummer's avatar
thomaskummer committed
570
571
    
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
572
573
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
thomaskummer's avatar
thomaskummer committed
574
575
576
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
577
        std::cout << "\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
578
579
    }

Simone Rossi's avatar
Simone Rossi committed
580
581
}

Simone Rossi's avatar
Simone Rossi committed
582
583


Simone Rossi's avatar
Simone Rossi committed
584
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
585
void
Thomas Kummer's avatar
Thomas Kummer committed
586
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
587
{
Simone Deparis's avatar
Simone Deparis committed
588
589
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
590
        std::cout << "\nEMSolver: setupElectroSolver ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
591
    }
Simone Rossi's avatar
Simone Rossi committed
592
593
594
595
596
	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
597
598
599

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
Simone Rossi's avatar
Simone Rossi committed
600
    M_electroSolverPtr->setParameters();
thomaskummer's avatar
thomaskummer committed
601
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Simone Rossi's avatar
Simone Rossi committed
602
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
thomaskummer's avatar
thomaskummer committed
603
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
604
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
605

Thomas Kummer's avatar
Thomas Kummer committed
606
607
608
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
609
610
611
612
613
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
614
//    }
thomaskummer's avatar
thomaskummer committed
615
    
Simone Deparis's avatar
Simone Deparis committed
616
617
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
618
        std::cout << "\nEMSolver: setupElectroSolver - done" << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
619
    }
Simone Rossi's avatar
Simone Rossi committed
620

Simone Rossi's avatar
Simone Rossi committed
621
622
623
}


Simone Rossi's avatar
Simone Rossi committed
624
625
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
626
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
627
void
Simone Rossi's avatar
Simone Rossi committed
628
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
629
{
Simone Deparis's avatar
Simone Deparis committed
630
631
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
632
        std::cout << "\nEMSolver: setupMechanicalSolver ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
633
    }
thomaskummer's avatar
thomaskummer committed
634
    
Simone Rossi's avatar
Simone Rossi committed
635
636
637
638
639
640
    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
641
642
643
644
645
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

thomaskummer's avatar
thomaskummer committed
646
    dFESpace->setQuadRule(quadRuleTetra4pt);
thomaskummer's avatar
thomaskummer committed
647
648
    //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
649
    
Simone Deparis's avatar
Simone Deparis committed
650
651
652
653
654
655
656
    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
657
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
thomaskummer's avatar
thomaskummer committed
658
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data, dataFile);
Simone Rossi's avatar
Simone Rossi committed
659

thomaskummer's avatar
thomaskummer committed
660
    M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
thomaskummer's avatar
thomaskummer committed
661
662
//    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
663
664
665
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
666
        std::cout << "\nEMSolver: setupMechanicalSolver - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
667
668
    }

Simone Rossi's avatar
Simone Rossi committed
669
670
671
672
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
673
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
674
void
Simone Rossi's avatar
Simone Rossi committed
675
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
676
677
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
678
{
Simone Deparis's avatar
Simone Deparis committed
679
680
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
681
        std::cout << "\nEMSolver: setupMechanicalBC ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
682
    }
thomaskummer's avatar
thomaskummer committed
683

Simone Deparis's avatar
Simone Deparis committed
684
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
685
686
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
687
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
thomaskummer's avatar
thomaskummer committed
688
689
690
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
691
        std::cout << "EMSolver: setupMechanicalBC - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
692
    }
Simone Rossi's avatar
Simone Rossi committed
693
694
}

695
696
697
698
699
700
701
    
/////////////////////
//restart hdf5
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::importHdf5 ()
{
thomaskummer's avatar
thomaskummer committed
702
703
    M_electroExporterPtr -> import (30);
    M_activationExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
704
    //M_activationTimeExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
705
    M_vonMisesStressExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
706
707
    //M_vonMisesStressExporterPtrP -> importHdf5 (30);
    //M_vonMisesStressExporterPtrA -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
708
    M_mechanicsExporterPtr -> importHdf5 (30);
709
710
}
                                                   
Simone Rossi's avatar
Simone Rossi committed
711
712
713

/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
714
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
715
void
thomaskummer's avatar
thomaskummer committed
716
717
718
EMSolver<Mesh, ElectroSolver>::setupExporters ( std::string problemFolder,
                                                std::string electroFileName,
                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
719
                                                //std::string activationTimeFileName,
thomaskummer's avatar
thomaskummer committed
720
                                                std::string mechanicsFileName,
thomaskummer's avatar
thomaskummer committed
721
                                                std::string vonMisesStressFileName)
Simone Rossi's avatar
Simone Rossi committed
722
{
Simone Deparis's avatar
Simone Deparis committed
723
724
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
725
        std::cout << "\nEMSolver: setupExporters ... " << '\r' << std::flush;
Simone Deparis's avatar
Simone Deparis committed
726
    }
thomaskummer's avatar
thomaskummer committed
727
728
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
729
730
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
731
732
733
734
735
736
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
737
    // Activation
Simone Deparis's avatar
Simone Deparis committed
738
739
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
thomaskummer's avatar
thomaskummer committed
740
741
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
thomaskummer's avatar
thomaskummer committed
742
                                             M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
743
744
                                             M_activationModelPtr -> fiberActivationPtr(),
                                             UInt (0) );
Simone Deparis's avatar
Simone Deparis committed
745

thomaskummer's avatar
thomaskummer committed
746
747
748
749
750
751
752
753
754
755
756
757
//    // 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
758
    
thomaskummer's avatar
thomaskummer committed
759
760
761
762
763
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
764
                                                "Von Mises Stress",
thomaskummer's avatar
thomaskummer committed
765
                                                M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
766
767
768
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
//    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
786

thomaskummer's avatar
thomaskummer committed
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
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
//    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
852
    
thomaskummer's avatar
thomaskummer committed
853
    // Mechanics
thomaskummer's avatar
thomaskummer committed
854
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
855
856
857
858
859
860
861
862
863
864
865
    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
866
867
868
869
870
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
thomaskummer's avatar
thomaskummer committed
871
872
873
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
874
        std::cout << "EMSolver: setupExporters - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
875
    }
Simone Rossi's avatar
Simone Rossi committed
876
877
}

thomaskummer's avatar
restart    
thomaskummer committed
878
879
880
881
882
883
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
884
    //M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
885
    M_vonMisesStressExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
886
887
    //M_vonMisesStressExporterPtrP -> setTimeIndex (time);
    //M_vonMisesStressExporterPtrA -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
888
889
890
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
891
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
892
void
thomaskummer's avatar
thomaskummer committed
893
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
894
{
thomaskummer's avatar
thomaskummer committed
895
896
897
898
    M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
    
    M_vonMisesStressExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
899
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
900
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
901
    //M_activationTimeExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
902
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
903
904
}

Simone Rossi's avatar
Simone Rossi committed
905
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
906
void
Simone Rossi's avatar
Simone Rossi committed
907
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
908
{
Simone Deparis's avatar
Simone Deparis committed
909
    M_electroExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
910
911
    M_activationExporterPtr -> closeFile();
    //M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
912
    M_mechanicsExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
913
    M_vonMisesStressExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
914
915
916
917
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
918
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
919
void
Simone Rossi's avatar
Simone Rossi committed
920
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
921
{
Simone Deparis's avatar
Simone Deparis committed
922
923
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
924
925
926
927
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
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>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
933
{
Simone Deparis's avatar
Simone Deparis committed
934
935
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
936
937
938
939
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
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>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
946
{
thomaskummer's avatar
thomaskummer committed
947
948
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
949
        std::cout << "\nEMSolver: solveElectrophysiology ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
950
951
    }
    
Simone Deparis's avatar
Simone Deparis committed
952
    setAppliedCurrent ( stimulus, time );
thomaskummer's avatar
thomaskummer committed
953

Simone Rossi's avatar
Simone Rossi committed
954
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
thomaskummer's avatar
thomaskummer committed
955
    //M_electroSolverPtr -> solveOneStepGatingVariablesRL();
thomaskummer's avatar
thomaskummer committed
956
    
Simone Rossi's avatar
Simone Rossi committed
957
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
958
    
thomaskummer's avatar
thomaskummer committed
959
    //M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
thomaskummer's avatar
thomaskummer committed
960
961
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
962
        std::cout << "EMSolver: solveElectrophysiology - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
963
    }
Simone Rossi's avatar
Simone Rossi committed
964
965
}

Simone Rossi's avatar
Simone Rossi committed
966
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
967
void
Simone Rossi's avatar
Simone Rossi committed
968
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
thomaskummer's avatar
thomaskummer committed
969
    
thomaskummer's avatar
thomaskummer committed
970
{
thomaskummer's avatar
thomaskummer committed
971
972
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
973
        std::cout << "\nEMSolver: solveActivation ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
974
975
    }
    
thomaskummer's avatar
thomaskummer committed
976
    computeI4f (M_activationModelPtr->I4f(), *M_EMStructuralOperatorPtr->EMMaterial()->fiberVectorPtr(), *M_EMStructuralOperatorPtr->displacementPtr(), M_EMStructuralOperatorPtr->dispFESpacePtr());
thomaskummer's avatar
thomaskummer committed
977

thomaskummer's avatar
thomaskummer committed
978
    M_activationModelPtr -> solveModelPathology ( dt, M_fullMeshPtr, M_EMStructuralOperatorPtr -> dispFESpacePtr() );
thomaskummer's avatar
thomaskummer committed
979
980
981
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
982
        std::cout << "EMSolver: solveActivation - done" << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
983
    }
Simone Rossi's avatar
Simone Rossi committed
984
985
}

thomaskummer's avatar
thomaskummer committed
986
987
template<typename Mesh , typename ElectroSolver>
void
thomaskummer's avatar
thomaskummer committed
988
EMSolver<Mesh, ElectroSolver>::computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr)
thomaskummer's avatar
thomaskummer committed
989
990
991
992
{
    VectorEpetra dUdx (disp);
    VectorEpetra dUdy (disp);
    VectorEpetra dUdz (disp);
thomaskummer's avatar
thomaskummer committed
993