FSIMonolithicGI.cpp 20.8 KB
Newer Older
1
2
3
//@HEADER
/*
*******************************************************************************
4

5
6
    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7

8
    This file is part of LifeV.
9

10
11
12
13
    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.
14

15
16
17
18
    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.
19

20
21
    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.
22

23
*******************************************************************************
24
*/
25
26
//@HEADER

27

28
#include <lifev/core/LifeV.hpp>
29

Radu Popescu's avatar
Radu Popescu committed
30
31
32
33
#include <lifev/fsi/solver/FSIMonolithicGI.hpp>
#include <lifev/fsi/solver/MonolithicBlockComposedDN.hpp>
#include <lifev/fsi/solver/MonolithicBlockComposedDND.hpp>
#include <lifev/fsi/solver/MonolithicBlockMatrixRN.hpp>
34

Alessio Fumagalli's avatar
Alessio Fumagalli committed
35
36
namespace LifeV
{
37

Cristiano Malossi's avatar
Cristiano Malossi committed
38
39
40
41
// ===================================================
//  Constructors and Descructor
// ===================================================
FSIMonolithicGI::FSIMonolithicGI() :
Alessio Fumagalli's avatar
Alessio Fumagalli committed
42
43
44
45
46
47
48
    super_Type              (),
    M_mapWithoutMesh        (),
    M_uk                    (),
    M_interface             (0),
    M_meshBlock             (),
    M_shapeDerivativesBlock (),
    M_solidDerBlock         ()
Cristiano Malossi's avatar
Cristiano Malossi committed
49
50
{
}
Paolo Crosetto's avatar
   
Paolo Crosetto committed
51

Cristiano Malossi's avatar
Cristiano Malossi committed
52
53
54
// ===================================================
//  Public Methods
// ===================================================
Alessio Fumagalli's avatar
Alessio Fumagalli committed
55
void FSIMonolithicGI::setup ( const GetPot& dataFile )
Cristiano Malossi's avatar
Cristiano Malossi committed
56
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
57
    super_Type::setup ( dataFile );
Cristiano Malossi's avatar
Cristiano Malossi committed
58
}
59

Alessio Fumagalli's avatar
Alessio Fumagalli committed
60
void FSIMonolithicGI::setupFluidSolid ( UInt const fluxes )
Cristiano Malossi's avatar
Cristiano Malossi committed
61
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
62
63
64
    super_Type::setupFluidSolid ( fluxes );
    UInt offset = M_monolithicMap->map ( Unique )->NumGlobalElements();
    M_mapWithoutMesh.reset ( new MapEpetra ( *M_monolithicMap ) );
65

66
    *M_monolithicMap += M_mmFESpace->map();
67

Alessio Fumagalli's avatar
Alessio Fumagalli committed
68
    M_interface = M_monolithicMatrix->interface();
69

Alessio Fumagalli's avatar
Alessio Fumagalli committed
70
71
72
73
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
    M_beta.reset ( new vector_Type (M_uFESpace->map() ) );
    M_rhs.reset (new vector_Type (*M_monolithicMap) );
    M_rhsFull.reset (new vector_Type (*M_monolithicMap) );
    if (M_data->dataFluid()->useShapeDerivatives() )
    {
        M_shapeDerivativesBlock.reset (new matrix_Type (*M_monolithicMap) );
    }
    M_uk.reset (new vector_Type (*M_monolithicMap) );

    M_meshMotion.reset (new meshMotion_Type (*M_mmFESpace,
                                             M_epetraComm,
                                             *M_monolithicMap,
                                             offset) );

    M_fluid.reset     (new fluid_Type (M_data->dataFluid(),
                                       *M_uFESpace,
                                       *M_pFESpace,
                                       *M_mmFESpace,
                                       M_epetraComm,
                                       *M_monolithicMap,
                                       fluxes) );
    M_solid.reset (new solid_Type() );

    M_solid->setup (M_data->dataSolid(),
                    M_dFESpace,
                    M_epetraComm,
                    M_monolithicMap,
                    M_offset
                   );
99

Cristiano Malossi's avatar
Cristiano Malossi committed
100
}
101

Paolo Crosetto's avatar
Paolo Crosetto committed
102
103
void
FSIMonolithicGI::buildSystem()
Cristiano Malossi's avatar
Cristiano Malossi committed
104
{
Cristiano Malossi's avatar
Cristiano Malossi committed
105
106
    super_Type::buildSystem();
    M_meshMotion->computeMatrix();
Cristiano Malossi's avatar
Cristiano Malossi committed
107
}
108

Paolo Crosetto's avatar
Paolo Crosetto committed
109
void
Alessio Fumagalli's avatar
Alessio Fumagalli committed
110
111
112
FSIMonolithicGI::evalResidual ( vector_Type&       res,
                                const vector_Type& disp,
                                const UInt          iter )
Cristiano Malossi's avatar
Cristiano Malossi committed
113
{
114
115
116
    // disp here is the current solution guess (u,p,ds,df)

    res = 0.;//this is important. Don't remove it!
117

Alessio Fumagalli's avatar
Alessio Fumagalli committed
118
    M_uk.reset (new vector_Type ( disp ) );
119

Alessio Fumagalli's avatar
Alessio Fumagalli committed
120
    UInt offset ( M_solidAndFluidDim + nDimensions * M_interface );
121

Alessio Fumagalli's avatar
Alessio Fumagalli committed
122
123
    vectorPtr_Type meshDisp ( new vector_Type (M_mmFESpace->map() ) );
    vectorPtr_Type mmRep ( new vector_Type (M_mmFESpace->map(), Repeated ) );
124

Alessio Fumagalli's avatar
Alessio Fumagalli committed
125
    meshDisp->subset (disp, offset); //if the conv. term is to be condidered implicitly
126
    mmRep = meshDisp;
Alessio Fumagalli's avatar
Alessio Fumagalli committed
127
    moveMesh ( *mmRep );
128

Paolo Tricerri's avatar
Paolo Tricerri committed
129
    //here should use extrapolationFirstDerivative instead of velocity
130
131
132
133
    if (iter == 0)
    {
        M_ALETimeAdvance->updateRHSFirstDerivative (M_data->dataFluid()->dataTime()->timeStep() );
    }
134
    vector_Type meshVelocityRepeated ( this->M_ALETimeAdvance->firstDerivative ( *meshDisp ), Repeated );
Alessio Fumagalli's avatar
Alessio Fumagalli committed
135
    vector_Type interpolatedMeshVelocity (this->M_uFESpace->map() );
Paolo Crosetto's avatar
Paolo Crosetto committed
136

Alessio Fumagalli's avatar
Alessio Fumagalli committed
137
    interpolateVelocity ( meshVelocityRepeated, interpolatedMeshVelocity );
Paolo Crosetto's avatar
Paolo Crosetto committed
138

Alessio Fumagalli's avatar
Alessio Fumagalli committed
139
140
    vectorPtr_Type fluid ( new vector_Type ( M_uFESpace->map() ) );
    M_beta->subset ( disp, 0 );
Paolo Tricerri's avatar
Paolo Tricerri committed
141
    *M_beta -= interpolatedMeshVelocity; // convective term, u^(n+1) - w^(n+1)
142

Alessio Fumagalli's avatar
Alessio Fumagalli committed
143
144
145
    assembleSolidBlock ( iter, disp );
    assembleFluidBlock ( iter, disp );
    assembleMeshBlock ( iter );
146

147
    *M_rhsFull = *M_rhs;
Cristiano Malossi's avatar
Cristiano Malossi committed
148

Alessio Fumagalli's avatar
Alessio Fumagalli committed
149
150
    M_monolithicMatrix->setRobin ( M_robinCoupling, M_rhsFull );
    M_precPtr->setRobin (M_robinCoupling, M_rhsFull);
151

Alessio Fumagalli's avatar
Alessio Fumagalli committed
152
    if (!M_monolithicMatrix->set() )
153
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
        M_BChs.push_back (M_BCh_d);
        M_BChs.push_back (M_BCh_u);
        M_FESpaces.push_back (M_dFESpace);
        M_FESpaces.push_back (M_uFESpace);

        M_BChs.push_back (M_BCh_mesh);
        M_FESpaces.push_back (M_mmFESpace);

        M_monolithicMatrix->push_back_matrix (M_solidBlockPrec, false);
        M_monolithicMatrix->push_back_matrix (M_fluidBlock, true);
        M_monolithicMatrix->push_back_matrix (M_meshBlock, false);
        M_monolithicMatrix->setConditions (M_BChs);
        M_monolithicMatrix->setSpaces (M_FESpaces);
        M_monolithicMatrix->setOffsets (3, M_offset, 0, M_solidAndFluidDim + nDimensions * M_interface);
        M_monolithicMatrix->coupler (M_monolithicMap, M_dofStructureToFluid->localDofMap(), M_numerationInterface, M_data->dataFluid()->dataTime()->timeStep(), M_solidTimeAdvance->coefficientFirstDerivative ( 0 ), M_solid->rescaleFactor() );
        M_monolithicMatrix->coupler ( M_monolithicMap, M_dofStructureToFluid->localDofMap(), M_numerationInterface, M_data->dataFluid()->dataTime()->timeStep(), M_solidTimeAdvance->coefficientFirstDerivative ( 0 ), M_solid->rescaleFactor(), 2);
170
171
172
    }
    else
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
173
174
175
        M_monolithicMatrix->replace_matrix (M_solidBlockPrec, 0);
        M_monolithicMatrix->replace_matrix (M_fluidBlock, 1);
        M_monolithicMatrix->replace_matrix (M_meshBlock, 2);
176
    }
177

178
    M_monolithicMatrix->blockAssembling();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
179
    super_Type::checkIfChangedFluxBC ( M_monolithicMatrix );
180
181
182



Alessio Fumagalli's avatar
Alessio Fumagalli committed
183
184
    if ( (M_data->dataSolid()->solidType().compare ("exponential") && M_data->dataSolid()->solidType().compare ("neoHookean") ) )
    {
185
        applyBoundaryConditions();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
186
    }
187

Paolo Crosetto's avatar
Paolo Crosetto committed
188
    M_monolithicMatrix->GlobalAssemble();
189

Alessio Fumagalli's avatar
Alessio Fumagalli committed
190
    super_Type::evalResidual ( disp, M_rhsFull, res, false );
Paolo Crosetto's avatar
   
Paolo Crosetto committed
191

Alessio Fumagalli's avatar
Alessio Fumagalli committed
192
    if ( ! ( M_data->dataSolid()->solidType().compare ( "exponential" ) && M_data->dataSolid()->solidType().compare ( "neoHookean" ) ) )
Cristiano Malossi's avatar
Cristiano Malossi committed
193
194
    {
        res += *M_meshBlock * disp;
Paolo Crosetto's avatar
   
Paolo Crosetto committed
195

Cristiano Malossi's avatar
Cristiano Malossi committed
196
        if ( !M_BCh_u->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
197
198
199
200
        {
            M_BCh_u->bcUpdate ( *M_uFESpace->mesh(), M_uFESpace->feBd(), M_uFESpace->dof() );
        }
        M_BCh_d->setOffset ( M_offset );
Cristiano Malossi's avatar
Cristiano Malossi committed
201
        if ( !M_BCh_d->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
202
203
204
205
        {
            M_BCh_d->bcUpdate ( *M_dFESpace->mesh(), M_dFESpace->feBd(), M_dFESpace->dof() );
        }
        M_BCh_mesh->setOffset ( M_solidAndFluidDim + nDimensions * M_interface );
Cristiano Malossi's avatar
Cristiano Malossi committed
206
        if ( !M_BCh_mesh->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
207
208
209
        {
            M_BCh_mesh->bcUpdate ( *M_mmFESpace->mesh(), M_mmFESpace->feBd(), M_mmFESpace->dof() );
        }
Paolo Crosetto's avatar
   
Paolo Crosetto committed
210

Alessio Fumagalli's avatar
Alessio Fumagalli committed
211
        M_monolithicMatrix->applyBoundaryConditions ( dataFluid()->dataTime()->time() /*, M_rhsFull*/);
Cristiano Malossi's avatar
Cristiano Malossi committed
212
        M_monolithicMatrix->GlobalAssemble();
213

Alessio Fumagalli's avatar
Alessio Fumagalli committed
214
215
        bcManageResidual ( res, *M_rhsFull, disp, *M_uFESpace->mesh(), M_uFESpace->dof(), *M_BCh_u,
                           M_uFESpace->feBd(), M_data->dataFluid()->dataTime()->time(), 1. );
216

Cristiano Malossi's avatar
Cristiano Malossi committed
217
        // below sol is repeated by BCManageResidual
Alessio Fumagalli's avatar
Alessio Fumagalli committed
218
219
220
221
        bcManageResidual ( res, *M_rhsFull, disp, *M_dFESpace->mesh(), M_dFESpace->dof(), *M_BCh_d,
                           M_dFESpace->feBd(), M_data->dataSolid()->dataTime()->time(), 1. );
        bcManageResidual ( res, *M_rhsFull, disp, *M_mmFESpace->mesh(), M_mmFESpace->dof(), *M_BCh_mesh,
                           M_mmFESpace->feBd(), M_data->dataFluid()->dataTime()->time(), 1. );
Cristiano Malossi's avatar
Cristiano Malossi committed
222
223
224
        res -= *M_rhsFull;
    }
}
Paolo Crosetto's avatar
   
Paolo Crosetto committed
225

Cristiano Malossi's avatar
Cristiano Malossi committed
226
227
void FSIMonolithicGI::applyBoundaryConditions()
{
228
    if ( !M_BCh_u->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
229
230
231
232
        M_BCh_u->bcUpdate ( *M_uFESpace->mesh(),
                            M_uFESpace->feBd(),
                            M_uFESpace->dof() );
    M_BCh_d->setOffset ( M_offset );
233
    if ( !M_BCh_d->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
234
235
236
237
        M_BCh_d->bcUpdate ( *M_dFESpace->mesh(),
                            M_dFESpace->feBd(),
                            M_dFESpace->dof() );
    M_BCh_mesh->setOffset ( M_solidAndFluidDim + nDimensions * M_interface );
238
    if ( !M_BCh_mesh->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
239
240
241
        M_BCh_mesh->bcUpdate ( *M_mmFESpace->mesh(),
                               M_mmFESpace->feBd(),
                               M_mmFESpace->dof() );
242

Alessio Fumagalli's avatar
Alessio Fumagalli committed
243
    M_monolithicMatrix->applyBoundaryConditions (dataFluid()->dataTime()->time(), M_rhsFull);
Cristiano Malossi's avatar
Cristiano Malossi committed
244
}
245

Alessio Fumagalli's avatar
Alessio Fumagalli committed
246
247
248
void FSIMonolithicGI::setALEVectorInStencil ( const vectorPtr_Type& fluidDisp,
                                              const UInt iter,
                                              const bool lastVector)
Paolo Tricerri's avatar
Paolo Tricerri committed
249
250
{

251
252
253
    //The fluid timeAdvance has size = orderBDF because it is seen as an equation frist order in time.
    //We need to add the solidVector to the fluidVector in the fluid TimeAdvance because we have the
    //extrapolation on it.
Alessio Fumagalli's avatar
Alessio Fumagalli committed
254
255
256
    if ( ( iter < M_fluidTimeAdvance->size() - 1 ) && !lastVector )
    {
        vectorPtr_Type harmonicSolutionRestartTime (new vector_Type ( *M_monolithicMap, Unique, Zero ) );
Paolo Tricerri's avatar
Paolo Tricerri committed
257

Alessio Fumagalli's avatar
Alessio Fumagalli committed
258
        *harmonicSolutionRestartTime *= 0.0;
Paolo Tricerri's avatar
Paolo Tricerri committed
259

Alessio Fumagalli's avatar
Alessio Fumagalli committed
260
261
        UInt givenOffset ( M_solidAndFluidDim + nDimensions * M_interface );
        harmonicSolutionRestartTime->subset (*fluidDisp, fluidDisp->map(), (UInt) 0, givenOffset );
Paolo Tricerri's avatar
Paolo Tricerri committed
262

Alessio Fumagalli's avatar
Alessio Fumagalli committed
263
264
265
        //We sum the vector in the first element of fluidtimeAdvance
        * ( M_fluidTimeAdvance->stencil() [ iter + 1 ] ) += *harmonicSolutionRestartTime;
    }
Paolo Tricerri's avatar
Paolo Tricerri committed
266

Alessio Fumagalli's avatar
Alessio Fumagalli committed
267
    if ( !lastVector )
268
269
270
    {
        //The shared_pointer for the vectors has to be trasformed into a pointer to VectorEpetra
        //That is the type of pointers that are used in TimeAdvance
Alessio Fumagalli's avatar
Alessio Fumagalli committed
271
272
        vector_Type* normalPointerToALEVector ( new vector_Type (*fluidDisp) );
        (M_ALETimeAdvance->stencil() ).push_back ( normalPointerToALEVector );
273
274
275
    }
    else
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
276
        vectorPtr_Type harmonicSolutionRestartTime (new vector_Type ( *M_monolithicMap, Unique, Zero ) );
277
278
279

        *harmonicSolutionRestartTime *= 0.0;

Alessio Fumagalli's avatar
Alessio Fumagalli committed
280
281
        UInt givenOffset ( M_solidAndFluidDim + nDimensions * M_interface );
        harmonicSolutionRestartTime->subset (*fluidDisp, fluidDisp->map(), (UInt) 0, givenOffset );
282
283

        //We sum the vector in the first element of fluidtimeAdvance
Alessio Fumagalli's avatar
Alessio Fumagalli committed
284
        * ( M_fluidTimeAdvance->stencil() [ 0 ] ) += *harmonicSolutionRestartTime;
285
    }
Paolo Tricerri's avatar
Paolo Tricerri committed
286
287
288
}


289

Cristiano Malossi's avatar
Cristiano Malossi committed
290
//============ Protected Methods ===================
291

Cristiano Malossi's avatar
Cristiano Malossi committed
292
293
void FSIMonolithicGI::setupBlockPrec()
{
paolo.crosetto's avatar
paolo.crosetto committed
294
295
296
297

    //The following part accounts for a possibly nonlinear structure model, should not be run when linear
    //elasticity is used

Alessio Fumagalli's avatar
Alessio Fumagalli committed
298
299
    if ( M_data->dataSolid()->getUseExactJacobian() && ( M_data->dataSolid()->solidType().compare ( "exponential" )
                                                         && M_data->dataSolid()->solidType().compare ( "neoHookean" ) ) )
Cristiano Malossi's avatar
Cristiano Malossi committed
300
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
301
302
303
304
305
306
        M_solid->material()->updateJacobianMatrix ( *M_uk * M_solid->rescaleFactor(),
                                                    dataSolid(),
                                                    M_solid->mapMarkersVolumes(),
                                                    M_solid->displayerPtr() ); // computing the derivatives if nonlinear (comment this for inexact Newton);
        M_solidBlockPrec.reset ( new matrix_Type ( *M_monolithicMap,
                                                   1 ) );
307
        *M_solidBlockPrec += *M_solid->massMatrix();
Paolo Tricerri's avatar
Paolo Tricerri committed
308
        *M_solidBlockPrec += *M_solid->material()->jacobian();
paolo.crosetto's avatar
paolo.crosetto committed
309
        M_solidBlockPrec->globalAssemble();
310
        *M_solidBlockPrec *= M_solid->rescaleFactor();
paolo.crosetto's avatar
paolo.crosetto committed
311

Alessio Fumagalli's avatar
Alessio Fumagalli committed
312
        M_monolithicMatrix->replace_matrix ( M_solidBlockPrec, 0 );
Cristiano Malossi's avatar
Cristiano Malossi committed
313
    }
314

Paolo Crosetto's avatar
   
Paolo Crosetto committed
315
    M_monolithicMatrix->blockAssembling();
316

Cristiano Malossi's avatar
Cristiano Malossi committed
317
318
    if ( M_data->dataFluid()->useShapeDerivatives() )
    {
crosetto's avatar
crosetto committed
319
        *M_shapeDerivativesBlock *= 0.;
Cristiano Malossi's avatar
Cristiano Malossi committed
320
        M_shapeDerivativesBlock->openCrsMatrix();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
321
        shapeDerivatives ( M_shapeDerivativesBlock );
Cristiano Malossi's avatar
Cristiano Malossi committed
322
        M_shapeDerivativesBlock->globalAssemble();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
323
        M_monolithicMatrix->addToGlobalMatrix ( M_shapeDerivativesBlock );
Cristiano Malossi's avatar
Cristiano Malossi committed
324
    }
crosetto's avatar
crosetto committed
325

paolo.crosetto's avatar
paolo.crosetto committed
326
    if ( M_data->dataFluid()->useShapeDerivatives() || M_data->dataSolid()->getUseExactJacobian() )
Cristiano Malossi's avatar
Cristiano Malossi committed
327
    {
paolo.crosetto's avatar
paolo.crosetto committed
328
        if ( !M_BCh_u->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
329
330
331
332
            M_BCh_u->bcUpdate ( *M_uFESpace->mesh(),
                                M_uFESpace->feBd(),
                                M_uFESpace->dof() );
        M_BCh_d->setOffset ( M_offset );
paolo.crosetto's avatar
paolo.crosetto committed
333
        if ( !M_BCh_d->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
334
335
336
337
            M_BCh_d->bcUpdate ( *M_dFESpace->mesh(),
                                M_dFESpace->feBd(),
                                M_dFESpace->dof() );
        M_BCh_mesh->setOffset ( M_solidAndFluidDim + nDimensions * M_interface );
paolo.crosetto's avatar
paolo.crosetto committed
338
        if ( !M_BCh_mesh->bcUpdateDone() )
Alessio Fumagalli's avatar
Alessio Fumagalli committed
339
340
341
            M_BCh_mesh->bcUpdate ( *M_mmFESpace->mesh(),
                                   M_mmFESpace->feBd(),
                                   M_mmFESpace->dof() );
paolo.crosetto's avatar
paolo.crosetto committed
342

Alessio Fumagalli's avatar
Alessio Fumagalli committed
343
        M_monolithicMatrix->applyBoundaryConditions ( dataFluid()->dataTime()->time() );
344
        M_monolithicMatrix->GlobalAssemble();
Cristiano Malossi's avatar
Cristiano Malossi committed
345
    }
346

Cristiano Malossi's avatar
Cristiano Malossi committed
347
    super_Type::setupBlockPrec();
paolo.crosetto's avatar
paolo.crosetto committed
348

Cristiano Malossi's avatar
Cristiano Malossi committed
349
350
    if ( M_precPtr->blockVector().size() < 3 )
    {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
351
352
353
354
355
356
357
358
        M_precPtr->push_back_matrix ( M_meshBlock,
                                      false );
        M_precPtr->setConditions ( M_BChs );
        M_precPtr->setSpaces ( M_FESpaces );
        M_precPtr->setOffsets ( 3, M_offset, 0,  M_solidAndFluidDim + nDimensions * M_interface );
        M_precPtr->coupler ( M_monolithicMap, M_dofStructureToFluid->localDofMap(), M_numerationInterface, M_data->dataFluid()->dataTime()->timeStep() , M_solidTimeAdvance->coefficientFirstDerivative ( 0 ), M_solid->rescaleFactor(), 2);

        if (M_data->dataFluid()->useShapeDerivatives() )
Cristiano Malossi's avatar
Cristiano Malossi committed
359
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
360
361
            boost::shared_ptr<MatrixEpetra<Real> > staticCast = boost::static_pointer_cast<MatrixEpetra<Real> > (M_shapeDerivativesBlock);
            M_precPtr->push_back_coupling ( staticCast );
Cristiano Malossi's avatar
Cristiano Malossi committed
362
363
        }
    }
uvilla's avatar
uvilla committed
364
    else
Cristiano Malossi's avatar
Cristiano Malossi committed
365
    {
366
367
        //M_precPtr->replace_matrix( M_solidBlockPrec, 0 );
        //M_precPtr->replace_matrix( M_fluidBlock, 1 );
Alessio Fumagalli's avatar
Alessio Fumagalli committed
368
        M_precPtr->replace_matrix ( M_meshBlock, 2 );
369

Cristiano Malossi's avatar
Cristiano Malossi committed
370
371
        if ( M_data->dataFluid()->useShapeDerivatives() )
        {
Alessio Fumagalli's avatar
Alessio Fumagalli committed
372
            M_precPtr->replace_coupling ( M_shapeDerivativesBlock, 2 );
Cristiano Malossi's avatar
Cristiano Malossi committed
373
374
375
        }
    }
}
376

Antonio Cervone's avatar
Antonio Cervone committed
377
void FSIMonolithicGI::shapeDerivatives ( FSIOperator::fluid_Type::matrixPtr_Type sdMatrix )
Cristiano Malossi's avatar
Cristiano Malossi committed
378
{
Alessio Fumagalli's avatar
Alessio Fumagalli committed
379
380
381
382
    Real alpha = M_fluidTimeAdvance->coefficientFirstDerivative ( 0 ) / M_data->dataFluid()->dataTime()->timeStep();
    vectorPtr_Type rhsNew (new vector_Type (*M_monolithicMap) );
    vector_Type un (M_uFESpace->map() );
    vector_Type uk (M_uFESpace->map() + M_pFESpace->map() );
383

Alessio Fumagalli's avatar
Alessio Fumagalli committed
384
    vectorPtr_Type meshVelRep ( new vector_Type ( M_mmFESpace->map(), Repeated ) );
385

386
    *meshVelRep = M_ALETimeAdvance->firstDerivative();
387

Paolo Tricerri's avatar
Paolo Tricerri committed
388
    //When this class is used, the convective term is used implictly
Alessio Fumagalli's avatar
Alessio Fumagalli committed
389
    un.subset ( *M_uk, 0 );
Cristiano Malossi's avatar
Cristiano Malossi committed
390

Alessio Fumagalli's avatar
Alessio Fumagalli committed
391
392
393
    uk.subset ( *M_uk, 0 );
    vector_Type veloFluidMesh ( M_uFESpace->map(), Repeated );
    this->transferMeshMotionOnFluid ( *meshVelRep, veloFluidMesh );
Cristiano Malossi's avatar
Cristiano Malossi committed
394

Paolo Tricerri's avatar
Paolo Tricerri committed
395
396
    //The last two flags are consistent with the currect interface.
    //When this class is used, they should not be changed.
grandper's avatar
grandper committed
397
398
399
400
401
402
403
404
    M_fluid->updateShapeDerivatives ( *sdMatrix, alpha,
                                      un,
                                      uk,
                                      veloFluidMesh, //(xk-xn)/dt (FI), or (xn-xn-1)/dt (CE)//Repeated
                                      M_solidAndFluidDim + M_interface * nDimensions,
                                      *M_mmFESpace,
                                      true /*This flag tells the method to consider the velocity of the domain implicitly*/,
                                      true /*This flag tells the method to consider the convective term implicitly */ );
Cristiano Malossi's avatar
Cristiano Malossi committed
405
}
406

grandper's avatar
grandper committed
407
void FSIMonolithicGI::assembleMeshBlock ( UInt /*iter*/ )
Cristiano Malossi's avatar
Cristiano Malossi committed
408
{
409

Alessio Fumagalli's avatar
Alessio Fumagalli committed
410
411
    M_meshBlock.reset ( new matrix_Type ( *M_monolithicMap ) );
    M_meshMotion->addSystemMatrixTo ( M_meshBlock );
412
    M_meshBlock->globalAssemble();
Alessio Fumagalli's avatar
Alessio Fumagalli committed
413
    UInt offset ( M_solidAndFluidDim + nDimensions * M_interface );
Cristiano Malossi's avatar
Cristiano Malossi committed
414
    std::map< ID, ID >::const_iterator ITrow;
Davide Forti's avatar
Davide Forti committed
415
    std::map< ID, ID > locdofmap ( M_dofStructureToFluid->localDofMap() );
416

crosetto's avatar
crosetto committed
417
    /******************alternative way************************/
Paolo Crosetto's avatar
   
Paolo Crosetto committed
418
419
420
421
    //     BCFunctionBase bcf(fZero);
    //     fluidBchandlerPtr_Type BCh(new fluidBchandler_Type() );
    //     BCh->addBC("Interface", 1, Essential, Full,
    //                bcf, 3);
crosetto's avatar
crosetto committed
422

Paolo Crosetto's avatar
   
Paolo Crosetto committed
423
    //     BCh->setOffset(M_solidAndFluidDim + nDimensions*M_interface);
crosetto's avatar
crosetto committed
424

Paolo Crosetto's avatar
   
Paolo Crosetto committed
425
426
    //     if ( !BCh->bcUpdateDone() )
    //         BCh->bcUpdate( *M_mmFESpace->mesh(), M_mmFESpace->feBd(), M_mmFESpace->dof() );
crosetto's avatar
crosetto committed
427

Paolo Crosetto's avatar
   
Paolo Crosetto committed
428
    //     bcManage( *M_meshBlock, *M_rhsFull, *M_mmFESpace->mesh(), M_mmFESpace->dof(), *BCh, M_mmFESpace->feBd(), 1., dataFluid()->dataTime()->time());
uvilla's avatar
uvilla committed
429
    /********************************************************/
crosetto's avatar
crosetto committed
430

Davide Forti's avatar
Davide Forti committed
431
432
    matrixPtr_Type diagFilter (new matrix_Type (*M_monolithicMap) );
    diagFilter->insertZeroDiagonal ();
Cristiano Malossi's avatar
Cristiano Malossi committed
433
    for ( ID dim = 0; dim < nDimensions; ++dim )
Davide Forti's avatar
Davide Forti committed
434
    {
Cristiano Malossi's avatar
Cristiano Malossi committed
435
        for ( ITrow = locdofmap.begin(); ITrow != locdofmap.end(); ++ITrow )
436
        {
Davide Forti's avatar
Davide Forti committed
437
438
            int i = ITrow->first;
            int iRow =  i + offset + dim * M_mmFESpace->dof().numTotalDof();
grandper's avatar
grandper committed
439
            if ( M_meshBlock->mapPtr()->map (Unique)->MyGID (iRow) )
Davide Forti's avatar
Davide Forti committed
440
441
442
443
444
445
446
447
448
            {
                M_meshBlock->diagonalize ( iRow, 1. );
            }
            else
            {
                double myValues[1];
                myValues[0] = -1;
                int myIndices[1];
                myIndices[0] = iRow;
grandper's avatar
grandper committed
449
                diagFilter->matrixPtr()->SumIntoGlobalValues ( iRow, 1, myValues, myIndices );
Davide Forti's avatar
Davide Forti committed
450
            }
451
        }
Davide Forti's avatar
Davide Forti committed
452
453
454
455
456
457
    }

    diagFilter->globalAssemble();
    // Processor informations
    Int  numLocalEntries = diagFilter->matrixPtr()->RowMap().NumMyElements();
    Int* globalEntries   = diagFilter->matrixPtr()->RowMap().MyGlobalElements();
grandper's avatar
grandper committed
458
    UInt globalRowIndex (0);
Davide Forti's avatar
Davide Forti committed
459
460

    // Source informations handlers
grandper's avatar
grandper committed
461
    double* diagValue;
Davide Forti's avatar
Davide Forti committed
462
463
    Int* indices;
    Int srcRow (0);
grandper's avatar
grandper committed
464
    Int controlValue (0); // This value should be one since it is a diagonal matrix
Davide Forti's avatar
Davide Forti committed
465
466
467
468
469
470
471
472
473
474

    for (Int i (0); i < numLocalEntries; ++i)
    {
        // Collecting the data from the source
        globalRowIndex = globalEntries[i];

        // Get the data of the row
        srcRow = diagFilter->matrixPtr()->LRID (globalRowIndex);
        diagFilter->matrixPtr()->ExtractMyRowView (srcRow, controlValue, diagValue, indices);

grandper's avatar
grandper committed
475
476
        ASSERT ( controlValue == 1, "Error: The matrix should be diagonal.");
        if (diagValue[0] < 0.0)
Davide Forti's avatar
Davide Forti committed
477
478
479
480
        {
            M_meshBlock->diagonalize ( globalRowIndex, 1. );
        }
    }
Cristiano Malossi's avatar
Cristiano Malossi committed
481
}
482

Cristiano Malossi's avatar
Cristiano Malossi committed
483
484
485
// ===================================================
//  Products registration
// ===================================================
Alessio Fumagalli's avatar
Alessio Fumagalli committed
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
bool FSIMonolithicGI::S_register = BlockPrecFactory::instance().registerProduct ( "AdditiveSchwarzGI",
                                   &MonolithicBlockMatrix::createAdditiveSchwarz )
                                   && BlockPrecFactory::instance().registerProduct ( "ComposedDNGI",
                                           &MonolithicBlockComposedDN::createComposedDNGI )
                                   && MonolithicBlockMatrix::Factory_Type::instance().registerProduct ( "AdditiveSchwarzGI",
                                           &MonolithicBlockMatrix::createAdditiveSchwarz )
                                   && MonolithicBlockMatrix::Factory_Type::instance().registerProduct ( "AdditiveSchwarzRNGI",
                                           &MonolithicBlockMatrixRN::createAdditiveSchwarzRN )
                                   && FSIOperator::FSIFactory_Type::instance().registerProduct ( "monolithicGI",
                                           &FSIMonolithicGI::instantiate )
                                   && BlockPrecFactory::instance().registerProduct ( "ComposedDNDGI",
                                           &MonolithicBlockComposedDND::createComposedDNDGI )
                                   && BlockPrecFactory::instance().registerProduct ( "ComposedDND2GI",
                                           &MonolithicBlockComposedDND::createComposedDND2GI )
                                   && BlockPrecFactory::instance().registerProduct ( "ComposedDN2GI",
                                           &MonolithicBlockComposedDN::createComposedDN2GI );
crosetto's avatar
crosetto committed
502

503
}