main.cpp 41.6 KB
Newer Older
thomaskummer's avatar
thomaskummer committed
1
2
3
4
5
6
7
8
//=======================================================
//  Simulating the human heart
//
//  Created by Thomas Kummer on 30.04.18.
//  Copyright © 2018 Thomas Kummer. All rights reserved.
//=======================================================


thomaskummer's avatar
thomaskummer committed
9
//============================================
thomaskummer's avatar
thomaskummer committed
10
// Includes
thomaskummer's avatar
thomaskummer committed
11
//============================================
thomaskummer's avatar
thomaskummer committed
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

#include <lifev/core/LifeV.hpp>

// Passive material
#include <lifev/electrophysiology/solver/ElectroETAMonodomainSolver.hpp>
#include <lifev/electrophysiology/util/HeartUtility.hpp>
#include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>
#include <lifev/em/solver/mechanics/EMStructuralOperator.hpp>
#include <lifev/em/solver/mechanics/EMStructuralConstitutiveLaw.hpp>
#include <lifev/bc_interface/3D/bc/BCInterface3D.hpp>

// Solver
#include <lifev/em/solver/EMSolver.hpp>

// Exporter
#include <lifev/core/filter/ExporterEnsight.hpp>
#include <lifev/core/filter/ExporterVTK.hpp>
#ifdef HAVE_HDF5
#include <lifev/core/filter/ExporterHDF5.hpp>
#endif
#include <lifev/core/filter/ExporterEmpty.hpp>

// Resize mesh
#include <lifev/core/mesh/MeshUtility.hpp>

// B.C. modification
#include <lifev/core/fem/BCVector.hpp>

// Circulation
#include <lifev/em/solver/circulation/Circulation.hpp>

// Volume computation
#include <lifev/em/solver/circulation/CirculationVolumeIntegrator.hpp>

thomaskummer's avatar
thomaskummer committed
46
47
// Heart solver
#include <lifev/em/solver/HeartSolver.hpp>
thomaskummer's avatar
thomaskummer committed
48

thomaskummer's avatar
thomaskummer committed
49
// PatchBC
thomaskummer's avatar
thomaskummer committed
50
#include <lifev/em/examples/example_EMHeart/EssentialPatchBCRectangle.hpp>
thomaskummer's avatar
thomaskummer committed
51
52
#include <lifev/em/examples/example_EMHeart/EssentialPatchBCCircular.hpp>
#include <lifev/em/examples/example_EMHeart/EssentialPatchBCCircularSmooth.hpp>
thomaskummer's avatar
thomaskummer committed
53
#include <lifev/em/examples/example_EMHeart/EssentialPatchBCEllipseSmooth.hpp>
thomaskummer's avatar
thomaskummer committed
54
#include <lifev/em/examples/example_EMHeart/EssentialPatchBCShell.hpp>
thomaskummer's avatar
thomaskummer committed
55

thomaskummer's avatar
thomaskummer committed
56
// Track nan
thomaskummer's avatar
thomaskummer committed
57
// #include <fenv.h>
thomaskummer's avatar
thomaskummer committed
58
59
60
61
62



using namespace LifeV;

thomaskummer's avatar
thomaskummer committed
63
#define mmHg * 1333.224
thomaskummer's avatar
thomaskummer committed
64
#define Pa * 10.0
thomaskummer's avatar
thomaskummer committed
65
66
67
68

int main (int argc, char** argv)
{

thomaskummer's avatar
thomaskummer committed
69
//    feenableexcept(FE_INVALID | FE_OVERFLOW);
thomaskummer's avatar
thomaskummer committed
70
    
thomaskummer's avatar
thomaskummer committed
71
    //============================================
thomaskummer's avatar
thomaskummer committed
72
    // Typedefs
thomaskummer's avatar
thomaskummer committed
73
    //============================================
thomaskummer's avatar
thomaskummer committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
    
    typedef RegionMesh<LinearTetra>                         mesh_Type;
    typedef boost::shared_ptr<mesh_Type>                    meshPtr_Type;
    
    typedef boost::function < Real (const Real & t,
                                    const Real &   x,
                                    const Real &   y,
                                    const Real & z,
                                    const ID&   /*i*/ ) >   function_Type;
    
    typedef VectorEpetra                                    vector_Type;
    typedef boost::shared_ptr<vector_Type>                  vectorPtr_Type;
    
    typedef BCVector                                        bcVector_Type;
    typedef boost::shared_ptr<bcVector_Type>                bcVectorPtr_Type;
    
    typedef EMMonodomainSolver<mesh_Type>                   monodomain_Type;

    typedef FESpace< mesh_Type, MapEpetra >                 solidFESpace_Type;
    typedef boost::shared_ptr<solidFESpace_Type>            solidFESpacePtr_Type;
    
    typedef ETFESpace< mesh_Type, MapEpetra, 3, 3 >         solidETFESpace_Type;
    typedef boost::shared_ptr<solidETFESpace_Type>          solidETFESpacePtr_Type;
    
    typedef StructuralConstitutiveLawData                   constitutiveLaw_Type;
    typedef boost::shared_ptr<constitutiveLaw_Type>         constitutiveLawPtr_Type;
    
    typedef BCHandler                                       bc_Type;
    typedef StructuralOperator< mesh_Type >                 physicalSolver_Type;
    typedef BCInterface3D< bc_Type, physicalSolver_Type >   bcInterface_Type;
    typedef boost::shared_ptr< bcInterface_Type >           bcInterfacePtr_Type;
    

thomaskummer's avatar
thomaskummer committed
107
    //============================================
thomaskummer's avatar
thomaskummer committed
108
    // Communicator and displayer
thomaskummer's avatar
thomaskummer committed
109
    //============================================
thomaskummer's avatar
thomaskummer committed
110
111
112
113
114
115
#ifdef HAVE_MPI
    MPI_Init ( &argc, &argv );
#endif

    boost::shared_ptr<Epetra_Comm>  comm ( new Epetra_MpiComm (MPI_COMM_WORLD) );
    Displayer displayer ( comm );
thomaskummer's avatar
thomaskummer committed
116
    
thomaskummer's avatar
thomaskummer committed
117
    displayer.leaderPrint("\nEMHeart running ...\n");
thomaskummer's avatar
thomaskummer committed
118
    
thomaskummer's avatar
thomaskummer committed
119
    //============================================
thomaskummer's avatar
thomaskummer committed
120
    // Read data file and create output folder
thomaskummer's avatar
thomaskummer committed
121
    //============================================
thomaskummer's avatar
thomaskummer committed
122
123
124
125
126
    GetPot command_line (argc, argv);
    const std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
    GetPot dataFile (data_file_name);
    std::string problemFolder = EMUtility::createOutputFolder (command_line, *comm);

thomaskummer's avatar
thomaskummer committed
127
    
thomaskummer's avatar
thomaskummer committed
128
    //============================================
thomaskummer's avatar
thomaskummer committed
129
    // Electromechanic solver
thomaskummer's avatar
thomaskummer committed
130
    //============================================
thomaskummer's avatar
thomaskummer committed
131
132
133
    EMSolver<mesh_Type, monodomain_Type> solver(comm);
    
    
thomaskummer's avatar
thomaskummer committed
134
    //============================================
thomaskummer's avatar
thomaskummer committed
135
    // Body circulation
thomaskummer's avatar
thomaskummer committed
136
    //============================================
thomaskummer's avatar
thomaskummer committed
137
138
139
140
141
142
143
144
145
146
    const std::string circulationInputFile = command_line.follow ("circulation", 2, "-cif", "--cifile");
    const std::string circulationOutputFile = command_line.follow ( (problemFolder + "solution.dat").c_str(), 2, "-cof", "--cofile");
    
    Circulation circulationSolver( circulationInputFile );
    
    // Flow rate between two vertices
    auto Q = [&circulationSolver] (const std::string& N1, const std::string& N2) { return circulationSolver.solution ( std::vector<std::string> {N1, N2} ); };
    auto p = [&circulationSolver] (const std::string& N1) { return circulationSolver.solution ( N1 ); };
    
    
thomaskummer's avatar
thomaskummer committed
147
    //============================================
thomaskummer's avatar
thomaskummer committed
148
    // Heart solver
thomaskummer's avatar
thomaskummer committed
149
    //============================================
thomaskummer's avatar
thomaskummer committed
150
    HeartSolver<EMSolver<mesh_Type, monodomain_Type> > heartSolver (solver, circulationSolver);
thomaskummer's avatar
thomaskummer committed
151
    
thomaskummer's avatar
thomaskummer committed
152
    
thomaskummer's avatar
thomaskummer committed
153
    //============================================
thomaskummer's avatar
thomaskummer committed
154
    // Setup material data
thomaskummer's avatar
thomaskummer committed
155
    //============================================
thomaskummer's avatar
thomaskummer committed
156
157
    EMData emdata;
    emdata.setup (dataFile);
thomaskummer's avatar
thomaskummer committed
158

thomaskummer's avatar
thomaskummer committed
159
    
thomaskummer's avatar
thomaskummer committed
160
    //============================================
thomaskummer's avatar
thomaskummer committed
161
    // Load mesh
thomaskummer's avatar
thomaskummer committed
162
    //============================================
thomaskummer's avatar
thomaskummer committed
163
    std::string meshType = dataFile("solid/space_discretization/mesh_type", ".mesh");
thomaskummer's avatar
thomaskummer committed
164
    std::string meshName = dataFile("solid/space_discretization/mesh_name", "cube4");
thomaskummer's avatar
thomaskummer committed
165
    std::string meshPath = dataFile("solid/space_discretization/mesh_dir", "mesh/humanHeart/") + meshName + "/";
thomaskummer's avatar
thomaskummer committed
166
    
thomaskummer's avatar
thomaskummer committed
167
    solver.loadMesh (meshName + meshType, meshPath);
thomaskummer's avatar
thomaskummer committed
168
169
    
    
thomaskummer's avatar
thomaskummer committed
170
    //============================================
thomaskummer's avatar
thomaskummer committed
171
    // Resize mesh
thomaskummer's avatar
thomaskummer committed
172
    //============================================
thomaskummer's avatar
thomaskummer committed
173
174
175
    Vector3D scale, rotate, translate;
    for ( UInt j (0); j < 3; ++j )
    {
thomaskummer's avatar
thomaskummer committed
176
177
178
        scale[j] = dataFile ( "solid/space_discretization/mesh_scaling", 1., j );
        rotate[j] = dataFile ( "solid/space_discretization/mesh_rotation", 0., j );
        translate[j] = dataFile ( "solid/space_discretization/mesh_translation", 0., j );
thomaskummer's avatar
thomaskummer committed
179
    }
thomaskummer's avatar
thomaskummer committed
180
181
182
    
    MeshUtility::MeshTransformer<mesh_Type> transformerFull (* (solver.fullMeshPtr() ) );
    MeshUtility::MeshTransformer<mesh_Type> transformerLocal (* (solver.localMeshPtr() ) );
thomaskummer's avatar
thomaskummer committed
183

thomaskummer's avatar
thomaskummer committed
184
185
186
    transformerFull.transformMesh (scale, rotate, translate);
    transformerLocal.transformMesh (scale, rotate, translate);
    
thomaskummer's avatar
thomaskummer committed
187
    if ( 0 == comm->MyPID() ) std::cout << "\nResizing mesh done" << std::endl;
thomaskummer's avatar
thomaskummer committed
188
    if ( 0 == comm->MyPID() ) solver.fullMeshPtr()->showMe();
thomaskummer's avatar
thomaskummer committed
189
190
    
    
thomaskummer's avatar
thomaskummer committed
191
    //============================================
thomaskummer's avatar
thomaskummer committed
192
    // Setup solver (including fe-spaces & b.c.)
thomaskummer's avatar
thomaskummer committed
193
    //============================================
thomaskummer's avatar
thomaskummer committed
194
    EMAssembler::quadRule.setQuadRule( dataFile ( "solid/space_discretization/quad_rule", "4pt") );
thomaskummer's avatar
thomaskummer committed
195
196
    solver.setup (dataFile);
    
thomaskummer's avatar
thomaskummer committed
197
198
199
    auto& disp = solver.structuralOperatorPtr() -> displacement();
    auto FESpace = solver.structuralOperatorPtr() -> dispFESpacePtr();
    auto dETFESpace = solver.electroSolverPtr() -> displacementETFESpacePtr();
thomaskummer's avatar
thomaskummer committed
200
    auto ETFESpace = solver.electroSolverPtr() -> ETFESpacePtr();
thomaskummer's avatar
thomaskummer committed
201
    
thomaskummer's avatar
thomaskummer committed
202
203
204
205
206
207
208
209
210
    
    //============================================
    // Create essential patch b.c.
    //============================================
    EssentialPatchBCHandler patchHandler ("listEssentialPatchBC", dataFile);
    patchHandler.addPatchBC(solver);
    
    if ( 0 == comm->MyPID() ) PRINT_FACTORY(EssentialPatchBC);
    
thomaskummer's avatar
thomaskummer committed
211

thomaskummer's avatar
thomaskummer committed
212
    //============================================
thomaskummer's avatar
thomaskummer committed
213
    // Setup anisotropy vectors
thomaskummer's avatar
thomaskummer committed
214
    //============================================
thomaskummer's avatar
thomaskummer committed
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
    bool anisotropy = dataFile ( "solid/space_discretization/anisotropic", false );

    if ( anisotropy )
    {
        std::string fiberFileName  =  dataFile ( "solid/space_discretization/fiber_name", "FiberDirection");
        std::string sheetFileName  =  dataFile ( "solid/space_discretization/sheet_name", "SheetsDirection");
        std::string fiberFieldName =  dataFile ( "solid/space_discretization/fiber_fieldname", "fibers");
        std::string sheetFieldName =  dataFile ( "solid/space_discretization/sheet_fieldname", "sheets");
        std::string fiberDir       =  meshPath; //dataFile ( "solid/space_discretization/fiber_dir", "./");
        std::string sheetDir       =  meshPath; //dataFile ( "solid/space_discretization/sheet_dir", "./");
        std::string elementOrder   =  dataFile ( "solid/space_discretization/order", "P1");

        solver.setupFiberVector ( fiberFileName, fiberFieldName, fiberDir, elementOrder );
        solver.setupMechanicalSheetVector ( sheetFileName, sheetFieldName, sheetDir, elementOrder );
    }
    else
    {
        solver.setupFiberVector (1., 0., 0.);
        solver.setupSheetVector (0., 1., 0.);
    }
    
thomaskummer's avatar
thomaskummer committed
236
    
thomaskummer's avatar
thomaskummer committed
237
    //============================================
thomaskummer's avatar
thomaskummer committed
238
    // Initialize electrophysiology
thomaskummer's avatar
thomaskummer committed
239
    //============================================
thomaskummer's avatar
thomaskummer committed
240
241
242
    solver.initialize();
    
    
thomaskummer's avatar
thomaskummer committed
243
    //============================================
thomaskummer's avatar
thomaskummer committed
244
    // Building Matrices
thomaskummer's avatar
thomaskummer committed
245
    //============================================
thomaskummer's avatar
thomaskummer committed
246
    std::string electroMechanicsCoupling = dataFile("electrophysiology/discretization/coupling", "one-way");
thomaskummer's avatar
thomaskummer committed
247
    if (electroMechanicsCoupling == "two-way") solver.twoWayCoupling();
thomaskummer's avatar
thomaskummer committed
248
249
    else solver.oneWayCoupling();
    
thomaskummer's avatar
thomaskummer committed
250
251
252
    solver.structuralOperatorPtr()->setNewtonParameters(dataFile);
    solver.buildSystem();
    
thomaskummer's avatar
thomaskummer committed
253
    if ( 0 == comm->MyPID() ) std::cout << "\n\nNode number: " << disp.size() / 3 << " -> dof: " << disp.size() << "\n\n";
thomaskummer's avatar
thomaskummer committed
254
    
thomaskummer's avatar
thomaskummer committed
255
    
thomaskummer's avatar
thomaskummer committed
256
    //============================================
thomaskummer's avatar
thomaskummer committed
257
    // Electric stimulus function
thomaskummer's avatar
thomaskummer committed
258
    //============================================
thomaskummer's avatar
thomaskummer committed
259
    function_Type stim = &HeartSolver<EMSolver<mesh_Type, monodomain_Type> >::Iapp;
thomaskummer's avatar
thomaskummer committed
260

thomaskummer's avatar
thomaskummer committed
261
    
thomaskummer's avatar
thomaskummer committed
262
263
    //============================================
    // Apply essential patch b.c.
thomaskummer's avatar
thomaskummer committed
264
265
    //============================================
    patchHandler.applyPatchBC(solver);
thomaskummer's avatar
thomaskummer committed
266
267
    heartSolver.setPatchDisplacementSumPtr(patchHandler.patchDisplacementSumPtr());
    heartSolver.setPatchLocationSumPtr(patchHandler.patchLocationSumPtr());
thomaskummer's avatar
thomaskummer committed
268
269

    
thomaskummer's avatar
thomaskummer committed
270
271
272
    //============================================
    // Setup exporters for EMSolver
    //============================================
thomaskummer's avatar
thomaskummer committed
273
274
    heartSolver.setupExporter(problemFolder);
    
thomaskummer's avatar
thomaskummer committed
275
    
thomaskummer's avatar
thomaskummer committed
276
    //============================================
thomaskummer's avatar
thomaskummer committed
277
    // Pressure b.c. on endocardia
thomaskummer's avatar
thomaskummer committed
278
    //============================================
thomaskummer's avatar
thomaskummer committed
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
    std::vector<vectorPtr_Type> pVecPtrs;
    std::vector<bcVectorPtr_Type> pBCVecPtrs;

    std::vector<ID> flagsBC;
    std::vector<UInt> ventIdx;

    // Endocardia
    UInt nVarBC = dataFile.vector_variable_size ( ( "solid/boundary_conditions/listVariableBC" ) );
    for ( UInt i (0) ; i < nVarBC ; ++i )
    {
        std::string varBCSection = dataFile ( ( "solid/boundary_conditions/listVariableBC" ), " ", i );
        flagsBC.push_back( dataFile ( ("solid/boundary_conditions/" + varBCSection + "/flag").c_str(), 0 ) );
        ventIdx.push_back ( dataFile ( ("solid/boundary_conditions/" + varBCSection + "/index").c_str(), 0 ) );
        
        pVecPtrs.push_back ( vectorPtr_Type ( new vector_Type ( solver.structuralOperatorPtr() -> displacement().map(), Repeated ) ) );
        pBCVecPtrs.push_back ( bcVectorPtr_Type( new bcVector_Type( *pVecPtrs[i], solver.structuralOperatorPtr() -> dispFESpacePtr() -> dof().numTotalDof(), 1 ) ) );
        solver.bcInterfacePtr() -> handler() -> addBC(varBCSection, flagsBC[i], Natural, Full, *pBCVecPtrs[i], 3);
    }
    
thomaskummer's avatar
thomaskummer committed
298
    // Functions to modify b.c.
thomaskummer's avatar
thomaskummer committed
299
    auto modifyPressureBC = [&] (const std::vector<Real>& bcValues)
thomaskummer's avatar
thomaskummer committed
300
301
302
    {
        for ( UInt i (0) ; i < nVarBC ; ++i )
        {
thomaskummer's avatar
thomaskummer committed
303
            *pVecPtrs[i] = - bcValues[ ventIdx[i] ] mmHg;
thomaskummer's avatar
thomaskummer committed
304
305
306
307
308
            pBCVecPtrs[i].reset ( ( new bcVector_Type (*pVecPtrs[i], solver.structuralOperatorPtr() -> dispFESpacePtr() -> dof().numTotalDof(), 1) ) );
            solver.bcInterfacePtr() -> handler() -> modifyBC(flagsBC[i], *pBCVecPtrs[i]);
        }
    };

thomaskummer's avatar
thomaskummer committed
309
    
thomaskummer's avatar
thomaskummer committed
310
    //============================================
thomaskummer's avatar
thomaskummer committed
311
    // Update and print bcHandler
thomaskummer's avatar
thomaskummer committed
312
    //============================================
thomaskummer's avatar
thomaskummer committed
313
    solver.bcInterfacePtr()->handler()->bcUpdate( *solver.structuralOperatorPtr()->dispFESpacePtr()->mesh(), solver.structuralOperatorPtr()->dispFESpacePtr()->feBd(), solver.structuralOperatorPtr()->dispFESpacePtr()->dof() );
thomaskummer's avatar
thomaskummer committed
314
    
thomaskummer's avatar
thomaskummer committed
315
316
317
    if ( 0 == comm->MyPID() ) solver.bcInterfacePtr() -> handler() -> showMe();
    
    
thomaskummer's avatar
thomaskummer committed
318
    //============================================
thomaskummer's avatar
thomaskummer committed
319
    // Volume integrators
thomaskummer's avatar
thomaskummer committed
320
    //============================================
thomaskummer's avatar
thomaskummer committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
    std::vector<int> LVFlags;
    std::vector<int> RVFlags;

    for ( UInt i (0) ; i < nVarBC ; ++i )
    {
        std::string varBCSection = dataFile ( ( "solid/boundary_conditions/listVariableBC" ), " ", i );
        ID flag  =  dataFile ( ("solid/boundary_conditions/" + varBCSection + "/flag").c_str(), 0 );
        ID index =  dataFile ( ("solid/boundary_conditions/" + varBCSection + "/index").c_str(), 0 );
        switch ( index )
        {
            case 0: LVFlags.push_back ( flag );
                break;
            case 1: RVFlags.push_back ( flag );
                break;
            default:
                break;
        }
    }
    
    VolumeIntegrator LV (LVFlags, "Left Ventricle", solver.fullMeshPtr(), solver.localMeshPtr(), ETFESpace, FESpace);
    VolumeIntegrator RV (RVFlags, "Right Ventricle", solver.fullMeshPtr(), solver.localMeshPtr(), ETFESpace, FESpace);
thomaskummer's avatar
thomaskummer committed
342
    if ( 0 == comm->MyPID() ) std::cout << "\n\n";
thomaskummer's avatar
thomaskummer committed
343
    
thomaskummer's avatar
thomaskummer committed
344
    //============================================
thomaskummer's avatar
thomaskummer committed
345
    // Set variables and functions
thomaskummer's avatar
thomaskummer committed
346
    //============================================
thomaskummer's avatar
thomaskummer committed
347
    const Real dt_activation = solver.data().electroParameter<Real>("timestep");
thomaskummer's avatar
thomaskummer committed
348
349
350
351
    const Real dt_loadstep =  dataFile ( "solid/time_discretization/dt_loadstep", 1.0 );
    const Real activationLimit_loadstep =  dataFile ( "solid/time_discretization/activation_limit_loadstep", 0.0 );
    const Real dt_mechanics = solver.data().solidParameter<Real>("timestep");
    const Real dt_save = dataFile ( "exporter/save", 10. );
thomaskummer's avatar
thomaskummer committed
352
    Real exportTime (0.);
thomaskummer's avatar
thomaskummer committed
353
354
355
356
357
358
359
360
361
362
    const Real endtime = solver.data().electroParameter<Real>("endtime");
    const UInt mechanicsLoadstepIter = static_cast<UInt>( dt_loadstep / dt_activation );
    const UInt mechanicsCouplingIter = static_cast<UInt>( dt_mechanics / dt_activation );
    const UInt maxiter = static_cast<UInt>( endtime / dt_activation ) ;
    
    const Real pPerturbationFe = dataFile ( "solid/coupling/pPerturbationFe", 1e-2 );
    const Real pPerturbationCirc = dataFile ( "solid/coupling/pPerturbationCirc", 1e-3 );
    const Real couplingError = dataFile ( "solid/coupling/couplingError", 1e-6 );
    const UInt couplingJFeSubIter = dataFile ( "solid/coupling/couplingJFeSubIter", 1 );
    const UInt couplingJFeSubStart = dataFile ( "solid/coupling/couplingJFeSubStart", 1 );
thomaskummer's avatar
thomaskummer committed
363
364
    
    const Real dpMax = dataFile ( "solid/coupling/dpMax", 0.1 );
thomaskummer's avatar
thomaskummer committed
365
    
thomaskummer's avatar
thomaskummer committed
366
367
368
    std::vector<std::vector<std::string> > bcNames { { "lv" , "p" } , { "rv" , "p" } };
    std::vector<double> bcValues { p ( "lv" ) , p ( "rv") };
    std::vector<double> bcValuesPre ( bcValues );
thomaskummer's avatar
thomaskummer committed
369
    std::vector<double> bcValues4thOAB ( bcValues );
thomaskummer's avatar
thomaskummer committed
370
371
372
373

    VectorSmall<2> VCirc, VCircNew, VCircPert, VFe, VFeNew, VFePert, R, dp;
    MatrixSmall<2,2> JFe, JCirc, JR;

thomaskummer's avatar
thomaskummer committed
374
375
    VectorSmall<2> AvgWorkVent;
    
thomaskummer's avatar
thomaskummer committed
376
377
378
379
380
    UInt iter (0);
    Real t (0);
    
    auto printCoupling = [&] ( std::string label ) { if ( 0 == comm->MyPID() )
    {
thomaskummer's avatar
thomaskummer committed
381
        std::cout.precision(5);
thomaskummer's avatar
thomaskummer committed
382
383
        std::cout << "\n=============================================================";
        std::cout << "\nCoupling iteration " << iter << " at time " << t << " (" << label << ")";
thomaskummer's avatar
thomaskummer committed
384
385
386
387
        std::cout << "\nPressure: \t\t" << bcValues[0] << " / " << bcValues[1];
        std::cout << "\nFE-Volume: \t\t" << VFeNew[0] << " / " << VFeNew[1];
        std::cout << "\nCirculation-Vol: \t" << VCircNew[0] << " / " << VCircNew[1];
        std::cout << "\nResidual: \t\t" << std::abs(VFeNew[0] - VCircNew[0]) << " / " << std::abs(VFeNew[1] - VCircNew[1]);
thomaskummer's avatar
thomaskummer committed
388
389
390
        //std::cout << "\nJFe   = " << JFe;
        //std::cout << "\nJCirc = " << JCirc;
        //std::cout << "\nJR    = " << JR;
thomaskummer's avatar
thomaskummer committed
391
        std::cout << "\n=============================================================\n"; }
thomaskummer's avatar
thomaskummer committed
392
    };
thomaskummer's avatar
thomaskummer committed
393
    
thomaskummer's avatar
thomaskummer committed
394

thomaskummer's avatar
thomaskummer committed
395
    //============================================
thomaskummer's avatar
thomaskummer committed
396
    // Load restart file
thomaskummer's avatar
thomaskummer committed
397
    //============================================
thomaskummer's avatar
thomaskummer committed
398
399
400
    std::string restartInput = command_line.follow ("noRestart", 2, "-r", "--restart");
    const bool restart ( restartInput != "noRestart" );

thomaskummer's avatar
b    
thomaskummer committed
401
402
    if ( restart )
    {
Thomas Kummer's avatar
Thomas Kummer committed
403
404
405
        LifeChrono chronoRestart;
        chronoRestart.start();
        
thomaskummer's avatar
b    
thomaskummer committed
406
        const std::string restartDir = command_line.follow (problemFolder.c_str(), 2, "-rd", "--restartDir");
thomaskummer's avatar
thomaskummer committed
407
408
        const std::string restoreAllPreviousTimestepsStr = command_line.follow ("no", 2, "-rsa", "--restoreAllPreviousTimesteps");
        const bool restoreAllPreviousTimesteps = ( restoreAllPreviousTimestepsStr != "no" );
thomaskummer's avatar
thomaskummer committed
409
        
thomaskummer's avatar
thomaskummer committed
410
        Real dtExport = dt_save; //5.;
thomaskummer's avatar
thomaskummer committed
411

thomaskummer's avatar
b    
thomaskummer committed
412
413
414
415
        // Set time variable
        const unsigned int restartInputStr = std::stoi(restartInput);
        const unsigned int nIter = (restartInputStr - 1) * dtExport / dt_mechanics;
        t = nIter * dt_mechanics;
thomaskummer's avatar
thomaskummer committed
416

thomaskummer's avatar
b    
thomaskummer committed
417
        // Set time exporter time index
thomaskummer's avatar
thomaskummer committed
418
        // if ( restoreAllPreviousTimesteps ) heartSolver.exporter()->setTimeIndex(restartInputStr); // + 1);
thomaskummer's avatar
b    
thomaskummer committed
419
420
421
422

        // Load restart solutions from output files
        std::string polynomialDegree = dataFile ( "solid/space_discretization/order", "P2");

thomaskummer's avatar
thomaskummer committed
423
        // Import and save initial conditions
thomaskummer's avatar
thomaskummer committed
424
        if ( restoreAllPreviousTimesteps )
thomaskummer's avatar
thomaskummer committed
425
        {
thomaskummer's avatar
thomaskummer committed
426
            if ( 0 == comm->MyPID() ) std::cout << "  TIME = " << "-1" << ": import frame " << "00000" << std::endl;
thomaskummer's avatar
thomaskummer committed
427
428
429
            
            heartSolver.exporter()->setTimeIndex(0);
            
thomaskummer's avatar
thomaskummer committed
430
431
432
433
434
435
436
437
438
439
            ElectrophysiologyUtility::importVectorField ( solver.structuralOperatorPtr() -> displacementPtr(), "humanHeartSolution" , "Displacement", solver.localMeshPtr(), restartDir, polynomialDegree, "00000" );
            
            for ( unsigned int i = 0; i < solver.electroSolverPtr()->ionicModelPtr()->Size() ; ++i )
            {
                ElectrophysiologyUtility::importScalarField (solver.electroSolverPtr()->globalSolution().at(i), "humanHeartSolution" , ("Ionic Variable " + std::to_string(i)), solver.localMeshPtr(), restartDir, polynomialDegree, "00000" );
            }
            
            ElectrophysiologyUtility::importScalarField (solver.activationModelPtr() -> fiberActivationPtr(), "humanHeartSolution" , "Activation", solver.localMeshPtr(), restartDir, polynomialDegree, "00000" );
            
            ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "humanHeartSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, "00000" );
thomaskummer's avatar
thomaskummer committed
440

thomaskummer's avatar
thomaskummer committed
441
442
443
444
445
446
447
448
449
            heartSolver.postProcess(-1.0);
            
            if ( 0 == comm->MyPID() )
            {
                //std::cout << "\n*****************************************************************";
                std::cout << "  Restart data at TIME = -1.0 imported in " << chronoRestart.diff() << " s";
                std::cout << "\n*****************************************************************\n";
            }
        }
thomaskummer's avatar
thomaskummer committed
450
451
452
453
        else
        {
            heartSolver.exporter()->setTimeIndex(restartInputStr);
        }
thomaskummer's avatar
thomaskummer committed
454
        
thomaskummer's avatar
thomaskummer committed
455
        // Import and save until desired restart frame
thomaskummer's avatar
thomaskummer committed
456
        Real t_ = ( restoreAllPreviousTimesteps ? 0. : t );
thomaskummer's avatar
thomaskummer committed
457

thomaskummer's avatar
thomaskummer committed
458
        for (t_ ; t_ <= t ; t_ += dtExport)
thomaskummer's avatar
b    
thomaskummer committed
459
        {
thomaskummer's avatar
thomaskummer committed
460
            solver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(t_);
thomaskummer's avatar
thomaskummer committed
461
            
thomaskummer's avatar
thomaskummer committed
462
            std::string importNumber = "00000" + std::to_string(int(t_ / dtExport + 1.0));
thomaskummer's avatar
thomaskummer committed
463
            importNumber = importNumber.substr(importNumber.length() - 5, importNumber.length());
thomaskummer's avatar
thomaskummer committed
464
            if ( 0 == comm->MyPID() ) std::cout << "TIME = " << t_ << ": import frame " << importNumber << std::endl;
thomaskummer's avatar
thomaskummer committed
465
466
            
            ElectrophysiologyUtility::importVectorField ( solver.structuralOperatorPtr() -> displacementPtr(), "humanHeartSolution" , "Displacement", solver.localMeshPtr(), restartDir, polynomialDegree, importNumber );
thomaskummer's avatar
b    
thomaskummer committed
467

thomaskummer's avatar
thomaskummer committed
468
469
            for ( unsigned int i = 0; i < solver.electroSolverPtr()->ionicModelPtr()->Size() ; ++i )
            {
thomaskummer's avatar
thomaskummer committed
470
                ElectrophysiologyUtility::importScalarField (solver.electroSolverPtr()->globalSolution().at(i), "humanHeartSolution" , ("Ionic Variable " + std::to_string(i)), solver.localMeshPtr(), restartDir, polynomialDegree, importNumber );
thomaskummer's avatar
thomaskummer committed
471
472
            }

thomaskummer's avatar
thomaskummer committed
473
            ElectrophysiologyUtility::importScalarField (solver.activationModelPtr() -> fiberActivationPtr(), "humanHeartSolution" , "Activation", solver.localMeshPtr(), restartDir, polynomialDegree, importNumber );
thomaskummer's avatar
thomaskummer committed
474

thomaskummer's avatar
thomaskummer committed
475
            ElectrophysiologyUtility::importScalarField (solver.activationTimePtr(), "humanHeartSolution" , "Activation Time", solver.localMeshPtr(), restartDir, polynomialDegree, importNumber );
thomaskummer's avatar
thomaskummer committed
476

thomaskummer's avatar
thomaskummer committed
477
            heartSolver.postProcess(t_);
thomaskummer's avatar
thomaskummer committed
478

Thomas Kummer's avatar
Thomas Kummer committed
479
480
            if ( 0 == comm->MyPID() )
            {
thomaskummer's avatar
thomaskummer committed
481
                //std::cout << "\n*****************************************************************";
thomaskummer's avatar
thomaskummer committed
482
                std::cout << "  Restart data at TIME = " << t_ << " imported after " << chronoRestart.diff() << " s";
thomaskummer's avatar
thomaskummer committed
483
                std::cout << "\n\n*****************************************************************\n";
Thomas Kummer's avatar
Thomas Kummer committed
484
485
            }
  
thomaskummer's avatar
thomaskummer committed
486
487
488
        }
        
        // Circulation export and/or restart
thomaskummer's avatar
thomaskummer committed
489
        for (t_ = 0 ; t_ <= t ; t_ += dtExport)
thomaskummer's avatar
thomaskummer committed
490
        {
thomaskummer's avatar
thomaskummer committed
491
            circulationSolver.restartFromFile ( restartDir + "solution.dat" , int(t_/dt_mechanics) );
thomaskummer's avatar
thomaskummer committed
492
            circulationSolver.exportSolution( circulationOutputFile );
thomaskummer's avatar
thomaskummer committed
493
            
thomaskummer's avatar
thomaskummer committed
494
495
            if ( 0 == comm->MyPID() ) std::cout << "  TIME = " << t_ << ": import circulation sub steps" << std::endl;
            
thomaskummer's avatar
thomaskummer committed
496
497
498
499
500
501
502
503
504
505
506
            if (t_ < t)
            {
                for (int nSub (1); nSub < dtExport/dt_mechanics; ++nSub)
                {
                    if ( 0 == comm->MyPID() ) std::cout << "  TIME = " << t_ + nSub*dt_mechanics << ": import circulation sub steps" << std::endl;
                    
                    circulationSolver.restartFromFile ( restartDir + "solution.dat" , int(t_/dt_mechanics) + nSub );
                    circulationSolver.exportSolution( circulationOutputFile );
                    
                }
            }
thomaskummer's avatar
thomaskummer committed
507

thomaskummer's avatar
thomaskummer committed
508
        }
thomaskummer's avatar
b    
thomaskummer committed
509
        
Thomas Kummer's avatar
Thomas Kummer committed
510
511
512
513
514
515
516
        if ( 0 == comm->MyPID() )
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nRestart data imported in " << chronoRestart.diff() << " s";
            std::cout << "\n*****************************************************************\n";
        }
        
thomaskummer's avatar
b    
thomaskummer committed
517
518
519
        // Set boundary mechanics conditions
        bcValues = { p ( "lv" ) , p ( "rv" ) };
        bcValuesPre = { p ( "lv" ) , p ( "rv" ) };
thomaskummer's avatar
thomaskummer committed
520
        bcValues4thOAB = { p ( "lv" ) , p ( "rv" ) };
thomaskummer's avatar
thomaskummer committed
521
        
thomaskummer's avatar
thomaskummer committed
522
523
        //modifyEssentialPatchBC(t);
        patchHandler.modifyPatchBC(solver, t);
thomaskummer's avatar
thomaskummer committed
524
        modifyPressureBC(bcValues);
thomaskummer's avatar
thomaskummer committed
525
526
        solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
        solver.solveMechanics();
thomaskummer's avatar
b    
thomaskummer committed
527
    }
thomaskummer's avatar
thomaskummer committed
528
529

    
thomaskummer's avatar
thomaskummer committed
530
    //============================================
thomaskummer's avatar
thomaskummer committed
531
    // Preload
thomaskummer's avatar
thomaskummer committed
532
    //============================================
thomaskummer's avatar
thomaskummer committed
533
    
thomaskummer's avatar
thomaskummer committed
534
535
536
537
538
539
540
    if ( ! restart )
    {
        solver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(0.0);
        
        const int preloadSteps = dataFile ( "solid/boundary_conditions/numPreloadSteps", 0);
        const bool exportPreload = dataFile ( "exporter/preload", false );
        const bool testPatchesAtPreload = dataFile ( "solid/patches/testAtPreload", false );
thomaskummer's avatar
thomaskummer committed
541
        const int testPatchesAtPreloadDt = dataFile ( "solid/patches/testAtPreloadDt", 1);
thomaskummer's avatar
thomaskummer committed
542

thomaskummer's avatar
thomaskummer committed
543
544
545
546
547
548
549
550
        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
551

thomaskummer's avatar
thomaskummer committed
552
553
        //solver.saveSolution (-1.0);
        heartSolver.postProcess( (exportPreload ? -preloadSteps : -1.0) );
thomaskummer's avatar
thomaskummer committed
554

thomaskummer's avatar
thomaskummer committed
555
556
557
        LV.volume(disp, dETFESpace, - 1);
        RV.volume(disp, dETFESpace, 1);
        
thomaskummer's avatar
thomaskummer committed
558
559
560
561
562
563
564
565
566
567
568
569
        if ( 0 == comm->MyPID() )
        {
            std::cout << "\n*****************************************************************";
            std::cout << "\nData stored in " << chronoSave.diff() << " s";
            std::cout << "\n*****************************************************************\n";
        }
        
        LifeChrono chronoPreload;
        chronoPreload.start();
        
        for (int i (1); i <= preloadSteps; i++)
        {
thomaskummer's avatar
thomaskummer committed
570
571
572
            if ( 0 == comm->MyPID() )
            {
                std::cout << "\n*****************************************************************";
thomaskummer's avatar
thomaskummer committed
573
                std::cout << "\nPreload step: " << i << " / " << preloadSteps;
thomaskummer's avatar
thomaskummer committed
574
575
                std::cout << "\n*****************************************************************\n";
            }
thomaskummer's avatar
thomaskummer committed
576

thomaskummer's avatar
thomaskummer committed
577
            // Update b.c.
thomaskummer's avatar
thomaskummer committed
578
            if ( !testPatchesAtPreload )
thomaskummer's avatar
thomaskummer committed
579
            {
thomaskummer's avatar
thomaskummer committed
580
                modifyPressureBC(preloadPressure(bcValues, i, preloadSteps));
thomaskummer's avatar
thomaskummer committed
581
            }
thomaskummer's avatar
thomaskummer committed
582
            else
thomaskummer's avatar
thomaskummer committed
583
            {
thomaskummer's avatar
thomaskummer committed
584
                //modifyEssentialPatchBC(t);
thomaskummer's avatar
thomaskummer committed
585
                patchHandler.modifyPatchBC(solver, i * testPatchesAtPreloadDt);
thomaskummer's avatar
thomaskummer committed
586
            }
thomaskummer's avatar
thomaskummer committed
587

thomaskummer's avatar
thomaskummer committed
588
589
590
591
            // Solve mechanics
            solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
            solver.solveMechanics();
            
thomaskummer's avatar
thomaskummer committed
592
            if (testPatchesAtPreload) heartSolver.postProcess(i-1);
thomaskummer's avatar
thomaskummer committed
593
        }
thomaskummer's avatar
thomaskummer committed
594
595
596
        
        if ( testPatchesAtPreload )
        {
thomaskummer's avatar
thomaskummer committed
597
598
            MPI_Barrier(MPI_COMM_WORLD);

thomaskummer's avatar
thomaskummer committed
599
600
601
602
603
            #ifdef HAVE_MPI
                MPI_Finalize();
            #endif
                return 0;
        }
thomaskummer's avatar
thomaskummer committed
604

thomaskummer's avatar
thomaskummer committed
605
606
        auto maxI4fValue ( solver.activationModelPtr()->I4f().maxValue() );
        auto minI4fValue ( solver.activationModelPtr()->I4f().minValue() );
thomaskummer's avatar
thomaskummer committed
607
        
thomaskummer's avatar
thomaskummer committed
608
        if ( 0 == comm->MyPID() )
thomaskummer's avatar
thomaskummer committed
609
        {
thomaskummer's avatar
thomaskummer committed
610
611
612
613
            std::cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
            std::cout << "\nI4fmax/I4fmin = " << maxI4fValue << "/" << minI4fValue << std::endl;
            std::cout << "\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
        }
thomaskummer's avatar
thomaskummer committed
614
        
thomaskummer's avatar
thomaskummer committed
615
        if ( 0 == comm->MyPID() )
thomaskummer's avatar
thomaskummer committed
616
        {
thomaskummer's avatar
thomaskummer committed
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
            std::cout << "\n*****************************************************************";
            std::cout << "\nPreload done in: " << chronoPreload.diff();
            std::cout << "\n*****************************************************************\n";
        }

    }
    

    //============================================
    // Time loop
    //============================================
    
    VFe[0] = LV.volume(disp, dETFESpace, - 1);
    VFe[1] = RV.volume(disp, dETFESpace, 1);
    VCirc = VFe;
    
    VectorEpetra dispCurrent ( disp );
thomaskummer's avatar
thomaskummer committed
634
635
    VectorEpetra dispPre ( disp );

thomaskummer's avatar
thomaskummer committed
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
    ID bdPowerFlag  =  dataFile ( ("solid/boundary_conditions/LVEndo/flag") , 0 );
    
    printCoupling("Initial values");
    
    auto perturbedPressure = [] (std::vector<double> p, const double& dp)
    {
        for (auto& i : p) {i += dp;}
        return p;
    };
    
    auto perturbedPressureComp = [] (std::vector<double> p, const double& dp, int comp)
    {
        p[comp] += dp;
        return p;
    };

    if ( ! restart )
    {
        heartSolver.postProcess(t);
        circulationSolver.exportSolution( circulationOutputFile );
    }

thomaskummer's avatar
thomaskummer committed
658
    LifeChrono chronoExport;
thomaskummer's avatar
thomaskummer committed
659
660
    chronoExport.start();

thomaskummer's avatar
thomaskummer committed
661
662
663
    for (int k (1); k <= maxiter; k++)
    {
        if ( 0 == comm->MyPID() )
thomaskummer's avatar
thomaskummer committed
664
        {
thomaskummer's avatar
thomaskummer committed
665
666
667
            std::cout << "\n*****************************************************************";
            std::cout << "\nTIME = " << t+dt_activation;
            std::cout << "\n*****************************************************************\n";
thomaskummer's avatar
thomaskummer committed
668
669
        }

thomaskummer's avatar
thomaskummer committed
670
        t = t + dt_activation;
thomaskummer's avatar
thomaskummer committed
671
        
thomaskummer's avatar
thomaskummer committed
672
673
674
        //============================================
        // Solve electrophysiology and activation
        //============================================
thomaskummer's avatar
thomaskummer committed
675

thomaskummer's avatar
thomaskummer committed
676
677
678
679
680
        auto maxI4fValue ( solver.activationModelPtr()->I4f().maxValue() );
        auto minI4fValue ( solver.activationModelPtr()->I4f().minValue() );
        
        solver.solveElectrophysiology (stim, t);
        solver.solveActivation (dt_activation);
thomaskummer's avatar
thomaskummer committed
681
682
        
        
thomaskummer's avatar
thomaskummer committed
683
684
685
        //============================================
        // Load steps mechanics (activation & b.c.)
        //============================================
thomaskummer's avatar
thomaskummer committed
686

thomaskummer's avatar
thomaskummer committed
687
        auto minActivationValue ( solver.activationModelPtr() -> fiberActivationPtr() -> minValue() );
thomaskummer's avatar
thomaskummer committed
688

thomaskummer's avatar
thomaskummer committed
689
690
691
        const bool activationBelowLoadstepThreshold (minActivationValue < activationLimit_loadstep);
        const bool makeLoadstep (k % mechanicsLoadstepIter == 0 && activationBelowLoadstepThreshold);
        const bool makeMechanicsCirculationCoupling (k % mechanicsCouplingIter == 0);
thomaskummer's avatar
thomaskummer committed
692

thomaskummer's avatar
thomaskummer committed
693
694
695
696
        if ( makeLoadstep && !makeMechanicsCirculationCoupling )
        {
            // Linear b.c. extrapolation
            auto bcValuesLoadstep ( bcValues );
thomaskummer's avatar
thomaskummer committed
697
698
            bcValuesLoadstep[0] = bcValues[0] + ( bcValues4thOAB[0] - bcValues[0] ) * ( k % mechanicsCouplingIter ) / mechanicsCouplingIter;
            bcValuesLoadstep[1] = bcValues[1] + ( bcValues4thOAB[1] - bcValues[1] ) * ( k % mechanicsCouplingIter ) / mechanicsCouplingIter;
thomaskummer's avatar
thomaskummer committed
699

thomaskummer's avatar
thomaskummer committed
700
            if ( 0 == comm->MyPID() )
thomaskummer's avatar
thomaskummer committed
701
            {
thomaskummer's avatar
thomaskummer committed
702
703
704
705
706
707
                std::cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
                std::cout << "\nLoad step at time = " << t;
                std::cout << "\nMinimal activation value = " << minActivationValue;
                std::cout << "\nLin. LV-Pressure extrapolation from " <<  bcValues[0] << " to " <<  bcValuesLoadstep[0];
                std::cout << "\nLin. RV-Pressure extrapolation from " <<  bcValues[1] << " to " <<  bcValuesLoadstep[1];
                std::cout << "\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
thomaskummer's avatar
thomaskummer committed
708
            }
thomaskummer's avatar
thomaskummer committed
709
710
711
712

            // Load step mechanics
            solver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(t);
            modifyPressureBC(bcValuesLoadstep);
thomaskummer's avatar
thomaskummer committed
713
714
            //modifyEssentialPatchBC(t);
            patchHandler.modifyPatchBC(solver, t);
thomaskummer's avatar
thomaskummer committed
715
716
717
718
719
720
721
722
723
724
725
726
727
728
            solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
            solver.solveMechanics();
        }
        
        
        //============================================
        // Iterate mechanics / circulation
        //============================================
        
        if ( makeMechanicsCirculationCoupling )
        {
            iter = 0;
            const double dt_circulation ( dt_mechanics / 1000 );
            solver.structuralOperatorPtr() -> data() -> dataTime() -> setTime(t);
thomaskummer's avatar
thomaskummer committed
729
            
thomaskummer's avatar
thomaskummer committed
730
            LifeChrono chronoCoupling;
thomaskummer's avatar
thomaskummer committed
731
            chronoCoupling.start();
thomaskummer's avatar
thomaskummer committed
732
            
thomaskummer's avatar
thomaskummer committed
733
            //============================================
thomaskummer's avatar
thomaskummer committed
734
            // 4th order Adam-Bashforth pressure extrapol.
thomaskummer's avatar
thomaskummer committed
735
            //============================================
thomaskummer's avatar
thomaskummer committed
736
            bcValues = bcValues4thOAB;
thomaskummer's avatar
thomaskummer committed
737
            
thomaskummer's avatar
thomaskummer committed
738
            if ( 0 == comm->MyPID() )
thomaskummer's avatar
thomaskummer committed
739
            {
thomaskummer's avatar
thomaskummer committed
740
741
742
743
744
745
746
747
748
749
750
                std::cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
                std::cout << "\nA.B. LV-Pressure extrapolation from " <<  bcValuesPre[0] << " to " <<  bcValues[0];
                std::cout << "\nA.B. RV-Pressure extrapolation from " <<  bcValuesPre[1] << " to " <<  bcValues[1];
                std::cout << "\nMinimal activation value = " << minActivationValue;
                std::cout << "\nI4fmax/I4fmin = " << maxI4fValue << "/" << minI4fValue;
                std::cout << "\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n";
            }
            
            //============================================
            // Solve mechanics
            //============================================
thomaskummer's avatar
thomaskummer committed
751
752
            //modifyEssentialPatchBC(t);
            patchHandler.modifyPatchBC(solver, t);
thomaskummer's avatar
thomaskummer committed
753
754
755
756
757
758
            modifyPressureBC(bcValues);
            solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
            solver.solveMechanics();
            
            VFeNew[0] = LV.volume(disp, dETFESpace, - 1);
            VFeNew[1] = RV.volume(disp, dETFESpace, 1);
thomaskummer's avatar
thomaskummer committed
759

thomaskummer's avatar
thomaskummer committed
760
761
762
763
764
765
            //============================================
            // Solve circlation
            //============================================
            circulationSolver.iterate(dt_circulation, bcNames, bcValues, iter);
            VCircNew[0] = VCirc[0] + dt_circulation * ( Q("la", "lv") - Q("lv", "sa") );
            VCircNew[1] = VCirc[1] + dt_circulation * ( Q("ra", "rv") - Q("rv", "pa") );
thomaskummer's avatar
thomaskummer committed
766

thomaskummer's avatar
thomaskummer committed
767
768
769
770
771
772
773
774
775
776
777
778
            //============================================
            // Residual computation
            //============================================
            R = VFeNew - VCircNew;
            printCoupling("Residual Computation");

            //============================================
            // Newton iterations
            //============================================
            while ( R.norm() > couplingError )
            {
                ++iter;
thomaskummer's avatar
thomaskummer committed
779

thomaskummer's avatar
thomaskummer committed
780
                //============================================
thomaskummer's avatar
thomaskummer committed
781
                // Jacobian circulation
thomaskummer's avatar
thomaskummer committed
782
                //============================================
thomaskummer's avatar
thomaskummer committed
783

thomaskummer's avatar
thomaskummer committed
784
785
786
787
                // Left ventricle
                circulationSolver.iterate(dt_circulation, bcNames, perturbedPressureComp(bcValues, pPerturbationCirc, 0), iter);
                VCircPert[0] = VCirc[0] + dt_circulation * ( Q("la", "lv") - Q("lv", "sa") );
                VCircPert[1] = VCirc[1] + dt_circulation * ( Q("ra", "rv") - Q("rv", "pa") );
thomaskummer's avatar
thomaskummer committed
788

thomaskummer's avatar
thomaskummer committed
789
790
                JCirc(0,0) = ( VCircPert[0] - VCircNew[0] ) / pPerturbationCirc;
                JCirc(1,0) = ( VCircPert[1] - VCircNew[1] ) / pPerturbationCirc;
thomaskummer's avatar
thomaskummer committed
791

thomaskummer's avatar
thomaskummer committed
792
793
794
795
                // Right ventricle
                circulationSolver.iterate(dt_circulation, bcNames, perturbedPressureComp(bcValues, pPerturbationCirc, 1), iter);
                VCircPert[0] = VCirc[0] + dt_circulation * ( Q("la", "lv") - Q("lv", "sa") );
                VCircPert[1] = VCirc[1] + dt_circulation * ( Q("ra", "rv") - Q("rv", "pa") );
thomaskummer's avatar
thomaskummer committed
796

thomaskummer's avatar
thomaskummer committed
797
798
                JCirc(0,1) = ( VCircPert[0] - VCircNew[0] ) / pPerturbationCirc;
                JCirc(1,1) = ( VCircPert[1] - VCircNew[1] ) / pPerturbationCirc;
thomaskummer's avatar
thomaskummer committed
799
800


thomaskummer's avatar
thomaskummer committed
801
802
803
                //============================================
                // Jacobian fe
                //============================================
thomaskummer's avatar
thomaskummer committed
804

thomaskummer's avatar
thomaskummer committed
805
806
                const bool jacobianFeSubIter ( ! ( (iter - couplingJFeSubStart) % couplingJFeSubIter) && iter >= couplingJFeSubStart );
                const bool jacobianFeEmpty ( JFe.norm() == 0 );
thomaskummer's avatar
thomaskummer committed
807

thomaskummer's avatar
thomaskummer committed
808
809
810
811
                if ( jacobianFeSubIter || jacobianFeEmpty )
                {
                    JFe *= 0.0;
                    dispCurrent = disp;
thomaskummer's avatar
thomaskummer committed
812

thomaskummer's avatar
thomaskummer committed
813
814
815
816
                    // Left ventricle
                    modifyPressureBC(perturbedPressureComp(bcValues, pPerturbationFe, 0));
                    solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
                    solver.solveMechanicsLin();
thomaskummer's avatar
thomaskummer committed
817

thomaskummer's avatar
thomaskummer committed
818
819
                    VFePert[0] = LV.volume(disp, dETFESpace, - 1);
                    VFePert[1] = RV.volume(disp, dETFESpace, 1);
thomaskummer's avatar
thomaskummer committed
820

thomaskummer's avatar
thomaskummer committed
821
822
                    JFe(0,0) = ( VFePert[0] - VFeNew[0] ) / pPerturbationFe;
                    JFe(1,0) = ( VFePert[1] - VFeNew[1] ) / pPerturbationFe;
thomaskummer's avatar
thomaskummer committed
823

thomaskummer's avatar
thomaskummer committed
824
                    disp = dispCurrent;
thomaskummer's avatar
thomaskummer committed
825

thomaskummer's avatar
thomaskummer committed
826
827
828
829
                    // Right ventricle
                    modifyPressureBC(perturbedPressureComp(bcValues, pPerturbationFe, 1));
                    solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
                    solver.solveMechanicsLin();
thomaskummer's avatar
thomaskummer committed
830

thomaskummer's avatar
thomaskummer committed
831
832
                    VFePert[0] = LV.volume(disp, dETFESpace, - 1);
                    VFePert[1] = RV.volume(disp, dETFESpace, 1);
thomaskummer's avatar
thomaskummer committed
833

thomaskummer's avatar
thomaskummer committed
834
835
                    JFe(0,1) = ( VFePert[0] - VFeNew[0] ) / pPerturbationFe;
                    JFe(1,1) = ( VFePert[1] - VFeNew[1] ) / pPerturbationFe;
thomaskummer's avatar
thomaskummer committed
836

thomaskummer's avatar
thomaskummer committed
837
838
                    disp = dispCurrent;
                }
thomaskummer's avatar
thomaskummer committed
839

thomaskummer's avatar
thomaskummer committed
840
841
842
843
                //============================================
                // Update pressure b.c.
                //============================================
                JR = JFe - JCirc;
thomaskummer's avatar
thomaskummer committed
844

thomaskummer's avatar
thomaskummer committed
845
846
847
848
849
850
851
852
                if ( JR.determinant() != 0 )
                {
                    dp = ( JR | R );
                    if ( iter > 5 ) dp *= 0.7;
                    if ( iter > 20 ) dp *= 0.5;
                    bcValues[0] -= std::min( std::max( dp(0) , - dpMax ) , dpMax );
                    bcValues[1] -= std::min( std::max( dp(1) , - dpMax ) , dpMax );
                }
thomaskummer's avatar
thomaskummer committed
853

thomaskummer's avatar
thomaskummer committed
854
                printCoupling("Pressure Update");
thomaskummer's avatar
thomaskummer committed
855

thomaskummer's avatar
thomaskummer committed
856
857
858
859
860
861
                //============================================
                // Solve circulation
                //============================================
                circulationSolver.iterate(dt_circulation, bcNames, bcValues, iter);
                VCircNew[0] = VCirc[0] + dt_circulation * ( Q("la", "lv") - Q("lv", "sa") );
                VCircNew[1] = VCirc[1] + dt_circulation * ( Q("ra", "rv") - Q("rv", "pa") );
thomaskummer's avatar
thomaskummer committed
862

thomaskummer's avatar
thomaskummer committed
863
                //============================================
thomaskummer's avatar
thomaskummer committed
864
                // Solve mechanics
thomaskummer's avatar
thomaskummer committed
865
                //============================================
thomaskummer's avatar
thomaskummer committed
866
867
868
869
870
871
872
                modifyPressureBC(bcValues);
                solver.bcInterfacePtr() -> updatePhysicalSolverVariables();
                solver.solveMechanics();

                VFeNew[0] = LV.volume(disp, dETFESpace, - 1);
                VFeNew[1] = RV.volume(disp, dETFESpace, 1);

thomaskummer's avatar
thomaskummer committed
873
                //============================================
thomaskummer's avatar
thomaskummer committed
874
                // Residual update
thomaskummer's avatar
thomaskummer committed
875
                //============================================
thomaskummer's avatar
thomaskummer committed
876
877
878
879
880
881
882
883
                R = VFeNew - VCircNew;
                printCoupling("Residual Update");
            }
 
            if ( 0 == comm->MyPID() )
            {
                std::cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
                std::cout << "\nCoupling converged after " << iter << " iteration" << ( iter > 1 ? "s" : "" );
thomaskummer's avatar
thomaskummer committed
884
                std::cout << "\nTime required: " << chronoCoupling.diff() << " s";
thomaskummer's avatar
thomaskummer committed
885
                std::cout << "\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n";
thomaskummer's avatar
thomaskummer committed
886
887
            }
            
thomaskummer's avatar
thomaskummer committed
888
            //============================================
thomaskummer's avatar
thomaskummer committed
889
            // Update volume variables
thomaskummer's avatar
thomaskummer committed
890
            //============================================
thomaskummer's avatar
thomaskummer committed
891
892
            VCirc = VCircNew;
            VFe = VFeNew;
thomaskummer's avatar
thomaskummer committed
893
            
thomaskummer's avatar
thomaskummer committed
894
895
896
897
898
899
            //============================================
            // 4th order Adam-Bashforth pressure extrapol.
            //============================================
            bcValues4thOAB = bcValues;
            heartSolver.extrapolate4thOrderAdamBashforth(bcValues4thOAB, bcValuesPre, dpMax);
            
thomaskummer's avatar
thomaskummer committed
900
901
902
903
904
            //============================================
            // Export circulation solution
            //============================================
            if ( 0 == comm->MyPID() ) circulationSolver.exportSolution( circulationOutputFile );
            
thomaskummer's avatar
thomaskummer committed
905
906
907
908
909
910
            
            //============================================
            // Power computations
            //============================================
            Real leftVentPower = heartSolver.externalPower(disp, dispPre, dETFESpace, p("lv"), dt_mechanics, 454);
            Real rightVentPower = heartSolver.externalPower(disp, dispPre, dETFESpace, p("rv"), dt_mechanics, 455);
thomaskummer's avatar
thomaskummer committed
911
912
            //Real patchPower1 = heartSolver.externalPower(disp, dispPre, dETFESpace, p("lv"), dt_mechanics, 900);
            //Real patchPower2 = heartSolver.externalPower(disp, dispPre, dETFESpace, p("rv"), dt_mechanics, 901);
thomaskummer's avatar
thomaskummer committed
913
            
thomaskummer's avatar
thomaskummer committed
914
915
916
            AvgWorkVent(0) += leftVentPower * dt_mechanics;
            AvgWorkVent(1) += rightVentPower * dt_mechanics;

thomaskummer's avatar
thomaskummer committed
917
918
919
            if ( 0 == comm->MyPID() )
            {
                std::cout << "\n******************************************";
thomaskummer's avatar
thomaskummer committed
920
921
                std::cout << "\nInstantaneous vent. power: \t" << leftVentPower << " / " << rightVentPower;
                std::cout << "\nAveraged vent. power: \t" << AvgWorkVent(0) / t << " / " << AvgWorkVent(1) / t;
thomaskummer's avatar
thomaskummer committed
922
923
924
925
926
                std::cout << "\n******************************************\n\n";
            }
            
            dispPre = disp;
            
thomaskummer's avatar
thomaskummer committed
927
928
        }
        
thomaskummer's avatar
thomaskummer committed
929
        
thomaskummer's avatar
thomaskummer committed
930
931
932
933
934
935
936
        //============================================
        // Export FE-solution
        //============================================
        bool save ( std::abs(std::remainder(t, dt_save)) < 0.01 );
        if ( save )
        {
            heartSolver.postProcess(t);
thomaskummer's avatar
thomaskummer committed
937
            
thomaskummer's avatar
thomaskummer committed
938
            Real chronoTimeNow = chronoExport.diff();
thomaskummer's avatar
thomaskummer committed
939
940
941
            Real chronoDiffToLastSave = chronoTimeNow - exportTime;
            exportTime = chronoTimeNow;
            
thomaskummer's avatar
thomaskummer committed
942
943
944
            if ( 0 == comm->MyPID() )
            {
                std::cout << "\n>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>";
thomaskummer's avatar
thomaskummer committed
945
                std::cout << "\nTotal simulation time at t = " << t << " ms is " << chronoTimeNow << " s";
thomaskummer's avatar
thomaskummer committed
946
                std::cout << "\nPrevious " << dt_save << " ms computed in " << chronoDiffToLastSave << " s";
thomaskummer's avatar
thomaskummer committed
947
948
                std::cout << "\n<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n\n";
            }
thomaskummer's avatar
thomaskummer committed
949
        }
thomaskummer's avatar
thomaskummer committed
950
951
        
    }
thomaskummer's avatar
thomaskummer committed
952
953

    
thomaskummer's avatar
thomaskummer committed
954
    //============================================
thomaskummer's avatar
thomaskummer committed
955
    // Close all exporters
thomaskummer's avatar
thomaskummer committed
956
    //============================================
thomaskummer's avatar
thomaskummer committed
957
    solver.closeExporters();
thomaskummer's avatar
thomaskummer committed
958
    heartSolver.exporter()->closeFile();
thomaskummer's avatar
thomaskummer committed
959
960
961
962
963
964
965
    

#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
}