LifeV
StabilizationSD.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 /*!
27  @file
28  @brief Streamline diffusion and SUPG stabilization.
29  @author M.A. Fernandez
30  @contributor Umberto Villa <uvilla@emory.edu>
31  @maintainer Umberto Villa <uvilla@emory.edu>
32  @date 01/05/2005
33 
34  This file contains a c++ class implementing the streamline diffusion (SD) and the SUPG
35  stabilization for the Navier-Stokes equations. Tested with P1/P1 and Q1/Q1
36 
37 */
38 
39 #ifndef _SDSTABILIZATION_HPP_
40 #define _SDSTABILIZATION_HPP_
41 
42 
43 #include <lifev/core/filter/GetPot.hpp>
44 #include <lifev/core/util/LifeDebug.hpp>
45 
46 namespace LifeV
47 {
48 
49 //! StabilizationSD Class
50 /*!
51  * @brief Streamline diffusion and SUPG for Navier-Stokes
52  * @author M.A. Fernandez
53  *
54  * Implementation of streamline diffusion (SD) and SUPG for the Navier-Stokes equations. <br>
55  * SUPG can be used only with equal order finite element for the discretization of velocity and pressure fields.
56  */
57 
58 template<typename MeshType, typename DofType>
60 {
61 public:
62 
63  //@name Public Types
64  //@{
65  typedef MeshType mesh_Type;
66  typedef DofType dof_Type;
67  //@}
68 
69  //! @name Constructor and Destructor
70  //@{
71  //! Constructor
72  /*!
73  * @param dataFile GetPot dataFile where to read the stabilization parameter
74  * @param mesh mesh_Type mesh
75  * @param dof dof_Type velocity field degree of freedom
76  * @param refFE refFE velocity field reference finite element
77  * @param quadRule QuadratureRule quadrature rule for the integration of the stabilization variational forms
78  * @param viscosity viscosity fluid viscosity @f$\nu@f$
79  */
80  StabilizationSD ( const GetPot& dataFile,
81  const mesh_Type& mesh,
82  const dof_Type& dof,
83  const ReferenceFE& refFE,
84  const QuadratureRule& quadRule,
85  Real viscosity);
86  //! ~Destructor
87  virtual ~StabilizationSD() {};
88  //@}
89 
90  //! @name Methods
91  //@{
92  //! compute the SUPG stabilization terms and add them into the monolithic N.S. matrix
93  /*!
94  * This function adds the following stabilization terms to the Navier-Stokes monolithic matrix:
95  * <ol>
96  * <li> Block(1,1): @f$c_\beta(\beta \nabla \mathbf{u} , \beta \nabla \mathbf{v}) + c_\beta( - \mu L \mathbf{u} , \beta \nabla \mathbf{v} )+ c_d( div \mathbf{u} , div \nabla \mathbf{v})@f$
97  * <li> Block(1,2): @f$c_\beta(\nabla p , \beta \nabla \mathbf{v})@f$
98  * <li> Block(2,1): @f$c_\beta(\beta \nabla \mathbf{u} , \nabla q) + c_\beta( - \mu L \mathbf{u} , \nabla q )@f$
99  * <li> Block(2,2): @f$c_\beta(\nabla p , \nabla q)@f$
100  * </ol>
101  * where @f$(\cdot, \cdot)@f$ represents the @f$L^2@f$ scalar product, and
102  * <ol>
103  * <li> @f$\displaystyle c_\beta = \frac{\gamma_\beta}{\sqrt{ 4/( dt^2)
104  + 4\|\beta\|_\infty/h^2
105  + 16*\nu/h^4 }} @f$, where @f$h@f$ is the diameter of the element;
106  * <li> @f$\displaystyle c_d = \gamma_d h \|\beta\|_\infty @f$.
107  * </ol>
108  * Both high Pechlet numbers and inf-sup incompatible FEM are stabilized.
109  *
110  * PREREQUISITE: The velocity and the pressure field should belong to the same finite element space
111  *
112  * Parameters are the followings:
113  * @param dt Real timestep (INPUT)
114  * @param matrix MatrixType where the stabilization terms are added into. (OUTPUT)
115  * @param state VectorType velocity field for the linearization of the stabilization (INPUT)
116  */
117  template<typename MatrixType, typename VectorType>
118  void applySUPG (const Real dt, MatrixType& matrix, const VectorType& state );
119 
120  //! compute the SD stabilization terms and add them into the Momentum matrix
121  /*!
122  * The following stabilization term is added:
123  * @f$c_\beta(\beta \nabla \mathbf{u} , \beta \nabla \mathbf{v})@f$
124  * @param dt Real timestep
125  * @param matrix MatrixType where the stabilization terms are added into
126  * @param state VectorType velocity field for the linearization of the stabilization
127  */
128  template<typename MatrixType, typename VectorType>
129  void applySD (const Real dt, MatrixType& matrix, const VectorType& state );
130 
131  //! compute the SUPG stabilization terms and add them into the right and side
132  /*!
133  * Add the following stabilization terms to rhs
134  * <ol>
135  * <li> Block(1,1): @f$(c_\beta f , \beta \nabla \mathbf{v})@f$;
136  * <li> Block(1,2): @f$(c_\beta f , \nabla q )@f$.
137  * </ol>
138  *
139  * @param dt Real timestep
140  * @param vector VectorType where the stabilization terms are added into
141  * @param state VectorType velocity field for the linearization of the stabilization
142  * @param source SourceType a functor f(time, x, y, z, ic) which represents the forcing term
143  * @param time REAL the actual time in which source should be evaluated
144  */
145  template <typename VectorType, typename SourceType >
146  void applyRHS (const Real dt, VectorType& vector, const VectorType& state,
147  const SourceType& source, const Real& time);
148 
149  //! Display class informations
150  void showMe (std::ostream& output = std::cout) const;
151  //@}
152 
153  //! @name Set Methods
154  //@{
155  //! Set the stabilization parameter @f$\gamma_\beta@f$ for @f$( \beta \nabla \mathbf{u} , \beta \nabla \mathbf{v})@f$
156  void setGammaBeta (const Real& gammaBeta)
157  {
158  M_gammaBeta = gammaBeta;
159  }
160  //! Set the stabilization parameter @f$\gamma_d@f$ for @f$( div \mathbf{u} , div \nabla \mathbf{v})@f$
161  void setGammaDiv (const Real& gammaDiv)
162  {
163  M_gammaDiv = gammaDiv;
164  }
165  //@}
166 private:
167 
168  //! @name Private Constructor
169  //@{
170  //! Default Constructor
171  StabilizationSD();
172  //! Copy Constructor
174  //@}
175 
176  //! @name Private Methods
177  //@{
178  // methods for elementary computations
179 
180  //! Compute the stabilization coefficients for each element
181  template <typename VectorType>
182  void computeParameters (const Real dt, const UInt iVol, const VectorType& state,
183  VectorElemental& beta, Real& coeffBeta, Real& coeffDiv) const;
184 
185  //!Evaluate the varf @f$(\beta \nabla \mathbf{u}, \beta \nabla \mathbf{v})@f$
186  void bgradu_bgradv (const Real& coef, VectorElemental& vel, ElemMat& elmat, const CurrentFE& fe,
187  UInt iblock, UInt jblock, UInt nb) const;
188 
189  //! Evaluate the varf @f$(\Delta \mathbf{u}, \beta \nabla \mathbf{v})@f$
190  void lapu_bgradv (const Real& coef, VectorElemental& vel, MatrixElemental& elmat, const CurrentFE& fe,
191  UInt iblock, UInt jblock, UInt nb) const;
192 
193  //! Evaluate the varf @f$(\nabla p, \beta \nabla \mathbf{v})@f$
194  void gradp_bgradv (const Real& coef, VectorElemental& vel, MatrixElemental& elmat, const CurrentFE& fe) const;
195 
196  //! Evaluate the varf @f$(\Delta \mathbf{u}, \nabla q)@f$
197  void lapu_gradq (const Real& coef, MatrixElemental& elmat, const CurrentFE& fe) const;
198 
199  //! Evaluate the varf @f$(f, \beta \nabla \mathbf{v})@f$
200  template <typename SourceType>
201  void f_bgradv (const Real& coef, SourceType& source, VectorElemental& vel,
202  VectorElemental& elvec, const CurrentFE& fe, UInt iblock, const Real& time) const;
203 
204  //! Evaluate the varf @f$(f, \nabla q)@f$
205  template<typename SourceType>
206  void f_gradq (const Real& coef, SourceType& source, VectorElemental& elvec,
207  const CurrentFE& fe, UInt iblock, const Real& time) const;
208  //@}
209 
210  //! @name Private Attributes
211  //@{
212  //! the mesh object
214  //! the dof object
215  const dof_Type& M_dof;
216  //! current fe for the assembling of stabilization terms.
218  //! fluid viscosity @f$\nu@f$
220  //! Stabilization coefficient of @f$(c(h,dt, |\beta|, nu) \beta \nabla \mathbf{u} , \beta \nabla \mathbf{v})@f$
222  //! Stabilization coefficient of @f$(c(h,dt) div \mathbf{u} , div \nabla \mathbf{v})@f$
224  //! Elementary Matrix for assembling the stabilization terms
226  //! Elementary Vector for assembling the stabilization terms
228  //@}
229 }; // class StabilizationSD
230 
231 //=============================================================================
232 // Constructor
233 //=============================================================================
234 
235 template<typename MeshType, typename DofType>
236 StabilizationSD<MeshType, DofType>::StabilizationSD ( const GetPot& dataFile,
237  const mesh_Type& mesh,
238  const dof_Type& dof,
239  const ReferenceFE& refFE,
240  const QuadratureRule& quadRule,
241  Real viscosity) :
242  M_mesh ( mesh ),
243  M_dof ( dof ),
245  M_viscosity ( viscosity ),
246  M_gammaBeta ( dataFile ( "fluid/sdstab/gammaBeta", 0. ) ),
247  M_gammaDiv ( dataFile ( "fluid/sdstab/gammaDiv", 0. ) ),
248  M_elMat ( M_fe.nbNode, nDimensions + 1, nDimensions + 1 ) ,
249  M_elVec ( M_fe.nbNode, nDimensions + 1 ) {}
250 
251 //=============================================================================
252 // Methods
253 //=============================================================================
254 
255 template<typename MeshType, typename DofType>
256 template <typename MatrixType, typename VectorType>
257 void StabilizationSD<MeshType, DofType>::applySUPG (const Real dt, MatrixType& matrix, const VectorType& state )
258 {
259  if ( M_gammaBeta == 0 && M_gammaDiv == 0)
260  {
261  return;
262  }
263 
264  Chrono chronoBeta;
265  Chrono chronoUpdate;
266  Chrono chronoElemComp;
267  Chrono chronoAssembly;
268 
269  // stabilization parameters
270  Real coeffBeta, coeffDiv;
271 
272  // local velocity
273  VectorElemental beta ( M_fe.nbNode, nDimensions );
274 
275  // loop on elements
276  for ( UInt iVol = 0; iVol < M_mesh.numVolumes(); iVol++ )
277  {
278  chronoUpdate.start();
279  // update current finite elements
280  M_fe.updateFirstSecondDeriv ( M_mesh.volumeList ( iVol ) );
281 
282  // stabilization parameters computation
283  chronoBeta.start();
284  this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
285  chronoBeta.stop();
286 
287  chronoElemComp.start();
288  M_elMat.zero();
289 
290  // coeffBeta (\beta \nabla \mathbf{u} , \beta \nabla \mathbf{v})
291  //
292  this->bgradu_bgradv (coeffBeta, beta, M_elMat, M_fe, 0, 0, nDimensions);
293 
294  // coeffBeta ( (\beta \nabla \mathbf{u} , \nabla q) + (\nabla p , \beta \nabla \mathbf{v}) )
295  //
296  this->gradp_bgradv (coeffBeta, beta, M_elMat, M_fe);
297 
298  // coeffBeta (\nabla p , \nabla q)
299  //
300  stiff ( coeffBeta, M_elMat, M_fe, nDimensions, nDimensions );
301 
302  // coeffBeta ( - \mu L \mathbf{u} , \beta \nabla \mathbf{v} )
303  //
304  this->lapu_bgradv (-coeffBeta * M_viscosity, beta, M_elMat, M_fe, 0, 0, nDimensions);
305 
306  // coeffBeta ( - \mu L \mathbf{u} , \nabla q )
307  //
308  this->lapu_gradq (-coeffBeta * M_viscosity, M_elMat, M_fe);
309 
310  // coeffDiv ( div \mathbf{u} , div \nabla \mathbf{v})
311  //
312  stiff_div (coeffDiv, M_elMat, M_fe );
313 
314  chronoElemComp.stop();
315 
316  chronoAssembly.start();
317  for ( UInt iComp = 0; iComp <= nDimensions; ++iComp )
318  for ( UInt jComp = 0; jComp <= nDimensions; ++jComp )
319  //assemb_mat( matrix, M_elMat, M_fe, M_dof, iComp, jComp );
320  assembleMatrix ( matrix, M_elMat, M_fe, M_dof,
321  iComp, jComp,
322  iComp * M_dof.numTotalDof(), jComp * M_dof.numTotalDof() );
323  chronoAssembly.stop();
324 
325 
326  }// loop on elements
327 
328  debugStream (7100) << std::endl;
329  debugStream (7100) << " Updating of element done in "
330  << chronoUpdate.diffCumul() << "s." << std::endl;
331  debugStream (7100) << " Determination of parameters done in "
332  << chronoBeta.diffCumul() << "s." << std::endl;
333  debugStream (7100) << " Element computations done in "
334  << chronoElemComp.diffCumul() << "s." << std::endl;
335  debugStream (7100) << " Assembly done in "
336  << chronoAssembly.diffCumul() << "s." << std::endl;
337 
338 
339 } // applySUPG(...)
340 
341 
342 template<typename MeshType, typename DofType>
343 template <typename MatrixType, typename VectorType>
344 void StabilizationSD<MeshType, DofType>::applySD (const Real dt, MatrixType& matrix, const VectorType& state )
345 {
346  if ( M_gammaBeta == 0 && M_gammaDiv == 0)
347  {
348  return;
349  }
350 
351  Chrono chronoBeta;
352  Chrono chronoUpdate;
353  Chrono chronoElemComp;
354  Chrono chronoAssembly;
355 
356  // Stabilization parameters
357  Real coeffBeta/*, coeffDiv*/;
358 
359  // local velocity
360  VectorElemental beta ( M_fe.nbNode, nDimensions );
361 
362  // loop on elements
363  for ( UInt iVol = 0; iVol < M_mesh.numVolumes(); iVol++ )
364  {
365  chronoUpdate.start();
366  // update current finite elements
367  M_fe.updateFirstSecondDeriv ( M_mesh.volumeList ( iVol ) );
368 
369  // stabilization parameters computation
370  chronoBeta.start();
371  this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
372  chronoBeta.stop();
373 
374  chronoElemComp.start();
375  M_elMat.zero();
376 
377  // coeffBeta (\beta \nabla \mathbf{u} , \beta \nabla \mathbf{v})
378  //
379  this->bgradu_bgradv (coeffBeta, beta, M_elMat, M_fe, 0, 0, nDimensions);
380 
381  // coeffBeta ( - \mu L \mathbf{u} , \beta \nabla \mathbf{v} )
382  //
383  //this->lapu_bgradv(-coeffBeta*M_viscosity, beta, M_elMat, M_fe, 0, 0, nDimensions);
384 
385  // coeffDiv ( div \mathbf{u} , div \nabla \mathbf{v})
386  //
387  //stiff_div(coeffDiv, M_elMat, M_fe );
388 
389  chronoElemComp.stop();
390 
391  chronoAssembly.start();
392  for ( UInt iComp = 0; iComp < nDimensions; ++iComp )
393  for ( UInt jComp = 0; jComp < nDimensions; ++jComp )
394  {
395  assemb_mat ( matrix, M_elMat, M_fe, M_dof, iComp, jComp );
396  }
397  chronoAssembly.stop();
398 
399 
400  }// loop on elements
401 
402  debugStream (7100) << std::endl;
403  debugStream (7100) << " Updating of element done in "
404  << chronoUpdate.diffCumul() << "s." << std::endl;
405  debugStream (7100) << " Determination of parameters done in "
406  << chronoBeta.diffCumul() << "s." << std::endl;
407  debugStream (7100) << " Element computations done in "
408  << chronoElemComp.diffCumul() << "s." << std::endl;
409  debugStream (7100) << " Assembly done in "
410  << chronoAssembly.diffCumul() << "s." << std::endl;
411 
412 
413 } // applySD(...)
414 
415 template<typename MeshType, typename DofType>
416 template <typename VectorType, typename SourceType >
417 void StabilizationSD<MeshType, DofType>::applyRHS (const Real dt, VectorType& vector, const VectorType& state,
418  const SourceType& source, const Real& time)
419 {
420  if ( M_gammaBeta == 0 )
421  {
422  return;
423  }
424 
425 
426  Chrono chronoBeta;
427  Chrono chronoUpdate;
428  Chrono chronoElemComp;
429  Chrono chronoAssembly;
430 
431  // stabilization parameters
432  Real coeffBeta, coeffDiv;
433 
434  // local velocity
435  VectorElemental beta ( M_fe.nbNode, nDimensions );
436 
437  // loop on elements
438  for ( UInt iVol = 0; iVol < M_mesh.numVolumes(); iVol++ )
439  {
440  chronoUpdate.start();
441  // update current finite elements
442  M_fe.updateFirstDeriv ( M_mesh.volumeList ( iVol ) );
443  // local mesh parameters
444  chronoUpdate.stop();
445 
446  chronoBeta.start();
447  this->computeParameters (dt, iVol, state, beta, coeffBeta, coeffDiv);
448  chronoBeta.stop();
449 
450  chronoElemComp.start();
451  M_elVec.zero();
452 
453  // coeffBeta ( f , \beta \nabla \mathbf{v})
454  //
455  this->f_bgradv (coeffBeta, source, beta, M_elVec, M_fe, 0, time);
456 
457  // coeffBeta ( f , \nabla q )
458  //
459  this->f_gradq (coeffBeta, source, M_elVec, M_fe, nDimensions, time);
460  chronoElemComp.stop();
461 
462  chronoAssembly.start();
463  for ( UInt iComp = 0; iComp <= nDimensions; ++iComp )
464  //assemb_vec( vector, M_elVec, M_fe, M_dof, iComp);
465  {
466  assembleVector ( vector, M_elVec, M_fe, M_dof, iComp, iComp * M_dof.numTotalDof() );
467  }
468  chronoAssembly.stop();
469 
470 
471  }// loop on elements
472 
473  debugStream (7100) << std::endl;
474  debugStream (7100) << " Updating of element done in "
475  << chronoUpdate.diffCumul() << "s." << std::endl;
476  debugStream (7100) << " Determination of parameters done in "
477  << chronoBeta.diffCumul() << "s." << std::endl;
478  debugStream (7100) << " Element computations done in "
479  << chronoElemComp.diffCumul() << "s." << std::endl;
480  debugStream (7100) << " Assembly done in "
481  << chronoAssembly.diffCumul() << "s." << std::endl;
482 
483 
484 } // applyRHS(...)
485 
486 template<typename MeshType, typename DofType>
487 void StabilizationSD<MeshType, DofType>::showMe (std::ostream& output) const
488 {
489  output << "StabilizationSD::showMe() " << std::endl;
490  output << "Fluid Viscosity: " << M_viscosity << std::endl;
491  output << "Stabilization coefficient SUPG/SD: " << M_gammaBeta << std::endl;
492  output << "Stabilization coefficient div grad: " << M_gammaDiv << std::endl;
493  M_mesh.showMe (output);
494  M_dof.showMe (output);
495 }
496 
497 //=============================================================================
498 // Private Method
499 //=============================================================================
500 template<typename MeshType, typename DofType>
501 template<typename VectorType>
502 void StabilizationSD<MeshType, DofType>::computeParameters (const Real dt, const UInt iVol, const VectorType& state,
503  VectorElemental& beta, Real& coeffBeta, Real& coeffDiv) const
504 {
505 
506  const UInt nDof = M_dof.numTotalDof();
507 
508  // square local mesh parameter in 1-norm
509  Real hK, hK2, hK4;
510 
511  hK = M_fe.diameter();
512  hK2 = hK * hK;
513  hK4 = hK2 * hK2;
514 
515  // determine bmax = ||\beta||_{0,\infty,K}
516  // first, get the local velocity into beta
517  for ( UInt iNode = 0; iNode < M_fe.nbNode; ++iNode )
518  {
519  for ( UInt iCoor = 0; iCoor < M_fe.nbLocalCoor(); ++iCoor )
520  {
521  UInt ig = M_dof.localToGlobalMap ( iVol, iNode ) + iCoor * nDof;
522  beta.vec() [ iCoor * M_fe.nbNode + iNode ] = state[ig];
523  }
524  }
525 
526  // second, calculate its max norm
527  Real bmax = std::fabs ( beta.vec() [ 0 ] );
528  for ( UInt l = 1; l < UInt ( M_fe.nbLocalCoor() *M_fe.nbNode ); ++l )
529  {
530  if ( bmax < std::fabs ( beta.vec() [ l ] ) )
531  {
532  bmax = std::fabs ( beta.vec() [ l ] );
533  }
534  }
535 
536  coeffBeta = M_gammaBeta / std::sqrt ( 4. / ( dt * dt)
537  + 4.*bmax * bmax / hK2
538  + 16.*M_viscosity * M_viscosity / hK4 );
539  coeffDiv = M_gammaDiv * bmax * hK;
540 
541 }
542 
543 
544 template<typename MeshType, typename DofType>
545 void StabilizationSD<MeshType, DofType>::gradp_bgradv (const Real& coef, VectorElemental& vel,
546  MatrixElemental& elmat, const CurrentFE& fe) const
547 {
548  ASSERT_PRE (fe.hasFirstDeriv(),
549  "advection_grad matrix needs at least the first derivatives");
550 
551  MatrixElemental::matrix_type v (fe.nbLocalCoor(), fe.nbQuadPt() );
552  Real s;
553 
554 
555  // local velocity at quadrature points
556  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
557  {
558  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
559  {
560  VectorElemental::vector_view velicoor = vel.block (icoor);
561  v (icoor, ig) = 0.;
562  for (UInt k (0); k < fe.nbNode; ++k)
563  {
564  v (icoor, ig) += velicoor (k) * fe.phi (k, ig); // velocity on the intgt point
565  }
566  }
567  }
568 
569  for (UInt ic (0); ic < fe.nbLocalCoor(); ++ic)
570  {
571  MatrixElemental::matrix_view mat_ic3 = elmat.block (ic, fe.nbLocalCoor() );
572  MatrixElemental::matrix_view mat_3ic = elmat.block (fe.nbLocalCoor(), ic);
573  for (UInt i = 0; i < fe.nbNode; ++i)
574  {
575  for (UInt j = 0; j < fe.nbNode; ++j)
576  {
577  s = 0.0;
578  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
579  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
580  {
581  s += fe.phiDer (j, ic, ig) * v (jcoor, ig) * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
582  }
583  mat_ic3 (i, j) += coef * s;
584  mat_3ic (j, i) += coef * s;
585  }
586  }
587  }
588 }
589 
590 
591 template<typename MeshType, typename DofType>
592 void StabilizationSD<MeshType, DofType>::bgradu_bgradv (const Real& coef, VectorElemental& vel, MatrixElemental& elmat, const CurrentFE& fe,
593  UInt iblock, UInt jblock, UInt nb) const
594 {
595  ASSERT_PRE (fe.hasFirstDeriv(),
596  "advection (vect) matrix needs at least the first derivatives");
597 
598 
599  MatrixElemental::matrix_type mat_tmp (fe.nbNode, fe.nbNode);
600  MatrixElemental::matrix_type v ( fe.nbLocalCoor(), fe.nbQuadPt() );
601  Real s;
602 
603 
604  // compute local vectors values
605  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
606  {
607  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
608  {
609  VectorElemental::vector_view velicoor = vel.block (icoor);
610  v (icoor, ig) = 0.;
611  for (UInt k (0); k < fe.nbNode; k++)
612  {
613  v (icoor, ig) += velicoor (k) * fe.phi (k, ig); // velocity on the intgt point
614  }
615  }
616  }
617 
618  for (UInt i (0); i < fe.nbNode; ++i)
619  {
620  for (UInt j (0); j < fe.nbNode; ++j)
621  {
622  s = 0.0;
623 
624  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
625  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
626  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
627  {
628  s += fe.phiDer (i, jcoor, ig) * v (jcoor, ig) * v (icoor, ig) * fe.phiDer (j, icoor, ig) * fe.weightDet (ig);
629  }
630  mat_tmp (i, j) = coef * s;
631  }
632  }
633 
634  // copy on the components
635  for (UInt icomp (0); icomp < nb; icomp++)
636  {
637  MatrixElemental::matrix_view mat_icomp = elmat.block (iblock + icomp, jblock + icomp);
638  for (UInt i (0); i < fe.nbDiag(); ++i)
639  {
640  for (UInt j (0); j < fe.nbDiag(); ++j)
641  {
642  mat_icomp (i, j) += mat_tmp (i, j);
643  }
644  }
645  }
646 }
647 
648 
649 
650 template<typename MeshType, typename DofType>
651 void StabilizationSD<MeshType, DofType>::lapu_bgradv (const Real& coef, VectorElemental& vel, MatrixElemental& elmat, const CurrentFE& fe,
652  UInt iblock, UInt jblock, UInt nb) const
653 {
654 
655 
656  ASSERT_PRE (fe.hasFirstDeriv(),
657  "lapu_bgradv matrix needs first derivatives");
658 
659  ASSERT_PRE (fe.hasSecondDeriv(),
660  "lapu_bgradv matrix needs second derivatives");
661 
662  MatrixElemental::matrix_type mat_tmp (fe.nbNode, fe.nbNode);
663  MatrixElemental::matrix_type v ( fe.nbLocalCoor(), fe.nbQuadPt() );
664  Real s;
665 
666 
667  // compute local vectors values at quadrature points
668  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
669  {
670  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
671  {
672  VectorElemental::vector_view velicoor = vel.block (icoor);
673  v (icoor, ig) = 0.;
674  for (UInt k (0); k < fe.nbNode; ++k)
675  {
676  v (icoor, ig) += velicoor (k) * fe.phi (k, ig); // velocity on the intgt point
677  }
678  }
679  }
680 
681 
682  // numerical integration
683  for (UInt i (0); i < fe.nbNode; ++i)
684  {
685  for (UInt j (0); j < fe.nbNode; ++j)
686  {
687  s = 0.0;
688  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
689  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
690  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
691  {
692  s += fe.phiDer2 (j, icoor, icoor, ig) * v (jcoor, ig) * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
693  }
694  mat_tmp (i, j) = coef * s;
695  }
696  }
697 
698  // copy on the components
699  for (UInt icomp (0); icomp < nb; ++icomp)
700  {
701  MatrixElemental::matrix_view mat_icomp = elmat.block (iblock + icomp, jblock + icomp);
702  for (UInt i (0); i < fe.nbDiag(); ++i)
703  {
704  for (UInt j (0); j < fe.nbDiag(); ++j)
705  {
706  mat_icomp (i, j) += mat_tmp (i, j);
707  }
708  }
709  }
710 }
711 
712 
713 template<typename MeshType, typename DofType>
714 void StabilizationSD<MeshType, DofType>::lapu_gradq (const Real& coef, MatrixElemental& elmat, const CurrentFE& fe) const
715 {
716 
717  ASSERT_PRE (fe.hasFirstDeriv(),
718  "lapu_gradq matrix needs first derivatives");
719 
720  ASSERT_PRE (fe.hasSecondDeriv(),
721  "lapu_gradq matrix needs second derivatives");
722 
723  Real s;
724 
725  for (UInt jc (0); jc < fe.nbLocalCoor(); ++jc) // loop on column blocks
726  {
727  MatrixElemental::matrix_view mat_view = elmat.block (fe.nbLocalCoor(), jc);
728  for (UInt i (0); i < fe.nbNode; ++i) // local rows
729  {
730  for (UInt j (0); j < fe.nbNode; ++j) // local columns
731  {
732  s = 0.0;
733  // quadrature formula
734  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
735  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor) // lap
736  {
737  s += fe.phiDer2 (j, jcoor, jcoor, ig) * fe.phiDer (i, jc, ig) * fe.weightDet (ig);
738  }
739  mat_view (i, j) += coef * s;
740  }
741  }
742  }
743 }
744 
745 
746 template<typename MeshType, typename DofType>
747 template<typename SourceType>
748 void StabilizationSD<MeshType, DofType>::f_bgradv (const Real& coef, SourceType& source, VectorElemental& vel,
749  VectorElemental& elvec, const CurrentFE& fe, UInt iblock, const Real& time) const
750 {
751 
752  ASSERT_PRE (fe.hasFirstDeriv(),
753  "f_bgradv vector needs at least the first derivatives");
754 
755  MatrixElemental::matrix_type v (fe.nbLocalCoor(), fe.nbQuadPt() );
756  Real s;
757 
758 
759  // local velocity at quadrature points
760  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
761  {
762  for (UInt icoor (0); icoor < fe.nbLocalCoor(); ++icoor)
763  {
764  VectorElemental::vector_view velicoor = vel.block (icoor);
765  v (icoor, ig) = 0.;
766  for (UInt k (0); k < fe.nbNode; ++k)
767  {
768  v (icoor, ig) += velicoor (k) * fe.phi (k, ig); // velocity on the intgt point
769  }
770  }
771  }
772 
773  // local vector per block
774  for (UInt ic (0); ic < fe.nbLocalCoor(); ++ic)
775  {
776  VectorElemental::vector_view vec_ic = elvec.block (ic + iblock);
777  for (UInt i (0); i < fe.nbNode; ++i)
778  {
779  s = 0.0;
780  for (UInt ig (0); ig < fe.nbQuadPt(); ++ig)
781  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
782  s += source (time, fe.quadPt (ig, 0), fe.quadPt (ig, 1), fe.quadPt (ig, 2), ic + 1)
783  * fe.phiDer (i, jcoor, ig) * v (jcoor, ig) * fe.weightDet (ig);
784  vec_ic (i) += coef * s;
785  }
786  }
787 }
788 
789 
790 
791 template<typename MeshType, typename DofType>
792 template<typename SourceType>
793 void StabilizationSD<MeshType, DofType>::f_gradq (const Real& coef, SourceType& source, VectorElemental& elvec, const CurrentFE& fe, UInt iblock, const Real& time) const
794 {
795 
796  ASSERT_PRE (fe.hasFirstDeriv(),
797  "f_gradq vector needs at least the first derivatives");
798 
799  Real s;
800 
801  VectorElemental::vector_view vec_ic = elvec.block (iblock);
802  for (UInt i (0); i < fe.nbNode; ++i)
803  {
804  s = 0.0;
805  for (UInt ig = 0; ig < fe.nbQuadPt(); ++ig)
806  for (UInt jcoor (0); jcoor < fe.nbLocalCoor(); ++jcoor)
807  s += source (time, fe.quadPt (ig, 0), fe.quadPt (ig, 1), fe.quadPt (ig, 2), jcoor + 1)
808  * fe.phiDer (i, jcoor, ig) * fe.weightDet (ig);
809  vec_ic (i) += coef * s;
810  }
811 }
812 
813 
814 
815 
816 
817 
818 
819 } // namespace LifeV
820 
821 #endif /* _SDSTABILIZATION_HPP_ */
StabilizationSD Class.
CurrentFE M_fe
current fe for the assembling of stabilization terms.
double operator()(const char *VarName, const double &Default) const
Definition: GetPot.hpp:2033
Real M_viscosity
fluid viscosity
Real M_gammaBeta
Stabilization coefficient of .
void bgradu_bgradv(const Real &coef, VectorElemental &vel, ElemMat &elmat, const CurrentFE &fe, UInt iblock, UInt jblock, UInt nb) const
Evaluate the varf .
void applySUPG(const Real dt, MatrixType &matrix, const VectorType &state)
compute the SUPG stabilization terms and add them into the monolithic N.S. matrix ...
StabilizationSD()
Default Constructor.
#define debugStream
Definition: LifeDebug.hpp:182
void f_gradq(const Real &coef, SourceType &source, VectorElemental &elvec, const CurrentFE &fe, UInt iblock, const Real &time) const
Evaluate the varf .
void lapu_gradq(const Real &coef, MatrixElemental &elmat, const CurrentFE &fe) const
Evaluate the varf .
void applySD(const Real dt, MatrixType &matrix, const VectorType &state)
compute the SD stabilization terms and add them into the Momentum matrix
void applyRHS(const Real dt, VectorType &vector, const VectorType &state, const SourceType &source, const Real &time)
compute the SUPG stabilization terms and add them into the right and side
void updateInverseJacobian(const UInt &iQuadPt)
StabilizationSD(const GetPot &dataFile, const mesh_Type &mesh, const dof_Type &dof, const ReferenceFE &refFE, const QuadratureRule &quadRule, Real viscosity)
Constructor.
void gradp_bgradv(const Real &coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe) const
Evaluate the varf .
Real M_gammaDiv
Stabilization coefficient of .
VectorElemental M_elVec
Elementary Vector for assembling the stabilization terms.
void setGammaBeta(const Real &gammaBeta)
Set the stabilization parameter for .
void lapu_bgradv(const Real &coef, VectorElemental &vel, MatrixElemental &elmat, const CurrentFE &fe, UInt iblock, UInt jblock, UInt nb) const
Evaluate the varf .
#define ASSERT_PRE(X, A)
Definition: LifeAssert.hpp:96
void computeParameters(const Real dt, const UInt iVol, const VectorType &state, VectorElemental &beta, Real &coeffBeta, Real &coeffDiv) const
Compute the stabilization coefficients for each element.
void setGammaDiv(const Real &gammaDiv)
Set the stabilization parameter for .
const dof_Type & M_dof
the dof object
StabilizationSD(const StabilizationSD< mesh_Type, dof_Type > &original)
Copy Constructor.
double Real
Generic real data.
Definition: LifeV.hpp:175
const UInt nDimensions(NDIM)
virtual ~StabilizationSD()
~Destructor
const mesh_Type & M_mesh
the mesh object
void showMe(std::ostream &output=std::cout) const
Display class informations.
void f_bgradv(const Real &coef, SourceType &source, VectorElemental &vel, VectorElemental &elvec, const CurrentFE &fe, UInt iblock, const Real &time) const
Evaluate the varf .
MatrixElemental M_elMat
Elementary Matrix for assembling the stabilization terms.
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191