LifeV
ALESolver.cpp
Go to the documentation of this file.
1 
2 #include <lifev/fsi_blocks/solver/ALESolver.hpp>
3 
4 namespace LifeV
5 {
6 
7 // ===================================================
8 // Constructors & Destructor
9 // ===================================================
10 
11 
12 ALESolver::ALESolver ( FESpace<mesh_Type, MapEpetra>& mmFESpace,
13  std::shared_ptr<Epetra_Comm> comm ) :
14  M_FESpace ( mmFESpace ),
15  M_localMap ( M_FESpace.map() ),
16  M_matrHE ( new matrix_Type (M_localMap ) ),
17  M_displayer ( comm ),
18  M_me ( comm->MyPID() ),
19  M_verbose ( M_me == 0 ),
20  M_elmat ( M_FESpace.fe().nbFEDof(), nDimensions, nDimensions ),
21  M_disp ( ),
22  M_secondRHS ( ),
23  M_diffusion ( 1. ),
24  M_offset (0),
25  M_linearElasticity (false),
26  M_young (0.0),
27  M_poisson (0.0)
28 {
29 }
30 
31 ALESolver::ALESolver ( FESpace<mesh_Type, MapEpetra>& mmFESpace,
32  std::shared_ptr<Epetra_Comm> comm ,
33  MapEpetra& localMap,
34  UInt offset) :
35  M_FESpace ( mmFESpace ),
36  M_localMap ( localMap),
37  M_matrHE ( new matrix_Type (M_localMap ) ),
38  M_displayer ( comm ),
39  M_me ( comm->MyPID() ),
40  M_verbose ( M_me == 0 ),
41  M_elmat ( M_FESpace.fe().nbFEDof(), nDimensions, nDimensions ),
42  M_secondRHS ( ),
43  M_diffusion ( 1. ),
44  M_offset (offset),
45  M_linearElasticity (false),
46  M_young (0.0),
47  M_poisson (0.0)
48 {
49 }
50 
51 // ===================================================
52 // Methods
53 // ===================================================
54 
55 
56 void ALESolver::setUp ( const GetPot& dataFile )
57 {
58  M_linearElasticity = dataFile ("ale/useLinearElasticity", false);
59 
60  M_diffusion = dataFile ("ale/diffusion", 1.0);
61 
62  M_young = dataFile ("solid/model/young", 0.0);
63 
64  M_poisson = dataFile ("solid/model/poisson", 0.0);
65 
66  computeMatrix( );
67 
68  M_secondRHS.reset (new vector_Type (M_FESpace.map() ) );
69  M_disp.reset (new vector_Type (M_FESpace.map() ) );
70 }
71 
72 void ALESolver::iterate ( BCHandler& BCh )
73 {
74  LifeChrono chrono;
75 
76  // matrix and vector assembling communication
77  M_displayer.leaderPrint (" ALE- Applying boundary conditions ... ");
78  chrono.start();
79 
80  *M_secondRHS *= 0.;
81  applyBoundaryConditions (*M_secondRHS, BCh);
82 
83  chrono.stop();
84  M_displayer.leaderPrintMax ("done in " , chrono.diff() );
85 }
86 
87 void ALESolver::applyBoundaryConditions (vector_Type& rhs, BCHandler& BCh)
88 {
89  // CHANGED BY S. QUINODOZ !
90  // "if" exchanged
91 
92  if ( ! BCh.bcUpdateDone() )
93  {
94  // BC boundary information update
95  BCh.bcUpdate ( *M_FESpace.mesh(), M_FESpace.feBd(), M_FESpace.dof() );
96  }
97 
98  if (M_offset) //mans that this is the fullMonolithic case
99  {
100  BCh.setOffset (M_offset);
101  }
102  else
103  {
104  bcManageRhs (rhs, *M_FESpace.mesh(), M_FESpace.dof(), BCh, M_FESpace.feBd(), 1., 0.0);
105  }
106 
107  bcManageMatrix ( *M_matrHE, *M_FESpace.mesh(), M_FESpace.dof(), BCh, M_FESpace.feBd(), 1.0, 0. );
108 }
109 
110 void ALESolver::applyBoundaryConditions (BCHandler& BCh)
111 {
112 
113  if ( ! BCh.bcUpdateDone() )
114  {
115  BCh.bcUpdate ( *M_FESpace.mesh(), M_FESpace.feBd(), M_FESpace.dof() );
116  }
117 
118  bcManageMatrix ( *M_matrHE_BC, *M_FESpace.mesh(), M_FESpace.dof(), BCh, M_FESpace.feBd(), 1.0, 0. );
119 }
120 
121 void ALESolver::computeMatrix( )
122 {
123  LifeChrono chrono;
124  chrono.start();
125  M_displayer.leaderPrint (" ALE- Computing constant matrices ... ");
126 
127  M_matrHE.reset ( new matrix_Type (M_localMap ) );
128 
129  if ( M_linearElasticity )
130  {
131  using namespace ExpressionAssembly;
132  M_aleFESpace_ETA.reset( new ETFESpaceAle_Type ( M_FESpace.mesh(), &(M_FESpace.refFE()), M_FESpace.map().commPtr() ) );
133 
134  Real lambda = ( M_young * M_poisson ) / ( ( 1.0 + M_poisson ) * ( 1.0 - 2.0 * M_poisson ) );
135  Real mu = M_young / ( 2.0 * ( 1.0 + M_poisson ) );
136 
137  integrate ( elements (M_aleFESpace_ETA->mesh() ),
138  M_FESpace.qr(),
139  M_aleFESpace_ETA,
140  M_aleFESpace_ETA,
141  pow( value(0.0228)/detJ, 0.8) * (
142  value( mu ) * dot ( grad(phi_i) , grad(phi_j) + transpose( grad(phi_j) ) ) +
143  value( lambda ) * div(phi_i) * div(phi_j)
144  )
145  ) >> M_matrHE;
146  }
147  else
148  {
149  UInt totalDof = M_FESpace.dof().numTotalDof();
150  // Loop on elements
151  for ( UInt i = 0; i < M_FESpace.mesh()->numVolumes(); ++i )
152  {
153  // Updating derivatives
154  M_FESpace.fe().updateFirstDerivQuadPt ( M_FESpace.mesh()->volumeList ( i ) );
155  M_elmat.zero();
156  AssemblyElemental::stiff ( M_diffusion, M_elmat, M_FESpace.fe(), 0, 0, 3 );
157  // Assembling
158  for ( UInt j = 0; j < M_FESpace.fieldDim(); ++j )
159  {
160  assembleMatrix ( *M_matrHE, M_elmat, M_FESpace.fe(), M_FESpace.dof(), j, j, j * totalDof + M_offset, j * totalDof + M_offset );
161  }
162  }
163  }
164 
165  M_matrHE->globalAssemble();
166 
167  M_matrHE_BC.reset ( new matrix_Type (M_localMap ) );
168  M_matrHE_BC->zero();
169  *M_matrHE_BC += *M_matrHE;
170 
171  chrono.stop();
172  M_displayer.leaderPrintMax ("done in " , chrono.diff() );
173 
174 }
175 
176 void ALESolver::updateShapeDerivatives ( Real& alpha,
177  const Real& density,
178  const Real& viscosity,
179  const vector_Type& un,
180  const vector_Type& uk,
181  //const vector_Type& disp,
182  const vector_Type& w,
183  FESpace<mesh_Type, MapEpetra>& velocityFESpace,
184  FESpace<mesh_Type, MapEpetra>& pressureFESpace,
185  bool wImplicit,
186  bool convectiveTermDerivative,
187  BCHandler& BCh )
188 
189 {
190  LifeChrono chrono;
191  chrono.start();
192  M_displayer.leaderPrint (" ALE- Computing shapeDerivatives ... ");
193 
194  UInt numVelocityComponent = nDimensions;
195 
196  // Velocity + Pressure. Probably needs also FESpaces for u and p
197  M_matrShapeDerVel.reset ( new matrix_Type ( velocityFESpace.map() ) );
198  M_matrShapeDerPressure.reset ( new matrix_Type ( pressureFESpace.map() ) );
199 
200  VectorElemental elementConvectionVelocity( velocityFESpace.fe().nbFEDof(), nDimensions );
201  VectorElemental elementMeshVelocity ( velocityFESpace.fe().nbFEDof(), nDimensions );
202  VectorElemental u_loc ( velocityFESpace.fe().nbFEDof(), nDimensions );
203  VectorElemental elementPressure ( pressureFESpace.fe().nbFEDof(), 1 );
204  VectorElemental elementVelocity ( velocityFESpace.fe().nbFEDof(), nDimensions );
205 
206 
207  vector_Type unRepeated ( un , Repeated );
208  vector_Type ukRepeated ( uk , Repeated );
209  // vector_Type dispRepeated( disp, Repeated );
210  vector_Type wRepeated ( w , Repeated );
211  // vector_Type dwRepeated ( dw , Repeated );
212 
213 
214  for ( UInt i = 0; i < this->M_FESpace.mesh()->numVolumes(); i++ )
215  {
216 
217  pressureFESpace.fe().update ( pressureFESpace.mesh()->volumeList ( i ) );
218  velocityFESpace.fe().updateFirstDerivQuadPt ( velocityFESpace.mesh()->volumeList ( i ) );
219  pressureFESpace.fe().updateFirstDerivQuadPt ( velocityFESpace.mesh()->volumeList ( i ) );
220 
221  // just to provide the id number in the assem_mat_mixed
222  //pressureFESpace.fe().updateFirstDeriv( velocityFESpace.mesh()->volumeList( i ) );
223  //as updateFirstDer
224  //velocityFESpace.fe().updateFirstDeriv( velocityFESpace.mesh()->volumeList( i ) );
225  M_FESpace.fe().updateFirstDerivQuadPt ( M_FESpace.mesh()->volumeList ( i ) );
226 
227  // initialization of elementary vectors
228  std::shared_ptr<MatrixElemental> elementMatrixPressure ( new MatrixElemental ( pressureFESpace.fe().nbFEDof(),
229  1,
230  0,
231  M_FESpace.fe().nbFEDof(),
232  0,
233  nDimensions ) );
234  std::shared_ptr<MatrixElemental> elementMatrixVelocity ( new MatrixElemental ( velocityFESpace.fe().nbFEDof(),
235  nDimensions,
236  0,
237  velocityFESpace.fe().nbFEDof(),
238  0,
239  nDimensions ) );
240  std::shared_ptr<MatrixElemental> elementMatrixConvective;
241 
242  if ( convectiveTermDerivative )
243  {
244  elementMatrixConvective.reset ( new MatrixElemental ( velocityFESpace.fe().nbFEDof(),
245  nDimensions,
246  0,
247  velocityFESpace.fe().nbFEDof(),
248  0,
249  nDimensions ) );
250  elementMatrixConvective->zero();
251  }
252 
253  elementMatrixPressure->zero();
254  elementMatrixVelocity->zero();
255 
256  for ( UInt iNode = 0 ; iNode < velocityFESpace.fe().nbFEDof() ; iNode++ )
257  {
258  UInt iLocal = velocityFESpace.fe().patternFirst ( iNode ); // iLocal = iNode
259 
260  for ( UInt iComponent = 0; iComponent < numVelocityComponent; ++iComponent )
261  {
262  UInt iGlobal = velocityFESpace.dof().localToGlobalMap ( i, iLocal ) + iComponent * velocityFESpace.dim();
263 
264  // if(!wImplicit)
265  // u^n - w^iNode local
266  elementConvectionVelocity.vec() [ iLocal + iComponent * velocityFESpace.fe().nbFEDof() ] = unRepeated (iGlobal)
267  - wRepeated ( iGlobal );
268  // else
269  // u^n - w^iNode local
270  // elementConvectionVelocity.vec() [ iLocal + iComponent*velocityFESpace.fe().nbFEDof() ] = ukRepeated(iGlobal)
271  // - wRepeated(iGlobal);
272  // w^iNode local
273  elementMeshVelocity.vec( ) [ iLocal + iComponent * velocityFESpace.fe().nbFEDof() ] = wRepeated ( iGlobal );
274  // u^iNode local
275  elementVelocity.vec( ) [ iLocal + iComponent * velocityFESpace.fe().nbFEDof() ] = ukRepeated ( iGlobal );
276  // dw local
277  //M_elementDisplacement.vec( ) [ iLocal + iComponent*velocityFESpace.fe().nbFEDof() ] = dispRepeated( iGlobal );
278  // dw local
279  //elementVelocityRightHandSide.vec( ) [ iLocal + iComponent*velocityFESpace.fe().nbFEDof() ] = dwRepeated( iGlobal );
280  // un local
281  u_loc.vec() [ iLocal + iComponent * velocityFESpace.fe().nbFEDof() ] = unRepeated ( iGlobal );
282  }
283  }
284  /*
285  std::cout << elementConvectionVelocity.vec() << std::endl;
286  std::cout << elementMeshVelocity.vec() << std::endl;
287  std::cout << elementVelocity.vec() << std::endl;
288  std::cout << M_elementDisplacement.vec() << std::endl;
289  std::cout << elementVelocityRightHandSide.vec() << std::endl;
290  std::cout << u_loc.vec() << std::endl;
291  */
292  for ( UInt iNode = 0 ; iNode < pressureFESpace.fe().nbFEDof() ; iNode++ )
293  {
294  // iLocal = iNode
295  UInt iLocal = pressureFESpace.fe().patternFirst ( iNode );
296  UInt iGlobal = pressureFESpace.dof().localToGlobalMap ( i, iLocal ) + numVelocityComponent * velocityFESpace.dim();
297  // p^iNode local
298  elementPressure[ iLocal ] = ukRepeated[ iGlobal ];
299  }
300 
301  AssemblyElemental::shape_terms ( //M_elementDisplacement,
302  density,
303  viscosity,
304  u_loc,
305  elementVelocity,
306  elementMeshVelocity,
307  elementConvectionVelocity,
308  elementPressure,
309  *elementMatrixVelocity,
310  M_FESpace.fe(),
311  pressureFESpace.fe(),
312  (ID) M_FESpace.fe().nbFEDof(),
313  *elementMatrixPressure,
314  0,
315  wImplicit,
316  alpha//,
317  //elementMatrixConvective
318  );
319 
320  //elementMatrixVelocity->showMe(std::cout);
321 
322  /*
323  source_mass2( density,
324  elementVelocity,
325  *M_elementMatrixConvective,
326  velocityFESpace.fe(),
327  alpha );
328  */
329 
330  AssemblyElemental::source_press ( 1.0,
331  elementVelocity,
332  *elementMatrixPressure,
333  M_FESpace.fe(),
334  pressureFESpace.fe(),
335  (ID) M_FESpace.fe().nbFEDof() );
336 
337  //derivative of the convective term
338  if ( convectiveTermDerivative )
339  AssemblyElemental::mass_gradu ( density,
340  elementVelocity,
341  *elementMatrixConvective,
342  velocityFESpace.fe() );
343  /*
344  std::cout << "source_press -> norm_inf( M_elementVectorVelocity )" << std::endl;
345  M_elementVectorPressure.showMe( std::cout );
346  */
347  //
348  // Assembling
349  //
350  /*
351  std::cout << "debut ====================" << std::endl;
352  M_elementVectorPressure.showMe( std::cout );
353  M_elementVectorVelocity.showMe( std::cout );
354  std::cout << "fin ====================" << std::endl;
355  */
356  UInt const velocityTotalDof ( velocityFESpace.dof().numTotalDof() );
357  UInt const meshTotalDof ( M_FESpace.dof().numTotalDof() );
358  for ( UInt iComponent = 0; iComponent < numVelocityComponent; ++iComponent )
359  {
360  for ( UInt jComponent = 0; jComponent < numVelocityComponent; ++jComponent )
361  {
362  assembleMatrix ( *M_matrShapeDerVel,
363  *elementMatrixVelocity,
364  velocityFESpace.fe(),
365  M_FESpace.fe(),
366  velocityFESpace.dof(),
367  M_FESpace.dof(),
368  iComponent,
369  jComponent,
370  iComponent * velocityTotalDof,
371  jComponent * meshTotalDof);
372 
373  //assembling the derivative of the convective term
374  if ( convectiveTermDerivative )
375  assembleMatrix ( *M_matrShapeDerVel,
376  *elementMatrixConvective,
377  velocityFESpace.fe(),
378  velocityFESpace.fe(),
379  velocityFESpace.dof(),
380  velocityFESpace.dof(),
381  iComponent,
382  jComponent,
383  iComponent * velocityTotalDof,
384  jComponent * velocityTotalDof );
385  }
386 
387  assembleMatrix ( *M_matrShapeDerPressure,
388  *elementMatrixPressure,
389  pressureFESpace.fe(),
390  M_FESpace.fe(),
391  pressureFESpace.dof(),
392  M_FESpace.dof(),
393  (UInt) 0,
394  iComponent,
395  (UInt) 0,
396  iComponent * meshTotalDof );
397  }
398  }
399 
400  if ( ! BCh.bcUpdateDone() )
401  {
402  // BC boundary information update
403  BCh.bcUpdate ( *velocityFESpace.mesh(), velocityFESpace.feBd(), velocityFESpace.dof() );
404  }
405 
406  bcManageMatrix ( *M_matrShapeDerVel, *velocityFESpace.mesh(), velocityFESpace.dof(), BCh, velocityFESpace.feBd(), 0.0, 0.0 );
407 
408  M_matrShapeDerVel->globalAssemble(M_FESpace.mapPtr(), velocityFESpace.mapPtr());
409  M_matrShapeDerPressure->globalAssemble(M_FESpace.mapPtr(), pressureFESpace.mapPtr());
410 
411 
412  chrono.stop();
413  M_displayer.leaderPrintMax ("done in ", chrono.diff() );
414 }
415 
416 
417 } // end namespace LifeV