EMSolver.hpp 34.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
55
#include <lifev/structure/solver/WallTensionEstimator.hpp>

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

Simone Rossi's avatar
Simone Rossi committed
59
60
61
62
63
64
65
66
67
68
69
//
//#include <lifev/em/solver/EMStructuralOperator.hpp>
//#include <lifev/em/solver/EMGeneralizedActiveHolzapfelOgdenMaterial.hpp>
//#include <lifev/em/solver/EMActiveStrainSolver.hpp>
//#include <lifev/core/interpolation/RBFlocallyRescaledVectorial.hpp>
//#include <lifev/core/interpolation/RBFlocallyRescaledScalar.hpp>
//#include <lifev/core/interpolation/RBFrescaledVectorial.hpp>
//#include <lifev/core/interpolation/RBFrescaledScalar.hpp>
////#include <lifev/core/interpolation/RBFscalar.hpp>
//#include <lifev/core/interpolation/RBFvectorial.hpp>
//
Simone Rossi's avatar
Simone Rossi committed
70
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>
Simone Rossi's avatar
Simone Rossi committed
71
72
73
//
//
//#include <lifev/em/solver/EMEvaluate.hpp>
Simone Rossi's avatar
Simone Rossi committed
74
75
76
77
78
79

namespace LifeV
{

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    typedef FESpace< RegionMesh<LinearTetra>, MapEpetra >      solidFESpace_Type;

    typedef boost::shared_ptr<solidFESpace_Type>                solidFESpacePtr_Type;

    typedef ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > solidETFESpace_Type;

    typedef boost::shared_ptr<solidETFESpace_Type>              solidETFESpacePtr_Type;

    typedef boost::function < Real (const Real& t,
Simone Deparis's avatar
Simone Deparis committed
135
136
                                    const Real&    x,
                                    const Real&    y,
Simone Rossi's avatar
Simone Rossi committed
137
138
139
140
141
                                    const Real& z,
                                    const ID&   /*i*/ ) >       function_Type;



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

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

Simone Rossi's avatar
Simone Rossi committed
146

Simone Rossi's avatar
Simone Rossi committed
147

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

Simone Deparis's avatar
Simone Deparis committed
163
    }
Simone Rossi's avatar
Simone Rossi committed
164

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

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

thomaskummer's avatar
thomaskummer committed
180
181
182
183
184
    void setupVonMisesStressExporter ( std::string problemFolder = "./", std::string outputFileName = "VonMisesStress" )
    {
        EMUtility::setupExporter<Mesh> (*M_vonMisesStressExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
    }

Simone Rossi's avatar
Simone Rossi committed
185
    void setupMechanicsExporter ( std::string problemFolder = "./", std::string outputFileName = "ElectroSolution" )
Simone Rossi's avatar
Simone Rossi committed
186
    {
Simone Deparis's avatar
Simone Deparis committed
187
188
189
190
        if (M_mechanicsExporterPtr)
        {
            EMUtility::setupExporter<Mesh> (*M_mechanicsExporterPtr, M_localMeshPtr, M_commPtr, outputFileName, problemFolder);
        }
Simone Rossi's avatar
Simone Rossi committed
191
192
    }

Simone Deparis's avatar
Simone Deparis committed
193
194
195
    void setupExporters (std::string problemFolder   = "./",
                         std::string electroFileName = "ElectroSolution",
                         std::string activationFileName  = "ActivationSolution",
thomaskummer's avatar
thomaskummer committed
196
                         std::string activationTimeFileName  = "ActivationTimeSolution",
thomaskummer's avatar
thomaskummer committed
197
198
                         std::string mechanicsFileName  = "MechanicalSolution",
                         std::string vonMisesStressFileName  = "VonMisesStress");
Simone Rossi's avatar
Simone Rossi committed
199

Thomas Kummer's avatar
Thomas Kummer committed
200
201
202
203
    void setupElectroSolver ( GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setupElectroSolver (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
204

Thomas Kummer's avatar
Thomas Kummer committed
205
206
    void setupElectroSolver ( GetPot& dataFile );
    
Simone Deparis's avatar
Simone Deparis committed
207
    void setupMechanicalSolver ( GetPot& dataFile);
Simone Rossi's avatar
Simone Rossi committed
208

Simone Deparis's avatar
Simone Deparis committed
209
210
211
    void setupMechanicalBC (std::string data_file_name,
                            std::string section,
                            solidFESpacePtr_Type dFESpace);
Simone Rossi's avatar
Simone Rossi committed
212

Simone Rossi's avatar
Simone Rossi committed
213
    void setupActivation (const MapEpetra& map)
Simone Deparis's avatar
Simone Deparis committed
214
215
216
217
218
    {
        if (M_commPtr -> MyPID() == 0)
        {
            std::cout << "EMS - setting up activation solver\n";
        }
Simone Rossi's avatar
Simone Rossi committed
219
220
221
222
        M_activationModelPtr.reset ( Activation::EMActivationFactory::instance().createObject ( M_data.activationParameter<std::string>( "ActivationModel" ) ) );
        M_activationModelPtr->setup(M_data, map);
        M_activationModelPtr->setVariablesPtr(*M_electroSolverPtr);
        M_activationModelPtr->setI4fPtr( M_EMStructuralOperatorPtr -> I4fPtr() );
Simone Deparis's avatar
Simone Deparis committed
223
    }
Simone Rossi's avatar
Simone Rossi committed
224

Thomas Kummer's avatar
Thomas Kummer committed
225
226
227
228
    void setup (GetPot& dataFile, Teuchos::ParameterList& list)
    {
        setup (dataFile);
    }
Simone Rossi's avatar
Simone Rossi committed
229

Thomas Kummer's avatar
Thomas Kummer committed
230
231
    void setup (GetPot& dataFile );
    
Simone Rossi's avatar
Simone Rossi committed
232
    void buildMechanicalSystem()
Simone Deparis's avatar
Simone Deparis committed
233
    {
Simone Rossi's avatar
Simone Rossi committed
234
235
    	//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
236
237
        M_EMStructuralOperatorPtr -> buildSystem (1.0);
    }
Simone Rossi's avatar
Simone Rossi committed
238

Simone Rossi's avatar
Simone Rossi committed
239
    void buildElectroSystem()
Simone Deparis's avatar
Simone Deparis committed
240
241
242
    {
        M_electroSolverPtr -> setupMatrices();
    }
Simone Rossi's avatar
Simone Rossi committed
243

Simone Rossi's avatar
Simone Rossi committed
244
     void buildSystem()
Simone Deparis's avatar
Simone Deparis committed
245
246
247
248
    {
        buildMechanicalSystem();
        buildElectroSystem();
    }
Simone Rossi's avatar
Simone Rossi committed
249

Simone Rossi's avatar
Simone Rossi committed
250
    void initializeElectroVariables()
Simone Deparis's avatar
Simone Deparis committed
251
252
253
    {
        M_electroSolverPtr -> setInitialConditions();
    }
Simone Rossi's avatar
Simone Rossi committed
254

Simone Rossi's avatar
Simone Rossi committed
255
    void initialize()
Simone Deparis's avatar
Simone Deparis committed
256
257
258
    {
        initializeElectroVariables();
    }
Simone Rossi's avatar
Simone Rossi committed
259

Simone Rossi's avatar
Simone Rossi committed
260
    bcInterfacePtr_Type bcInterfacePtr()
Simone Deparis's avatar
Simone Deparis committed
261
262
263
    {
        return M_bcInterfacePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
264
265


266
267
268
269
270
271
272
273
274
275
276
277
278
279
    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
280
        
281
        ElectrophysiologyUtility::importVectorField (getMechanicsFibers(),  fileName,  fieldName, M_localMeshPtr, postDir, polynomialDegree );
thomaskummer's avatar
thomaskummer committed
282
283
284
285
286
287
288
289
290
        
//        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) );
//        }

291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
    }

    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
311
    void setupMechanicalFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
312
313
314
    {
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector ( fx, fy, fz);
    }
Simone Rossi's avatar
Simone Rossi committed
315

Simone Rossi's avatar
Simone Rossi committed
316
    void setupElectroFiberVector ( VectorSmall<3>& fibers)
Simone Deparis's avatar
Simone Deparis committed
317
318
319
    {
        M_electroSolverPtr -> setupFibers (fibers);
    }
Simone Rossi's avatar
Simone Rossi committed
320
321


Simone Rossi's avatar
Simone Rossi committed
322
    void setupFiberVector ( Real fx, Real fy, Real fz )
Simone Deparis's avatar
Simone Deparis committed
323
324
325
326
327
328
    {
        VectorSmall<3> f;
        f[0] = fx;
        f[1] = fy;
        f[2] = fz;
        setupElectroFiberVector (f);
Simone Rossi's avatar
Simone Rossi committed
329
        M_EMStructuralOperatorPtr -> EMMaterial() -> setupFiberVector(fx, fy, fz);
Simone Deparis's avatar
Simone Deparis committed
330
    }
Simone Rossi's avatar
Simone Rossi committed
331

Simone Rossi's avatar
Simone Rossi committed
332
333
334
335
336
337
    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
338
339
340
    {
        return M_electroSolverPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
341

Simone Rossi's avatar
Simone Rossi committed
342
    structuralOperatorPtr_Type structuralOperatorPtr()
Simone Deparis's avatar
Simone Deparis committed
343
344
345
    {
        return M_EMStructuralOperatorPtr;
    }
Simone Rossi's avatar
Simone Rossi committed
346

Simone Rossi's avatar
Simone Rossi committed
347
    activationModelPtr_Type  activationModelPtr()
Simone Deparis's avatar
Simone Deparis committed
348
349
350
    {
        return M_activationModelPtr;
    }
thomaskummer's avatar
thomaskummer committed
351
352
353
354
355
    
    vectorPtr_Type activationTimePtr()
    {
        return M_activationTimePtr;
    }
Simone Rossi's avatar
Simone Rossi committed
356

thomaskummer's avatar
thomaskummer committed
357
    void saveSolution (Real time, const bool& restart = 0);
thomaskummer's avatar
restart    
thomaskummer committed
358
359
360
    
    void setTimeIndex (const UInt& time);

Simone Rossi's avatar
Simone Rossi committed
361

Simone Deparis's avatar
Simone Deparis committed
362
    void closeExporters();
Simone Rossi's avatar
Simone Rossi committed
363

Simone Deparis's avatar
Simone Deparis committed
364
    void oneWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
365

Simone Deparis's avatar
Simone Deparis committed
366
    void twoWayCoupling();
Simone Rossi's avatar
Simone Rossi committed
367

Simone Rossi's avatar
Simone Rossi committed
368
    void setAppliedCurrent (function_Type& stimulus, Real time = 0.0)
Simone Deparis's avatar
Simone Deparis committed
369
370
371
    {
        M_electroSolverPtr -> setAppliedCurrentFromFunction (stimulus, time);
    }
Simone Rossi's avatar
Simone Rossi committed
372
373


Simone Rossi's avatar
Simone Rossi committed
374
    void solveMechanics()
Simone Deparis's avatar
Simone Deparis committed
375
376
377
    {
        M_EMStructuralOperatorPtr -> iterate ( M_bcInterfacePtr -> handler() );
    }
Simone Rossi's avatar
Simone Rossi committed
378

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

Simone Rossi's avatar
Simone Rossi committed
381

Simone Rossi's avatar
Simone Rossi committed
382
    void solveActivation (Real dt);
Simone Rossi's avatar
Simone Rossi committed
383

Simone Rossi's avatar
Simone Rossi committed
384
385
386
387
388
389
390
391
392
393
394
395
396
    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
397
398
399
    //  {
    //      return *M_bcInterfacePtr;
    //  }
Simone Rossi's avatar
Simone Rossi committed
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
    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
415

416
    void showMe() const {}
417
418
419
420
421
422
423
424
425
426
    
    meshPtr_Type fullMeshPtr()
    {
        return M_fullMeshPtr;
    }
    
    meshPtr_Type localMeshPtr()
    {
        return M_localMeshPtr;
    }
427

428
    
Simone Rossi's avatar
Simone Rossi committed
429
protected:
Simone Rossi's avatar
Simone Rossi committed
430
public:
Simone Deparis's avatar
Simone Deparis committed
431
432
433
434
    electroSolverPtr_Type                M_electroSolverPtr;
    activationModelPtr_Type              M_activationModelPtr;
    bcInterfacePtr_Type                  M_bcInterfacePtr;
    structuralOperatorPtr_Type           M_EMStructuralOperatorPtr;
Simone Rossi's avatar
Simone Rossi committed
435

Simone Deparis's avatar
Simone Deparis committed
436
437
438
    exporterPtr_Type                     M_electroExporterPtr;
    exporterPtr_Type                     M_activationExporterPtr;
    exporterPtr_Type                     M_mechanicsExporterPtr;
thomaskummer's avatar
thomaskummer committed
439
    exporterPtr_Type                     M_activationTimeExporterPtr;
thomaskummer's avatar
thomaskummer committed
440
441
    exporterPtr_Type                     M_vonMisesStressExporterPtr;
    
Simone Deparis's avatar
Simone Deparis committed
442
443
    meshPtr_Type                         M_localMeshPtr;
    meshPtr_Type                         M_fullMeshPtr;
thomaskummer's avatar
thomaskummer committed
444
445
    
    vectorPtr_Type                       M_activationTimePtr;
thomaskummer's avatar
thomaskummer committed
446
    vectorPtr_Type                       M_vonMisesStressPtr;
Simone Rossi's avatar
Simone Rossi committed
447
448

    bool                                 M_oneWayCoupling;
thomaskummer's avatar
thomaskummer committed
449
450
    
    WallTensionEstimator<RegionMesh<LinearTetra> > M_wte;
Simone Rossi's avatar
Simone Rossi committed
451

Simone Deparis's avatar
Simone Deparis committed
452
    commPtr_Type                         M_commPtr;
Simone Rossi's avatar
Simone Rossi committed
453

Simone Rossi's avatar
Simone Rossi committed
454
455
    EMData                               M_data;

Simone Rossi's avatar
Simone Rossi committed
456

Simone Rossi's avatar
Simone Rossi committed
457
};
Simone Rossi's avatar
Simone Rossi committed
458
459
460

/////////////////////
// CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
461
462
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver(commPtr_Type comm) :
Simone Rossi's avatar
Simone Rossi committed
463
464
465
466
467
468
469
470
471
    M_electroSolverPtr       ( ),
    M_activationModelPtr    ( ),
    M_bcInterfacePtr        ( ),
    M_EMStructuralOperatorPtr(),
    M_electroExporterPtr ( ),
    M_activationExporterPtr ( ),
    M_mechanicsExporterPtr  ( ),
    M_localMeshPtr      ( ),
    M_fullMeshPtr      ( ),
thomaskummer's avatar
thomaskummer committed
472
    M_activationTimePtr     ( ),
Simone Rossi's avatar
Simone Rossi committed
473
    M_oneWayCoupling     (true),
thomaskummer's avatar
thomaskummer committed
474
    M_wte ( ),
Simone Rossi's avatar
Simone Rossi committed
475
476
    M_commPtr               (comm),
    M_data                    ()
Simone Rossi's avatar
Simone Rossi committed
477
478
479
480
481
482
{
}


/////////////////////
// COPY CONSTRUCTORS
Simone Rossi's avatar
Simone Rossi committed
483
484
template<typename Mesh , typename ElectroSolver>
EMSolver<Mesh, ElectroSolver>::EMSolver (const EMSolver& solver) :
Simone Deparis's avatar
Simone Deparis committed
485
486
    M_electroSolverPtr (solver.M_electroSolverPtr),
    M_activationModelPtr (solver.M_activationModelPtr),
Simone Rossi's avatar
Simone Rossi committed
487
    M_bcInterfacePtr        ( solver.M_bcInterfacePtr),
Simone Deparis's avatar
Simone Deparis committed
488
    M_EMStructuralOperatorPtr (solver.M_EMStructuralOperatorPtr),
Simone Rossi's avatar
Simone Rossi committed
489
490
491
492
493
    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
494
    M_activationTimePtr     ( solver.M_activationTimePtr),
Simone Rossi's avatar
Simone Rossi committed
495
    M_oneWayCoupling     ( solver.M_oneWayCoupling),
thomaskummer's avatar
thomaskummer committed
496
    M_wte                   (solver.M_wte),
Simone Rossi's avatar
Simone Rossi committed
497
498
    M_commPtr               ( solver.M_commPtr),
    M_data                   (solver.M_data)
Simone Rossi's avatar
Simone Rossi committed
499
500
501
502
503
504
505

{
}


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

Simone Rossi's avatar
Simone Rossi committed
507
508
509
510


/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
511
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
512
void
Thomas Kummer's avatar
Thomas Kummer committed
513
EMSolver<Mesh, ElectroSolver>::setup ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
514
{
Simone Rossi's avatar
Simone Rossi committed
515
516
    M_data.setup (dataFile);
    std::cout << "\nEMSolver - endtime = " << M_data.activationParameter<Real>("endtime");
Thomas Kummer's avatar
Thomas Kummer committed
517
    setupElectroSolver ( dataFile );
Simone Deparis's avatar
Simone Deparis committed
518
519
    if (M_commPtr -> MyPID() == 0)
    {
Simone Rossi's avatar
Simone Rossi committed
520
        std::cout << "\nEMS - electro solver setup done! ";
Simone Deparis's avatar
Simone Deparis committed
521
522
523
    }
    setupMechanicalSolver ( dataFile );
    setupActivation ( M_electroSolverPtr -> potentialPtr() ->map() );
Simone Rossi's avatar
Simone Rossi committed
524

Simone Rossi's avatar
Simone Rossi committed
525
526
}

Simone Rossi's avatar
Simone Rossi committed
527
528


Simone Rossi's avatar
Simone Rossi committed
529
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
530
void
Thomas Kummer's avatar
Thomas Kummer committed
531
EMSolver<Mesh, ElectroSolver>::setupElectroSolver ( GetPot& dataFile )
Simone Rossi's avatar
Simone Rossi committed
532
{
Simone Deparis's avatar
Simone Deparis committed
533
534
535
536
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - creating ionic model ";
    }
Simone Rossi's avatar
Simone Rossi committed
537
538
539
540
541
	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
542
543
544
545
546
547
548
549

    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up electrophysiology solver ";
    }

    M_electroSolverPtr.reset ( new ElectroSolver() );
    M_electroSolverPtr -> setIonicModelPtr (ionicModelPtr);
Simone Rossi's avatar
Simone Rossi committed
550
551
552
553
    M_electroSolverPtr->setParameters();
    M_electroSolverPtr ->showParameters();
	M_electroSolverPtr -> setParametersFromEMData ( M_data );
    M_electroSolverPtr ->showParameters();
Thomas Kummer's avatar
Thomas Kummer committed
554
    M_electroSolverPtr->init (M_localMeshPtr); //(M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
555

Thomas Kummer's avatar
Thomas Kummer committed
556
557
558
//    if (M_localMeshPtr)
//    {
//        M_electroSolverPtr -> setLocalMeshPtr (M_localMeshPtr);
Simone Deparis's avatar
Simone Deparis committed
559
560
561
562
563
        if (M_fullMeshPtr)
        {
            M_electroSolverPtr -> setFullMeshPtr (M_fullMeshPtr);
        }
        M_electroSolverPtr ->  setup (dataFile, ionicModelPtr->Size() );
Thomas Kummer's avatar
Thomas Kummer committed
564
//    }
Simone Deparis's avatar
Simone Deparis committed
565
566
567
568
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "... `Done\n";
    }
Simone Rossi's avatar
Simone Rossi committed
569

Simone Rossi's avatar
Simone Rossi committed
570
571
572
}


Simone Rossi's avatar
Simone Rossi committed
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
Simone Rossi's avatar
Simone Rossi committed
577
EMSolver<Mesh, ElectroSolver>::setupMechanicalSolver ( GetPot& dataFile)
Simone Rossi's avatar
Simone Rossi committed
578
{
Simone Deparis's avatar
Simone Deparis committed
579
580
581
582
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up mechanical solver\n";
    }
Simone Rossi's avatar
Simone Rossi committed
583
584
585
586
587
588
    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
589
590
591
592
593
594
595
596
597
598
599
600
                                                                  & (dFESpace->refFE() ),
                                                                  & (dFESpace->fe().geoMap() ),
                                                                  M_commPtr) );
    std::string data_file_name = dataFile.get (0, "NO_DATA_FILENAME_FOUND");

    setupMechanicalBC (data_file_name, "solid",  dFESpace);
    M_EMStructuralOperatorPtr.reset (new structuralOperator_Type() );
    M_EMStructuralOperatorPtr->setup ( dataStructure,
                                       dFESpace,
                                       dETFESpace,
                                       M_bcInterfacePtr->handler(),
                                       M_commPtr);
Simone Rossi's avatar
Simone Rossi committed
601
    M_EMStructuralOperatorPtr->setDataFromGetPot (dataFile);
Simone Rossi's avatar
Simone Rossi committed
602
603
    M_EMStructuralOperatorPtr->EMMaterial()->setParameters(M_data);

thomaskummer's avatar
thomaskummer committed
604
    M_wte.setup(dataStructure, dFESpace, dETFESpace, M_commPtr, 0);
Simone Rossi's avatar
Simone Rossi committed
605
606
607
608
}

/////////////////////
// Setting up the electrophysiology solver
Simone Rossi's avatar
Simone Rossi committed
609
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
610
void
Simone Rossi's avatar
Simone Rossi committed
611
EMSolver<Mesh, ElectroSolver>::setupMechanicalBC (std::string data_file_name,
Simone Deparis's avatar
Simone Deparis committed
612
613
                                                                   std::string section,
                                                                   solidFESpacePtr_Type dFESpace)
Simone Rossi's avatar
Simone Rossi committed
614
{
Simone Deparis's avatar
Simone Deparis committed
615
616
617
618
619
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up bc interface\n";
    }
    M_bcInterfacePtr.reset (new bcInterface_Type() );
Simone Rossi's avatar
Simone Rossi committed
620
621
    M_bcInterfacePtr->createHandler();
    M_bcInterfacePtr->fillHandler ( data_file_name, "solid" );
622
    // M_bcInterfacePtr->handler()->bcUpdate ( *dFESpace->mesh(), dFESpace->feBd(), dFESpace->dof() );
Simone Rossi's avatar
Simone Rossi committed
623
624
625
626
627
}


/////////////////////
//Setup exporters
Simone Rossi's avatar
Simone Rossi committed
628
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
629
void
Simone Rossi's avatar
Simone Rossi committed
630
EMSolver<Mesh, ElectroSolver>::setupExporters (std::string problemFolder,
Simone Deparis's avatar
Simone Deparis committed
631
632
                                                                std::string electroFileName,
                                                                std::string activationFileName,
thomaskummer's avatar
thomaskummer committed
633
                                                                std::string activationTimeFileName,
thomaskummer's avatar
thomaskummer committed
634
635
                                                                std::string mechanicsFileName,
                                                                std::string vonMisesStressFileName)
Simone Rossi's avatar
Simone Rossi committed
636
{
Simone Deparis's avatar
Simone Deparis committed
637
638
639
640
    if (M_commPtr -> MyPID() == 0)
    {
        std::cout << "EMS - setting up exporters\n";
    }
thomaskummer's avatar
thomaskummer committed
641
642
    
    // Electrophysiology
Simone Deparis's avatar
Simone Deparis committed
643
644
    M_electroExporterPtr.reset (new exporter_Type() );
    setupElectroExporter (problemFolder, electroFileName);
Simone Rossi's avatar
Simone Rossi committed
645
646
647
648
649
650
    M_electroExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "fibers",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_electroSolverPtr -> fiberPtr(),
                                            UInt (0) );

thomaskummer's avatar
thomaskummer committed
651
    // Activation
Simone Deparis's avatar
Simone Deparis committed
652
653
654
655
656
    M_activationExporterPtr.reset (new exporter_Type() );
    setupActivationExporter (problemFolder, activationFileName );
    M_activationExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                             "Activation",
                                             M_electroSolverPtr -> feSpacePtr(),
Simone Rossi's avatar
Simone Rossi committed
657
                                             M_activationModelPtr -> fiberActivationPtr(),
Simone Deparis's avatar
Simone Deparis committed
658
659
660
                                             UInt (0) );
    M_mechanicsExporterPtr.reset (new exporter_Type() );

thomaskummer's avatar
thomaskummer committed
661
662
663
664
665
666
    // Activation time
    M_activationTimeExporterPtr.reset (new exporter_Type() );
    setupActivationTimeExporter (problemFolder, activationTimeFileName );
    
    M_activationTimePtr.reset (new vector_Type ( M_electroSolverPtr->potentialPtr() -> map() ) );
    *M_activationTimePtr = -1.0;
thomaskummer's avatar
thomaskummer committed
667

thomaskummer's avatar
thomaskummer committed
668
    M_activationTimeExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
thomaskummer's avatar
thomaskummer committed
669
670
671
672
673
                                            "Activation Time",
                                            M_electroSolverPtr -> feSpacePtr(),
                                            M_activationTimePtr,
                                            UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
674
675
676
677
678
679
680
681
682
683
684
    // Von Mises stress
    M_vonMisesStressExporterPtr.reset (new exporter_Type() );
    setupVonMisesStressExporter (problemFolder, vonMisesStressFileName );
    
    M_vonMisesStressPtr.reset (new vector_Type ( M_electroSolverPtr->potentialPtr() -> map() ) );
    M_vonMisesStressExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField,
                                                "Von Mises Stress",
                                                M_electroSolverPtr -> feSpacePtr(),
                                                M_wte.vonMisesStressPtr(),
                                                UInt (0) );
    
thomaskummer's avatar
thomaskummer committed
685
    // Mechanics
thomaskummer's avatar
thomaskummer committed
686
    M_mechanicsExporterPtr.reset (new exporter_Type() );
Simone Deparis's avatar
Simone Deparis committed
687
688
689
690
691
692
693
694
695
696
697
    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
698
699
700
701
702
    M_mechanicsExporterPtr -> addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
                                            "sheets",
                                            M_EMStructuralOperatorPtr -> dispFESpacePtr(),
                                            M_EMStructuralOperatorPtr -> EMMaterial() -> sheetVectorPtr(),
                                            UInt (0) );
Simone Rossi's avatar
Simone Rossi committed
703
704
}

thomaskummer's avatar
restart    
thomaskummer committed
705
706
707
708
709
710
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
711
    M_activationTimeExporterPtr -> setTimeIndex (time);
thomaskummer's avatar
restart    
thomaskummer committed
712
713
714
    M_mechanicsExporterPtr -> setTimeIndex (time);
}
    
Simone Rossi's avatar
Simone Rossi committed
715
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
716
void
thomaskummer's avatar
thomaskummer committed
717
EMSolver<Mesh, ElectroSolver>::saveSolution (Real time, const bool& restart)
Simone Rossi's avatar
Simone Rossi committed
718
{
thomaskummer's avatar
thomaskummer committed
719
720
721
    M_wte.setDisplacement ( M_EMStructuralOperatorPtr -> displacement() );
    M_wte.analyzeTensionsRecoveryVonMisesStress();
    
thomaskummer's avatar
thomaskummer committed
722
    M_electroExporterPtr -> postProcess (time);//, restart);
thomaskummer's avatar
thomaskummer committed
723
724
    //if(M_activationExporterPtr) std::cout << "\nActivation exporter available.";
    //if(M_activationModelPtr -> fiberActivationPtr()) std::cout << "\nFiber activation exporter available.";
thomaskummer's avatar
thomaskummer committed
725
    M_activationExporterPtr -> postProcess (time);//, restart );
thomaskummer's avatar
thomaskummer committed
726
    M_activationTimeExporterPtr -> postProcess (time);
thomaskummer's avatar
thomaskummer committed
727
728
    M_vonMisesStressExporterPtr -> postProcess (time);

thomaskummer's avatar
thomaskummer committed
729
    //if(M_mechanicsExporterPtr) std::cout << "\nMechanics exporter available.";
thomaskummer's avatar
thomaskummer committed
730
    M_mechanicsExporterPtr -> postProcess (time);//, restart);
Simone Rossi's avatar
Simone Rossi committed
731
732
}

Simone Rossi's avatar
Simone Rossi committed
733
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
734
void
Simone Rossi's avatar
Simone Rossi committed
735
EMSolver<Mesh, ElectroSolver>::closeExporters()
Simone Rossi's avatar
Simone Rossi committed
736
{
Simone Deparis's avatar
Simone Deparis committed
737
738
    M_electroExporterPtr -> closeFile();
    M_activationExporterPtr -> closeFile();
thomaskummer's avatar
thomaskummer committed
739
    M_activationTimeExporterPtr -> closeFile();
Simone Deparis's avatar
Simone Deparis committed
740
    M_mechanicsExporterPtr -> closeFile();
Simone Rossi's avatar
Simone Rossi committed
741
742
743
744
}


///////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
745
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
746
void
Simone Rossi's avatar
Simone Rossi committed
747
EMSolver<Mesh, ElectroSolver>::oneWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
748
{
Simone Deparis's avatar
Simone Deparis committed
749
750
    M_electroSolverPtr -> setMechanicsModifiesConductivity (false);
    M_electroSolverPtr -> displacementPtr().reset();
Simone Rossi's avatar
Simone Rossi committed
751
752
753
754
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
755
756
}

Simone Rossi's avatar
Simone Rossi committed
757
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
758
void
Simone Rossi's avatar
Simone Rossi committed
759
EMSolver<Mesh, ElectroSolver>::twoWayCoupling()
Simone Rossi's avatar
Simone Rossi committed
760
{
Simone Deparis's avatar
Simone Deparis committed
761
762
    M_electroSolverPtr -> setMechanicsModifiesConductivity (true);
    M_electroSolverPtr -> setDisplacementPtr ( M_EMStructuralOperatorPtr -> displacementPtr() );
Simone Rossi's avatar
Simone Rossi committed
763
764
765
766
    M_activationModelPtr -> setupActivationPtrs( M_EMStructuralOperatorPtr -> EMMaterial() -> fiberActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> sheetActivationPtr(),
												 M_EMStructuralOperatorPtr -> EMMaterial() -> normalActivationPtr()
												 );
Simone Rossi's avatar
Simone Rossi committed
767
768
769
}

////////////////////////////
Simone Rossi's avatar
Simone Rossi committed
770
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
771
void
Simone Rossi's avatar
Simone Rossi committed
772
EMSolver<Mesh, ElectroSolver>::solveElectrophysiology (function_Type& stimulus, Real time )
Simone Rossi's avatar
Simone Rossi committed
773
{
Simone Deparis's avatar
Simone Deparis committed
774
    setAppliedCurrent ( stimulus, time );
Simone Rossi's avatar
Simone Rossi committed
775
776
    auto v = M_electroSolverPtr -> diffusionTensor();
    std::cout << "\nEMS - " << v[0] << ", " << v[1] << ", " << v[2];
Simone Rossi's avatar
Simone Rossi committed
777
    M_electroSolverPtr -> solveOneStepGatingVariablesFE();
Thomas Kummer's avatar
Thomas Kummer committed
778
//    M_electroSolverPtr -> solveOneStepGatingVariablesRL();
Simone Rossi's avatar
Simone Rossi committed
779
    M_electroSolverPtr -> solveOneICIStep();
thomaskummer's avatar
thomaskummer committed
780
    M_electroSolverPtr -> registerActivationTime (*M_activationTimePtr, time, 0.9);
Simone Rossi's avatar
Simone Rossi committed
781
782
}

Simone Rossi's avatar
Simone Rossi committed
783
template<typename Mesh , typename ElectroSolver>
Simone Rossi's avatar
Simone Rossi committed
784
void
Simone Rossi's avatar
Simone Rossi committed
785
EMSolver<Mesh, ElectroSolver>::solveActivation (Real dt)
Simone Rossi's avatar
Simone Rossi committed
786
{
Simone Rossi's avatar
Simone Rossi committed
787
    M_activationModelPtr -> solveModel ( dt);
Simone Rossi's avatar
Simone Rossi committed
788
789
}

Simone Rossi's avatar
Simone Rossi committed
790
791
792

} // namespace LifeV

thomaskummer's avatar
thomaskummer committed
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
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
//void constructGlobalStressVector()
//{
//    //Chrono
//    LifeChrono chrono;
//    chrono.start();
//    
//    // Reset stress components
//    boost::shared_ptr<Epetra_SerialDenseMatrix> sigma ( new Epetra_SerialDenseMatrix ( M_FESpace->fieldDim(), M_FESpace->fieldDim() ) );
//    boost::shared_ptr<VectorEpetra> sigmaX ( new VectorEpetra (*M_FESpace->mapPtr() ) );
//    boost::shared_ptr<VectorEpetra> sigmaY ( new VectorEpetra (*M_FESpace->mapPtr() ) );
//    boost::shared_ptr<VectorEpetra> sigmaZ ( new VectorEpetra (*M_FESpace->mapPtr() ) );
////    *sigmaX *= 0.;
////    *sigmaY *= 0.;
////    *sigmaZ *= 0.;
//    
//    //Constructing the patch area vector for reconstruction purposes
//    VectorEpetra patchArea (*M_displacement, Unique, Add);
//    patchArea *= 0.0;
//    
//    constructPatchAreaVector ( patchArea );
//    
//    //Before assembling the reconstruction process is done
//    solutionVect_Type patchAreaR (patchArea, Repeated);
//    
//    //Compute the area of reference element
//    Real refElemArea (0);
//    for (UInt iq (0); iq < M_FESpace->qr().nbQuadPt(); ++iq)
//    {
//        refElemArea += M_FESpace->qr().weight (iq);
//    }
//    
//    //Setting the quadrature Points = DOFs of the element and weight = 1
//    Real wQuad (refElemArea / M_FESpace->refFE().nbDof() );
//    std::vector<Real> weights (M_FESpace->fe().nbFEDof(), wQuad);
//    std::vector<GeoVector> coords = M_FESpace->refFE().refCoor();
//    
//    QuadratureRule fakeQuadratureRule;
//    fakeQuadratureRule.setDimensionShape ( shapeDimension (M_FESpace->refFE().shape() ), M_FESpace->refFE().shape() );
//    fakeQuadratureRule.setPoints (coords, weights);
//    
//    //Creating a copy of the FESpace
//    feSpace_Type fakeFESpace ( M_FESpace->mesh(), M_FESpace->refFE(), M_FESpace->qr(), M_FESpace->bdQr(), 3, M_FESpace->map().commPtr() );
//    
//    //Set the new quadrature rule
//    fakeFESpace.setQuadRule (fakeQuadratureRule);
//    
//    //Preliminary variables
//    UInt totalDof = fakeFESpace.dof().numTotalDof();
//    VectorElemental dk_loc (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//    
//    //Vectors for the deformation tensor
//    (*M_deformationF).Scale (0.0);
//    std::vector<matrix_Type> vectorDeformationF (fakeFESpace.fe().nbFEDof(), *M_deformationF);
//    
//    //Copying the displacement field into a vector with repeated map for parallel computations
//    solutionVect_Type dRep (*M_displacement, Repeated);
//    
//    //Creating the local stress tensors
//    VectorElemental elVecSigmaX (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//    VectorElemental elVecSigmaY (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//    VectorElemental elVecSigmaZ (fakeFESpace.fe().nbFEDof(), fakeFESpace.fieldDim() );
//    
//    //Loop on each volume
//    for ( UInt i (0); i < fakeFESpace.mesh()->numVolumes(); ++i )
//    {
//        //fakeFESpace.fe().update ( fakeFESpace.mesh()->volumeList ( i ), UPDATE_DPHI | UPDATE_WDET );
//        fakeFESpace.fe().updateFirstDerivQuadPt ( fakeFESpace.mesh()->volumeList ( i ) );
//        
//        elVecSigmaX.zero();
//        elVecSigmaY.zero();
//        elVecSigmaZ.zero();
//        
//        M_marker = fakeFESpace.mesh()->volumeList ( i ).markerID();
//        
//        UInt eleID = fakeFESpace.fe().currentLocalId();
//        
//        //Extracting the local displacement
//        for ( UInt iNode (0); iNode < ( UInt ) fakeFESpace.fe().nbFEDof(); ++iNode )
//        {
//            UInt  iloc = fakeFESpace.fe().patternFirst ( iNode );
//            
//            for ( UInt iComp = 0; iComp < fakeFESpace.fieldDim(); ++iComp )
//            {
//                UInt ig = fakeFESpace.dof().localToGlobalMap ( eleID, iloc ) + iComp * fakeFESpace.dim() + M_offset;
//                dk_loc[iloc + iComp * fakeFESpace.fe().nbFEDof() ] = dRep[ig];
//            }
//        }
//        
//        //Compute the element tensor F
//        AssemblyElementalStructure::computeLocalDeformationGradient ( dk_loc, vectorDeformationF, fakeFESpace.fe() );
//        
//        //Compute the local vector of the principal stresses
//        for ( UInt nDOF (0); nDOF < ( UInt ) fakeFESpace.fe().nbFEDof(); ++nDOF )
//        {
//            UInt  iloc = fakeFESpace.fe().patternFirst ( nDOF );
//            
//            sigma->Scale (0.0);
//            M_firstPiola->Scale (0.0);
//            M_cofactorF->Scale (0.0);
//            
//            //Compute the rightCauchyC tensor
//            AssemblyElementalStructure::computeInvariantsRightCauchyGreenTensor (M_invariants, vectorDeformationF[nDOF], *M_cofactorF);
//            
//            //Compute the first Piola-Kirchhoff tensor
//            M_material->computeLocalFirstPiolaKirchhoffTensor (*M_firstPiola, vectorDeformationF[nDOF], *M_cofactorF, M_invariants, M_marker);
//            
//            //Compute the Cauchy tensor
//            AssemblyElementalStructure::computeCauchyStressTensor (*sigma, *M_firstPiola, M_invariants[3], vectorDeformationF[nDOF]);
//            
//            //Assembling the local vectors for local tensions Component X, Y, Z
//            for ( UInt coor (0); coor < fakeFESpace.fieldDim(); ++coor )
//            {
//                (elVecSigmaX) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 0);
//                (elVecSigmaY) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 1);
//                (elVecSigmaZ) [iloc + coor * fakeFESpace.fe().nbFEDof() ] = (*sigma) (coor, 2);
//            }
//        }
//        
//        reconstructElementaryVector ( elVecSigmaX, patchAreaR, fakeFESpace );
//        reconstructElementaryVector ( elVecSigmaY, patchAreaR, fakeFESpace );
//        reconstructElementaryVector ( elVecSigmaZ, patchAreaR, fakeFESpace );
//        
//        //Assembling the three elemental vector in the three global
//        for ( UInt ic = 0; ic < fakeFESpace.fieldDim(); ++ic )
//        {
//            assembleVector (*sigmaX, elVecSigmaX, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset +  ic * totalDof );
//            assembleVector (*sigmaY, elVecSigmaY, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset +  ic * totalDof );
//            assembleVector (*sigmaZ, elVecSigmaZ, fakeFESpace.fe(), fakeFESpace.dof(), ic, M_offset +  ic * totalDof );
//        }
//    }
//    
//    sigmaX->globalAssemble();
//    sigmaY->globalAssemble();
//    sigmaZ->globalAssemble();
//    
//    //Chrono
//    chrono.stop();
//    M_displayer->leaderPrint ("  S-  Cauchy stresses recovered in:             ", chrono.globalDiff ( *M_displayer->comm() ), " s\n" );
//}
Simone Rossi's avatar
Simone Rossi committed
932
933
934



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