LifeV
OseenSolverImplementation.cpp
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 /*!
28  @file
29  @brief File containing the implementation of the file OseenSolver.hpp
30 
31  @author Davide Forti <davide.forti@epfl.ch>
32 
33  @date 28-01-2014
34 
35  */
36 
37 #include <lifev/core/LifeV.hpp>
38 
39 namespace LifeV
40 {
41 
42 class buildVector
43 {
44 public:
45  typedef VectorSmall<3> return_Type;
46 
47  inline return_Type operator() (Real a, Real b, Real c)
48  {
49  VectorSmall<3> o;
50  o[0] = a;
51  o[1] = b;
52  o[2] = c;
53  return o;
54  }
55 
56  buildVector() {}
57  buildVector (const buildVector&) {}
58  ~buildVector() {}
59 };
60 
61 // ===================================================
62 // Constructors & Destructor
63 // ===================================================
64 
65 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
66 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
67 OseenSolver ( boost::shared_ptr<data_Type> dataType,
68  FESpace<mesh_Type, MapEpetra>& velocityFESpace,
69  FESpace<mesh_Type, MapEpetra>& pressureFESpace,
70  boost::shared_ptr<Epetra_Comm>& communicator,
71  const Int lagrangeMultiplier ) :
72  M_oseenData ( dataType ),
73  M_velocityFESpace ( velocityFESpace ),
74  M_pressureFESpace ( pressureFESpace ),
75  M_Displayer ( communicator ),
76  M_fluxMap ( lagrangeMultiplier, communicator),
77  M_localMap ( M_velocityFESpace.map() + M_pressureFESpace.map() + lagrangeMultiplier),
78  M_velocityMatrixMass ( ),
79  M_pressureMatrixMass ( ),
80  M_matrixStokes ( ),
81  M_matrixNoBC ( ),
82  M_matrixStabilization ( ),
83  M_rightHandSideNoBC ( ),
84  M_solution ( new vector_Type ( M_localMap ) ),
85  M_residual ( new vector_Type (M_localMap ) ),
86  M_linearSolver ( new linearSolver_Type (communicator) ),
87  M_steady ( ),
88  M_postProcessing ( new PostProcessingBoundary<mesh_Type> ( M_velocityFESpace.mesh(),
89  &M_velocityFESpace.feBd(),
90  &M_velocityFESpace.dof(),
91  &M_pressureFESpace.feBd(),
92  &M_pressureFESpace.dof(),
93  M_localMap ) ),
94  M_stabilization ( false ),
95  M_reuseStabilization ( false ),
96  M_resetStabilization ( false ),
97  M_iterReuseStabilization ( -1 ),
98  // M_ipStabilization ( M_velocityFESpace.mesh(),
99  // M_velocityFESpace.dof(),
100  // M_velocityFESpace.refFE(),
101  // M_velocityFESpace.feBd(),
102  // M_velocityFESpace.qr(),
103  // 0., 0., 0.,
104  // M_oseenData->viscosity() ),
105  M_betaFunction ( 0 ),
106  M_divBetaUv ( false ),
107  M_stiffStrain ( false ),
108  M_diagonalize ( false ),
109  M_count ( 0 ),
110  M_recomputeMatrix ( false ),
111  M_elementMatrixStiff ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim(), velocityFESpace.fieldDim() ),
112  M_elementMatrixMass ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim(), velocityFESpace.fieldDim() ),
113  M_elementMatrixPreconditioner ( M_pressureFESpace.fe().nbFEDof(), 1, 1 ),
114  M_elementMatrixDivergence ( M_pressureFESpace.fe().nbFEDof(), 1, 0,
115  M_velocityFESpace.fe().nbFEDof(), 0, velocityFESpace.fieldDim() ),
116  M_elementMatrixGradient ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim(), 0,
117  M_pressureFESpace.fe().nbFEDof(), 0, 1 ),
118  M_elementRightHandSide ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
119  M_wLoc ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
120  M_uLoc ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
121  M_un ( new vector_Type (M_localMap) ),
122  M_velocityPreviousTimestep (new vector_Type (M_velocityFESpace.map()) ),
123  M_pressurePreviousTimestep (new vector_Type (M_pressureFESpace.map()) ),
124  M_pressureExtrapolated ( new vector_Type (M_pressureFESpace.map() ) ),
125  M_fespaceUETA ( new ETFESpace_velocity(M_velocityFESpace.mesh(), &(M_velocityFESpace.refFE()), communicator)),
126  M_fespacePETA ( new ETFESpace_pressure(M_pressureFESpace.mesh(), &(M_pressureFESpace.refFE()), communicator)),
127  M_supgStabilization (new StabilizationSUPG<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
128  M_supgVmsStabilization (new StabilizationSUPGVMS<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
129  M_VMSLESStabilization (new StabilizationVMSLES<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace))
130 {
131  *M_solution *= 0;
132  // if(M_stabilization = ( &M_velocityFESpace.refFE() == &M_pressureFESpace.refFE() ))
133  {
134  M_ipStabilization.setFeSpaceVelocity (M_velocityFESpace);
135  M_ipStabilization.setViscosity (M_oseenData->viscosity() );
136  }
137 }
138 
139 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
140 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
141 OseenSolver ( boost::shared_ptr<data_Type> dataType,
142  FESpace<mesh_Type, MapEpetra>& velocityFESpace,
143  FESpace<mesh_Type, MapEpetra>& pressureFESpace,
144  boost::shared_ptr<Epetra_Comm>& communicator,
145  MapEpetra monolithicMap,
146  UInt offset ) :
147  M_oseenData ( dataType ),
148  M_velocityFESpace ( velocityFESpace ),
149  M_pressureFESpace ( pressureFESpace ),
150  M_Displayer ( communicator ),
151  M_fluxMap ( offset, communicator),
152  M_localMap ( monolithicMap ),
153  M_velocityMatrixMass ( ),
154  M_matrixStokes ( ),
155  M_matrixNoBC ( ),
156  M_matrixStabilization ( ),
157  M_rightHandSideNoBC ( ),
158  M_solution ( ),
159  M_residual ( ),
160  M_linearSolver ( ),
161  M_postProcessing ( new PostProcessingBoundary<mesh_Type> (M_velocityFESpace.mesh(),
162  &M_velocityFESpace.feBd(),
163  &M_velocityFESpace.dof(),
164  &M_pressureFESpace.feBd(),
165  &M_pressureFESpace.dof(),
166  M_localMap ) ),
167  M_stabilization ( false ),
168  M_reuseStabilization ( false ),
169  M_resetStabilization ( false ),
170  M_iterReuseStabilization ( -1 ),
171  M_betaFunction ( 0 ),
172  M_divBetaUv ( false ),
173  M_stiffStrain ( false ),
174  M_diagonalize ( false ),
175  M_count ( 0 ),
176  M_recomputeMatrix ( false ),
177  M_elementMatrixStiff ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), M_velocityFESpace.fieldDim() ),
178  M_elementMatrixMass ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), M_velocityFESpace.fieldDim() ),
179  M_elementMatrixPreconditioner ( M_pressureFESpace.fe().nbFEDof(), 1, 1 ),
180  M_elementMatrixDivergence ( M_pressureFESpace.fe().nbFEDof(), 1, 0,
181  M_velocityFESpace.fe().nbFEDof(), 0, M_velocityFESpace.fieldDim() ),
182  M_elementMatrixGradient ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), 0,
183  M_pressureFESpace.fe().nbFEDof(), 0, 1 ),
184  M_elementRightHandSide ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim() ),
185  M_wLoc ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim() ),
186  M_uLoc ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim() ),
187  M_un ( /*new vector_Type(M_localMap)*/ ),
188  M_velocityPreviousTimestep (new vector_Type (M_velocityFESpace.map()) ),
189  M_pressurePreviousTimestep (new vector_Type (M_pressureFESpace.map()) ),
190  M_pressureExtrapolated ( new vector_Type (M_pressureFESpace.map() ) ),
191  M_fespaceUETA ( new ETFESpace_velocity(M_velocityFESpace.mesh(), &(M_velocityFESpace.refFE()), communicator)),
192  M_fespacePETA ( new ETFESpace_pressure(M_pressureFESpace.mesh(), &(M_pressureFESpace.refFE()), communicator)),
193  M_supgStabilization (new StabilizationSUPG<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
194  M_supgVmsStabilization (new StabilizationSUPGVMS<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
195  M_VMSLESStabilization (new StabilizationVMSLES<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace))
196 {
197  ASSERT(0!=0,"ENTERING IN THE MONOLITHIC CONSTRUCTOR");
198  // if(M_stabilization = ( &M_velocityFESpace.refFE() == &M_pressureFESpace.refFE() ))
199  {
200  M_ipStabilization.setFeSpaceVelocity (M_velocityFESpace);
201  M_ipStabilization.setViscosity (M_oseenData->viscosity() );
202  }
203 }
204 
205 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
206 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
207 OseenSolver ( boost::shared_ptr<data_Type> dataType,
208  FESpace<mesh_Type, MapEpetra>& velocityFESpace,
209  FESpace<mesh_Type, MapEpetra>& pressureFESpace,
210  const std::vector<Int>& lagrangeMultipliers,
211  boost::shared_ptr<Epetra_Comm>& communicator ) :
212  M_oseenData ( dataType ),
213  M_velocityFESpace ( velocityFESpace ),
214  M_pressureFESpace ( pressureFESpace ),
215  M_Displayer ( communicator ),
216  M_fluxMap ( lagrangeMultipliers, communicator),
217  M_localMap ( M_velocityFESpace.map() + M_pressureFESpace.map() + lagrangeMultiplier),
218  M_velocityMatrixMass ( ),
219  M_matrixStokes ( ),
220  M_matrixNoBC ( ),
221  M_matrixStabilization ( ),
222  M_rightHandSideNoBC ( ),
223  M_solution ( new vector_Type ( M_localMap ) ),
224  M_residual ( ),
225  M_linearSolver ( new linearSolver_Type (communicator) ),
226  M_postProcessing ( new PostProcessingBoundary<mesh_Type> (M_velocityFESpace.mesh(),
227  &M_velocityFESpace.feBd(),
228  &M_velocityFESpace.dof(),
229  &M_pressureFESpace.feBd(),
230  &M_pressureFESpace.dof(),
231  M_localMap ) ),
232  M_stabilization ( false ),
233  M_reuseStabilization ( false ),
234  M_resetStabilization ( false ),
235  M_iterReuseStabilization ( -1 ),
236  M_betaFunction ( 0 ),
237  M_divBetaUv ( false ),
238  M_stiffStrain ( false ),
239  M_diagonalize ( false ),
240  M_count ( 0 ),
241  M_recomputeMatrix ( false ),
242  M_elementMatrixStiff ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), M_velocityFESpace.fieldDim() ),
243  M_elementMatrixMass ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), M_velocityFESpace.fieldDim() ),
244  M_elementMatrixPreconditioner ( M_pressureFESpace.fe().nbFEDof(), 1, 1 ),
245  M_elementMatrixDivergence ( M_pressureFESpace.fe().nbFEDof(), 1, 0,
246  M_velocityFESpace.fe().nbFEDof(), 0, M_velocityFESpace.fieldDim() ),
247  M_elementMatrixGradient ( M_velocityFESpace.fe().nbFEDof(), M_velocityFESpace.fieldDim(), 0,
248  M_pressureFESpace.fe().nbFEDof(), 0, 1 ),
249  M_elementRightHandSide ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
250  M_wLoc ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
251  M_uLoc ( M_velocityFESpace.fe().nbFEDof(), velocityFESpace.fieldDim() ),
252  M_un ( new vector_Type (M_localMap) ),
253  M_velocityPreviousTimestep (new vector_Type (M_velocityFESpace.map()) ),
254  M_pressurePreviousTimestep (new vector_Type (M_pressureFESpace.map()) ),
255  M_pressureExtrapolated ( new vector_Type (M_pressureFESpace.map() ) ),
256  M_fespaceUETA ( new ETFESpace_velocity(M_velocityFESpace.mesh(), &(M_velocityFESpace.refFE()), communicator)),
257  M_fespacePETA ( new ETFESpace_pressure(M_pressureFESpace.mesh(), &(M_pressureFESpace.refFE()), communicator)),
258  M_supgStabilization (new StabilizationSUPG<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
259  M_supgVmsStabilization (new StabilizationSUPGVMS<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace)),
260  M_VMSLESStabilization (new StabilizationVMSLES<mesh_Type, MapEpetra, SpaceDim>(velocityFESpace, pressureFESpace))
261 {
262  *M_solution *= 0;
263  // if(M_stabilization = ( &M_velocityFESpace.refFE() == &M_pressureFESpace.refFE() ))
264  {
265  M_ipStabilization.setFeSpaceVelocity (M_velocityFESpace);
266  M_ipStabilization.setViscosity (M_oseenData->viscosity() );
267  }
268 }
269 
270 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
271 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
272 ~OseenSolver()
273 {
274 
275 }
276 
277 
278 // ===================================================
279 // Methods
280 // ===================================================
281 
282 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
283 void
284 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::setUp ( const GetPot& dataFile )
285 {
286  if (M_linearSolver.get() )
287  {
288  M_linearSolver->setupPreconditioner ( dataFile, "fluid/prec" );
289  M_linearSolver->setDataFromGetPot ( dataFile, "fluid/solver" );
290  }
291 
292  M_stabilization = dataFile ( "fluid/stabilization/use", /*(&M_velocityFESpace.refFE() == &M_pressureFESpace.refFE() ) ,*/ false);
293 
294  // If using P1-P1 the use of the stabilization is necessary
295  //if(&M_velocityFESpace.refFE() == &M_pressureFESpace.refFE())
296  // M_stabilization = true;
297 
298  M_steady = dataFile ( "fluid/miscellaneous/steady", 0 );
299 
300  Real bdfOrder = dataFile ( "fluid/time_discretization/BDF_order", 1.0 );
301 
302  if (M_stabilization)
303  {
304  if(M_oseenData->stabilizationType() == "IP")
305  {
306  M_gammaBeta = dataFile ( "fluid/ipstab/gammaBeta", 0. );
307  M_gammaDiv = dataFile ( "fluid/ipstab/gammaDiv", 0. );
308  M_gammaPress = dataFile ( "fluid/ipstab/gammaPress", 0. );
309  M_reuseStabilization = dataFile ( "fluid/ipstab/reuse", false );
310 
311  M_ipStabilization.setGammaBeta ( M_gammaBeta );
312  M_ipStabilization.setGammaDiv ( M_gammaDiv );
313  M_ipStabilization.setGammaPress ( M_gammaPress );
314 
315  if (M_linearSolver.get() )
316  M_iterReuseStabilization = dataFile ( "fluid/ipstab/max_iter_reuse", static_cast<Int> ( M_linearSolver->maxNumIterations() * 8. / 10. ) );
317  }
318  else if (M_oseenData->stabilizationType() == "SUPG")
319  {
320  int vel_order = dataFile ( "fluid/space_discretization/vel_order", 1 );
321  M_supgStabilization->setConstant ( vel_order );
322  M_supgStabilization->setETvelocitySpace(M_fespaceUETA);
323  M_supgStabilization->setETpressureSpace(M_fespacePETA);
324  M_supgStabilization->setCommunicator(M_velocityFESpace.map().commPtr());
325  M_supgStabilization->setDensity(M_oseenData->density());
326  M_supgStabilization->setBDForder(bdfOrder);
327  M_supgStabilization->setViscosity(M_oseenData->viscosity());
328  M_supgStabilization->setTimeStep(M_oseenData->dataTime()->timeStep());
329 
330  }
331  else if (M_oseenData->stabilizationType() == "SUPGVMS")
332  {
333  int vel_order = dataFile ( "fluid/space_discretization/vel_order", 1 );
334  M_supgVmsStabilization->setConstant ( vel_order );
335  M_supgVmsStabilization->setETvelocitySpace(M_fespaceUETA);
336  M_supgVmsStabilization->setETpressureSpace(M_fespacePETA);
337  M_supgVmsStabilization->setCommunicator(M_velocityFESpace.map().commPtr());
338  M_supgVmsStabilization->setDensity(M_oseenData->density());
339  M_supgVmsStabilization->setBDForder(bdfOrder);
340  M_supgVmsStabilization->setViscosity(M_oseenData->viscosity());
341  M_supgVmsStabilization->setTimeStep(M_oseenData->dataTime()->timeStep());
342 
343  }
344  else if (M_oseenData->stabilizationType() == "VMSLES")
345  {
346  *M_velocityPreviousTimestep *= 0;
347  *M_pressurePreviousTimestep *= 0;
348  int vel_order = dataFile ( "fluid/space_discretization/vel_order", 1 );
349  M_VMSLESStabilization->setConstant ( vel_order );
350  M_VMSLESStabilization->setETvelocitySpace(M_fespaceUETA);
351  M_VMSLESStabilization->setETpressureSpace(M_fespacePETA);
352  M_VMSLESStabilization->setCommunicator(M_velocityFESpace.map().commPtr());
353  M_VMSLESStabilization->setDensity(M_oseenData->density());
354  M_VMSLESStabilization->setBDForder(bdfOrder);
355  M_VMSLESStabilization->setViscosity(M_oseenData->viscosity());
356  M_VMSLESStabilization->setTimeStep(M_oseenData->dataTime()->timeStep());
357  }
358  }
359  // Energetic stabilization term
360  M_divBetaUv = dataFile ( "fluid/space_discretization/div_beta_u_v", false);
361  // Enable grad( u )^T in stress tensor
362  M_stiffStrain = dataFile ( "fluid/space_discretization/stiff_strain", false);
363  M_diagonalize = dataFile ( "fluid/space_discretization/diagonalize", 1. );
364 
365  // M_linearSolver.setAztecooPreconditioner( dataFile, "fluid/solver" );
366 
367 }
368 
369 
370 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
371 void
372 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
373 initialize ( const function_Type& velocityFunction, const function_Type& pressureFunction )
374 {
375  vector_Type velocityInitialGuess ( M_velocityFESpace.map() );
376  M_velocityFESpace.interpolate ( velocityFunction,
377  velocityInitialGuess,
378  M_oseenData->dataTime()->time() );
379 
380  vector_Type pressureInitialGuess ( M_pressureFESpace.map() );
381  M_pressureFESpace.interpolate ( pressureFunction,
382  pressureInitialGuess,
383  M_oseenData->dataTime()->time() );
384 
385  initialize ( velocityInitialGuess, pressureInitialGuess );
386 }
387 
388 
389 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
390 void
391 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
392 initialize ( const vector_Type& velocityInitialGuess, const vector_Type& pressureInitialGuess )
393 {
394 
395  *M_solution = velocityInitialGuess;
396  *M_un = velocityInitialGuess;
397  M_solution->add ( pressureInitialGuess, M_velocityFESpace.fieldDim() * M_velocityFESpace.dof().numTotalDof() );
398  M_un->add ( pressureInitialGuess, M_velocityFESpace.fieldDim() * M_velocityFESpace.dof().numTotalDof() );
399 
400 }
401 
402 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
403 void
404 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
405 initialize ( const vector_Type& velocityAndPressure )
406 {
407 
408  *M_un = velocityAndPressure;
409  if ( M_solution.get() )
410  {
411  *M_solution = velocityAndPressure;
412  }
413 
414 }
415 
416 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
417 void
418 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::preprocessing( BCHandler& bcloads,
419  const GetPot& dataFile, BCHandler& bcDrag,
420  BCHandler& bcHLift,
421  vector_Type& solDrag,
422  vector_Type& solLift )
423 {
424  matrixPtr_Type laplacian;
425  laplacian.reset( new matrix_Type ( M_velocityFESpace.map() ) );
426 
427  {
428  using namespace ExpressionAssembly;
429 
430  integrate(
431  elements(M_fespaceUETA->mesh()),
432  M_velocityFESpace.qr(),
433  M_fespaceUETA,
434  M_fespaceUETA,
435  dot( grad(phi_i) + transpose(grad(phi_i)) , grad(phi_j) + transpose(grad(phi_j)) )
436  ) >> laplacian;
437  }
438 
439  vector_Type rhs_drag( M_velocityFESpace.map(), Unique );
440  rhs_drag *= 0;
441 
442  vector_Type rhs_lift( M_velocityFESpace.map(), Unique );
443  rhs_lift *= 0;
444 
445  vectorPtr_Type solution_drag;
446  solution_drag.reset( new vector_Type( M_velocityFESpace.map(), Unique ) );
447 
448  vectorPtr_Type solution_lift;
449  solution_lift.reset( new vector_Type( M_velocityFESpace.map(), Unique ) );
450 
451  bcloads.bcUpdate ( *M_velocityFESpace.mesh(), M_velocityFESpace.feBd(), M_velocityFESpace.dof() );
452 
453  bcManage ( *laplacian, rhs_drag, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcloads, M_velocityFESpace.feBd(), 1.0, 0.0 );
454 
455  bcManageRhs ( rhs_drag, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcDrag, M_velocityFESpace.feBd(), 1., 0.);
456  bcManageRhs ( rhs_lift, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHLift, M_velocityFESpace.feBd(), 1., 0.);
457 
458  laplacian->globalAssemble();
459  rhs_drag.globalAssemble();
460  rhs_lift.globalAssemble();
461 
462  // drag
463  linearSolverPtr_Type linearSolver_drag;
464  linearSolver_drag.reset( new linearSolver_Type (M_velocityFESpace.map().commPtr()) );
465 
466  linearSolver_drag->setupPreconditioner ( dataFile, "preprocessing/prec" );
467  linearSolver_drag->setDataFromGetPot ( dataFile, "preprocessing/solver" );
468 
469  linearSolver_drag->setMatrix ( *laplacian );
470 
471  boost::shared_ptr<MatrixEpetra<Real> > staticCast_drag = boost::static_pointer_cast<MatrixEpetra<Real> > (laplacian);
472  Int numIter_drag = linearSolver_drag->solveSystem ( rhs_drag, *solution_drag, staticCast_drag );
473 
474  solDrag = *solution_drag;
475 
476  // lift
477  linearSolverPtr_Type linearSolver_lift;
478  linearSolver_lift.reset( new linearSolver_Type (M_velocityFESpace.map().commPtr()) );
479 
480  linearSolver_lift->setupPreconditioner ( dataFile, "preprocessing/prec" );
481  linearSolver_lift->setDataFromGetPot ( dataFile, "preprocessing/solver" );
482 
483  linearSolver_lift->setMatrix ( *laplacian );
484 
485  boost::shared_ptr<MatrixEpetra<Real> > staticCast_lift = boost::static_pointer_cast<MatrixEpetra<Real> > (laplacian);
486  Int numIter_lift = linearSolver_lift->solveSystem ( rhs_lift, *solution_lift, staticCast_lift );
487 
488  solLift = *solution_lift;
489 }
490 
491 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
492 void
493 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::buildSystem()
494 {
495  // New ETA PART
496  /*
497  * TESTING THE ASSEMBLY USING ETA
498  *
499  */
500 
501  M_Displayer.leaderPrint ( " F- Computing constant matrices ... " );
502  LifeChrono chrono;
503  chrono.start();
504 
505  M_velocityMatrixMass.reset ( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
506  M_matrixStokes.reset(new matrix_block_Type( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ));
507 
508  *M_velocityMatrixMass *= 0;
509  *M_matrixStokes *= 0;
510 
511  {
512  using namespace ExpressionAssembly;
513 
514  if ( !M_steady )
515  {
516  integrate(
517  elements(M_fespaceUETA->mesh()),
518  M_velocityFESpace.qr(),
519  M_fespaceUETA,
520  M_fespaceUETA,
521  M_oseenData->density() * dot(phi_i, phi_j)
522  ) >> M_velocityMatrixMass->block(0,0);
523 
524  M_velocityMatrixMass->globalAssemble();
525 
526  M_matrixMass.reset(new matrix_Type(M_localMap));
527  *M_matrixMass += *M_velocityMatrixMass;
528  M_matrixMass->globalAssemble();
529  }
530 
531  if ( M_stiffStrain )
532  {
533  integrate(
534  elements(M_fespaceUETA->mesh()),
535  M_velocityFESpace.qr(),
536  M_fespaceUETA,
537  M_fespaceUETA,
538  value( 0.5 * M_oseenData->viscosity() ) * dot( grad(phi_i) + transpose(grad(phi_i)) , grad(phi_j) + transpose(grad(phi_j)) )
539  ) >> M_matrixStokes->block(0,0);
540  }
541  else
542  {
543  integrate(
544  elements(M_fespaceUETA->mesh()),
545  M_velocityFESpace.qr(),
546  M_fespaceUETA,
547  M_fespaceUETA,
548  M_oseenData->viscosity() * dot( grad(phi_i) , grad(phi_j) + transpose(grad(phi_j)) )
549  ) >> M_matrixStokes->block(0,0);
550  }
551 
552 
553  integrate(
554  elements(M_fespaceUETA->mesh()),
555  M_velocityFESpace.qr(),
556  M_fespaceUETA,
557  M_fespacePETA,
558  value(-1.0) * phi_j * div(phi_i)
559  ) >> M_matrixStokes->block(0,1);
560 
561  integrate(
562  elements(M_fespaceUETA->mesh()),
563  M_pressureFESpace.qr(),
564  M_fespacePETA,
565  M_fespaceUETA,
566  phi_i * div(phi_j)
567  ) >> M_matrixStokes->block(1,0);
568  }
569 
570  M_matrixStokes->globalAssemble();
571 
572  chrono.stop();
573  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
574 
575 }
576 
577 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
578 void
579 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
580 updateSystem ( const Real alphaOverTimestep,
581  const vector_Type& u_star,
582  const vector_Type& rightHandSide )
583 {
584  if ( M_matrixNoBC.get() )
585  {
586  M_matrixNoBC.reset ( new matrix_Type ( M_localMap, M_matrixNoBC->meanNumEntries() ) );
587  }
588  else
589  {
590  M_matrixNoBC.reset ( new matrix_Type ( M_localMap ) );
591  }
592 
593  updateSystem ( alphaOverTimestep, u_star, rightHandSide, M_matrixNoBC, *M_un );
594  if ( alphaOverTimestep != 0. )
595  {
596  M_matrixNoBC->globalAssemble();
597  }
598 
599 }
600 
601 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
602 void
603 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::
604 updateSystem ( const Real alphaOverTimestep,
605  const vector_Type& u_star,
606  const vector_Type& rightHandSide,
607  matrixPtr_Type matrixNoBC,
608  const vector_Type& un )
609 {
610  M_Displayer.leaderPrint ( " F- Setting the right hand side ..." );
611  LifeChrono chrono;
612  chrono.start();
613 
614  updateRightHandSide ( rightHandSide );
615 
616  chrono.stop();
617  M_Displayer.leaderPrintMax ( " done in ", chrono.diff() );
618 
619  M_matrixNoBC_block.reset(new matrix_block_Type(M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap));
620 
621  if ( M_recomputeMatrix )
622  buildSystem();
623 
624  //! managing the convective term : semi-implicit approximation of the convective term
625 
626  Real normInf;
627  u_star.normInf ( &normInf );
628 
629  MatrixSmall<3, 3> Eye;
630  Eye *= 0.0;
631  Eye[0][0] = 1;
632  Eye[1][1] = 1;
633  Eye[2][2] = 1;
634 
635  // u_star: extrapolation of the velocity
636  // NOTE: for ALE formulation it has to be already: extrapolation of fluid velocity - extrapolation of fluid mesh velocity
637  vector_Type u_starRepeated ( u_star, Repeated );
638 
639  M_convectiveMatrix.reset ( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
640  *M_convectiveMatrix *= 0;
641 
642  // only if the extrapolation of the velocity differs from zero we update the convective term
643  if ( normInf != 0. )
644  {
645  M_convectiveMatrix.reset ( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
646  *M_convectiveMatrix *= 0;
647 
648  if(M_oseenData->conservativeFormulation())
649  {
650  using namespace ExpressionAssembly;
651  integrate(
652  elements(M_fespaceUETA->mesh()), // Mesh
653  M_velocityFESpace.qr(), // QR
654  M_fespaceUETA,
655  M_fespaceUETA,
656  dot(grad(phi_j) * M_oseenData->density() * value(M_fespaceUETA, u_starRepeated), phi_i) // semi-implicit convective term
657  + 0.5 * M_oseenData->density() * dot ( value ( Eye ) , grad(M_fespaceUETA, u_starRepeated) ) * dot( phi_j , phi_i ) // consistency term
658  // - M_oseenData->density() * dot ( value ( Eye ) , grad(M_fespaceUETA, wRepeated) ) * dot( phi_j , phi_i ) // conservative formulation
659  )
660  >> M_convectiveMatrix->block(0,0);
661  }
662  else
663  {
664  using namespace ExpressionAssembly;
665  integrate(
666  elements(M_fespaceUETA->mesh()), // Mesh
667  M_velocityFESpace.qr(), // Quadrature rule
668  M_fespaceUETA,
669  M_fespaceUETA,
670  value(M_oseenData->density()) * dot(phi_i, value(M_fespaceUETA, u_starRepeated)*grad(phi_j)) // semi-implicit treatment of the convective term
671  // value( -1.0*M_oseenData->density() )* dot( grad(phi_i), outerProduct( value(M_fespaceUETA, u_starRepeated), phi_j ) )
672  )
673  >> M_convectiveMatrix->block(0,0);
674  }
675 
676  // Stabilising term: -rho div(u^n) u v
677  if ( M_divBetaUv )
678  {
679  ASSERT(0!=0,"If really needed, please code it using ET");
680  // OLD-STYLE assembly:
681  // mass_divw ( -0.5 * M_oseenData->density(),
682  // M_uLoc,
683  // M_elementMatrixStiff,
684  // M_velocityFESpace.fe(), 0, 0, numVelocityComponent );
685  }
686  }
687 
688  if(M_stabilization)
689  computeStabilization(u_starRepeated, alphaOverTimestep);
690 
691  if ( alphaOverTimestep != 0. )
692  {
693  *M_matrixNoBC_block += (*M_velocityMatrixMass) * alphaOverTimestep;
694  }
695 
696  *M_convectiveMatrix += *M_matrixStokes;
697  M_convectiveMatrix->globalAssemble();
698 
699  *M_matrixNoBC_block += *M_convectiveMatrix;
700  M_matrixNoBC_block->globalAssemble();
701  MapEpetra fullMap ( M_velocityFESpace.map() + M_pressureFESpace.map() + M_fluxMap );
702  M_matrixNoBC.reset(new matrix_Type(fullMap));
703  *M_matrixNoBC += *M_matrixNoBC_block;
704  M_matrixNoBC->globalAssemble();
705 
706 }
707 
708 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
709 void
710 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::computeStabilization ( const vector_Type& u_star, const Real& alpha )
711 {
712  Real normInf;
713  u_star.normInf ( &normInf );
714 
715  LifeChrono chrono;
716 
717  if ( normInf != 0. )
718  {
719  if( M_oseenData->stabilizationType() == "IP" /*&& ( M_resetStabilization || !M_reuseStabilization || ( M_matrixStabilization.get() == 0 ) )*/ )
720  {
721  vector_Type u_starRepeated ( u_star, Repeated );
722  M_Displayer.leaderPrint ( " F- Updating the IP stabilization terms ... " );
723  chrono.start();
724 
725  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
726  M_ipStabilization.apply ( *M_matrixStabilization, u_starRepeated, false );
727  M_matrixStabilization->globalAssemble();
728  M_resetStabilization = false;
729  chrono.stop();
730  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
731  }
732  else if(M_oseenData->stabilizationType() == "SUPG")
733  {
734 
735  M_Displayer.leaderPrint ( " F- Updating the SUPG stabilization terms ... " );
736  chrono.start();
737 
738  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
739  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
740  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
741  *M_matrixStabilizationET *= 0;
742  M_supgStabilization->applySUPG_Matrix_semi_implicit(M_matrixStabilizationET, u_star, alpha);
743  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
744  *M_matrixStabilization += *M_matrixStabilizationET;
745  M_matrixStabilization->globalAssemble();
746 
747  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
748  *M_rhsStabilization *= 0;
749  M_supgStabilization->applySUPG_RHS_semi_implicit(M_rhsStabilization, u_star, *M_velocityRhs);
750 
751  chrono.stop();
752  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
753 
754  }
755  else if(M_oseenData->stabilizationType() == "SUPGVMS")
756  {
757 
758  M_Displayer.leaderPrint ( " F- Updating the SUPG-VMS stabilization terms ... " );
759  chrono.start();
760 
761  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
762  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
763  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
764  *M_matrixStabilizationET *= 0;
765  M_supgVmsStabilization->applySUPGVMS_Matrix_semi_implicit(M_matrixStabilizationET, u_star, alpha);
766  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
767  *M_matrixStabilization += *M_matrixStabilizationET;
768  M_matrixStabilization->globalAssemble();
769 
770  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
771  *M_rhsStabilization *= 0;
772  M_supgVmsStabilization->applySUPGVMS_RHS_semi_implicit(M_rhsStabilization, u_star, *M_velocityRhs);
773 
774  chrono.stop();
775  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
776 
777  }
778  else if(M_oseenData->stabilizationType() == "VMSLES")
779  {
780 
781  M_Displayer.leaderPrint ( " F- Updating the VMSLES stabilization terms ... " );
782  chrono.start();
783 
784  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
785  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
786  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
787  *M_matrixStabilizationET *= 0;
788  M_VMSLESStabilization->applyVMSLES_Matrix_semi_implicit(M_matrixStabilizationET,
789  u_star,
790  alpha,
791  *M_pressureExtrapolated,
792  *M_velocityRhs);
793 
794  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
795  *M_matrixStabilization += *M_matrixStabilizationET;
796  M_matrixStabilization->globalAssemble();
797 
798  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
799  *M_rhsStabilization *= 0;
800  M_VMSLESStabilization->applyVMSLES_RHS_semi_implicit(M_rhsStabilization,
801  u_star,
802  *M_velocityRhs,
803  alpha,
804  *M_pressureExtrapolated);
805 
806  chrono.stop();
807  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
808 
809  }
810  }
811  else
812  {
813  if (M_oseenData->stabilizationType() == "IP")
814  {
815  M_Displayer.leaderPrint ( " F- Updating the IP stabilization terms ... " );
816  chrono.start();
817 
818  if ( M_resetStabilization || !M_reuseStabilization || ( M_matrixStabilization.get() == 0 ) )
819  {
820  vector_Type u_starUnique ( u_star, Unique );
821  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
822  M_ipStabilization.apply ( *M_matrixStabilization, u_starUnique, false );
823  M_resetStabilization = false;
824  chrono.stop();
825  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
826  }
827  else
828  {
829  M_Displayer.leaderPrint ( "reusing\n" );
830  }
831  }
832  else if(M_oseenData->stabilizationType() == "SUPG")
833  {
834 
835  M_Displayer.leaderPrint ( " F- Updating the SUPG stabilization terms ... " );
836  chrono.start();
837 
838  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
839  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
840  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
841  *M_matrixStabilizationET *= 0;
842  M_supgStabilization->applySUPG_Matrix_semi_implicit(M_matrixStabilizationET, u_star, alpha);
843 
844  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
845  *M_matrixStabilization += *M_matrixStabilizationET;
846  M_matrixStabilization->globalAssemble();
847 
848  // comment because if steady simulation this is not needed
849  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
850  *M_rhsStabilization *= 0;
851  M_supgStabilization->applySUPG_RHS_semi_implicit(M_rhsStabilization, u_star, *M_velocityRhs);
852 
853  chrono.stop();
854  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
855 
856  }
857  else if(M_oseenData->stabilizationType() == "SUPGVMS")
858  {
859 
860  M_Displayer.leaderPrint ( " F- Updating the SUPG-VMS stabilization terms ... " );
861  chrono.start();
862 
863  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
864  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
865  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
866  *M_matrixStabilizationET *= 0;
867  M_supgVmsStabilization->applySUPGVMS_Matrix_semi_implicit(M_matrixStabilizationET, u_star, alpha);
868 
869  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
870  *M_matrixStabilization += *M_matrixStabilizationET;
871  M_matrixStabilization->globalAssemble();
872 
873  // comment because if steady simulation this is not needed
874  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
875  *M_rhsStabilization *= 0;
876  M_supgVmsStabilization->applySUPGVMS_RHS_semi_implicit(M_rhsStabilization, u_star, *M_velocityRhs);
877 
878  chrono.stop();
879  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
880 
881  }
882  else if(M_oseenData->stabilizationType() == "VMSLES")
883  {
884 
885  M_Displayer.leaderPrint ( " F- Updating the VMSLES stabilization terms ... " );
886  chrono.start();
887 
888  // TIpically here alpha is already divided by the timestep, but I want to use the actual alfa, so I multiply
889  Real alfa = alpha*M_oseenData->dataTime()->timeStep();
890  M_matrixStabilizationET.reset( new matrix_block_Type ( M_fespaceUETA->map() | M_fespacePETA->map() | M_fluxMap ) );
891  *M_matrixStabilizationET *= 0;
892  M_VMSLESStabilization->applyVMSLES_Matrix_semi_implicit(M_matrixStabilizationET,
893  u_star,
894  alpha,
895  *M_pressureExtrapolated,
896  *M_velocityRhs);
897 
898  M_matrixStabilization.reset ( new matrix_Type ( M_localMap ) );
899  *M_matrixStabilization += *M_matrixStabilizationET;
900  M_matrixStabilization->globalAssemble();
901 
902  // comment because if steady simulation this is not needed
903  M_rhsStabilization.reset(new vector_block_Type( M_velocityFESpace.map() | M_pressureFESpace.map() | M_fluxMap ));
904  *M_rhsStabilization *= 0;
905  M_VMSLESStabilization->applyVMSLES_RHS_semi_implicit( M_rhsStabilization,
906  u_star,
907  *M_velocityRhs,
908  alpha,
909  *M_pressureExtrapolated );
910 
911  chrono.stop();
912  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
913 
914  }
915  }
916 }
917 
918 
919 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
920 void
921 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::updateStabilization ( matrix_Type& matrixFull )
922 {
923  if ( M_stabilization )
924  {
925  M_matrixStabilization->globalAssemble();
926  matrixFull += *M_matrixStabilization;
927  }
928 
929 }
930 
931 template <typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
932 void OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::updateSourceTerm ( source_Type const& source )
933 {
934  vector_Type rhs ( M_velocityFESpace.map() + M_pressureFESpace.map() + M_fluxMap );
935 
936  VectorElemental M_elvec (M_velocityFESpace->fe().nbFEDof(), nDimensions);
937  UInt nc = nDimensions;
938 
939  // loop on volumes: assembling source term
940  for ( UInt i = 1; i <= M_velocityFESpace->mesh()->numVolumes(); ++i )
941  {
942 
943  M_velocityFESpace->fe().updateFirstDerivQuadPt ( M_velocityFESpace->mesh()->volumeList ( i ) );
944 
945  M_elvec.zero();
946 
947  for ( UInt ic = 0; ic < nc; ++ic )
948  {
949  compute_vec ( source, M_elvec, M_velocityFESpace->fe(), M_oseenData->dataTime()->time(), ic ); // compute local vector
950  assembleVector ( *rhs, M_elvec, M_velocityFESpace->fe(), M_velocityFESpace->dof(), ic, ic * M_velocityFESpace->getDim() ); // assemble local vector into global one
951  }
952  }
953  M_rightHandSideNoBC *= 0;
954  M_rightHandSideNoBC += rhs;
955 }
956 
957 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
958 void
959 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::iterate ( bcHandler_Type& bcHandler )
960 {
961 
962  LifeChrono chrono;
963 
964  // matrix and vector assembling communication
965  M_Displayer.leaderPrint ( " F- Updating the boundary conditions ... " );
966 
967  // HERE I SHOULD APPLY THE STABILIZATION ON THE RHS
968 
969  chrono.start();
970 
971  M_matrixNoBC->globalAssemble();
972 
973  matrixPtr_Type matrixFull ( new matrix_Type ( M_velocityFESpace.map() + M_pressureFESpace.map() + M_fluxMap ) );
974 
975  updateStabilization ( *matrixFull );
976 
977  getFluidMatrix ( *matrixFull );
978 
979  vector_Type rightHandSideFull ( *M_rightHandSideNoBC );
980 
981  if(M_stabilization)
982  {
983  if(M_oseenData->stabilizationType() == "SUPG" || M_oseenData->stabilizationType() == "SUPGVMS" || M_oseenData->stabilizationType() == "VMSLES" )
984  {
985  rightHandSideFull += *M_rhsStabilization;
986  }
987  }
988 
989  chrono.stop();
990 
991  M_Displayer.leaderPrintMax ( "done in ", chrono.diff() );
992 
993  // boundary conditions update
994  M_Displayer.leaderPrint (" F- Applying boundary conditions ... ");
995 
996  chrono.start();
997 
998  applyBoundaryConditions ( *matrixFull, rightHandSideFull, bcHandler );
999 
1000  matrixFull->globalAssemble();
1001 
1002  chrono.stop();
1003 
1004  M_Displayer.leaderPrintMax ( "done in " , chrono.diff() );
1005 
1006  // solving the system
1007  M_linearSolver->setMatrix ( *matrixFull );
1008 
1009  boost::shared_ptr<MatrixEpetra<Real> > staticCast = boost::static_pointer_cast<MatrixEpetra<Real> > (matrixFull);
1010 
1011  Int numIter = M_linearSolver->solveSystem ( rightHandSideFull, *M_solution, staticCast );
1012 
1013  // if the preconditioner has been rese the stab terms are to be updated
1014  if ( numIter < 0 || numIter > M_iterReuseStabilization )
1015  {
1016  resetStabilization();
1017  }
1018 
1019  if (M_stabilization && M_oseenData->stabilizationType() == "VMSLES" )
1020  {
1021  *M_velocityPreviousTimestep *= 0;
1022  *M_pressurePreviousTimestep *= 0;
1023 
1024  M_velocityPreviousTimestep->subset(*M_solution, M_velocityFESpace.map(), 0, 0);
1025  M_pressurePreviousTimestep->subset(*M_solution, M_pressureFESpace.map(), 0, 0);
1026  }
1027 
1028  // Computation of the residual for the coefficients
1029 
1030  matrixPtr_Type matrixNoBC ( new matrix_Type ( M_velocityFESpace.map() + M_pressureFESpace.map() + M_fluxMap ) );
1031  *matrixNoBC *= 0;
1032  updateStabilization ( *matrixNoBC );
1033  getFluidMatrix ( *matrixNoBC );
1034  matrixNoBC->globalAssemble();
1035 
1036  vector_Type rightHandSideNoBC ( *M_rightHandSideNoBC );
1037  rightHandSideNoBC.globalAssemble();
1038 
1039  *M_residual *= 0;
1040 
1041  *M_residual = rightHandSideNoBC;
1042  *M_residual -= (*matrixNoBC) * (*M_solution);
1043 
1044  //M_residual.spy("residual");
1045 } // iterate()
1046 
1047 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1048 void
1049 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::h1normVelocity(Real& uh1error )
1050 {
1051  boost::shared_ptr<uExactFunctor> uExactFct( new uExactFunctor );
1052  boost::shared_ptr<vector_Type> uExactVec(new vector_Type(M_velocityFESpace.map(),Unique));
1053  M_velocityFESpace.interpolate(RossEthierSteinmanUnsteadyDec::uexact, *uExactVec, 0.0);
1054 
1055  vector_Type uComputed ( M_fespaceUETA->map() , Unique );
1056  uComputed.subset( *M_solution );
1057 
1058  Real errorH1SquaredLocal( 0.0 );
1059  Real errorH1Squared( 0.0 );
1060 
1061  boost::shared_ptr<gradUExactFunctor> gradUExactFct( new gradUExactFunctor );
1062 
1063  {
1064  using namespace ExpressionAssembly;
1065  integrate (
1066  elements (M_fespaceUETA->mesh() ), // Mesh
1067  M_velocityFESpace.qr(), // QR
1068  dot ( ( eval ( gradUExactFct, value(M_oseenData->dataTime()->time()), X) + (-1) * grad( M_fespaceUETA , uComputed ) ) ,
1069  ( eval ( gradUExactFct, value(M_oseenData->dataTime()->time()), X) + (-1) * grad( M_fespaceUETA , uComputed ) ) )
1070  )
1071  >> errorH1SquaredLocal;
1072 
1073  }
1074 
1075  M_fespaceUETA->mesh()->comm()->Barrier();
1076  M_fespaceUETA->mesh()->comm()->SumAll (&errorH1SquaredLocal, &errorH1Squared, 1);
1077  uh1error= std::sqrt(errorH1Squared);
1078 }
1079 
1080 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1081 void
1082 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::reduceSolution ( Vector& velocityVector, Vector& pressureVector )
1083 {
1084  vector_Type solution ( *M_solution, 0 );
1085 
1086  if ( false /*S_verbose*/ )
1087  {
1088  for ( UInt iDof = 0; iDof < M_velocityFESpace.fieldDim() * dimVelocity(); ++iDof )
1089  {
1090  velocityVector[ iDof ] = solution[ iDof ];
1091  }
1092 
1093  for ( UInt iDof = 0; iDof < dimPressure(); ++iDof )
1094  {
1095  pressureVector[ iDof ] = solution[ iDof + M_velocityFESpace.fieldDim() * dimVelocity() ];
1096  }
1097  }
1098 
1099 }
1100 
1101 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1102 void
1103 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::reduceResidual ( Vector& residualVector )
1104 {
1105  vector_Type residual ( *M_residual, 0 );
1106 
1107  if ( false /*S_verbose*/ )
1108  {
1109  for ( UInt iDof = 0; iDof < M_velocityFESpace.fieldDim() * dimVelocity(); ++iDof )
1110  {
1111  residualVector[ iDof ] = residual[ iDof ];
1112  }
1113 
1114  }
1115 }
1116 
1117 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1118 void
1119 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::getFluidMatrix ( matrix_Type& matrixFull )
1120 {
1121  M_matrixNoBC->globalAssemble();
1122  matrixFull += *M_matrixNoBC;
1123 }
1124 
1125 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1126 void
1127 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::postProcessingSetArea()
1128 {
1129  M_postProcessing->set_area();
1130 }
1131 
1132 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1133 void
1134 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::postProcessingSetNormal()
1135 {
1136  M_postProcessing->set_normal();
1137 }
1138 
1139 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1140 void
1141 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::postProcessingSetPhi()
1142 {
1143  M_postProcessing->set_phi();
1144 }
1145 
1146 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1147 Real
1148 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::flux ( const markerID_Type& flag )
1149 {
1150  return flux ( flag, *M_solution );
1151 }
1152 
1153 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1154 Real
1155 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::flux ( const markerID_Type& flag,
1156  const vector_Type& solution )
1157 {
1158  vector_Type velocityAndPressure ( solution, Repeated );
1159  vector_Type velocity ( this->M_velocityFESpace.map(), Repeated );
1160  velocity.subset ( velocityAndPressure );
1161 
1162  return M_postProcessing->flux ( velocity, flag );
1163 }
1164 
1165 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1166 Real
1167 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::kineticNormalStress ( const markerID_Type& flag )
1168 {
1169  return kineticNormalStress ( flag, *M_solution );
1170 }
1171 
1172 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1173 Real
1174 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::kineticNormalStress ( const markerID_Type& flag,
1175  const vector_Type& solution )
1176 {
1177  vector_Type velocityAndPressure ( solution, Repeated );
1178  vector_Type velocity ( this->M_velocityFESpace.map(), Repeated );
1179  velocity.subset ( velocityAndPressure );
1180 
1181  return M_postProcessing->kineticNormalStress ( velocity, M_oseenData->density(), flag );
1182 }
1183 
1184 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1185 Real
1186 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::area ( const markerID_Type& flag )
1187 {
1188  return M_postProcessing->measure ( flag );
1189 }
1190 
1191 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1192 Vector
1193 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::normal ( const markerID_Type& flag )
1194 {
1195  return M_postProcessing->normal ( flag );
1196 }
1197 
1198 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1199 Vector
1200 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::geometricCenter ( const markerID_Type& flag )
1201 {
1202  return M_postProcessing->geometricCenter ( flag );
1203 }
1204 
1205 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1206 Real
1207 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::pressure ( const markerID_Type& flag )
1208 {
1209  return pressure ( flag, *M_solution );
1210 }
1211 
1212 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1213 Real
1214 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::pressure (const markerID_Type& flag,
1215  const vector_Type& solution)
1216 {
1217  vector_Type velocityAndPressure ( solution, Repeated );
1218  vector_Type pressure ( this->M_pressureFESpace.map(), Repeated );
1219  pressure.subset ( velocityAndPressure,
1220  this->M_velocityFESpace.dim() *this->M_velocityFESpace.fieldDim() );
1221 
1222  // third argument is 1, to use the pressure finite element space (see PostProcessingBoundary docs)
1223  return M_postProcessing->average ( pressure, flag, 1 ) [0];
1224 }
1225 
1226 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1227 Real
1228 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::meanNormalStress ( const markerID_Type& flag, bcHandler_Type& bcHandler )
1229 {
1230  return meanNormalStress ( flag, bcHandler, *M_solution );
1231 }
1232 
1233 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1234 Real
1235 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::meanNormalStress (const markerID_Type& flag, bcHandler_Type& bcHandler, const vector_Type& solution )
1236 {
1237  if ( bcHandler.findBCWithFlag ( flag ).type() == Flux )
1238  {
1239  return -lagrangeMultiplier ( flag, bcHandler, solution );
1240  }
1241  else
1242  {
1243 #ifdef HAVE_LIFEV_DEBUG
1244  M_Displayer.leaderPrint ( " !!! WARNING - OseenSolver::meanNormalStress( flag, bcHandler, solution) is returning an approximation \n" );
1245 #endif
1246  return -pressure ( flag, solution ); // TODO: This is an approximation of the mean normal stress as the pressure.
1247  // A proper method should be coded in the PostprocessingBoundary class
1248  // to compute the exact mean normal stress.
1249  }
1250 }
1251 
1252 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1253 Real
1254 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::meanTotalNormalStress ( const markerID_Type& flag, bcHandler_Type& bcHandler )
1255 {
1256  return meanTotalNormalStress ( flag, bcHandler, *M_solution );
1257 }
1258 
1259 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1260 Real
1261 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::meanTotalNormalStress (const markerID_Type& flag, bcHandler_Type& bcHandler, const vector_Type& solution )
1262 {
1263  return meanNormalStress ( flag, bcHandler, solution ) - kineticNormalStress ( flag, solution );
1264 }
1265 
1266 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1267 Real
1268 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::lagrangeMultiplier ( const markerID_Type& flag, bcHandler_Type& bcHandler )
1269 {
1270  return lagrangeMultiplier ( flag, bcHandler, *M_solution );
1271 }
1272 
1273 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1274 Real
1275 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::lagrangeMultiplier ( const markerID_Type& flag,
1276  bcHandler_Type& bcHandler,
1277  const vector_Type& solution )
1278 {
1279  // Create a list of Flux bcName_Type ??
1280  std::vector< bcName_Type > fluxBCVector = bcHandler.findAllBCWithType ( Flux );
1281  bcName_Type fluxbcName_Type = bcHandler.findBCWithFlag ( flag ).name();
1282 
1283  // Create a Repeated vector for looking to the lambda
1284  vector_Type velocityPressureLambda ( solution, Repeated );
1285 
1286  // Find the index associated to the correct Lagrange multiplier
1287  for ( UInt lmIndex = 0; lmIndex < static_cast <UInt> ( fluxBCVector.size() ); ++lmIndex )
1288  if ( fluxbcName_Type.compare ( fluxBCVector[ lmIndex ] ) == 0 )
1289  return velocityPressureLambda[M_velocityFESpace.fieldDim() * M_velocityFESpace.dof().numTotalDof()
1290  + M_pressureFESpace.dof().numTotalDof() + lmIndex];
1291 
1292  // If lmIndex has not been found a warning message is printed
1293  M_Displayer.leaderPrint ( "!!! Warning - Lagrange multiplier for Flux BC not found!\n" );
1294  return 0;
1295 }
1296 
1297 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1298 Real
1299 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::energy ( )
1300 {
1301  Real energy = 0.0;
1302 
1303  vector_Type tmp(M_localMap);
1304  tmp *= 0;
1305  tmp = (*M_matrixMass) * (*M_solution);
1306  energy = M_solution->dot(tmp);
1307  energy *= 0.5;
1308 
1309  return energy;
1310 }
1311 
1312 
1313 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1314 VectorSmall<2>
1315 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::computeForces ( bcHandler_Type& bcHandlerDrag,
1316  bcHandler_Type& bcHandlerLift )
1317 {
1318  vector_Type onesOnBodyDrag(M_localMap, Unique);
1319  onesOnBodyDrag *= 0;
1320 
1321  vector_Type onesOnBodyLift(M_localMap, Unique);
1322  onesOnBodyLift *= 0;
1323 
1324  bcManageRhs ( onesOnBodyDrag, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerDrag, M_velocityFESpace.feBd(), 1., 0.);
1325  bcManageRhs ( onesOnBodyLift, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerLift, M_velocityFESpace.feBd(), 1., 0.);
1326 
1327  Real drag (0.0);
1328  Real lift (0.0);
1329 
1330  drag = M_residual->dot(onesOnBodyDrag);
1331  lift = M_residual->dot(onesOnBodyLift);
1332 
1333  M_Displayer.leaderPrint ( " F- Value of the drag: ", drag );
1334  M_Displayer.leaderPrint ( " F- Value of the lift: ", lift );
1335 
1336  VectorSmall<2> Forces;
1337  Forces[0] = drag;
1338  Forces[1] = lift;
1339 
1340  return Forces;
1341 }
1342 
1343 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1344 VectorSmall<2>
1345 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::computeForcesNewTest ( const vector_Type& drag_vector,
1346  const vector_Type& lift_vector )
1347 {
1348  Real drag (0.0);
1349  Real lift (0.0);
1350 
1351  drag = M_residual->dot(drag_vector);
1352  lift = M_residual->dot(lift_vector);
1353 
1354  M_Displayer.leaderPrint ( " F- Value of the drag: ", drag );
1355  M_Displayer.leaderPrint ( " F- Value of the lift: ", lift );
1356 
1357  VectorSmall<2> Forces;
1358  Forces[0] = drag;
1359  Forces[1] = lift;
1360 
1361  return Forces;
1362 }
1363 
1364 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1365 VectorSmall<2>
1366 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::computeForces ( const markerID_Type& flag,
1367  bcHandler_Type& bcHandlerDrag,
1368  bcHandler_Type& bcHandlerLift)
1369 {
1370  // 1) flag che dice in quale componente calcolare la resistenza (0 in x, 1 in y, 2 in z)
1371  // 2) creare vettore pieno di zeri lungo quanto il vettore soluzione -> vettore soluzione u organizzato come ( ux | uy | uz | p)
1372  // 3) mettere degli 1 nei nodi corrispondenti al bordo nel blocco del vettore, ovvero:
1373  // se la resistenza é in direzione x: vettore da creare ( 1suBody | 0 | 0 | 0 )
1374  // 4) prodotto scalare del vettore del residuo con questo vettore definito al punto 3)
1375 
1376 
1377  vector_Type onesOnBodyDrag(M_localMap, Unique);
1378  onesOnBodyDrag *= 0;
1379 
1380  vector_Type onesOnBodyLift(M_localMap, Unique);
1381  onesOnBodyLift *= 0;
1382 
1383  bcManageRhs ( onesOnBodyDrag, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerDrag, M_velocityFESpace.feBd(), 1., 0.);
1384  bcManageRhs ( onesOnBodyLift, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerLift, M_velocityFESpace.feBd(), 1., 0.);
1385 
1386  Real drag (0.0);
1387  Real lift (0.0);
1388 
1389  drag = M_residual->dot(onesOnBodyDrag);
1390  lift = M_residual->dot(onesOnBodyLift);
1391 
1392  M_Displayer.leaderPrint ( " F- Value of the drag: ", drag );
1393  M_Displayer.leaderPrint ( " F- Value of the lift: ", lift );
1394 
1395  VectorSmall<2> Forces;
1396  Forces[0] = drag;
1397  Forces[1] = lift;
1398 
1399  return Forces;
1400 }
1401 
1402 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1403 VectorSmall<2>
1404 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::computeDrag ( const markerID_Type& flag,
1405  bcHandler_Type& bcHandlerDrag,
1406  bcHandler_Type& bcHandlerLift,
1407  const Real& velocityInfty,
1408  const Real& Area)
1409 {
1410  // 1) flag che dice in quale componente calcolare la resistenza (0 in x, 1 in y, 2 in z)
1411  // 2) creare vettore pieno di zeri lungo quanto il vettore soluzione -> vettore soluzione u organizzato come ( ux | uy | uz | p)
1412  // 3) mettere degli 1 nei nodi corrispondenti al bordo nel blocco del vettore, ovvero:
1413  // se la resistenza é in direzione x: vettore da creare ( 1suBody | 0 | 0 | 0 )
1414  // 4) prodotto scalare del vettore del residuo con questo vettore definito al punto 3)
1415 
1416 
1417  vector_Type onesOnBodyDrag(M_localMap, Unique);
1418  onesOnBodyDrag *= 0;
1419 
1420  vector_Type onesOnBodyLift(M_localMap, Unique);
1421  onesOnBodyLift *= 0;
1422 
1423  bcManageRhs ( onesOnBodyDrag, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerDrag, M_velocityFESpace.feBd(), 1., 0.);
1424  bcManageRhs ( onesOnBodyLift, *M_velocityFESpace.mesh(), M_velocityFESpace.dof(), bcHandlerLift, M_velocityFESpace.feBd(), 1., 0.);
1425 
1426  Real drag (0.0);
1427  Real lift (0.0);
1428 
1429  drag = M_residual->dot(onesOnBodyDrag);
1430  lift = M_residual->dot(onesOnBodyLift);
1431 
1432  M_Displayer.leaderPrint ( " F- Value of the drag: ", drag );
1433  M_Displayer.leaderPrint ( " F- Value of the lift: ", lift );
1434 
1435  M_Displayer.leaderPrint ( "\n\n" );
1436 
1437  drag /= (0.5*M_oseenData->density()*velocityInfty*velocityInfty*Area);
1438  lift /= (0.5*M_oseenData->density()*velocityInfty*velocityInfty*Area);
1439 
1440  M_Displayer.leaderPrint ( " F- Value of the drag coefficient: ", drag );
1441  M_Displayer.leaderPrint ( " F- Value of the lift coefficient: ", lift );
1442 
1443  VectorSmall<2> Coefficients;
1444  Coefficients[0] = drag;
1445  Coefficients[1] = lift;
1446 
1447  return Coefficients;
1448 }
1449 
1450 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1451 Real
1452 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::removeMean ( vector_Type& x )
1453 {
1454 
1455  LifeChrono chrono;
1456  chrono.start();
1457 
1458  const UInt numVelocityComponent ( M_velocityFESpace.fieldDim() );
1459  const UInt velocityTotalDof ( M_velocityFESpace.dof().numTotalDof() );
1460 
1461 
1462  if ( M_pressureMatrixMass.get() == 0 )
1463  {
1464  M_pressureMatrixMass.reset ( new matrix_Type ( M_localMap ) );
1465  }
1466 
1467  for ( UInt iElement = 0; iElement < M_velocityFESpace.mesh()->numElements(); iElement++ )
1468  {
1469  chrono.start();
1470  // just to provide the id number in the assem_mat_mixed
1471  M_pressureFESpace.fe().update ( M_pressureFESpace.mesh()->element ( iElement ) );
1472 
1473  M_elementMatrixPreconditioner.zero();
1474  // mass
1475  chrono.start();
1476  mass ( 1, M_elementMatrixPreconditioner, M_pressureFESpace.fe(), 0, 0, M_velocityFESpace.fieldDim() );
1477  chrono.stop();
1478 
1479  chrono.start();
1480  assembleMatrix ( *M_pressureMatrixMass,
1481  M_elementMatrixPreconditioner,
1482  M_pressureFESpace.fe(),
1483  M_pressureFESpace.fe(),
1484  M_pressureFESpace.dof(),
1485  M_pressureFESpace.dof(),
1486  numVelocityComponent,
1487  numVelocityComponent,
1488  numVelocityComponent * velocityTotalDof,
1489  numVelocityComponent * velocityTotalDof );
1490  chrono.stop();
1491  }
1492 
1493  M_pressureMatrixMass->globalAssemble();
1494 
1495  vector_Type ones ( *M_solution );
1496  ones = 1.0;
1497 
1498  Real mean;
1499  mean = ones * ( M_pressureMatrixMass * x );
1500  x += ( -mean );
1501 
1502  ASSERT ( std::fabs ( ones * ( M_pressureMatrixMass * x ) ) < 1e-9 , "after removeMean the mean pressure should be zero!");
1503 
1504  return mean;
1505 
1506 } // removeMean()
1507 
1508 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1509 void
1510 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::applyBoundaryConditions ( matrix_Type& matrix,
1511  vector_Type& rightHandSide,
1512  bcHandler_Type& bcHandler )
1513 {
1514  // M_rightHandSideFull = M_rightHandSideNoBC;
1515 
1516  // BC manage for the velocity
1517  if ( !bcHandler.bcUpdateDone() || M_recomputeMatrix )
1518  {
1519  bcHandler.bcUpdate ( *M_velocityFESpace.mesh(),
1520  M_velocityFESpace.feBd(),
1521  M_velocityFESpace.dof() );
1522  }
1523 
1524  // ignoring non-local entries, Otherwise they are summed up lately
1525  //vector_Type rightHandSideFull( rightHandSide, Repeated, Zero );
1526  // ignoring non-local entries, Otherwise they are summed up lately
1527  vector_Type rightHandSideFull ( rightHandSide, Unique );
1528 
1529  bcManage ( matrix, rightHandSideFull,
1530  *M_velocityFESpace.mesh(),
1531  M_velocityFESpace.dof(),
1532  bcHandler,
1533  M_velocityFESpace.feBd(),
1534  1.,
1535  M_oseenData->dataTime()->time() );
1536 
1537  rightHandSide = rightHandSideFull;
1538 
1539  if ( bcHandler.hasOnlyEssential() && M_diagonalize )
1540  {
1541  matrix.diagonalize ( M_velocityFESpace.fieldDim() *dimVelocity(),
1542  M_diagonalize,
1543  rightHandSide,
1544  0. );
1545  }
1546 
1547 } // applyBoundaryCondition
1548 
1549 // ===================================================
1550 // Set Methods
1551 // ===================================================
1552 
1553 
1554 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1555 void
1556 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::setupPostProc( )
1557 {
1558  M_postProcessing.reset ( new PostProcessingBoundary<mesh_Type> ( M_velocityFESpace.mesh(),
1559  &M_velocityFESpace.feBd(),
1560  &M_velocityFESpace.dof(),
1561  &M_pressureFESpace.feBd(),
1562  &M_pressureFESpace.dof(),
1563  M_localMap ) );
1564 }
1565 
1566 template<typename MeshType, typename SolverType, typename MapType , UInt SpaceDim, UInt FieldDim>
1567 void
1568 OseenSolver<MeshType, SolverType, MapType , SpaceDim, FieldDim>::setTolMaxIteration ( const Real& tolerance, const Int& maxIteration )
1569 {
1570  M_linearSolver->setTolerance ( tolerance );
1571  M_linearSolver->setMaxNumIterations ( maxIteration );
1572 }
1573 
1574 } //end namespace LifeV
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::buildVector::buildVector
buildVector(const buildVector &)
Definition: OseenSolverImplementation.cpp:57
LifeV::buildVector::operator()
return_Type operator()(Real a, Real b, Real c)
Definition: OseenSolverImplementation.cpp:47
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191