FSIMonolithic.hpp 22.4 KB
Newer Older
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
//@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
26

27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
/**
 * @file monolithic.hpp
 * @author Paolo Crosetto
 * @date 13-09-2008
 * Class handling the monolithic solver for FSI problems. The block structure of the matrix can be
 \f$\left(\begin{array}{cc}
 C&B\\
 D&N
 \end{array}\right)\f$
 if the time discretization at hand is the Geometry-Explicit one (implemented in monolithicGE.hpp), or
 \f$\left(\begin{array}{ccc}
 C&B&D\\
 D&N&0\\
 0&I&H
 \end{array}\right)\f$
 if the time discretization at hand is the Geometry-Implicit one (implemented in monolithicGI.hpp)
*/
#ifndef _MONOLITHIC_HPP
#define _MONOLITHIC_HPP

#include <EpetraExt_MatrixMatrix.h>
//#include <EpetraExt_Reindex_MultiVector.h>
//#include <EpetraExt_Reindex_CrsMatrix.h>

51
52
#include <lifev/core/util/LifeChrono.hpp>
#include <lifev/core/fem/FESpace.hpp>
Radu Popescu's avatar
Radu Popescu committed
53
#include <lifev/fsi/solver/FSIOperator.hpp>
54

55
56
#include <lifev/core/algorithm/PreconditionerComposed.hpp>
#include <lifev/core/algorithm/ComposedOperator.hpp>
57
#ifdef HAVE_TRILINOS_ANASAZI
58
#include <lifev/core/algorithm/EigenSolver.hpp>
59
60
#endif

Radu Popescu's avatar
Radu Popescu committed
61
#include <lifev/fsi/solver/MonolithicBlockMatrix.hpp>
Paolo Crosetto's avatar
   
Paolo Crosetto committed
62
#include <lifev/fsi/solver/MonolithicBlockComposed.hpp>
63
64
65
#ifdef HAVE_NS_PREC
#include <life/lifealg/PreconditionerPCD.hpp>
#endif
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93

namespace LifeV
{

class WRONG_PREC_EXCEPTION;

//! FSIMonolithic.hpp pure virtual class containing the core methods of the FSIMonolithic FSI solver
/*!
 * Class handling the monolithic solver for FSI problems. The block structure of the matrix can be
 \f$\left(\begin{array}{cc}
 C&B\\
 D&N
 \end{array}\right)\f$
 if the time discretization at hand is the Geometry-Explicit one (implemented in monolithicGE.hpp), or
 \f$\left(\begin{array}{ccc}
 C&B&D\\
 D&N&0\\
 0&I&H
 \end{array}\right)\f$
 if the time discretization at hand is the Geometry-Implicit one (implemented in monolithicGI.hpp),
 * where \f$N\f$ represents the solid block, \f$C\f$ the fluid block, \f$H\f$ the harmonic extension block, while the extra
 * diagonal blocks represent the coupling. The implementation of the stress continuity coupling condition
 * is obtained by means of an augmented formulation.
 * Different possible preconditioners are implemented.
 * The flag semiImplicit in the data file is used to distinguish between the GCE and CE (with quasi Newton) time discretizations.
 * Exact Newton method and full implicit time discretization are implemented in the FSIMonolithicGI class.
 */

94
class FSIMonolithic : public FSIOperator
95
96
97
98
99
{
public:

    //!@name Typedefs
    //@{
100
101

    typedef FSIOperator                                               super_Type;
Antonio Cervone's avatar
Antonio Cervone committed
102
    typedef FSIOperator::fluid_Type::matrix_Type                      matrix_Type;
103
    typedef boost::shared_ptr<matrix_Type>                            matrixPtr_Type;
104
105
    typedef super_Type::solution_Type                                 solution_Type;
    typedef super_Type::solutionPtr_Type                              solutionPtr_Type;
106
107
108
109
110
111
112
    typedef MonolithicBlock                                           prec_Type;
    typedef boost::shared_ptr<prec_Type>                              precPtr_Type;
    typedef MonolithicBlockMatrix                                     blockMatrix_Type;
    typedef boost::shared_ptr<blockMatrix_Type>                       blockMatrixPtr_Type;
    typedef FactorySingleton< Factory< FSIMonolithic, std::string > > factory_Type;
    typedef SolverAztecOO                                             solver_Type;

113
114
115
116
117
118
119
    //@}

    // constructors

    //!@name Constructors, Destructor
    //!@{
    FSIMonolithic();
120

121
122
123
124
125
126
127
128
    ~FSIMonolithic();
    //!@}

    //!@name Public Setup Methods
    //!@{

    /**
       create FEspace
Cristiano Malossi's avatar
Cristiano Malossi committed
129
     */
130
131
132
133
134
135
    void setupFEspace();


    /**
       sets the interface map between the fluid and solid meshes
       (non scalable, do not use for massively parallel simulations)
Cristiano Malossi's avatar
Cristiano Malossi committed
136
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
137
    virtual void setupDOF ( void );
138

139
#ifdef HAVE_HDF5
140
141
    /**
       reads the interface map between the fluid and solid meshes from file.
Cristiano Malossi's avatar
Cristiano Malossi committed
142
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
143
    void setupDOF ( meshFilter_Type& filterMesh );
144
#endif
145
146
147
148

    //! sets the parameters from data file
    /**
       Calls the setup of the fluid problem and the setUp method.
Cristiano Malossi's avatar
Cristiano Malossi committed
149
     */
150
151
152
153
154
155
156
157
    virtual void setupSystem( );

    //     /** stores the data file into a member */
    //     virtual void setDataFile( const GetPot& dataFile );

    //! setup method for the FSIMonolithic solver
    /**
       sets some parameters specific to the FSIMonolithic class
Cristiano Malossi's avatar
Cristiano Malossi committed
158
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
159
    virtual void setup ( const GetPot& dataFile );
160
161
162
163

    //! builds the global Epetra map
    /**
       assigns each mesh partition to the corresponding processor, builds the monolithic map
Cristiano Malossi's avatar
Cristiano Malossi committed
164
     */
165
166
167
168
169
170
    virtual void setupFluidSolid();

    //! builds the global Epetra map
    /**
       assigns each mesh partition to the corresponding processor, builds the monolithic map with a number of fluxes
       specified from input
Cristiano Malossi's avatar
Cristiano Malossi committed
171
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
172
    virtual void setupFluidSolid ( UInt const fluxes );
173
174
175
176
177
178

    //!@}

    //!@name Public Methods
    //!@{

179
    //! Transfers a vector to the interface
180
181
    /**
       restricts a vector with a monolithic map on the solid interface map
182
183
       \param lambdaSolid: vector on the solid interface
       \param disp: monolithic vector
Cristiano Malossi's avatar
Cristiano Malossi committed
184
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
185
    void monolithicToInterface (    vector_Type& lambdaSolid, const vector_Type& sol) ;
186
187


188
    //!Transfers a vector to a subdomain
189
190
191
    /**
       restricts a vector with a monolithic map on another map that must
       have a sequential numbering (not the interface map)
192
193
194
195
       \param disp: monolithic vector
       \param dispFluid: vector on the fluid domain
       \param map: MapEpetra of the part of vector that we want to transfer
       \param offset: offset for the monolithic vector (also alpplied to the input map)
Cristiano Malossi's avatar
Cristiano Malossi committed
196
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
197
    void monolithicToX (const vector_Type& disp, vector_Type& dispFluid, MapEpetra& map, UInt offset = (UInt) 0);
198
199
200
201


    /**
       builds the constant part of the monolithic matrix
Cristiano Malossi's avatar
Cristiano Malossi committed
202
     */
203
204
205
    void buildSystem();


206
207
208
209
    //! Merges the flux boundary conditions into the fluid BCHandler
    /*!two separate BCHandlers are initially created for the flux-type boundary conditions, these are later merged with the fluid BCHandler
      automatically, using this method
     */
210
211
    void mergeBCHandlers()
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
212
213
214
215
        if (M_BCh_u.get() )
        {
            M_BCh_u->merge (*M_BCh_flux);
        }
Cristiano Malossi's avatar
Cristiano Malossi committed
216
        else // in this case only fluxes are imposed on the fluid
Alessio Fumagalli's avatar
Alessio Fumagalli committed
217
218
219
        {
            M_BCh_u = M_BCh_flux;
        }
Cristiano Malossi's avatar
Cristiano Malossi committed
220
        M_BCh_flux.reset();
221
222
223
224
225
    }

#ifdef HAVE_TRILINOS_ANASAZI
    //! Computes the maximum singular value
    /**
226
       \small Computes the maximum singular value of the preconditioned system \f$P^{-1}A\f$ where \f$P\f$ is an
227
       instance of ComposedOperator and \f$A\f$ is the system matrix.
Cristiano Malossi's avatar
Cristiano Malossi committed
228
     */
229
230
231
232
233
234
235
    LifeV::Real& computeMaxSingularValue();
#endif
    //! Computes the normals to the fluid domain
    /**
       \small Computes the normals to the fluid domain in the nodes. It is an example of how
       to use the boundary conditions methods to compute the normal field on
       a surface.
Cristiano Malossi's avatar
Cristiano Malossi committed
236
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
237
    void computeFluidNormals ( vector_Type& normals);
238
239
240
241
242
243
244


    //!Evaluates the nonlinear residual
    /**
       This class is pure virtual, it depends on which type of monolithic solver is used
       \param res: output
       \param _sol: monolithic solution
245
       \param iter: current NonLinearRichardson (Newton) iteration
Cristiano Malossi's avatar
Cristiano Malossi committed
246
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
247
248
249
    virtual void   evalResidual (vector_Type&        res,
                                 const vector_Type& sol,
                                 const UInt         iter) = 0 ;
250
251
252

    /**
       solves the Jacobian system
253
254
255
       \param muk: output, solution at the current Newton step
       \param res: nonlinear residual
       \param linearRelTol: not used
256
257
258
259
260

       \small The preconditioner type is usually an algebraic additive Schwarz. The following values
       assigned to the field DDBlockPrec in the data file correspond to different variants:

       Only for the FSIMonolithic Geometry Explicit:
261
       - DDBlockPrec = AdditiveSchwarz is AAS on a the system matrix
262
263
264
265
266
267
       - DDBlockPrec = MonolithicBlockComposedDN is AAS on a Dirichlet-Neumann preconditioner using the ComposedOperator strategy
       - DDBlockPrec = ComposedDN2 is AAS on an alternative Dirichlet-Neumann preconditioner using the ComposedOperator strategy
       - DDBlockPrec = MonolithicBlockComposedNN is AAS on a Neumann-Neumann preconditioner using the ComposedOperator strategy
       - DDBlockPrec = MonolithicBlockComposedDNND is AAS on a Dirichler-Neumamm -- Neumann-Dirichlet preconditioner using the ComposedOperator strategy

       Only for the Geometry Implicit:
268
269
270
271
272
273
274
       - DDBlockPrec = AdditiveSchwarzGI is AAS on the whole matrix.
       - DDBlockPrec = MonolithicBlockComposedDNDGI is AAS on the quasi-newton matrix obtained with the ComposedOperator strategy. Split in 3 factors.
       - DDBlockPrec = ComposedDND2GI is AAS on the quasi-newton matrix obtained with the ComposedOperator strategy. Split in 3 factors.
       - DDBlockPrec = MonolithicBlockComposedDNGI is AAS on the full Jacobian matrix obtained with the ComposedOperator strategy by neglecting part
       of the fluid--structure coupling, split in 3 factors
       - DDBlockPrec = MonolithicBlockComposedDN2GI is AAS on an alternative matrix obtained with the previous strategy
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
275
276
277
    virtual void   solveJac (vector_Type&       muk,
                             const vector_Type& res,
                             const Real       linearRelTol);
278
279
280
281

    /**
       updates the meshmotion, advances of a time step
       \param _sol: solution
Cristiano Malossi's avatar
Cristiano Malossi committed
282
     */
283
284
    virtual void updateSystem();

paolo.crosetto's avatar
paolo.crosetto committed
285
286
287
288
    //! activates the computation of the wall stress on the boundary with a specified flag.
    /**
       Notice that the specified flag must be in the coupling fluid-structure interface
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
289
    void enableStressComputation (UInt  flag);
paolo.crosetto's avatar
paolo.crosetto committed
290
291
292
293
294
295
296

    /**
       Computes the stress on the coupling boundary (the traction vector)
     */
    vectorPtr_Type computeStress();


297
298
299
300
301
302
    //@}

    //!@name Set Methods
    //@{

    //! returns a non-const pointer to the preconditioner. Can be used either as a setter or a getter.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
303
304
305
306
    precPtr_Type& precPtrView()
    {
        return M_precPtr;
    }
307
308

    //! returns a non-const pointer to the preconditioner. Can be used either as a setter or a getter.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
309
310
311
312
    blockMatrixPtr_Type& operatorPtrView()
    {
        return M_monolithicMatrix;
    }
313
314

    /**
315
       \small sets the solid BCHandle
Cristiano Malossi's avatar
Cristiano Malossi committed
316
     */
317
318
    virtual void setSolidBC     ( const fluidBchandlerPtr_Type& bc_solid )
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
319
        super_Type::setSolidBC (bc_solid);
320
321
322
323
324
325
        //bc_solid->merge(*M_BCh_Robin);
    }

    //! initializes the solution by reference (through a shared_ptr)
    /*!
      \param sol: input pointer
Cristiano Malossi's avatar
Cristiano Malossi committed
326
     */
327

328
329
    void setFluidBC     ( const fluidBchandlerPtr_Type& bc_fluid )
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
330
331
        super_Type::setFluidBC (bc_fluid);
        if (M_BChs.size() )
332
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
333
            UInt nfluxes (M_BChs[1]->numberOfBCWithType (Flux) );
334
            //M_substituteFlux.reset(new std::vector<bool>(nfluxes))
Alessio Fumagalli's avatar
Alessio Fumagalli committed
335
336
337
            M_fluxOffset.resize (nfluxes);
            M_BCFluxNames = M_BChs[1]->findAllBCWithType (Flux);
            for (UInt i = 0; i < nfluxes; ++i)
Cristiano Malossi's avatar
Cristiano Malossi committed
338
            {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
339
340
                const BCBase* bc = M_BChs[1]->findBCWithName (M_BCFluxNames[i]);
                M_fluxOffset[i] = bc->offset();
Cristiano Malossi's avatar
Cristiano Malossi committed
341
            }
Alessio Fumagalli's avatar
Alessio Fumagalli committed
342
343
344
            M_BChs[1] = bc_fluid;
            M_monolithicMatrix->setConditions (M_BChs);
            M_precPtr->setConditions (M_BChs);
345
346
347
        }
    }

348
349

    //!get the total dimension of the FS interface
Alessio Fumagalli's avatar
Alessio Fumagalli committed
350
351
352
353
    UInt dimInterface() const
    {
        return nDimensions * M_monolithicMatrix->interface();
    }
354
355
356
357
358

    //! Returns true if CE of FI methods are used, false otherwise (GCE)
    //bool const isFullMonolithic(){return M_fullMonolithic;}

    //! Returns the offset assigned to the solid block
Alessio Fumagalli's avatar
Alessio Fumagalli committed
359
360
361
362
    UInt  offset() const
    {
        return M_offset;
    }
363

364
    //!Get the solid displacement from the solution
365
    /*!
366
      \param solidDisplacement: input vector
Cristiano Malossi's avatar
Cristiano Malossi committed
367
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
368
    void exportSolidDisplacement ( vector_Type& solidDisplacement )
369
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
370
        solidDisplacement.subset ( M_fluidTimeAdvance->singleElement (0), M_offset );
371
        solidDisplacement *= M_solid->rescaleFactor();
372
373
374
375
    }

    //!Get the solid velocity
    /*!
376
      fills an input vector with the solid displacement from the solution.
377
      \param solidVelocity: input vector (output solid velocity)
Cristiano Malossi's avatar
Cristiano Malossi committed
378
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
379
    void exportSolidVelocity ( vector_Type& solidVelocity )
380
    {
381
        solidVelocity.subset ( M_solidTimeAdvance->firstDerivative(), M_offset );
382
        solidVelocity *= M_solid->rescaleFactor();
383
    }
384

385
    //!Get the solid accelration
386
387
388
    /*!
      fills an input vector with the solid displacement from the solution.
      \param solidVelocity: input vector (output solid acceleration)
Cristiano Malossi's avatar
Cristiano Malossi committed
389
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
390
    void exportSolidAcceleration ( vector_Type& solidAcceleration )
391
    {
392
        solidAcceleration.subset ( M_solidTimeAdvance->secondDerivative(), M_offset );
Cristiano Malossi's avatar
Cristiano Malossi committed
393
        solidAcceleration *= M_solid->rescaleFactor();
394
    }
395

396
397
398
399
    //! Export the fluid velocity by copying it to an external vector
    /*!
     * @param fluidVelocity vector to be filled with the fluid velocity
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
400
401
402
403
    void exportFluidVelocity ( vector_Type& fluidVelocity )
    {
        fluidVelocity.subset ( M_fluidTimeAdvance->singleElement (0), 0 );
    }
404
405
406
407
408

    //! Export the fluid pressure by copying it to an external vector
    /*!
     * @param fluidPressure vector to be filled with the fluid pressure
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
409
410
411
412
    void exportFluidPressure ( vector_Type& fluidPressure )
    {
        fluidPressure.subset ( M_fluidTimeAdvance->singleElement (0), static_cast<UInt> (3 * M_uFESpace->dof().numTotalDof() ) );
    }
413

414
415
    //! Gets the fluid and pressure
    /**
416
417
       fills an input vector with the fluid and pressure from the solution M_un.
       It performs a trilinos import. Thus it works also for the velocity, depending on the map of the input vector
418
       \param fluidVelocityandPressure: input vector
Cristiano Malossi's avatar
Cristiano Malossi committed
419
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
420
    void exportFluidVelocityAndPressure ( vector_Type& fluidVelocityAndPressure )
421
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
422
        fluidVelocityAndPressure = M_fluidTimeAdvance->singleElement (0);
423
424
425
    }

    //! Returns the monolithic map
Alessio Fumagalli's avatar
Alessio Fumagalli committed
426
427
428
429
    virtual boost::shared_ptr<MapEpetra>& couplingVariableMap()
    {
        return M_monolithicMap;
    }
430
431
432

    //! get the solution vector
    virtual const vector_Type& solution() const = 0;
433
434
435
436
437

    //! Update the solution after NonLinearRichardson is called.
    /*!
     *  Here it is used also to update the velocity for the post-processing.
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
438
    virtual void updateSolution ( const vector_Type& solution )
439
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
440
441
442
443
444
445
        this->M_fluidTimeAdvance->shiftRight (solution);
        if (M_data->dataFluid()->conservativeFormulation() )
        {
            this->M_fluidMassTimeAdvance->shiftRight (M_fluid->matrixMass() *solution);
        }
        this->M_solidTimeAdvance->shiftRight (solution);
446
447
448
449
450
451
452
453
454
455
456
457
458
    }

    //! Updates the right hand side
    /**
       Adds to the rhs the fluid time discretization terms
       \todo this should be handled externally
     */
    void updateRHS();

    //! Set vectors for restart
    /*!
     *  Set vectors for restart
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
459
460
461
462
463
    void setVectorInStencils ( const vectorPtr_Type& vel,
                               const vectorPtr_Type& pressure,
                               const vectorPtr_Type& solidDisp,
                               //const vectorPtr_Type& fluidDisp,
                               const UInt iter);
464

Alessio Fumagalli's avatar
Alessio Fumagalli committed
465
    void setFluidVectorInStencil ( const vectorPtr_Type& vel, const vectorPtr_Type& pressure, const UInt iter);
466

Alessio Fumagalli's avatar
Alessio Fumagalli committed
467
    void setSolidVectorInStencil ( const vectorPtr_Type& solidDisp, const UInt iter);
468

Alessio Fumagalli's avatar
Alessio Fumagalli committed
469
    virtual void setALEVectorInStencil ( const vectorPtr_Type& fluidDisp, const UInt iter, const bool lastVector) = 0;
470
471
472

    void finalizeRestart();

473
474
475
476
477
478
479
480
481
482
    //@}


protected:


    //!@name Protected methods
    //!@{

    //! pure virtual: creates the operator (either of type FSIMonolithicGI or FSIMonolithicGE)
Alessio Fumagalli's avatar
Alessio Fumagalli committed
483
    virtual void createOperator ( std::string& operType ) = 0;
484
485
486

    /**
       \small solves the monolithic system, once a solver, a preconditioner and a rhs have been defined.
Cristiano Malossi's avatar
Cristiano Malossi committed
487
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
488
    void iterateMonolithic (const vector_Type& rhs, vector_Type& step);
489
490
491
492
493

    /**
       \small adds the part due to coupling to the rhs
       \param rhs: right hand side
       \param un: current solution
Cristiano Malossi's avatar
Cristiano Malossi committed
494
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
495
    void couplingRhs ( vectorPtr_Type rhs);
496
497
498
499
500
501
502


    /**\small evaluates the linear residual
       \param sol: solution vector
       \param rhs: right-hand side
       \param res: the output residual
       \param diagonalScaling: flag stating wether to perform diagonal scaling
Cristiano Malossi's avatar
Cristiano Malossi committed
503
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
504
    void evalResidual ( const vector_Type& sol, const vectorPtr_Type& rhs,  vector_Type& res, bool diagonalScaling = false);
505
506

    //!\small says if the preconditioner will be recomputed
Alessio Fumagalli's avatar
Alessio Fumagalli committed
507
508
509
510
    bool recomputePrec()
    {
        return (!M_reusePrec || M_resetPrec);
    }
511
512

    //!\small updates the rhs of the solid block.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
513
    void updateSolidSystem (vectorPtr_Type& rhsFluidCoupling);
514
515
516
517

    /**\small scales matrix and rhs
       \param rhs: the output rhs
       \param matrFull: the output matrix*/
Alessio Fumagalli's avatar
Alessio Fumagalli committed
518
    void diagonalScale (vector_Type& rhs, matrixPtr_Type matrFull);
519

520

521
522
523
    //! Constructs the solid FESpace
    /**
       Creates the solid FESpace with an unpartitioned mesh, necessary step to create the dof interconnections
524
       at the interface. The solid FESpace will be reset in variablesInit using the partitioned mesh.export
525
526
       If the interface map is created offline this method is never called.
       \param dOrder: discretization order
Cristiano Malossi's avatar
Cristiano Malossi committed
527
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
528
    void solidInit (std::string const& dOrder);
529
530
531
532
533

    //! Constructs the solid FESpace and initializes the coupling variables at the interface
    /**
       If the mesh is partitioned online the previous FESpace constructed with the unpartitioned mesh is discarded and
       replaced with one using a partitioned mesh.
Cristiano Malossi's avatar
Cristiano Malossi committed
534
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
535
    void variablesInit (std::string const& dOrder);
536
537

    //!
538
    virtual void setupBlockPrec();
539
540
541
542
543


    //! assembles the solid problem (the matrix and the rhs due to the time derivative)
    /*
      \param iter: current nonlinear iteration: used as flag to distinguish the first nonlinear iteration from the others
544
      \param solution: current solution, used for the time advance implementation, and thus for the update of the right hand side
Cristiano Malossi's avatar
Cristiano Malossi committed
545
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
546
    void assembleSolidBlock (UInt iter, const vector_Type& solution);
547
548
549
550
551


    //! assembles the fluid problem (the matrix and the rhs due to the time derivative)
    /*
      \param iter: current nonlinear iteration: used as flag to distinguish the first nonlinear iteration from the others
552
      \param solution: current solution, used for the time advance implementation, and thus for the update of the right hand side
Cristiano Malossi's avatar
Cristiano Malossi committed
553
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
554
    void assembleFluidBlock (UInt iter, const vector_Type& solution);
555

paolo.crosetto's avatar
paolo.crosetto committed
556
557
558
559
560
    //! Checks if the flux bcs changed during the simulation, e.g. if a flux b.c. has been changed to Natural
    //! (this can be useful when modeling valves)
    /**
       When the fluxes bcs changed a '1' is added in the line corresponding to the Lagrange multiplier. This method must
       be called for both operator and preconditioner
Cristiano Malossi's avatar
Cristiano Malossi committed
561
     */
Alessio Fumagalli's avatar
Alessio Fumagalli committed
562
    void checkIfChangedFluxBC ( precPtr_Type oper );
paolo.crosetto's avatar
paolo.crosetto committed
563

564
565
566
567
568
    //!@}


    //!@name Protected attributes
    //@{
569
570
    boost::shared_ptr<MapEpetra>                      M_monolithicMap;
    boost::shared_ptr< MapEpetra >                    M_interfaceMap;///the solid interface map
571
572
573
574
575
    boost::shared_ptr<vector_Type>                    M_beta;
    boost::shared_ptr<MonolithicBlockMatrix >                   M_monolithicMatrix;
    precPtr_Type                                      M_precPtr;
    boost::shared_ptr<vector_Type>                    M_rhsFull;

576
577
    fluidBchandlerPtr_Type                            M_BCh_flux;
    solidBchandlerPtr_Type                            M_BChWS;
paolo.crosetto's avatar
paolo.crosetto committed
578
    BCFunctionRobin                                   M_bcfWs;
579
580
    UInt                                              M_offset;
    UInt                                              M_solidAndFluidDim;
Antonio Cervone's avatar
Antonio Cervone committed
581
    FSIOperator::fluid_Type::matrixPtr_Type           M_fluidBlock;
582
583
    matrixPtr_Type                                    M_solidBlockPrec;
    matrixPtr_Type                                    M_robinCoupling; //uninitialized if not needed
584
    matrixPtr_Type                                    M_boundaryMass;
585
586
587
    boost::shared_ptr<solver_Type>                    M_linearSolver;
    boost::shared_ptr<vector_Type>                    M_numerationInterface;
    std::vector<fluidBchandlerPtr_Type>                 M_BChs;
588
    std::vector<boost::shared_ptr<FESpace<mesh_Type, MapEpetra> > >          M_FESpaces;
589
590
591
    bool                                              M_diagonalScale;
    bool                                              M_reusePrec;//!\todo to move to private
    bool                                              M_resetPrec;//!\todo to move to private
592
    Int                                               M_maxIterSolver;//!\todo to move to private
593
    bool                                              M_restarts;
594
595
596
597
598
599
600
    //@}

private:
    //!@name Private attributes
    //!@{
    //! operator \f$P^{-1}AA^TP^{-T}\f$, where P is the preconditioner and A is the monolithic matrix
    boost::shared_ptr<ComposedOperator<Epetra_Operator> > M_preconditionedSymmetrizedMatrix;
paolo.crosetto's avatar
paolo.crosetto committed
601
    boost::shared_ptr<vector_Type>                    M_stress;
602
    UInt                                              M_fluxes;
paolo.crosetto's avatar
paolo.crosetto committed
603
604
    std::vector<bcName_Type>                          M_BCFluxNames;
    std::vector<UInt>                                 M_fluxOffset;
paolo.crosetto's avatar
paolo.crosetto committed
605

606
607
608
609
610
611
612
613
    //@}
};

class WRONG_PREC_EXCEPTION
{
public:

    //! Exception thrown when a wrong preconditioning strategy is set
Alessio Fumagalli's avatar
Alessio Fumagalli committed
614
615
    WRONG_PREC_EXCEPTION() {}
    virtual ~WRONG_PREC_EXCEPTION() {}
616
617
618
619

};

}
620
621


622
#endif