HeartSolver.hpp 13 KB
Newer Older
thomaskummer's avatar
thomaskummer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
//
//  HeartSolver.cpp
//  LifeV
//
//  Created by Thomas Kummer on 03.05.16.
//  Copyright © 2016 Thomas Kummer. All rights reserved.
//

#ifndef _HEARTDATA_H_
#define _HEARTDATA_H_


#include <stdio.h>
thomaskummer's avatar
thomaskummer committed
14
#include <lifev/em/solver/EMSolver.hpp>
thomaskummer's avatar
thomaskummer committed
15
#include <lifev/em/solver/circulation/Circulation.hpp>
thomaskummer's avatar
thomaskummer committed
16

thomaskummer's avatar
thomaskummer committed
17
#include <lifev/em/solver/HeartData.hpp>
thomaskummer's avatar
thomaskummer committed
18

thomaskummer's avatar
thomaskummer committed
19
20
21
22

namespace LifeV
{

thomaskummer's avatar
thomaskummer committed
23
//#define PI 3.14159265359
thomaskummer's avatar
thomaskummer committed
24

thomaskummer's avatar
thomaskummer committed
25
    
thomaskummer's avatar
thomaskummer committed
26
27
28
29
template <class EmSolver>
class HeartSolver {
   
public:
thomaskummer's avatar
thomaskummer committed
30
    
thomaskummer's avatar
thomaskummer committed
31
32
    HeartSolver(EmSolver& emSolver,  Circulation& circulationSolver) :
        M_emSolver          (emSolver),
thomaskummer's avatar
thomaskummer committed
33
34
        M_circulationSolver (circulationSolver),
        M_heartData         (HeartData())
thomaskummer's avatar
thomaskummer committed
35
    {}
thomaskummer's avatar
thomaskummer committed
36
    
thomaskummer's avatar
thomaskummer committed
37
    virtual ~HeartSolver() {}
thomaskummer's avatar
thomaskummer committed
38
39
40
41
42
43
44
45
46
47

    EmSolver& emSolver()
    {
        return M_emSolver;
    }
    
    Circulation& circulation()
    {
        return M_heartData;
    }
thomaskummer's avatar
thomaskummer committed
48
    
thomaskummer's avatar
thomaskummer committed
49
    const HeartData& data() const
thomaskummer's avatar
thomaskummer committed
50
51
52
53
    {
        return M_heartData;
    }
    
thomaskummer's avatar
thomaskummer committed
54
    void setup(const GetPot& datafile)
thomaskummer's avatar
thomaskummer committed
55
    {
thomaskummer's avatar
thomaskummer committed
56
        M_heartData.setup(datafile);
thomaskummer's avatar
thomaskummer committed
57
    }
thomaskummer's avatar
thomaskummer committed
58
    
thomaskummer's avatar
thomaskummer committed
59
    template <class lambda>
thomaskummer's avatar
thomaskummer committed
60
    void preload(const lambda& modifyFeBC, const std::vector<Real>& bcValues)
thomaskummer's avatar
thomaskummer committed
61
    {
thomaskummer's avatar
thomaskummer committed
62
        M_emSolver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(0.0);
thomaskummer's avatar
thomaskummer committed
63
64
65
66
67
68
69
70
71
72
        
        auto preloadPressure = [] (std::vector<double> p, const int& step, const int& steps)
        {
            for (auto& i : p) {i *= double(step) / double(steps);}
            return p;
        };
        
        LifeChrono chronoSave;
        chronoSave.start();
        
thomaskummer's avatar
thomaskummer committed
73
        M_emSolver.saveSolution (-1.0);
thomaskummer's avatar
thomaskummer committed
74
        
thomaskummer's avatar
thomaskummer committed
75
        if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
76
77
78
79
80
81
82
83
84
85
86
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nData stored in " << chronoSave.diff() << " s";
            std::cout << "\n*****************************************************************\n";
        }
        
        LifeChrono chronoPreload;
        chronoPreload.start();
        
        for (int i (1); i <= data().preloadSteps(); i++)
        {
thomaskummer's avatar
thomaskummer committed
87
            if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
88
89
90
91
92
93
94
            {
                std::cout << "\n*****************************************************************";
                std::cout << "\nPreload step: " << i << " / " << data().preloadSteps();
                std::cout << "\n*****************************************************************\n";
            }
            
            // Update pressure b.c.
thomaskummer's avatar
thomaskummer committed
95
            modifyFeBC(preloadPressure(bcValues, i, data().preloadSteps() ));
thomaskummer's avatar
thomaskummer committed
96
97
            
            // Solve mechanics
thomaskummer's avatar
thomaskummer committed
98
99
            M_emSolver.bcInterfacePtr() -> updatePhysicalSolverVariables();
            M_emSolver.solveMechanics();
thomaskummer's avatar
thomaskummer committed
100
101
            
            // Safe preload steps
thomaskummer's avatar
thomaskummer committed
102
            if ( data().safePreload() ) M_emSolver.saveSolution (i-1);
thomaskummer's avatar
thomaskummer committed
103
104
        }
        
thomaskummer's avatar
thomaskummer committed
105
        if ( 0 == M_emSolver.comm()->MyPID() )
thomaskummer's avatar
thomaskummer committed
106
107
108
109
110
111
112
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nPreload done in: " << chronoPreload.diff();
            std::cout << "\n*****************************************************************\n";
        }

    }
thomaskummer's avatar
thomaskummer committed
113
114
    
    
thomaskummer's avatar
thomaskummer committed
115
    void restart(std::string& restartInput, const GetPot& command_line, Real& t)
thomaskummer's avatar
thomaskummer committed
116
    {
thomaskummer's avatar
thomaskummer committed
117
        const std::string restartDir = ""; //command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
thomaskummer's avatar
thomaskummer committed
118
119
120
121
122
123
124
125
126
127
128
        
        Real dtExport = 10.;
        
        // Set time variable
        const unsigned int restartInputStr = std::stoi(restartInput);
        const unsigned int nIter = (restartInputStr - 1) * dtExport / data().dt_mechanics();
        t = nIter * data().dt_mechanics();
        
        // Set time exporter time index
        M_emSolver.setTimeIndex(restartInputStr + 1);
        //solver.importHdf5();
thomaskummer's avatar
thomaskummer committed
129

thomaskummer's avatar
thomaskummer committed
130
        // Load restart solutions from output files
thomaskummer's avatar
thomaskummer committed
131
        std::string polynomialDegree = data().elementOrder();
thomaskummer's avatar
thomaskummer committed
132

thomaskummer's avatar
thomaskummer committed
133
        ElectrophysiologyUtility::importVectorField ( M_emSolver.structuralOperatorPtr() -> displacementPtr(), "MechanicalSolution" , "displacement", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
thomaskummer's avatar
thomaskummer committed
134
135
136
        ElectrophysiologyUtility::importScalarField (M_emSolver.activationModelPtr() -> fiberActivationPtr(), "ActivationSolution" , "Activation", M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        //ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "ActivationTimeSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        
thomaskummer's avatar
thomaskummer committed
137
138
139
140
        for ( unsigned int i = 0; i < M_emSolver.electroSolverPtr()->globalSolution().size() ; ++i )
        {
            ElectrophysiologyUtility::importScalarField (M_emSolver.electroSolverPtr()->globalSolution().at(i), "ElectroSolution" , ("Variable" + std::to_string(i)), M_emSolver.localMeshPtr(), restartDir, polynomialDegree, restartInput );
        }
thomaskummer's avatar
thomaskummer committed
141

thomaskummer's avatar
thomaskummer committed
142
143
144
145
        if ( 0 == M_emSolver.comm()->MyPID() )
        {
            std::cout << "\nLoad from restart: " << restartInput << ",  nIterCirculation = " << nIter << ",  time = " << t << std::endl;
        }
thomaskummer's avatar
thomaskummer committed
146
        
thomaskummer's avatar
thomaskummer committed
147
        M_circulationSolver.restartFromFile ( restartDir + "solution.dat" , nIter );
thomaskummer's avatar
thomaskummer committed
148
149
    }
    
thomaskummer's avatar
thomaskummer committed
150
    
thomaskummer's avatar
thomaskummer committed
151
    void createPatch (EMSolver<RegionMesh<LinearTetra>, EMMonodomainSolver<RegionMesh<LinearTetra> > >& solver, const Vector3D& center, const Real& radius, const int& currentFlag, const int& newFlag)
thomaskummer's avatar
thomaskummer committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
    {
        for (auto& mesh : solver.mesh())
        {
            for (int j(0); j < mesh->numBoundaryFacets(); j++)
            {
                auto& face = mesh->boundaryFacet(j);
                auto faceFlag = face.markerID();
                
                if (faceFlag == currentFlag || faceFlag == 470 || faceFlag == 471)
                {
                    int numPointsInsidePatch (0);
                    
                    for (int k(0); k < 3; ++k)
                    {
                        auto coord = face.point(k).coordinates();
                        bool pointInPatch = (coord - center).norm() < radius;
                        
                        if (pointInPatch)
                        {
                            ++numPointsInsidePatch;
                        }
                    }
                    
                    if (numPointsInsidePatch > 2)
                    {
                        face.setMarkerID(newFlag);
                    }
                    
                }
            }
        }
    }
    
thomaskummer's avatar
thomaskummer committed
185
186
187
188
189
    boost::shared_ptr<VectorEpetra> directionalVectorField (const boost::shared_ptr<FESpace<RegionMesh<LinearTetra>, MapEpetra >> dFeSpace, Vector3D& direction, const Real& disp)
    {
        boost::shared_ptr<VectorEpetra> vectorField (new VectorEpetra( dFeSpace->map(), Repeated ));
        auto nCompLocalDof = vectorField->epetraVector().MyLength() / 3;
        
thomaskummer's avatar
thomaskummer committed
190
        direction.normalize();
thomaskummer's avatar
thomaskummer committed
191
192
193
194
195
196
197
        direction *= disp;
        
        for (int j (0); j < nCompLocalDof; ++j)
        {
            UInt iGID = vectorField->blockMap().GID (j);
            UInt jGID = vectorField->blockMap().GID (j + nCompLocalDof);
            UInt kGID = vectorField->blockMap().GID (j + 2 * nCompLocalDof);
thomaskummer's avatar
thomaskummer committed
198

thomaskummer's avatar
thomaskummer committed
199
200
201
202
203
204
205
206
            (*vectorField)[iGID] = direction[0];
            (*vectorField)[jGID] = direction[1];
            (*vectorField)[kGID] = direction[2];
        }
        
        return vectorField;
    }
    
thomaskummer's avatar
thomaskummer committed
207
    static Real Iapp (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/)
thomaskummer's avatar
thomaskummer committed
208
209
210
211
212
213
214
    {
        bool coords ( Y < -7. );
        //bool coords ( Y > 4. ); //( Y > 1.5 && Y < 3 );
        bool time ( fmod(t, 800.) < 4 && fmod(t, 800.) > 2);
        return ( coords && time ? 30 : 0 );
    }
    
thomaskummer's avatar
thomaskummer committed
215
    Real sinSquared (const Real& time, const Real& Tmax, const Real& tmax, const Real& tduration)
thomaskummer's avatar
thomaskummer committed
216
    {
thomaskummer's avatar
thomaskummer committed
217
        Real timeInPeriod = fmod(time-tmax+0.5*tduration, 800.);
thomaskummer's avatar
thomaskummer committed
218
        bool inPeriod ( timeInPeriod < tduration && timeInPeriod > 0);
thomaskummer's avatar
thomaskummer committed
219
        Real sinusSquared = std::pow( std::sin(timeInPeriod*3.14159265359/tduration) , 2 ) * Tmax;
thomaskummer's avatar
thomaskummer committed
220
        return ( inPeriod ? sinusSquared : 0 );
thomaskummer's avatar
thomaskummer committed
221
222
    }
    
thomaskummer's avatar
thomaskummer committed
223
    template<class bcVectorType>
thomaskummer's avatar
thomaskummer committed
224
    void extrapolate4thOrderAdamBashforth(bcVectorType& bcValues, bcVectorType& bcValuesPre, const Real& dpMax)
thomaskummer's avatar
thomaskummer committed
225
226
227
228
229
230
    {
        VectorSmall<4> ABcoef;
        ABcoef (0) = 55/24; ABcoef (1) = -59/24; ABcoef (2) = 37/24; ABcoef (3) = -3/8;
        
        for ( unsigned int i = ABcoef.size() - 1; i > 0; --i )
        {
thomaskummer's avatar
thomaskummer committed
231
232
            m_ABdplv(i) = m_ABdplv(i-1);
            m_ABdprv(i) = m_ABdprv(i-1);
thomaskummer's avatar
thomaskummer committed
233
234
        }
        
thomaskummer's avatar
thomaskummer committed
235
236
        m_ABdplv(0) = bcValues[0] - bcValuesPre[0];
        m_ABdprv(0) = bcValues[1] - bcValuesPre[1];
thomaskummer's avatar
thomaskummer committed
237
238
239
        
        bcValuesPre = bcValues;
        
thomaskummer's avatar
thomaskummer committed
240
241
        bcValues[0] += std::min( std::max( ABcoef.dot( m_ABdplv ) , - dpMax ) , dpMax );
        bcValues[1] += std::min( std::max( ABcoef.dot( m_ABdprv ) , - dpMax ) , dpMax );
thomaskummer's avatar
thomaskummer committed
242
    }
thomaskummer's avatar
thomaskummer committed
243
    
thomaskummer's avatar
thomaskummer committed
244
    
thomaskummer's avatar
thomaskummer committed
245
protected:
thomaskummer's avatar
thomaskummer committed
246
    
thomaskummer's avatar
thomaskummer committed
247
    
thomaskummer's avatar
thomaskummer committed
248
249
    EmSolver& M_emSolver;
    Circulation& M_circulationSolver;
thomaskummer's avatar
thomaskummer committed
250
    
thomaskummer's avatar
thomaskummer committed
251
    HeartData M_heartData;
thomaskummer's avatar
thomaskummer committed
252
    
thomaskummer's avatar
thomaskummer committed
253
    
thomaskummer's avatar
thomaskummer committed
254
255
256
    VectorSmall<2> M_pressure;
    VectorSmall<2> M_volume;

thomaskummer's avatar
thomaskummer committed
257
    VectorSmall<4> m_ABdplv, m_ABdprv;
thomaskummer's avatar
thomaskummer committed
258

thomaskummer's avatar
thomaskummer committed
259
    
thomaskummer's avatar
thomaskummer committed
260
    
thomaskummer's avatar
thomaskummer committed
261
262
263
    
    
    
thomaskummer's avatar
thomaskummer committed
264
265
266
267
268
269
270
271
272
273
274
    std::string pipeToString ( const char* command )
    {
        FILE* file = popen( command, "r" ) ;
        std::ostringstream stm ;
        char line[6] ;
        fgets( line, 6, file );
        stm << line;
        pclose(file) ;
        return stm.str() ;
    };
    
thomaskummer's avatar
thomaskummer committed
275
276
277


    Real patchDispFunNormal (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& i)
thomaskummer's avatar
thomaskummer committed
278
    {
thomaskummer's avatar
thomaskummer committed
279
        return (-0.000 - 0.00005*t);// sinSquared(t, 0.1, 50, 100)); // -0.001;// (t * 1e-5);
thomaskummer's avatar
thomaskummer committed
280
281
    }
    
thomaskummer's avatar
thomaskummer committed
282
    Real patchFunction (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/)
thomaskummer's avatar
thomaskummer committed
283
    {
thomaskummer's avatar
thomaskummer committed
284
        Real disp = std::pow( std::sin(fmod(t, 800.) * 3.14159265359/300) , 2 )*15;
thomaskummer's avatar
thomaskummer committed
285
286
287
        return disp;
    }
    
thomaskummer's avatar
thomaskummer committed
288
    Real potentialMultiplyerFcn (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& /*i*/)
thomaskummer's avatar
thomaskummer committed
289
    {
thomaskummer's avatar
thomaskummer committed
290
291
        bool time ( fmod(t, 800.) < 4 && fmod(t, 800.) > 2);
        return 1.4 * time; // ( Y < 2.5 && Y > 0.5 ? 1.0 : 0.0 );
thomaskummer's avatar
thomaskummer committed
292
293
    }
    
thomaskummer's avatar
thomaskummer committed
294
    Real patchDispFun1 (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& i)
thomaskummer's avatar
thomaskummer committed
295
    {
thomaskummer's avatar
thomaskummer committed
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
        switch (i)
        {
            case 0:
                return (t);
                break;
            case 1:
                return 0;
                break;
            case 2:
                return (t);
                break;
            default:
                ERROR_MSG ("This entry is not allowed");
                return 0.;
                break;
        }
    }
    
    Real patchDispFun2 (const Real& t, const Real&  X, const Real& Y, const Real& Z, const ID& i)
    {
        switch (i)
        {
            case 0:
                return (-t);
                break;
            case 1:
                return 0;
                break;
            case 2:
                return (-t);
                break;
            default:
                ERROR_MSG ("This entry is not allowed");
                return 0.;
                break;
        }
    }
    
    Real normalDirection ( const Real& /*t*/, const Real& x , const Real& y, const Real& z, const ID& i)
    {
        Real nnorm = std::sqrt(x * x + y * y + z * z);
        switch (i)
        {
            case 0:
                return 1/std::sqrt(2);
                break;
            case 1:
                return 0;
                break;
            case 2:
                return 1/std::sqrt(2);
                break;
            default:
                ERROR_MSG ("This entry is not allowed");
                return 0.;
                break;
        }
thomaskummer's avatar
thomaskummer committed
353
    }
thomaskummer's avatar
thomaskummer committed
354
355

    
thomaskummer's avatar
thomaskummer committed
356
    
thomaskummer's avatar
thomaskummer committed
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
    Real
    externalPower ( const VectorEpetra& dispCurrent,
                    const VectorEpetra& dispPrevious,
                    const boost::shared_ptr<ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 > > dispETFESpace,
                    Real pressure,
                    Real dt,
                    const unsigned int bdFlag) const
    {
        VectorEpetra traction ( dispCurrent.map() );
        VectorEpetra velocity ( (dispCurrent - dispPrevious) / dt );
        
        MatrixSmall<3,3> Id;
        Id(0,0) = 1.; Id(0,1) = 0., Id(0,2) = 0.;
        Id(1,0) = 0.; Id(1,1) = 1., Id(1,2) = 0.;
        Id(2,0) = 0.; Id(2,1) = 0., Id(2,2) = 1.;
        
        {
            using namespace ExpressionAssembly;
            
            auto I = value(Id);
            auto Grad_u = grad( dispETFESpace, dispCurrent, 0);
            auto F =  Grad_u + I;
            auto FmT = minusT(F);
            auto J = det(F);
            auto p = value(pressure);
            
            QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria7pt) );
            
            integrate ( boundary ( dispETFESpace->mesh(), bdFlag),
                       myBDQR,
                       dispETFESpace,
                       p * J * dot( FmT * Nface,  phi_i)
                       //p * J * dot( FmT * Nface,  phi_i)
                       //value(-1.0) * J * dot (vE1, FmT * Nface) * phi_i) >> intergral
                       ) >> traction;
            
            traction.globalAssemble();
        }
        
        return traction.dot(velocity);
    }
    
    
thomaskummer's avatar
thomaskummer committed
400
401
402
403
};

}

thomaskummer's avatar
thomaskummer committed
404
#endif