AdvectionDiffusionReactionSolver.hpp 44.8 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
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
51
52
53
54
55
56
57
/* -*- mode: c++ -*-

 This file is part of the LifeV library

 Author(s): Miguel A. Fernandez <miguel.fernandez@inria.fr>
            Christoph Winkelmann <christoph.winkelmann@epfl.ch>
      Date: 2003-06-09

 Copyright (C) 2004 EPFL

 This library 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 2.1 of the License, or (at your option) any later version.

 This library 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 this library; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/
/**
  \file Oseen.hpp
  \author G. Fourestey
  \date 2/2007

  \brief This file contains a Advection Diffusion Reaction equation solver class
*/
#ifndef _ADR_H_
#define _ADR_H_

#include <life/lifearray/elemMat.hpp>
#include <life/lifearray/elemVec.hpp>
#include <life/lifefem/elemOper.hpp>
#include <life/lifefem/values.hpp>
#include <life/lifearray/pattern.hpp>
#include <life/lifefem/assemb.hpp>
#include <life/lifefem/bcManage.hpp>
#include <life/lifefilters/medit_wrtrs.hpp>

#include <life/lifealg/SolverTrilinos.hpp>
#include <life/lifealg/EpetraMap.hpp>
#include <life/lifearray/EpetraMatrix.hpp>
#include <life/lifearray/EpetraVector.hpp>
//
#include <life/lifefem/bcHandler.hpp>
#include <life/lifecore/chrono.hpp>
#include <life/lifefem/sobolevNorms.hpp>
#include <life/lifefem/geoMap.hpp>
#include <life/lifesolver/nsipterms.hpp>
#include <life/lifesolver/dataADR.hpp>
//
#include <boost/shared_ptr.hpp>

58
#include <life/lifefem/FESpace.hpp>
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

#include <list>

namespace LifeV
{
/*!
  \class ADRSolver

  This class implements a NavierStokes solver.
  The resulting linear systems are solved by GMRES on the full
  matrix ( u and p coupled ).

*/



template< typename Mesh,
76
          typename SolverType = LifeV::SolverTrilinos >
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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
class ADRSolver
//     :
//     public NavierStokesHandler<Mesh>, EpetraHandler
{

public:

    typedef DataADR<Mesh>  data_type;

    typedef Real ( *Function ) ( const Real&, const Real&, const Real&,const Real&, const ID& );
    typedef boost::function<Real ( Real const&, Real const&, Real const&,
                                   Real const&, ID const& )> source_type;

    typedef Mesh mesh_type;

    typedef BCHandler                             bchandler_raw_type;
    typedef boost::shared_ptr<bchandler_raw_type> bchandler_type;

//    typedef SolverType                    solver_type;
    typedef typename SolverType::matrix_type      matrix_type;
    typedef boost::shared_ptr<matrix_type>        matrix_ptrtype;
    typedef typename SolverType::vector_type      vector_type;

    typedef typename SolverType::prec_raw_type    prec_raw_type;
    typedef typename SolverType::prec_type        prec_type;


    //! Constructor
    /*!
      \param dataType
      \param localMap
      \param velocity FE space
      \param pressure FE space
      \param bcHandler boundary conditions for the velocity
    */
    ADRSolver( const data_type&          dataType,
               FESpace<Mesh, EpetraMap>& FESpace,
               FESpace<Mesh, EpetraMap>& betaFESpace,
               BCHandler&                bcHandler,
               Epetra_Comm&              comm );

    /*!
      \param dataType
      \param localMap
      \param velocity FE space
      \param pressure FE space
    */
    ADRSolver( const data_type&      dataType,
               FESpace<Mesh, EpetraMap>& FESpace,
               FESpace<Mesh, EpetraMap>& betaFESpace,
               Epetra_Comm&              comm );

    //! virtual destructor

    virtual ~ADRSolver();

    //! Update convective term, bc treatment and solve the linearized ns system

    virtual void iterate( bchandler_raw_type& bch );

    virtual void setUp        ( const GetPot& dataFile );

    virtual void buildSystem();

    virtual void updateRHS(vector_type const& sourceVec)
    {
        M_rhsNoBC = sourceVec;
        M_rhsNoBC.GlobalAssemble();
    }

147
    virtual void updateSystem(Real       alpha,
148
149
150
151
                              vector_type& betaVec,
                              vector_type& sourceVec
                              );

152
//     void updateLinearSystem(Real       alpha,
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
//                             vector_type& betaVec,
//                             vector_type& sourceVec
//                             );

    void initialize( const Function& );
    void initialize( const vector_type& );

    //! returns the local solution vector

    const vector_type& solution() const {return M_sol;}

    //! returns the local residual vector

    const vector_type& residual() const {return M_residual;}

    //! reduce the local solution solution in global vectors

    void reduceSolution( Vector& u,
                         Vector& p );

    void reduceResidual( Vector& res);


    FESpace<Mesh, EpetraMap>& velFESpace()   {return M_FESpace;}


    //! Bounday Conditions
    /*const*/ bool BCset() const {return M_setBC;}
    //! set the BCs
    void setBC(BCHandler &BCh)
        {
            M_BCh = &BCh; M_setBC = true;
        }

//     BCHandler& bcHandler()
//         {
//             return *M_BCh;
//         }

    //! set the source term functor
    void setSourceTerm( source_type __s )
        {
            M_source = __s;
        }

    //! get the source term functor
    source_type sourceTerm() const
        {
            return M_source;
        }

    // compute area on a face with given flag
    Real area(const EntityFlag& flag);

    // compute flux on a face with given flag
    Real flux(const EntityFlag& flag);

    //! Postprocessing
    void postProcess(bool _writeMesh = false);

213
    void resetPrec(bool reset = true) { if (reset) M_linearSolver.precReset(); }
214
215
216
217
218
219
220
221
222
223
224

    EpetraMap const& getMap() const { return M_localMap; }

    const Epetra_Comm& comm() const {return *M_comm;}

    bool isLeader() const
    {
        assert( M_comm != 0);
        return comm().MyPID() == 0;
    }

225
    void leaderPrint   (string const message, Real const number) const;
226
    void leaderPrint   (string const message) const;
227
    void leaderPrintMax(string const message, Real const number) const;
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314


    void recomputeMatrix(bool const recomp){M_recomputeMatrix = recomp;}

    matrix_type& matrNoBC()
        {
            return *M_matrNoBC;
        }
    matrix_type& matrMass()
        {
            return *M_matrMass;
        }


protected:

    UInt dim() const           { return M_FESpace.dim(); }

    void applyBoundaryConditions(  matrix_type&        matrix,
                                   vector_type&        rhs,
                                   bchandler_raw_type& BCh);

    void echo(std::string message);

    //! removes mean of component comp of vector x
    Real removeMean( vector_type& x );

    //private members

    //! data for NS solvers
    const data_type&               M_data;

    // FE spaces
    FESpace<Mesh, EpetraMap>&      M_FESpace;

    // beta FE spaces
    FESpace<Mesh, EpetraMap>&      M_betaFESpace;

    //! MPI communicator
    Epetra_Comm*                   M_comm;
    int                            M_me;

    //! Bondary Conditions Handler
    BCHandler*                     M_BCh;
    bool                           M_setBC;

    EpetraMap                      M_localMap;

    //! mass matrix

    matrix_ptrtype                 M_matrMass;

    //! Stiffness Matrix mu*stiff
    matrix_ptrtype                 M_matrStiff;

    //! Advection Matrix
    matrix_ptrtype                 M_matrAdv;

    //! matrix to be solved
//    matrix_ptrtype                 M_matrFull;

    //! matrix without boundary conditions

    matrix_ptrtype                 M_matrNoBC;

    //! stabilization matrix
    matrix_ptrtype                 M_matrStab;

    //! source term for NS
    source_type                    M_source;


    //! Right hand side for the velocity
    vector_type                    M_rhsNoBC;

    //! Right hand side global
    vector_type                    M_rhsFull;

    //! Global solution _u and _p
    vector_type                    M_sol;

    //! residual

    vector_type                    M_residual;

    SolverType                     M_linearSolver;

gfourestey's avatar
gfourestey committed
315
     boost::shared_ptr<EpetraPreconditioner> M_prec;
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330

    bool                           M_steady;

    //! Stabilization
    std::string                    M_stab;
    Real                           M_gammaBeta;

//     details::
//    IPStabilization<Mesh, Dof>     M_ipStab;

    const Function*                M_betaFct;

    bool                           M_divBetaUv;

    //
331
    Real                         M_diagonalize;
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398

    UInt                           M_count;

    //! boolean that indicates if output is sent to cout

    bool                           M_verbose;

    //! boolean that indicates if the matrix is updated for the current iteration

    bool                           M_updated;

    //! boolean that indicates if le precond has to be recomputed

    bool                           M_reusePrec;

    //! interger storing the max number of solver iteration with prec recomputing

    int                            M_maxIterSolver;

    //!

    bool                           M_recomputeMatrix;

private:

    //! Elementary matrices and vectors
    ElemMat                        M_elmatStiff;      // stiffness
    ElemMat                        M_elmatAdv;        // advection
    ElemMat                        M_elmatMass;       // mass
    ElemMat                        M_elmatStab;
//    ElemVec                        M_elvec;           // Elementary right hand side
    ElemVec                        M_elvec_u;         // Elementary right hand side


}; // class ADRSolver



//
// IMPLEMENTATION
//


template<typename Mesh, typename SolverType>
ADRSolver<Mesh, SolverType>::
ADRSolver( const data_type&          dataType,
           FESpace<Mesh, EpetraMap>& FESpace,
           FESpace<Mesh, EpetraMap>& betaFESpace,
           BCHandler&                BCh,
           Epetra_Comm&              comm ):
    M_data                   ( dataType ),
    M_FESpace                ( FESpace ),
    M_betaFESpace            ( betaFESpace ),
    M_comm                   ( &comm ),
    M_me                     ( M_comm->MyPID() ),
    M_BCh                    ( &BCh ),
    M_setBC                  ( true ),
    M_localMap               ( M_FESpace.map() ),
    M_matrMass               ( ),
    M_matrStiff              ( ),
//    M_matrFull               ( ),
    M_matrNoBC               ( ),
    M_rhsNoBC                ( M_localMap ),
    M_rhsFull                ( M_localMap ),
    M_sol                    ( M_localMap ),
    M_residual               ( M_localMap ),
    M_linearSolver           ( ),
399
    M_prec                   ( ),
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
//     M_ipStab                 ( M_FESpace.mesh(),
//                                M_FESpace.dof(), M_FESpace.refFE(),
//                                M_FESpace.feBd(), M_FESpace.qr(),
//                                0., 0., 0.,
//                                M_data.diffusivity() ),
    M_betaFct                ( 0 ),
    M_count                  ( 0 ),
    M_verbose                ( M_me == 0),
    M_updated                ( false ),
    M_reusePrec              ( true ),
    M_maxIterSolver          ( -1 ),
    M_recomputeMatrix        ( false ),
    M_elmatStiff             ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatMass              ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatAdv               ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatStab              ( M_FESpace.fe().nbNode, 1, 1 ),
//    M_elvec                  ( M_FESpace.fe().nbNode, nDimensions ),
quinodoz's avatar
   
quinodoz committed
417
    M_elvec_u                ( M_betaFESpace.fe().nbNode, nDimensions ) //SQ: from M_FESpace to M_betaFESpace
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
{
}



template<typename Mesh, typename SolverType>
ADRSolver<Mesh, SolverType>::
ADRSolver( const data_type&          dataType,
           FESpace<Mesh, EpetraMap>& FESpace,
           FESpace<Mesh, EpetraMap>& betaFESpace,
           Epetra_Comm&              comm ):
    M_data                   ( dataType ),
    M_FESpace                ( FESpace ),
    M_betaFESpace            ( betaFESpace ),
    M_comm                   ( &comm ),
    M_me                     ( M_comm->MyPID() ),
    M_setBC                  ( false ),
    M_localMap               ( M_FESpace.map() ),
    M_matrMass               ( ),
gfourestey's avatar
gfourestey committed
437
    M_matrStiff              ( ),
438
439
440
441
442
443
    M_matrNoBC               ( ),
    M_rhsNoBC                ( M_localMap ),
    M_rhsFull                ( M_localMap ),
    M_sol                    ( M_localMap ),
    M_residual               ( M_localMap ),
    M_linearSolver           ( ),
444
    M_prec                   ( ),
445
446
447
448
449
450
451
452
453
454
455
    M_betaFct                ( 0 ),
    M_count                  ( 0 ),
    M_verbose                ( M_me == 0),
    M_updated                ( false ),
    M_reusePrec              ( true ),
    M_maxIterSolver          ( -1 ),
    M_recomputeMatrix        ( false ),
    M_elmatStiff             ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatMass              ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatAdv               ( M_FESpace.fe().nbNode, 1, 1 ),
    M_elmatStab              ( M_FESpace.fe().nbNode, 1, 1 ),
quinodoz's avatar
   
quinodoz committed
456
    M_elvec_u                ( M_betaFESpace.fe().nbNode, nDimensions ) //SQ: from M_FESpace to M_betaFESpace
457
458
459
460
461
462
463
464
465
466
467
468
469
470
{
}



template<typename Mesh, typename SolverType>
ADRSolver<Mesh, SolverType>::
~ADRSolver()
{

}

template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
471
leaderPrint(string const message, Real const number) const
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
{
  if ( isLeader() )
    std::cout << message << number << std::endl;

}

template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
leaderPrint(string const message) const
{
  if ( isLeader() )
    std::cout << message << std::flush;

}

template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
489
leaderPrintMax(string const message, Real const number) const
490
{
491
492
  Real num(number);
  Real globalMax;
493
494
495
496
497
498
499
500
501
502
503
504
505
506
  M_comm->MaxAll(&num, &globalMax, 1);

  leaderPrint( message , globalMax );

}



template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::setUp( const GetPot& dataFile )
{
    M_steady      = dataFile( "adr/miscellaneous/steady",        0  );
    M_stab        = dataFile( "adr/stab/type",                   "ip");
    M_gammaBeta   = dataFile( "adr/stab/gammaBeta",              0. );
507
    M_diagonalize = dataFile( "adr/space_discretization/diagonalize",  0. );
508
509
510


    M_linearSolver.setDataFromGetPot( dataFile, "adr/solver" );
crosetto's avatar
bug fix    
crosetto committed
511
    M_linearSolver.setUpPrec( dataFile, "adr/prec" );
512
513
514
515

    M_maxIterSolver   = dataFile( "adr/solver/max_iter", -1);
    M_reusePrec       = dataFile( "adr/prec/reuse", true);

516
517
518
519
520
    std::string precType = dataFile( "adr/prec/prectype", "Ifpack");

    M_prec.reset( PRECFactory::instance().createObject( precType ) );
    ASSERT(M_prec.get() != 0, "AdvectionDiffusionSolver : Preconditioner not set");

521
522
523
524
525
526
527
528
529
530
531
532
}

template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::buildSystem()
{


    M_matrMass.reset  ( new matrix_type(M_localMap) );
    M_matrStiff.reset( new matrix_type(M_localMap) );

//    M_comm->Barrier();

533
    leaderPrint("  adr-  Computing constant matrices ...          ");
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559

    Chrono chrono;

    Chrono chronoDer;
    Chrono chronoStiff;
    Chrono chronoMass;
    Chrono chronoGrad;

    Chrono chronoStiffAssemble;
    Chrono chronoMassAssemble;
    Chrono chronoGradAssemble;
    Chrono chronoDivAssemble;
    Chrono chronoStab;
    Chrono chronoZero;

    // Number of velocity components
    UInt nbCompU = 1;

    // Elementary computation and matrix assembling
    // Loop on elements

    UInt velTotalDof   = M_FESpace.dof().numTotalDof();


    chrono.start();

560
    for ( UInt iVol = 1; iVol <= M_FESpace.mesh()->numElements(); iVol++ )
561
562
    {
        chronoDer.start();
563
        M_FESpace.fe().updateFirstDeriv( M_FESpace.mesh()->element( iVol ) );
564
565
566
567
568
569
570
571
572
573
574
575

        chronoDer.stop();

        chronoZero.start();
        M_elmatStiff.zero();
        M_elmatMass.zero();
        chronoZero.stop();


        // stiffness strain
        chronoStiff.start();
        //stiff_strain( 2.0*M_data.viscosity(), M_elmatStiff, M_FESpace.fe() );
576
        stiff( M_data.diffusivity(), M_elmatStiff,  M_FESpace.fe(), 0, 0 );
577
578
579
580
581
582
583
        //stiff_div( 0.5*M_FESpace.fe().diameter(), M_elmatStiff, M_FESpace.fe() );
        chronoStiff.stop();

        // mass
        if ( !M_steady )
        {
            chronoMass.start();
584
            mass( 1., M_elmatMass, M_FESpace.fe(), 0, 0);
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
            chronoMass.stop();
        }

        // stiffness

        chronoStiffAssemble.start();
        assembleMatrix( *M_matrStiff,
                        M_elmatStiff,
                        M_FESpace.fe(),
                        M_FESpace.fe(),
                        M_FESpace.dof(),
                        M_FESpace.dof(),
                        0, 0,
                        0, 0);
        chronoStiffAssemble.stop();

        if ( !M_steady )
        {
            chronoMassAssemble.start();
            assembleMatrix( *M_matrMass,
                            M_elmatMass,
                            M_FESpace.fe(),
                            M_FESpace.fe(),
                            M_FESpace.dof(),
                            M_FESpace.dof(),
                            0, 0, 0, 0);
            chronoMassAssemble.stop();
        }

    }



    M_comm->Barrier();

    chrono.stop();
    leaderPrintMax( "done in " , chrono.diff());


624
    leaderPrint( "  adr-  Finalizing the matrices ...              ");
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
    chrono.start();

    M_matrStiff->GlobalAssemble();
    M_matrMass->GlobalAssemble();

//     M_matrStiff->spy("stiff");
//     M_matrMass->spy("mass");

    chrono.stop();
    leaderPrintMax("done in " , chrono.diff() );;

    if (false)
        std::cout << "partial times:  \n"
                  << " Der            " << chronoDer.diff_cumul() << " s.\n"
                  << " Zero           " << chronoZero.diff_cumul() << " s.\n"
                  << " Stiff          " << chronoStiff.diff_cumul() << " s.\n"
                  << " Stiff Assemble " << chronoStiffAssemble.diff_cumul() << " s.\n"
                  << " Mass           " << chronoMass.diff_cumul() << " s.\n"
                  << " Mass Assemble  " << chronoMassAssemble.diff_cumul() << " s.\n"
                  << " Grad           " << chronoGrad.diff_cumul() << " s.\n"
                  << " Grad Assemble  " << chronoGradAssemble.diff_cumul() << " s.\n"
                  << " Div Assemble   " << chronoDivAssemble.diff_cumul() << " s.\n"
                  << std::endl;

}


template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
initialize( const Function& u0 )
{
     vector_type u(M_FESpace.map());
     M_FESpace.interpolate(u0, u, 0.0);
quinodoz's avatar
   
quinodoz committed
658
     M_sol = u;
659
660
661
662
663
664
665
666
667
668
669
670
671

}


template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
initialize( const vector_type& u0 )
{
    M_sol = u0;
}

template<typename Mesh, typename SolverType>
void ADRSolver<Mesh, SolverType>::
672
updateSystem( Real       alpha,
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
              vector_type& betaVec,
              vector_type& sourceVec
              )
{

    Chrono chrono;

    leaderPrint("  adr-  Updating mass term on right hand side... ");

    chrono.start();

    UInt velTotalDof   = M_FESpace.dof().numTotalDof();

    // Right hand side for the velocity at time

    updateRHS(sourceVec);

    chrono.stop();

    leaderPrintMax("done in ", chrono.diff());


    M_updated = false;

//

    if (M_recomputeMatrix)
        buildSystem();

    leaderPrint( "  adr-  Copying the matrices ...                 ");

    chrono.start();

     if (M_matrNoBC)
         M_matrNoBC.reset(new matrix_type(M_localMap, M_matrNoBC->getMeanNumEntries() ));
     else
         M_matrNoBC.reset(new matrix_type(M_localMap, M_matrStiff->getMeanNumEntries() ));

     M_matrStab.reset( new matrix_type(M_localMap) );

    if ( M_data.diffusivity() != 0. )
        *M_matrNoBC += *M_matrStiff;

    if (alpha != 0. )
    {
        *M_matrNoBC += *M_matrMass*alpha;
    }


    chrono.stop();
    leaderPrintMax( "done in " , chrono.diff() );

//    UInt nbCompU       = nDimensions;

    //! managing the convective term

729
    Real normInf;
730
731
732
733
734
735
736
737
738
739
740
741
742
    betaVec.NormInf(&normInf);

    if (normInf != 0.)
    {

        leaderPrint("  adr-  Updating the convective terms ...        ");

        // vector with repeated nodes over the processors

        vector_type betaVecRep(betaVec, Repeated );
        int uDim = betaVec.size()/nDimensions;
        chrono.start();

743
        for ( UInt iVol = 1; iVol <= M_FESpace.mesh()->numElements(); ++iVol )
744
745
        {

746
747
            M_FESpace.fe().updateFirstDeriv( M_FESpace.mesh()->element( iVol ) ); //as updateFirstDer
            M_betaFESpace.fe().updateFirstDeriv( M_FESpace.mesh()->element( iVol ) ); //as updateFirstDer
748
749
750

            M_elmatAdv.zero();

751
            UInt eleID = M_betaFESpace.fe().currentLocalId();
752
753
            // Non linear term, Semi-implicit approach
            // M_elvec contains the velocity values in the nodes
754
            for ( UInt iNode = 0 ; iNode < ( UInt ) M_betaFESpace.fe().nbNode ; iNode++ )
755
            {
quinodoz's avatar
   
quinodoz committed
756
	         UInt  iloc = M_betaFESpace.fe().patternFirst( iNode );
757
758
                for ( UInt iComp = 0; iComp < nDimensions; ++iComp )
                {
quinodoz's avatar
   
quinodoz committed
759
		    UInt ig = M_betaFESpace.dof().localToGlobal( eleID, iloc + 1 ) + iComp*uDim;
760
                    M_elvec_u.vec()[ iloc + iComp*M_betaFESpace.fe().nbNode ] = betaVecRep[ig]; // BASEINDEX + 1
761
762
763
764
765
                }
            }


            // compute local convective term and assembling
766
767
            for(UInt iComp=0; iComp<nDimensions; iComp++)
            	grad( iComp, M_elvec_u, M_elmatAdv, M_FESpace.fe(), M_FESpace.fe(), M_betaFESpace.fe());
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792




            assembleMatrix( *M_matrNoBC,
                            M_elmatAdv,
                            M_FESpace.fe(),
                            M_FESpace.fe(),
                            M_FESpace.dof(),
                            M_FESpace.dof(),
                            0, 0,
                            0, 0
                            );


            // Streamline Diffusion ( only stab for now )


            if (M_stab == "sd")
            {
                M_elmatStab.zero();
                Real VLoc_infty = 0.;
                Real VLoc_mean  = 0.;
                Real VLoc_c     = 0.;

793
                for ( UInt ih_c = 0 ; ih_c < ( UInt ) this->M_betaFESpace.fe().nbNode ; ih_c++ )
794
                {
795
                    UInt iloc = this->M_betaFESpace.fe().patternFirst( ih_c );
796
797
                    for ( UInt iComp = 0; iComp < nDimensions; ++iComp)
                    {
798
799
                        UInt ig = M_betaFESpace.dof().localToGlobal( eleID, iloc + 1 ) + iComp*uDim;
                        M_elvec_u.vec()[ iloc + iComp * this->M_betaFESpace.fe().nbNode ] = betaVecRep[ ig ];
800
801
802
803
804
805
806
807
808
809
                        VLoc_c += betaVecRep[ ig ] * betaVecRep[ ig ];
                    }

                    VLoc_c     = sqrt( VLoc_c );
                    VLoc_mean += VLoc_c;

                    if ( VLoc_c > VLoc_infty )
                        VLoc_infty = VLoc_c;
                }

810
                VLoc_mean = VLoc_mean / this->M_betaFESpace.fe().nbNode;
811
812

                Real coef_stab, Pe_loc = 0;
813
                coef_stab=M_gammaBeta*this->M_betaFESpace.fe().diameter()*VLoc_infty; // Alessandro - method
814

815
                stiff_sd( coef_stab / ( VLoc_mean*VLoc_mean ), M_elvec_u, M_elmatStab, this->M_FESpace.fe(), this->M_betaFESpace.fe() );
816
817
818
819
820
821
822
823

                assembleMatrix( *M_matrStab,
                                M_elmatStab,
                                M_FESpace.fe(),
                                M_FESpace.fe(),
                                M_FESpace.dof(),
                                M_FESpace.dof(),
                                0, 0, 0, 0 );
824
				}
825
826
        }

827
828
//TODO: check both the 2D and 3D implementation of ip-stabilization (on the Laplacian test doesn't work properly)

829
830
831

        if (M_stab == "ip")
        {
quinodoz's avatar
   
quinodoz committed
832
//	      leaderPrint("   adr- IP stab");
833
//            if ( M_resetStab )
834
835
836
837
838
839
840
841
842
843
844
845
846
            {
                const UInt nDof = M_betaFESpace.dof().numTotalDof();

                CurrentFE fe1(M_FESpace.refFE(),
                              getGeoMap(*M_FESpace.mesh()),
                              M_FESpace.qr());
                CurrentFE fe2(M_FESpace.refFE(),
                              getGeoMap(*M_FESpace.mesh()),
                              M_FESpace.qr());
                CurrentFE fe3(M_betaFESpace.refFE(),
                              getGeoMap(*M_betaFESpace.mesh()),
                              M_betaFESpace.qr());

847
848
849
850
851

#ifdef TWODIM
                typedef ID ( *ETOP )( ID const localFace, ID const point );
                ETOP  eToP;
                switch( M_FESpace.fe().refFE.type )
852
                {
853
854
                    case FE_P1_2D:
                    	eToP = LinearTriangle::eToP;
855
                        break;
856
857
                    case FE_P2_2D:
                    	eToP = QuadraticTriangle::eToP;
858
                        break;
859
860
                    case FE_Q1_2D:
                        eToP = LinearQuad::eToP;
861
                        break;
862
863
                    case FE_Q2_2D:
                        eToP = QuadraticQuad::eToP;
864
865
                        break;
                    default:
866
867
                    	eToP=0;
                        ERROR_MSG( "This refFE is not allowed with IP stabilization" );
868
869
                        break;
                }
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
                for ( UInt iEdge = M_FESpace.mesh()->numBEdges() + 1; iEdge <= M_FESpace.mesh()->numEdges();
                      ++iEdge )
                {
                    const UInt iElAd1 = M_FESpace.mesh()->edge( iEdge ).ad_first();
                    const UInt iElAd2 = M_FESpace.mesh()->edge( iEdge ).ad_second();

                    if ( iElAd1 == iElAd2 || iElAd1 == 0 || iElAd2 == 0)
                    {
                        continue;
                    }

                    M_elmatStab.zero();

                    M_betaFESpace.feBd().updateMeas( M_betaFESpace.mesh()->edge( iEdge ) );
                    const Real hK2  = std::pow(M_betaFESpace.feBd().measure(), 2.);

                    M_betaFESpace.feBd().updateMeasNormal( M_betaFESpace.mesh()->edge( iEdge ) );
                    KNM<Real>& normal = M_betaFESpace.feBd().normal;

                    fe1.updateFirstDeriv( M_FESpace.mesh()->element( iElAd1 ) );
                    fe2.updateFirstDeriv( M_FESpace.mesh()->element( iElAd2 ) );

                   ElemVec beta(M_betaFESpace.feBd().nbNode, nDimensions);

                    // first, get the local trace of the velocity into beta
                    // local id of the face in its adjacent element

                    UInt iEdEl = M_betaFESpace.mesh()->edge( iEdge ).pos_first();
                    for ( int iNode = 0; iNode < M_betaFESpace.feBd().nbNode; ++iNode )
                    {
                        UInt iloc = eToP( iEdEl, iNode+1 );
                        for ( int iCoor = 0; iCoor < fe1.nbCoor; ++iCoor )
                        {
                            UInt ig = M_betaFESpace.dof().localToGlobal( iElAd1, iloc + 1 ) - 1 +iCoor*nDof;
                            if (betaVecRep.BlockMap().LID(ig + 1) >= 0)
                                beta.vec()[ iCoor*M_betaFESpace.feBd().nbNode + iNode ] = betaVecRep( ig + 1); // BASEINDEX + 1
                        }
                    }

#elif defined THREEDIM
                typedef ID ( *FTOP )( ID const localFace, ID const point );
                 FTOP  fToP;
                switch( M_FESpace.fe().refFE.type )
                {
                case FE_P1_3D:
                case FE_P1bubble_3D:
                	fToP = LinearTetra::fToP;
                	break;
                case FE_P2_3D:
                	fToP = QuadraticTetra::fToP;
                	break;
                case FE_Q1_3D:
					fToP = LinearHexa::fToP;
					break;
                case FE_Q2_3D:
                	fToP = QuadraticHexa::fToP;
                	break;
				default:
					fToP = 0;
					ERROR_MSG( "This refFE is not allowed with IP stabilisation" );
					break;
                }

933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953

                for ( UInt iFace = M_FESpace.mesh()->numBFaces() + 1; iFace <= M_FESpace.mesh()->numFaces();
                      ++iFace )
                {

                    const UInt iElAd1 = M_FESpace.mesh()->face( iFace ).ad_first();
                    const UInt iElAd2 = M_FESpace.mesh()->face( iFace ).ad_second();

                    if ( iElAd1 == iElAd2 || iElAd1 == 0 || iElAd2 == 0)
                    {
                        continue;
                    }

                    M_elmatStab.zero();

                    M_betaFESpace.feBd().updateMeas( M_betaFESpace.mesh()->face( iFace ) );
                    const Real hK2  = std::pow(M_betaFESpace.feBd().measure(), 2.);

                    M_betaFESpace.feBd().updateMeasNormal( M_betaFESpace.mesh()->face( iFace ) );
                    KNM<Real>& normal = M_betaFESpace.feBd().normal;

954
955
                    fe1.updateFirstDeriv( M_FESpace.mesh()->element( iElAd1 ) );
                    fe2.updateFirstDeriv( M_FESpace.mesh()->element( iElAd2 ) );
956
957
958
959

                    Real bn   = 0;
                    Real bmax = 0;

quinodoz's avatar
   
quinodoz committed
960
961
		    // Old version, removed by SQ
		    /*
962
963
964
965
966
967
968
969
970
971
972
973
974
                    ElemVec beta( M_betaFESpace.feBd().nbNode, nDimensions);

                    // first, get the local trace of the velocity into beta
                    // local id of the face in its adjacent element

                    UInt iFaEl = M_betaFESpace.mesh()->face( iFace ).pos_first();
                    for ( int iNode = 0; iNode < M_betaFESpace.feBd().nbNode; ++iNode )
                    {
                        UInt iloc = fToP( iFaEl, iNode+1 );
                        for ( int iCoor = 0; iCoor < fe1.nbCoor; ++iCoor )
                        {
                            UInt ig = M_betaFESpace.dof().localToGlobal( iElAd1, iloc + 1 ) - 1 +iCoor*nDof;
                            if (betaVecRep.BlockMap().LID(ig + 1) >= 0)
quinodoz's avatar
   
quinodoz committed
975
976
977
978
			      {
				Real value( betaVecRep( ig + 1) );
                                beta.vec()[ iCoor*M_betaFESpace.feBd().nbNode + iNode ] = value; // BASEINDEX + 1
			      };
979
                        }
quinodoz's avatar
   
quinodoz committed
980
		    }*/
981

quinodoz's avatar
   
quinodoz committed
982
983
984
		    // New version, added by SQ : beta on the domain, not the boundary!
		    // See elemOper.cpp for justification of this usage
                    ElemVec beta( M_betaFESpace.fe().nbNode, nDimensions);
985

quinodoz's avatar
   
quinodoz committed
986
987
988
989
990
991
992
                    for ( int iNode = 0; iNode < M_betaFESpace.fe().nbNode; ++iNode )
                    {
		        UInt  iloc = M_betaFESpace.fe().patternFirst( iNode );
                        for ( int iCoor = 0; iCoor < fe1.nbCoor; ++iCoor )
                        {
                            UInt ig = M_betaFESpace.dof().localToGlobal( iElAd1, iloc + 1 ) + iCoor*nDof;
			    beta.vec()[ iloc + iCoor*M_betaFESpace.fe().nbNode ] = betaVecRep[ig]; // BASEINDEX + 1
993

quinodoz's avatar
   
quinodoz committed
994
995
                        }
		    }
996

quinodoz's avatar
   
quinodoz committed
997
998
999
1000
                    // second, calculate its max norm
                    for ( int l = 0; l < int( M_betaFESpace.fe().nbCoor*M_betaFESpace.fe().nbNode ); ++l ) // SQ: feBd->fe
                    {
		      if ( bmax < fabs( beta.vec()[ l ] ) )
For faster browsing, not all history is shown. View entire blame