LifeV
StabilizationSUPG_semi_implicit_ale.cpp
Go to the documentation of this file.
1 #include <lifev/navier_stokes_blocks/solver/StabilizationSUPG_semi_implicit_ale.hpp>
2 
3 // MACRO TO DEFINE TAU_M
4 #define TAU_M value(1.0)/( eval(squareroot,TAU_M_DEN) )
5 #define TAU_M_DEN TAU_M_DEN_DT + TAU_M_DEN_VEL + TAU_M_DEN_VISC
6 #define TAU_M_DEN_DT value(M_density*M_density)*value(M_bdfOrder*M_bdfOrder)/value(M_timestep * M_timestep)
7 #define TAU_M_DEN_VEL value(M_density*M_density)*dot(value(M_fespaceUETA, beta_rep), value(M_fespaceUETA, beta_rep))/(h_K*h_K)
8 #define TAU_M_DEN_VISC value(M_C_I)*value(M_viscosity*M_viscosity)/(h_K*h_K*h_K*h_K)
9 
10 #define TAU_M_TILDE value(1.0)/( TAU_M_DEN_DT_NOSQUARED + eval(squareroot,TAU_M_DEN_VEL+TAU_M_DEN_VISC) )
11 #define TAU_M_DEN_DT_NOSQUARED value(M_density)*value(M_alpha)/value(M_timestep)
12 
13 // MACRO TO DEFINE TAU_C
14 #define TAU_C (h_K*h_K)/(TAU_M)
15 
16 // MACRO BELOW FOR FINE SCALE
17 #define TAU_M_NO_DT value(1.0)/( eval( squareroot, TAU_M_DEN_VEL + TAU_M_DEN_VISC ) )
18 #define TAU_C_NO_DT (h_K*h_K)/(TAU_M_NO_DT)
19 
20 // MACRO FOR EVALUATION OF THE FINE SCALE VELOCITY
21 #define TAU_M_UPRIME value(1.0)/( TAU_M_DEN_DT_NOSQUARED + eval(squareroot,TAU_M_DEN_UPRIME) )
22 #define TAU_M_DEN_UPRIME TAU_M_DEN_VEL_UPRIME + TAU_M_DEN_VISC_UPRIME
23 #define TAU_M_DEN_VEL_UPRIME value(M_density*M_density)*dot(value(M_fespaceUETA, *velocity_repeated), value(M_fespaceUETA, *velocity_repeated))/(h_K*h_K)
24 #define TAU_M_DEN_VISC_UPRIME value(M_C_I)*value(M_viscosity*M_viscosity)/(h_K*h_K*h_K*h_K)
25 
26 // MACRO FOR EVALUATION OF THE FINE SCALE PRESSURE
27 #define TAU_M_NO_DT_UPRIME value(1.0)/( eval( squareroot, TAU_M_DEN_VEL_UPRIME + TAU_M_DEN_VISC_UPRIME ) )
28 #define TAU_C_NO_DT_PPRIME (h_K*h_K)/(TAU_M_NO_DT_UPRIME)
29 
30 namespace LifeV
31 {
32 
33 //=============================================================================
34 // Constructor
35 //=============================================================================
36 
37 StabilizationSUPG_semi_implicit_ale::StabilizationSUPG_semi_implicit_ale():
38  M_label("SUPG_semi_implicit_ale"),
39  M_useODEfineScale ( false )
40 {
41 }
42 
43 //=============================================================================
44 // Methods
45 //=============================================================================
46 
47 void StabilizationSUPG_semi_implicit_ale::setConstant(const int & value)
48 {
49  if ( value == 1 )
50  M_C_I = 30;
51  else if ( value == 2 )
52  M_C_I = 60;
53  else
54  ASSERT(0!=0, "Please implement a suitable value for M_C_I for your velocity FE order");
55 }
56 
57 void StabilizationSUPG_semi_implicit_ale::setUseODEfineScale ( const bool& useODEfineScale )
58 {
59  M_useODEfineScale = useODEfineScale;
60  setupODEfineScale();
61 }
62 
63 void StabilizationSUPG_semi_implicit_ale::setupODEfineScale ( )
64 {
65  int numVolumes = M_uFESpace->mesh()->numVolumes();
66  UInt numQuadraturePointsVelocity = M_uFESpace->qr().nbQuadPt();
67 
68  std::vector<std::vector<VectorSmall<3>>> velocity_fine;
69  velocity_fine.resize(numVolumes);
70 
71  M_fineScaleVelocity.resize(numVolumes);
72  M_fineScalePressure.resize(numVolumes);
73  M_fineScaleVelocityRhs.resize(numVolumes);
74 
75  for ( int i = 0; i < numVolumes; ++i )
76  {
77  velocity_fine[i].resize( numQuadraturePointsVelocity );
78  M_fineScaleVelocityRhs[i].resize( numQuadraturePointsVelocity );
79  M_fineScaleVelocity[i].resize( numQuadraturePointsVelocity );
80  M_fineScalePressure[i].resize( numQuadraturePointsVelocity );
81 
82  // Here we assume that same quadrature rule
83  // is used for pressure and velocity.
84  for ( int j = 0; j < numQuadraturePointsVelocity; ++j )
85  {
86  velocity_fine[i][j](0) = 0.0;
87  velocity_fine[i][j](1) = 0.0;
88  velocity_fine[i][j](2) = 0.0;
89 
90  M_fineScaleVelocity[i][j](0) = 0.0;
91  M_fineScaleVelocity[i][j](1) = 0.0;
92  M_fineScaleVelocity[i][j](2) = 0.0;
93 
94  M_fineScalePressure[i][j](0) = 0.0;
95 
96  M_fineScaleVelocityRhs[i][j](0) = 0.0;
97  M_fineScaleVelocityRhs[i][j](1) = 0.0;
98  M_fineScaleVelocityRhs[i][j](2) = 0.0;
99  }
100  }
101 
102  M_handlerFineScaleVelocity.reset( new TimeAndExtrapolationHandlerQuadPts<3> ( ) );
103  M_handlerFineScaleVelocity->setBDForder(M_bdfOrder);
104  M_handlerFineScaleVelocity->setTimeStep(M_timestep);
105  std::vector<std::vector<std::vector<VectorSmall<3>>>> initial_state_velocity;
106  for ( int i = 0; i < M_bdfOrder; ++i )
107  {
108  initial_state_velocity.push_back ( velocity_fine );
109  }
110  M_handlerFineScaleVelocity->initialize(initial_state_velocity);
111 
112 }
113 
114 void StabilizationSUPG_semi_implicit_ale::updateODEfineScale ( const vectorPtr_Type& velocity, const vectorPtr_Type& pressure )
115 {
116  computeFineScales(velocity, pressure );
117  computeFineScalesForVisualization ( velocity, pressure );
118  M_handlerFineScaleVelocity->shift(M_fineScaleVelocity);
119  M_handlerFineScaleVelocity->rhsContribution(M_fineScaleVelocityRhs);
120 }
121 
122 void StabilizationSUPG_semi_implicit_ale::computeFineScales ( const vectorPtr_Type& velocity, const vectorPtr_Type& pressure )
123 {
124  vectorPtr_Type velocity_repeated( new vector_Type( *velocity, Repeated ) );
125  vectorPtr_Type pressure_repeated( new vector_Type( *pressure, Repeated ) );
126  vectorPtr_Type velocity_rhs_repeated( new vector_Type( *M_rhsVelocity, Repeated) );
127 
128  using namespace ExpressionAssembly;
129 
130  std::shared_ptr<SquareRoot_supg_semi_implicit_ale> squareroot(new SquareRoot_supg_semi_implicit_ale());
131 
132  EvaluateAtQuadrature ( elements ( M_uFESpace->mesh() ),
133  M_uFESpace->qr(),
134  M_fespaceUETA,
135  TAU_M_UPRIME * (
136  value(-M_density*M_alpha/M_timestep)*value(M_fespaceUETA, *velocity_repeated)
137  +value(M_density)*value(M_fespaceUETA, *velocity_rhs_repeated)
138  -value(M_density)*value(M_fespaceUETA, *velocity_repeated)*grad(M_fespaceUETA, *velocity_repeated)
139  -grad(M_fespacePETA, *pressure_repeated)
140  +value(M_viscosity)*laplacian(M_fespaceUETA, *velocity_repeated)
141  +value(M_density)*quadpts(M_fespaceUETA, M_fineScaleVelocityRhs )
142  )
143  ) >> M_fineScaleVelocity;
144 
145  EvaluateAtQuadrature ( elements ( M_uFESpace->mesh() ),
146  M_uFESpace->qr(),
147  M_fespacePETA,
148  value(-1.0) * TAU_C_NO_DT_PPRIME * trace( grad ( M_fespaceUETA, *velocity_repeated ) )
149  ) >> M_fineScalePressure;
150 }
151 
152 //=============================================================================================
153 void StabilizationSUPG_semi_implicit_ale::computeFineScalesForVisualization ( const vectorPtr_Type& velocity, const vectorPtr_Type& pressure )
154 {
155  QuadratureRule qr ( quadRuleTetra1pt );
156 
157  using namespace ExpressionAssembly;
158 
159  vectorPtr_Type velocity_repeated( new vector_Type( *velocity, Repeated ) );
160  vectorPtr_Type pressure_repeated( new vector_Type( *pressure, Repeated ) );
161  vectorPtr_Type velocity_rhs_repeated( new vector_Type( *M_rhsVelocity, Repeated) );
162 
163  std::shared_ptr<SquareRoot_supg_semi_implicit_ale> squareroot(new SquareRoot_supg_semi_implicit_ale());
164 
165  ComputeFineScaleVel ( elements ( M_fespaceUETA->mesh() ),
166  qr,
167  M_fespaceUETA,
168  TAU_M_UPRIME * (
169  value(-M_density*M_alpha/M_timestep)*value(M_fespaceUETA, *velocity_repeated)
170  +value(M_density)*value(M_fespaceUETA, *velocity_rhs_repeated)
171  -value(M_density)*value(M_fespaceUETA, *velocity_repeated)*grad(M_fespaceUETA, *velocity_repeated)
172  -grad(M_fespacePETA, *pressure_repeated)
173  +value(M_viscosity)*laplacian(M_fespaceUETA, *velocity_repeated)
174  +value(M_density)*quadpts(M_fespaceUETA, M_fineScaleVelocityRhs )
175  )
176  ) >> *M_fineVelocity;
177 
178  ComputeFineScalePres ( elements ( M_pFESpace->mesh() ),
179  qr,
180  M_fespacePETA,
181  value(-1.0) * TAU_C_NO_DT_PPRIME * trace( grad ( M_fespaceUETA, *velocity_repeated ) )
182  ) >> *M_finePressure;
183 
184 }
185 
186 void StabilizationSUPG_semi_implicit_ale::setExportFineScaleVelocity ( ExporterHDF5<mesh_Type> & exporter, const int& numElementsTotal )
187 {
188  int numVolumes = M_uFESpace->mesh()->numVolumes();
189  std::vector<int> id_elem;
190  std::vector<int> id_elem_scalar;
191  for ( int i = 0; i < numVolumes; ++i )
192  {
193  id_elem_scalar.push_back ( M_uFESpace->mesh()->element(i).id() );
194 
195  for ( int j = 0; j < 3; ++j )
196  {
197  id_elem.push_back ( M_uFESpace->mesh()->element(i).id() + numElementsTotal * j );
198  }
199  }
200 
201  int* pointerToDofs_scalar (0);
202  pointerToDofs_scalar = &id_elem_scalar[0];
203  std::shared_ptr<MapEpetra> map_scalar ( new MapEpetra ( -1, static_cast<int> (id_elem_scalar.size() ), pointerToDofs_scalar, M_uFESpace->map().commPtr() ) );
204 
205  int* pointerToDofs (0);
206  pointerToDofs = &id_elem[0];
207  std::shared_ptr<MapEpetra> map ( new MapEpetra ( -1, static_cast<int> (id_elem.size() ), pointerToDofs, M_uFESpace->map().commPtr() ) );
208 
209  M_fineVelocity.reset ( new vector_Type (*map, Unique ) );
210  M_fineVelocity->zero();
211 
212  M_finePressure.reset ( new vector_Type (*map, Unique ) );
213  M_finePressure->zero();
214 
215  exporter.addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "fine scale velocity", M_uFESpace, M_fineVelocity, UInt (0),
216  ExporterData< mesh_Type >::UnsteadyRegime, ExporterData< mesh_Type >::Cell );
217  exporter.addVariable ( ExporterData<RegionMesh<LinearTetra> >::ScalarField, "fine scale pressure", M_pFESpace, M_finePressure, UInt (0),
218  ExporterData< mesh_Type >::UnsteadyRegime, ExporterData< mesh_Type >::Cell );
219 }
220 
221 void StabilizationSUPG_semi_implicit_ale::buildGraphs()
222 {
223  /*
224  vector_Type velocity_extrapolated_rep( M_uFESpace->map(), Repeated);
225 
226  velocity_extrapolated_rep.zero();
227 
228  velocity_extrapolated_rep += 1;
229 
230  std::shared_ptr<SquareRoot_supg_semi_implicit_ale> squareroot(new SquareRoot_supg_semi_implicit_ale());
231 
232  {
233  using namespace ExpressionAssembly;
234 
235  // Graph for block 00
236  M_graph_block00.reset (new Epetra_FECrsGraph (Copy, * (M_uFESpace->map().map (Unique) ), 0) );
237  buildGraph ( elements ( M_uFESpace->mesh() ),
238  quadRuleTetra4pt,
239  M_fespaceUETA,
240  M_fespaceUETA,
241 
242  TAU_M*value(M_density*M_density)*value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_i), phi_j )
243  +TAU_M*value(M_density*M_density) * dot( value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_i), value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_j) )
244  +TAU_C*div(phi_i)*div(phi_j)
245  -TAU_M*value(M_density*M_viscosity)*dot( value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_i), laplacian(phi_j) )
246 
247  ) >> M_graph_block00;
248 
249  M_graph_block00->GlobalAssemble();
250  M_graph_block00->OptimizeStorage();
251 
252  // Graph for block 10
253  M_graph_block10.reset (new Epetra_FECrsGraph (Copy, * (M_pFESpace->map().map (Unique) ), 0) );
254  buildGraph ( elements (M_pFESpace->mesh() ),
255  quadRuleTetra4pt,
256  M_fespacePETA,
257  M_fespaceUETA,
258 
259  TAU_M*value(M_density*M_alpha/M_timestep)*dot( grad(phi_i), phi_j )
260  +TAU_M*value(M_density) * dot( grad(phi_i), value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_j) )
261  -TAU_M*value(M_viscosity)*dot(grad(phi_i), laplacian(phi_j))
262 
263  ) >> M_graph_block10;
264 
265  M_graph_block10->GlobalAssemble ( * (M_uFESpace->map().map (Unique) ), * (M_pFESpace->map().map (Unique) ) );
266  M_graph_block10->OptimizeStorage();
267 
268  // Graph for block 01
269  M_graph_block01.reset (new Epetra_FECrsGraph (Copy, * (M_uFESpace->map().map (Unique) ), 0) );
270  buildGraph ( elements (M_pFESpace->mesh() ),
271  quadRuleTetra4pt,
272  M_fespaceUETA,
273  M_fespacePETA,
274 
275  TAU_M*value(M_density)*dot( value(M_fespaceUETA, velocity_extrapolated_rep)*grad(phi_i), grad(phi_j) )
276 
277  ) >> M_graph_block01;
278 
279  M_graph_block01->GlobalAssemble ( * (M_pFESpace->map().map (Unique) ), * (M_uFESpace->map().map (Unique) ) );
280  M_graph_block01->OptimizeStorage();
281 
282  // Graph for block 11
283  M_graph_block11.reset (new Epetra_FECrsGraph (Copy, * (M_pFESpace->map().map (Unique) ), 0) );
284  buildGraph ( elements (M_pFESpace->mesh() ),
285  quadRuleTetra4pt,
286  M_fespacePETA,
287  M_fespacePETA,
288 
289  TAU_M*dot( grad(phi_i), grad(phi_j) )
290 
291  ) >> M_graph_block11;
292 
293  M_graph_block11->GlobalAssemble ( );
294  M_graph_block11->OptimizeStorage();
295 
296  }
297 
298  M_block_00.reset (new matrix_Type ( M_uFESpace->map(), *M_graph_block00 ) );
299  M_block_00->zero();
300  M_block_00->globalAssemble();
301 
302  M_block_01.reset (new matrix_Type ( M_uFESpace->map(), *M_graph_block01 ) );
303  M_block_01->zero();
304  M_block_01->globalAssemble( M_pFESpace->mapPtr(), M_uFESpace->mapPtr() );
305 
306  M_block_10.reset (new matrix_Type ( M_pFESpace->map(), *M_graph_block10 ) );
307  M_block_10->zero();
308  M_block_10->globalAssemble( M_uFESpace->mapPtr(), M_pFESpace->mapPtr() );
309 
310  M_block_11.reset (new matrix_Type ( M_pFESpace->map(), *M_graph_block11 ) );
311  M_block_11->zero();
312  M_block_11->globalAssemble( );
313  */
314 }
315 
316 // This will be applicied to the system matrix
317 void StabilizationSUPG_semi_implicit_ale::apply_matrix( const vector_Type& velocityExtrapolated, const vector_Type& velocityALE )
318 {
319  if ( !M_useGraph )
320  {
321  M_block_00.reset (new matrix_Type ( M_uFESpace->map() ) );
322 
323  M_block_01.reset (new matrix_Type ( M_uFESpace->map() ) );
324 
325  M_block_10.reset (new matrix_Type ( M_pFESpace->map() ) );
326 
327  M_block_11.reset (new matrix_Type ( M_pFESpace->map() ) );
328  }
329 
330  // missing force
331 
332  M_block_00->zero();
333  M_block_01->zero();
334  M_block_10->zero();
335  M_block_11->zero();
336 
337  vector_Type beta( velocityExtrapolated - velocityALE );
338  vector_Type beta_rep ( beta, Repeated );
339 
340  vector_Type velocity_extrapolated_rep( velocityExtrapolated, Repeated);
341  vector_Type velocity_ALE_rep( velocityALE, Repeated);
342 
343  std::shared_ptr<SquareRoot_supg_semi_implicit_ale> squareroot(new SquareRoot_supg_semi_implicit_ale());
344 
345  using namespace ExpressionAssembly;
346 
347  if ( M_useODEfineScale )
348  {
349  integrate(
350  elements(M_uFESpace->mesh()),
351  M_uFESpace->qr(),
352  M_fespaceUETA, // test w -> phi_i
353  M_fespaceUETA, // trial u^{n+1} -> phi_j
354 
355  TAU_M_TILDE*value(M_density*M_density)*value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), phi_j )
356  +TAU_M_TILDE*value(M_density*M_density) * dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), value(M_fespaceUETA, beta_rep)*grad(phi_j) )
357  +TAU_C_NO_DT*div(phi_i)*div(phi_j)
358  -TAU_M_TILDE*value(M_density*M_viscosity)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), laplacian(phi_j) )
359 
360  ) >> M_block_00;
361  }
362  else
363  {
364  integrate(
365  elements(M_uFESpace->mesh()),
366  M_uFESpace->qr(),
367  M_fespaceUETA, // test w -> phi_i
368  M_fespaceUETA, // trial u^{n+1} -> phi_j
369 
370  TAU_M * (
371 
372  dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), value(M_density*M_density*M_alpha/M_timestep) * phi_j
373  +value(M_density*M_density)*value(M_fespaceUETA, beta_rep)*grad(phi_j)
374  -value(M_density*M_viscosity)*laplacian(phi_j)
375  )
376 
377  )
378 
379  +TAU_C*div(phi_i)*div(phi_j)
380 
381  /*
382  TAU_M*value(M_density*M_density)*value(M_alpha/M_timestep) * dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), phi_j )
383  +TAU_M*value(M_density*M_density) * dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), value(M_fespaceUETA, beta_rep)*grad(phi_j) )
384  +TAU_C*div(phi_i)*div(phi_j)
385  -TAU_M*value(M_density*M_viscosity)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), laplacian(phi_j) )
386  */
387  ) >> M_block_00;
388  }
389 
390  M_block_00->globalAssemble();
391 
392  if ( M_useODEfineScale )
393  {
394  integrate(
395  elements(M_uFESpace->mesh()),
396  M_pFESpace->qr(),
397  M_fespacePETA, // test q -> phi_i
398  M_fespaceUETA, // trial u^{n+1} -> phi_j
399 
400  TAU_M_TILDE*value(M_density*M_alpha/M_timestep)*dot( grad(phi_i), phi_j )
401  +TAU_M_TILDE*value(M_density) * dot( grad(phi_i), value(M_fespaceUETA, beta_rep)*grad(phi_j) )
402  -TAU_M_TILDE*value(M_viscosity)*dot(grad(phi_i), laplacian(phi_j))
403 
404  ) >> M_block_10;
405  }
406  else
407  {
408  integrate(
409  elements(M_uFESpace->mesh()),
410  M_pFESpace->qr(),
411  M_fespacePETA, // test q -> phi_i
412  M_fespaceUETA, // trial u^{n+1} -> phi_j
413 
414  TAU_M * (
415  dot( grad(phi_i), value(M_density*M_alpha/M_timestep)*phi_j
416  +value(M_density)*value(M_fespaceUETA, beta_rep)*grad(phi_j)
417  -value(M_viscosity)*laplacian(phi_j)
418  )
419 
420  )
421 
422  /*
423  TAU_M*value(M_density*M_alpha/M_timestep)*dot( grad(phi_i), phi_j )
424  +TAU_M*value(M_density) * dot( grad(phi_i), value(M_fespaceUETA, beta_rep)*grad(phi_j) )
425  -TAU_M*value(M_viscosity)*dot(grad(phi_i), laplacian(phi_j))
426  */
427 
428  ) >> M_block_10;
429  }
430  M_block_10->globalAssemble( M_uFESpace->mapPtr(), M_pFESpace->mapPtr() );
431 
432  if ( M_useODEfineScale )
433  {
434  integrate(
435  elements(M_uFESpace->mesh()),
436  M_uFESpace->qr(),
437  M_fespaceUETA, // test w -> phi_i
438  M_fespacePETA, // trial p^{n+1} -> phi_j
439 
440  TAU_M_TILDE*value(M_density)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), grad(phi_j) )
441 
442  ) >> M_block_01;
443  }
444  else
445  {
446  integrate(
447  elements(M_uFESpace->mesh()),
448  M_uFESpace->qr(),
449  M_fespaceUETA, // test w -> phi_i
450  M_fespacePETA, // trial p^{n+1} -> phi_j
451 
452  TAU_M*value(M_density)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), grad(phi_j) )
453 
454  ) >> M_block_01;
455  }
456  M_block_01->globalAssemble( M_pFESpace->mapPtr(), M_uFESpace->mapPtr() );
457 
458  if ( M_useODEfineScale )
459  {
460  integrate(
461  elements(M_uFESpace->mesh()),
462  M_pFESpace->qr(),
463  M_fespacePETA, // test q -> phi_i
464  M_fespacePETA, // trial p^{n+1} -> phi_j
465 
466  TAU_M_TILDE*dot( grad(phi_j), grad(phi_i) )
467 
468  ) >> M_block_11;
469  }
470  else
471  {
472  integrate(
473  elements(M_uFESpace->mesh()),
474  M_pFESpace->qr(),
475  M_fespacePETA, // test q -> phi_i
476  M_fespacePETA, // trial p^{n+1} -> phi_j
477 
478  TAU_M*dot( grad(phi_j), grad(phi_i) )
479 
480  ) >> M_block_11;
481  }
482  M_block_11->globalAssemble();
483 
484 }
485 
486 void StabilizationSUPG_semi_implicit_ale::apply_vector( vectorPtr_Type& rhs_velocity,
487  vectorPtr_Type& rhs_pressure,
488  const vector_Type& velocityExtrapolated,
489  const vector_Type& velocityALE,
490  const vector_Type& velocity_rhs)
491 {
492  M_rhsVelocity.reset( new vector_Type ( velocity_rhs, Unique ) );
493 
494  // missing force, terms 11 and 12
495 
496  vector_Type beta( velocityExtrapolated - velocityALE );
497  vector_Type beta_rep ( beta, Repeated );
498 
499  vector_Type velocity_rhs_rep( velocity_rhs, Repeated);
500  vector_Type velocity_extrapolated_rep( velocityExtrapolated, Repeated);
501 
502  std::shared_ptr<SquareRoot_supg_semi_implicit_ale> squareroot(new SquareRoot_supg_semi_implicit_ale());
503 
504  using namespace ExpressionAssembly;
505 
506  if ( M_useODEfineScale )
507  {
508  integrate(
509  elements(M_uFESpace->mesh()),
510  M_uFESpace->qr(),
511  M_fespaceUETA,
512 
513  TAU_M_TILDE*value(M_density*M_density)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), value(M_fespaceUETA, velocity_rhs_rep) )
514  +TAU_M_TILDE*value(M_density*M_density)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), quadpts(M_fespaceUETA, M_fineScaleVelocityRhs) )
515 
516  ) >> rhs_velocity;
517 
518  integrate(
519  elements(M_uFESpace->mesh()),
520  M_pFESpace->qr(),
521  M_fespacePETA,
522 
523  TAU_M_TILDE*value(M_density)*dot( grad(phi_i), value(M_fespaceUETA, velocity_rhs_rep))
524  +TAU_M_TILDE*value(M_density)*dot( grad(phi_i), quadpts(M_fespaceUETA, M_fineScaleVelocityRhs))
525 
526  ) >> rhs_pressure;
527 
528 
529  }
530  else
531  {
532  integrate(
533  elements(M_uFESpace->mesh()),
534  M_uFESpace->qr(),
535  M_fespaceUETA,
536 
537  TAU_M*value(M_density*M_density)*dot( value(M_fespaceUETA, beta_rep)*grad(phi_i), value(M_fespaceUETA, velocity_rhs_rep) )
538 
539  ) >> rhs_velocity;
540 
541  integrate(
542  elements(M_uFESpace->mesh()),
543  M_pFESpace->qr(),
544  M_fespacePETA,
545 
546  TAU_M*value(M_density)*dot( grad(phi_i), value(M_fespaceUETA, velocity_rhs_rep))
547 
548  ) >> rhs_pressure;
549  }
550 
551 
552 }
553 
554 } // namespace LifeV
TAU_M_DEN_UPRIME
#define TAU_M_DEN_UPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:22
TAU_C_NO_DT
#define TAU_C_NO_DT
Definition: StabilizationSUPG_semi_implicit_ale.cpp:18
TAU_M_DEN
#define TAU_M_DEN
Definition: StabilizationSUPG.cpp:5
TAU_M_DEN_VISC_UPRIME
#define TAU_M_DEN_VISC_UPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:24
TAU_M_DEN_DT_NOSQUARED
#define TAU_M_DEN_DT_NOSQUARED
Definition: StabilizationSUPG_semi_implicit_ale.cpp:11
TAU_M_DEN_VEL
#define TAU_M_DEN_VEL
Definition: StabilizationSUPG.cpp:7
TAU_M_DEN_VEL_UPRIME
#define TAU_M_DEN_VEL_UPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:23
TAU_M_UPRIME
#define TAU_M_UPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:21
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
TAU_M_DEN_VISC
#define TAU_M_DEN_VISC
Definition: StabilizationSUPG.cpp:8
TAU_C
#define TAU_C
Definition: StabilizationSUPG.cpp:11
TAU_M
#define TAU_M
Definition: StabilizationSUPG.cpp:4
TAU_M_DEN_DT
#define TAU_M_DEN_DT
Definition: StabilizationSUPG.cpp:6
TAU_C_NO_DT_PPRIME
#define TAU_C_NO_DT_PPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:28
TAU_M_NO_DT
#define TAU_M_NO_DT
Definition: StabilizationSUPG_semi_implicit_ale.cpp:17
TAU_M_NO_DT_UPRIME
#define TAU_M_NO_DT_UPRIME
Definition: StabilizationSUPG_semi_implicit_ale.cpp:27
TAU_M_TILDE
#define TAU_M_TILDE
Definition: StabilizationSUPG_semi_implicit_ale.cpp:10