EMSolver.hpp 39.8 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
#include <lifev/structure/solver/WallTensionEstimator.hpp>
Simone Rossi's avatar
Simone Rossi committed
55
#include <lifev/em/solver/activation/ActivationModelsList.hpp>
Simone Rossi's avatar
Simone Rossi committed
56
57

#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
thomaskummer's avatar
thomaskummer committed
58

Simone Rossi's avatar
Simone Rossi committed
59
60
61
62
63
64

namespace LifeV
{

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

Simone Rossi's avatar
Simone Rossi committed
65
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
66
67
class EMSolver
{
Simone Rossi's avatar
Simone Rossi committed
68
public:
Simone Deparis's avatar
Simone Deparis committed
69
    typedef Mesh                                              mesh_Type;
Simone Rossi's avatar
Simone Rossi committed
70

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

Simone Deparis's avatar
Simone Deparis committed
73
    typedef Epetra_Comm                                       comm_Type;
Simone Rossi's avatar
Simone Rossi committed
74

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

Simone Deparis's avatar
Simone Deparis committed
77
    typedef VectorEpetra                                      vector_Type;
thomaskummer's avatar
thomaskummer committed
78
    
Simone Deparis's avatar
Simone Deparis committed
79
    typedef StructuralConstitutiveLawData                     structureData_Type;
Simone Rossi's avatar
Simone Rossi committed
80

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

Simone Deparis's avatar
Simone Deparis committed
83
    typedef EMStructuralOperator< mesh_Type >                 structuralOperator_Type;
Simone Rossi's avatar
Simone Rossi committed
84

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

Simone Deparis's avatar
Simone Deparis committed
87
    typedef BCHandler                                          bc_Type;
Simone Rossi's avatar
Simone Rossi committed
88

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

Simone Deparis's avatar
Simone Deparis committed
91
    typedef StructuralOperator< mesh_Type >                   physicalSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
92

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

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

Simone Rossi's avatar
Simone Rossi committed
97
    typedef boost::shared_ptr<Activation>                     activationModelPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
98

Simone Deparis's avatar
Simone Deparis committed
99
    typedef ElectroSolver                                      electroSolver_Type;
Simone Rossi's avatar
Simone Rossi committed
100

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

Simone Deparis's avatar
Simone Deparis committed
103
    typedef ElectroIonicModel                                  ionicModel_Type;
Simone Rossi's avatar
Simone Rossi committed
104

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

Simone Deparis's avatar
Simone Deparis committed
107
    typedef ExporterHDF5< Mesh >                               exporter_Type;
Simone Rossi's avatar
Simone Rossi committed
108

Simone Deparis's avatar
Simone Deparis committed
109
    typedef boost::shared_ptr<ExporterHDF5< Mesh > >           exporterPtr_Type;
Simone Rossi's avatar
Simone Rossi committed
110
111
112
113
114
115
116
117
118
119

    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
120
121
                                    const Real&    x,
                                    const Real&    y,
Simone Rossi's avatar
Simone Rossi committed
122
123
124
125
126
                                    const Real& z,
                                    const ID&   /*i*/ ) >       function_Type;



Simone Rossi's avatar
Simone Rossi committed
127
    EMSolver(commPtr_Type comm);
Simone Rossi's avatar
Simone Rossi committed
128

Simone Deparis's avatar
Simone Deparis committed
129
    EMSolver (const EMSolver& solver);
Simone Rossi's avatar
Simone Rossi committed
130

Simone Rossi's avatar
Simone Rossi committed
131

Simone Rossi's avatar
Simone Rossi committed
132

Simone Rossi's avatar
Simone Rossi committed
133
    void loadMesh (std::string meshName, std::string meshPath)
Simone Deparis's avatar
Simone Deparis committed
134
    {
thomaskummer's avatar
thomaskummer committed
135
136
137
138
//        if (M_commPtr -> MyPID() == 0)
//        {
//            std::cout << "\n\nEMSolver: loadMesh ... " << '\r' << std::flush;
//        }
thomaskummer's avatar
thomaskummer committed
139

Simone Rossi's avatar
Simone Rossi committed
140
        M_fullMeshPtr.reset( new Mesh() );
Simone Deparis's avatar
Simone Deparis committed
141
        MeshUtility::loadMesh (M_localMeshPtr, M_fullMeshPtr, meshName, meshPath);
Simone Rossi's avatar
Simone Rossi committed
142
143
144
145
146
147
148
149
150
        if(M_commPtr)
        {
			M_localMeshPtr->setComm(M_commPtr);
			M_fullMeshPtr->setComm(M_commPtr);
        }
        else
        {
        	M_commPtr = M_localMeshPtr -> comm();
        }
thomaskummer's avatar
thomaskummer committed
151
152
153
        
        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
154
            std::cout << "EMSolver: loadMesh - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
155
        }
thomaskummer's avatar
thomaskummer committed
156

Simone Deparis's avatar
Simone Deparis committed
157
    }
Simone Rossi's avatar
Simone Rossi committed
158

thomaskummer's avatar
thomaskummer committed
159
    void setupElectroExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Deparis's avatar
Simone Deparis committed
160
161
162
    {
        M_electroSolverPtr -> setupExporter (*M_electroExporterPtr, outputFileName, problemFolder);
    }
Simone Rossi's avatar
Simone Rossi committed
163

Simone Rossi's avatar
Simone Rossi committed
164
    void setupActivationExporter ( std::string problemFolder = "./", std::string outputFileName = "ActivationSolution" )
Simone Rossi's avatar
Simone Rossi committed
165
    {
Simone Deparis's avatar
Simone Deparis committed
166
        EMUtility::setupExporter<Mesh> (*M_activationExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
Simone Rossi's avatar
Simone Rossi committed
167
    }
thomaskummer's avatar
thomaskummer committed
168
    
thomaskummer's avatar
thomaskummer committed
169
170
171
172
//    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
173

thomaskummer's avatar
thomaskummer committed
174
175
176
177
    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
178
    
thomaskummer's avatar
thomaskummer committed
179
180
181
182
183
184
185
186
187
//    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
188

thomaskummer's avatar
thomaskummer committed
189
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "MechanicalSolution" )
Simone Rossi's avatar
Simone Rossi committed
190
    {
Simone Deparis's avatar
Simone Deparis committed
191
192
193
194
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
195
196
    }

197
198
    void importHdf5 ();

Simone Deparis's avatar
Simone Deparis committed
199
200
201
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
thomaskummer's avatar
thomaskummer committed
202
                         //std::string activationTimeFileName  = "ActivationTimeSolution",
thomaskummer's avatar
thomaskummer committed
203
                         std::string mechanicsFileName  = "MechanicalSolution",
thomaskummer's avatar
thomaskummer committed
204
                         std::string vonMisesStressFileName  = "VonMisesStress");
Simone Rossi's avatar
Simone Rossi committed
205

Thomas Kummer's avatar
Thomas Kummer committed
206
207
208
209
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
210

Thomas Kummer's avatar
Thomas Kummer committed
211
212
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
213
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
214

Simone Deparis's avatar
Simone Deparis committed
215
216
217
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
218

Simone Rossi's avatar
Simone Rossi committed
219
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
220
    {
thomaskummer's avatar
thomaskummer committed
221
222
223
224
//        if (M_commPtr -> MyPID() == 0)
//        {
//            std::cout << "\nEMSolver: setupActivation ... " << '\r' << std::flush;
//        }
thomaskummer's avatar
thomaskummer committed
225
        
Simone Rossi's avatar
Simone Rossi committed
226
227
228
229
        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
230
231
232

        if (M_commPtr -> MyPID() == 0)
        {
thomaskummer's avatar
thomaskummer committed
233
            std::cout << "EMSolver: setupActivation - done " << std::flush;
thomaskummer's avatar
thomaskummer committed
234
235
        }

Simone Deparis's avatar
Simone Deparis committed
236
    }
Simone Rossi's avatar
Simone Rossi committed
237

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

Thomas Kummer's avatar
Thomas Kummer committed
243
244
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
245
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
246
    {
thomaskummer's avatar
thomaskummer committed
247
//        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: buildMechanicalSystem ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
248
        
Simone Rossi's avatar
Simone Rossi committed
249
250
    	//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
251
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
thomaskummer's avatar
thomaskummer committed
252
        
thomaskummer's avatar
thomaskummer committed
253
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: buildMechanicalSystem - done" << std::flush;
Simone Deparis's avatar
Simone Deparis committed
254
    }
Simone Rossi's avatar
Simone Rossi committed
255

Simone Rossi's avatar
Simone Rossi committed
256
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
257
    {
thomaskummer's avatar
thomaskummer committed
258
//        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: buildElectroSystem ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
259

Simone Deparis's avatar
Simone Deparis committed
260
        M_electroSolverPtr -> setupMatrices();
thomaskummer's avatar
thomaskummer committed
261
        
thomaskummer's avatar
thomaskummer committed
262
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: buildElectroSystem - done" << std::flush;
Simone Deparis's avatar
Simone Deparis committed
263
    }
Simone Rossi's avatar
Simone Rossi committed
264

Simone Rossi's avatar
Simone Rossi committed
265
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
266
267
268
269
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
270

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

Simone Deparis's avatar
Simone Deparis committed
275
        M_electroSolverPtr -> setInitialConditions();
thomaskummer's avatar
thomaskummer committed
276
        
thomaskummer's avatar
thomaskummer committed
277
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: initializeElectroVariables - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
278

Simone Deparis's avatar
Simone Deparis committed
279
    }
Simone Rossi's avatar
Simone Rossi committed
280

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

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


292
293
294
295
296
    void setupFiberVector( const std::string& fileName,
						   const std::string& fieldName,
						   const std::string& postDir = "./",
						   const std::string& polynomialDegree = "P1"  )
    {
thomaskummer's avatar
thomaskummer committed
297
//        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: setupFiberVector ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
298

299
300
    	setupMechanicalFiberVector(fileName, fieldName, postDir, polynomialDegree);
    	M_electroSolverPtr->setFiberPtr(getMechanicsFibers());
thomaskummer's avatar
thomaskummer committed
301
        
thomaskummer's avatar
thomaskummer committed
302
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: setupFiberVector - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
303

304
305
306
307
308
309
310
    }

    void setupMechanicalFiberVector( const std::string& fileName,
    		                         const std::string& fieldName,
									 const std::string& postDir = "./",
                                     const std::string& polynomialDegree = "P1"  )
    {
thomaskummer's avatar
thomaskummer committed
311
//        if (M_commPtr -> MyPID() == 0) std::cout << "\nEMSolver: setupMechanicalFiberVector ... " << '\r' << std::flush;
thomaskummer's avatar
thomaskummer committed
312
        
313
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
314
315
316
317
318
319
320
321
        
//        if ( polynomialDegree == "P1" )
//        {
//            FESpace<RegionMesh<LinearTetra> , MapEpetra > p2FESpace ( fullMeshPtr(), "P2", 1, fullMeshPtr() -> comm() );
//            vector_Type p2FibersRep (*getMechanicsFibers(), Repeated);
//            vector_Type p1FibersRep = M_EMStructuralOperatorPtr -> dispFESpacePtr() -> feToFEInterpolate(p2FESpace, p2FibersRep);
//            getMechanicsFibers().reset ( new vector_Type (p1FibersRep, Unique) );
//        }
thomaskummer's avatar
thomaskummer committed
322
        if (M_commPtr -> MyPID() == 0) std::cout << "EMSolver: setupMechanicalFiberVector - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
323

324
325
326
327
328
329
330
    }

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

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

336
337
338
339
340
341
342
343
344
345
346
347
    }

    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
348
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
349
350
351
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
352

Simone Rossi's avatar
Simone Rossi committed
353
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
354
355
356
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
357
358


Simone Rossi's avatar
Simone Rossi committed
359
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
360
361
362
363
364
365
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
366
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
367
    }
Simone Rossi's avatar
Simone Rossi committed
368

Simone Rossi's avatar
Simone Rossi committed
369
370
371
372
373
374
    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
375
376
377
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
378

Simone Rossi's avatar
Simone Rossi committed
379
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
380
381
382
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
383

Simone Rossi's avatar
Simone Rossi committed
384
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
385
386
387
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
388
    
thomaskummer's avatar
thomaskummer committed
389
390
391
392
//    vectorPtr_Type activationTimePtr()
//    {
//        return M_activationTimePtr;
//    }
Simone Rossi's avatar
Simone Rossi committed
393

thomaskummer's avatar
thomaskummer committed
394
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
395
396
397
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
398

Simone Deparis's avatar
Simone Deparis committed
399
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
400

Simone Deparis's avatar
Simone Deparis committed
401
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
402

Simone Deparis's avatar
Simone Deparis committed
403
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
404

Simone Rossi's avatar
Simone Rossi committed
405
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
406
407
408
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
409
410


Simone Rossi's avatar
Simone Rossi committed
411
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
412
413
414
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
415

416
417
418
419
420
    void solveMechanicsLin()
    {
        M_EMStructuralOperatorPtr -> solveLin ();
    }

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

Simone Rossi's avatar
Simone Rossi committed
423

Simone Rossi's avatar
Simone Rossi committed
424
    void solveActivation (Real dt);
thomaskummer's avatar
thomaskummer committed
425
    
thomaskummer's avatar
thomaskummer committed
426
    void computeI4f (VectorEpetra& i4f, VectorEpetra& f0_, VectorEpetra& disp, solidFESpacePtr_Type feSpacePtr);
Simone Rossi's avatar
Simone Rossi committed
427

Simone Rossi's avatar
Simone Rossi committed
428
429
430
431
432
433
434
435
436
437
438
439
440
    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
441
442
443
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
    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
459

460
    void showMe() const {}
461
    
thomaskummer's avatar
thomaskummer committed
462
463
464
465
466
    WallTensionEstimator<RegionMesh<LinearTetra> >& tensionEstimator()
    {
        return M_wteTotal;
    }
    
467
468
469
470
471
472
473
474
475
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
476

thomaskummer's avatar
thomaskummer committed
477
    std::vector<meshPtr_Type> mesh()
thomaskummer's avatar
thomaskummer committed
478
    {
thomaskummer's avatar
thomaskummer committed
479
        std::vector<meshPtr_Type> meshVector;
thomaskummer's avatar
thomaskummer committed
480
481
        meshVector.push_back(M_fullMeshPtr);
        meshVector.push_back(M_localMeshPtr);
thomaskummer's avatar
thomaskummer committed
482
        return meshVector;
thomaskummer's avatar
thomaskummer committed
483
    }
484
    
Simone Rossi's avatar
Simone Rossi committed
485
protected:
Simone Rossi's avatar
Simone Rossi committed
486
public:
Simone Deparis's avatar
Simone Deparis committed
487
488
489
490
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
491

Simone Deparis's avatar
Simone Deparis committed
492
493
494
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
495
//    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
496
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
thomaskummer's avatar
thomaskummer committed
497
498
//    exporterPtr_Type                     M_vonMisesStressExporterPtrP;
//    exporterPtr_Type                     M_vonMisesStressExporterPtrA;
thomaskummer's avatar
thomaskummer committed
499
    
Simone Deparis's avatar
Simone Deparis committed
500
501
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
502
    
thomaskummer's avatar
thomaskummer committed
503
    //vectorPtr_Type                       M_activationTimePtr;
Simone Rossi's avatar
Simone Rossi committed
504
505

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
506
    
thomaskummer's avatar
thomaskummer committed
507
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteTotal;
thomaskummer's avatar
thomaskummer committed
508
509
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wtePassive;
//    WallTensionEstimator<RegionMesh<LinearTetra> > M_wteActive;
thomaskummer's avatar
thomaskummer committed
510

Simone Rossi's avatar
Simone Rossi committed
511

Simone Deparis's avatar
Simone Deparis committed
512
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
513

Simone Rossi's avatar
Simone Rossi committed
514
515
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
516

Simone Rossi's avatar
Simone Rossi committed
517
};
Simone Rossi's avatar
Simone Rossi committed
518
519
520

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
521
522
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
523
524
525
526
527
528
529
530
531
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
532
    //M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
533
    M_oneWayCoupling     (true),
thomaskummer's avatar
thomaskummer committed
534
    M_wteTotal ( ),
thomaskummer's avatar
thomaskummer committed
535
536
//    M_wtePassive ( ),
//    M_wteActive ( ),
Simone Rossi's avatar
Simone Rossi committed
537
538
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
539
540
541
542
543
544
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
545
546
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
547
548
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
549
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
550
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
551
552
553
554
555
    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
556
    //M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
557
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
thomaskummer's avatar
thomaskummer committed
558
    M_wteTotal                   (solver.M_wteTotal),
thomaskummer's avatar
thomaskummer committed
559
560
//    M_wtePassive                   (solver.M_wtePassive),
//    M_wteActive                   (solver.M_wteActive),
Simone Rossi's avatar
Simone Rossi committed
561
562
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
563
564
565
566
567
568
569

{
}


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

Simone Rossi's avatar
Simone Rossi committed
571
572
573
574


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
575
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
576
void
Thomas Kummer's avatar
Thomas Kummer committed
577
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
578
{
Simone Rossi's avatar
Simone Rossi committed
579
    M_data.setup (dataFile);
thomaskummer's avatar
thomaskummer committed
580
    
thomaskummer's avatar
thomaskummer committed
581
582
583
584
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\n\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " ...";
//    }
thomaskummer's avatar
thomaskummer committed
585
586
    
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
587
588
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
thomaskummer's avatar
thomaskummer committed
589
590
591
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
592
        std::cout << "\nEMSolver: setup (simulation endtime = " << M_data.activationParameter<Real>("endtime") << ")" << " - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
593
594
    }

Simone Rossi's avatar
Simone Rossi committed
595
596
}

Simone Rossi's avatar
Simone Rossi committed
597
598


Simone Rossi's avatar
Simone Rossi committed
599
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
600
void
Thomas Kummer's avatar
Thomas Kummer committed
601
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
602
{
thomaskummer's avatar
thomaskummer committed
603
604
605
606
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\nEMSolver: setupElectroSolver ... " << '\r' << std::flush;
//    }
Simone Rossi's avatar
Simone Rossi committed
607
608
609
610
611
	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
612
613
614

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
Simone Rossi's avatar
Simone Rossi committed
615
    M_electroSolverPtr->setParameters();
thomaskummer's avatar
thomaskummer committed
616
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Simone Rossi's avatar
Simone Rossi committed
617
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
thomaskummer's avatar
thomaskummer committed
618
    //if (M_commPtr -> MyPID() == 0) M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
619
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
620

Thomas Kummer's avatar
Thomas Kummer committed
621
622
623
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
624
625
626
627
628
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
629
//    }
thomaskummer's avatar
thomaskummer committed
630
    
Simone Deparis's avatar
Simone Deparis committed
631
632
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
633
        std::cout << "\nEMSolver: setupElectroSolver - done" << std::flush;
Simone Deparis's avatar
Simone Deparis committed
634
    }
Simone Rossi's avatar
Simone Rossi committed
635

Simone Rossi's avatar
Simone Rossi committed
636
637
638
}


Simone Rossi's avatar
Simone Rossi committed
639
640
/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
641
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
642
void
Simone Rossi's avatar
Simone Rossi committed
643
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
644
{
thomaskummer's avatar
thomaskummer committed
645
646
647
648
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\nEMSolver: setupMechanicalSolver ... " << '\r' << std::flush;
//    }
thomaskummer's avatar
thomaskummer committed
649
    
Simone Rossi's avatar
Simone Rossi committed
650
651
652
653
654
655
    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
656
657
658
659
660
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

thomaskummer's avatar
thomaskummer committed
661
    dFESpace->setQuadRule(quadRuleTetra4pt);
thomaskummer's avatar
thomaskummer committed
662
663
    //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
664
    
Simone Deparis's avatar
Simone Deparis committed
665
666
667
668
669
670
671
    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
672
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
thomaskummer's avatar
thomaskummer committed
673
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data, dataFile);
Simone Rossi's avatar
Simone Rossi committed
674

thomaskummer's avatar
thomaskummer committed
675
    M_wteTotal.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0, M_EMStructuralOperatorPtr->EMMaterial());
thomaskummer's avatar
thomaskummer committed
676
677
//    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
678
679
680
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
681
        std::cout << "\nEMSolver: setupMechanicalSolver - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
682
683
    }

Simone Rossi's avatar
Simone Rossi committed
684
685
686
687
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
688
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
689
void
Simone Rossi's avatar
Simone Rossi committed
690
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
691
692
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
693
{
thomaskummer's avatar
thomaskummer committed
694
695
696
697
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\nEMSolver: setupMechanicalBC ... " << '\r' << std::flush;
//    }
thomaskummer's avatar
thomaskummer committed
698

Simone Deparis's avatar
Simone Deparis committed
699
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
700
701
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
702
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
thomaskummer's avatar
thomaskummer committed
703
704
705
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
706
        std::cout << "EMSolver: setupMechanicalBC - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
707
    }
Simone Rossi's avatar
Simone Rossi committed
708
709
}

710
711
712
713
714
715
716
    
/////////////////////
//restart hdf5
template<typename Mesh , typename ElectroSolver>
void
EMSolver<Mesh, ElectroSolver>::importHdf5 ()
{
thomaskummer's avatar
thomaskummer committed
717
718
    M_electroExporterPtr -> import (30);
    M_activationExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
719
    //M_activationTimeExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
720
    M_vonMisesStressExporterPtr -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
721
722
    //M_vonMisesStressExporterPtrP -> importHdf5 (30);
    //M_vonMisesStressExporterPtrA -> importHdf5 (30);
thomaskummer's avatar
thomaskummer committed
723
    M_mechanicsExporterPtr -> importHdf5 (30);
724
725
}
                                                   
Simone Rossi's avatar
Simone Rossi committed
726
727
728

/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
729
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
730
void
thomaskummer's avatar
thomaskummer committed
731
732
733
EMSolver<Mesh, ElectroSolver>::setupExporters ( std::string problemFolder,
                                                std::string electroFileName,
                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
734
                                                //std::string activationTimeFileName,
thomaskummer's avatar
thomaskummer committed
735
                                                std::string mechanicsFileName,
thomaskummer's avatar
thomaskummer committed
736
                                                std::string vonMisesStressFileName)
Simone Rossi's avatar
Simone Rossi committed
737
{
thomaskummer's avatar
thomaskummer committed
738
739
740
741
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\nEMSolver: setupExporters ... " << '\r' << std::flush;
//    }
thomaskummer's avatar
thomaskummer committed
742
743
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
744
745
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
746
747
748
749
750
751
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
752
    // Activation
Simone Deparis's avatar
Simone Deparis committed
753
754
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
thomaskummer's avatar
thomaskummer committed
755
756
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
thomaskummer's avatar
thomaskummer committed
757
                                             M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
758
759
                                             M_activationModelPtr -> fiberActivationPtr(),
                                             UInt (0) );
Simone Deparis's avatar
Simone Deparis committed
760

thomaskummer's avatar
thomaskummer committed
761
762
763
764
765
766
767
768
769
770
771
772
//    // 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
773
    
thomaskummer's avatar
thomaskummer committed
774
775
776
777
778
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
779
                                                "Von Mises Stress",
thomaskummer's avatar
thomaskummer committed
780
                                                M_electroSolverPtr -> feSpacePtr(),
thomaskummer's avatar
thomaskummer committed
781
782
783
                                                M_wteTotal.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
//    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
801

thomaskummer's avatar
thomaskummer committed
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
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
//    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
867
    
thomaskummer's avatar
thomaskummer committed
868
    // Mechanics
thomaskummer's avatar
thomaskummer committed
869
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
870
871
872
873
874
875
876
877
878
879
880
    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
881
882
883
884
885
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
thomaskummer's avatar
thomaskummer committed
886
887
888
    
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
889
        std::cout << "EMSolver: setupExporters - done" << std::flush;
thomaskummer's avatar
thomaskummer committed
890
        std::cout << "\n\n";
thomaskummer's avatar
thomaskummer committed
891
    }
Simone Rossi's avatar
Simone Rossi committed
892
893
}

thomaskummer's avatar
restart    
thomaskummer committed
894
895
896
897
898
899
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
900
    //M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
901
    M_vonMisesStressExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
thomaskummer committed
902
903
    //M_vonMisesStressExporterPtrP -> setTimeIndex (time);
    //M_vonMisesStressExporterPtrA -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
904
905
906
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
907
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
908
void
thomaskummer's avatar
thomaskummer committed
909
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
910
{
thomaskummer's avatar
thomaskummer committed
911
912
913
914
    M_wteTotal.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
    M_wteTotal.analyzeTensionsRecoveryVonMisesStress();
    
    M_vonMisesStressExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
915
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
916
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
917
    //M_activationTimeExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
918
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
919
920
}

Simone Rossi's avatar
Simone Rossi committed
921
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
922
void
Simone Rossi's avatar
Simone Rossi committed
923
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
924
{
Simone Deparis's avatar
Simone Deparis committed
925
    M_electroExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
926
927
    M_activationExporterPtr -> closeFile();
    //M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
928
    M_mechanicsExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
929
    M_vonMisesStressExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
930
931
932
933
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
934
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
935
void
Simone Rossi's avatar
Simone Rossi committed
936
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
937
{
Simone Deparis's avatar
Simone Deparis committed
938
939
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
940
941
942
943
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
944
945
}

Simone Rossi's avatar
Simone Rossi committed
946
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
947
void
Simone Rossi's avatar
Simone Rossi committed
948
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
949
{
Simone Deparis's avatar
Simone Deparis committed
950
951
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
952
953
954
955
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
956
957
958
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
959
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
960
void
Simone Rossi's avatar
Simone Rossi committed
961
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
962
{
thomaskummer's avatar
thomaskummer committed
963
964
965
966
//    if (M_commPtr -> MyPID() == 0)
//    {
//        std::cout << "\nEMSolver: solveElectrophysiology ... " << '\r' << std::flush;
//    }
thomaskummer's avatar
thomaskummer committed
967
    
Simone Deparis's avatar
Simone Deparis committed
968
    setAppliedCurrent ( stimulus, time );
thomaskummer's avatar
thomaskummer committed
969

Simone Rossi's avatar
Simone Rossi committed
970
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
thomaskummer's avatar
thomaskummer committed
971
    //M_electroSolverPtr -> solveOneStepGatingVariablesRL();
thomaskummer's avatar
thomaskummer committed
972
    
Simone Rossi's avatar
Simone Rossi committed
973
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
974
    
thomaskummer's avatar
thomaskummer committed
975
    //M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
thomaskummer's avatar
thomaskummer committed
976
977
    if (M_commPtr -> MyPID() == 0)
    {
thomaskummer's avatar
thomaskummer committed
978
        std::cout << "EMSolver: solveElectrophysiology - done" << std::flush;