LifeV
lifev/fsi/examples/challenge_VPH/main.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  * @include fluidstructure.dox
29  * @file
30  * @brief for testing the benchmark contained in the
31  *
32  * @date 2009-04-09
33  * @author Paolo Crosetto <crosetto@iacspc70.epfl.ch>
34  *
35  * @maintainer Paolo Crosetto <crosetto@iacspc70.epfl.ch>
36  *
37  * Monolithic problem. Features:
38  * - fullMonolithic (CE):
39  * -# solution with exact Newton (semiImplicit = false, useShapeDerivatives = true, domainVelImplicit = false, convectiveTermDer = false)
40  * -# solution with quasi Newton (semiImplicit = false, useShapeDerivatives = false, domainVelImplicit = false, convectiveTermDer = false)
41  * -# preconditioner choice: see the classes Monolithic and fullMonolithic
42  * - Monolithic (GCE):
43  * -# solution extrapolating the fluid domain (semiImplicit = true, useShapeDerivatives = false, domainVelImplicit = false, convectiveTermDer = false)
44  * -# preconditioner choice: see the classes Monolithic and fullMonolithic
45  * - boundary conditions:
46  * -# Neumann
47  * -# Dirichlet
48  * -# Robin
49  * -# Fluxes (defective)
50  * -# absorbing \cite BadiaNobileVergara2008 :
51  * through the class flowConditions.
52  * - optional: computation of wall shear stress (not properly tested in parallel)
53  * - optional: computation of the largest singular values of the preconditioned matrix
54  *
55  * \b Features:
56  * This test by default solves the FSI probem discretized in time using the GCE or CE methods, implemented respectively
57  * in the files monolithicGE.hpp and monolithicGI.hpp . The geometry is that of a tube (benchmark test introduced in \cite Gerbeau2003).
58  * In this test the boundary conditions assigned are of type:
59  * - flux (defective b.c.) at the inlet
60  * - absorbing (see \cite BadiaNobileVergara2008) at the outlet
61  * - Robin b.c. on the solid external wall
62  * - Dirichlet homogeneous at the solid rings on the inlet-outlet (clamped tube).
63  *
64  * The output is written at every timestep, in HDF5 (if available) or ensight format.
65  * This test implements an inlet flux bundary condition for the first three time steps, then at the fourth time step
66  * the inlet boundary condition is replaced by a Neumann one (this mechanism is useful to implement rudimental valves).
67  * The outflow boundary condition is of absorbing type. At the outer wall for the structure a Robin condition is imposed.
68  */
69 
70 
71 #include <cassert>
72 #include <cstdlib>
73 
74 #include <boost/timer.hpp>
75 
76 #include <Epetra_ConfigDefs.h>
77 #ifdef EPETRA_MPI
78 #include <mpi.h>
79 #include <Epetra_MpiComm.h>
80 #else
81 #include <Epetra_SerialComm.h>
82 #endif
83 
84 
85 #include "ud_functions.hpp"
86 // LifeV includes
87 #include <lifev/core/fem/FESpace.hpp>
88 #include <lifev/core/array/MapEpetra.hpp>
89 #include <lifev/core/fem/BCHandler.hpp>
90 #include <lifev/core/LifeV.hpp>
91 
92 #include <lifev/core/algorithm/PreconditionerIfpack.hpp>
93 #include <lifev/core/algorithm/PreconditionerML.hpp>
94 
95 #include <lifev/fsi/solver/FSISolver.hpp>
96 #include <lifev/structure/solver/StructuralOperator.hpp>
97 #include <lifev/fsi/solver/FSIMonolithicGI.hpp>
98 
99 #include <lifev/core/filter/ExporterEnsight.hpp>
100 #include <lifev/core/filter/ExporterEmpty.hpp>
101 #ifdef HAVE_HDF5
102 #include <lifev/core/filter/ExporterHDF5.hpp>
103 #endif
104 
105 #include "ud_functions.hpp"
106 #include "boundaryConditions.hpp"
107 
108 class Problem
109 {
110 public:
111 
112  typedef std::shared_ptr<LifeV::FSISolver> FSISolverPtr_Type;
113 
114  typedef LifeV::FSIOperator::data_Type data_Type;
115  typedef LifeV::FSIOperator::dataPtr_Type dataPtr_Type;
116 
117  typedef LifeV::FSIOperator::vector_Type vector_Type;
118  typedef LifeV::FSIOperator::vectorPtr_Type vectorPtr_Type;
119 
120  typedef std::shared_ptr< LifeV::Exporter<LifeV::RegionMesh<LifeV::LinearTetra> > > filterPtr_Type;
121 
122  typedef LifeV::ExporterEnsight<LifeV::FSIOperator::mesh_Type> ensightFilter_Type;
123  typedef std::shared_ptr<ensightFilter_Type> ensightFilterPtr_Type;
124 #ifdef HAVE_HDF5
125  typedef LifeV::ExporterHDF5<LifeV::FSIOperator::mesh_Type> hdf5Filter_Type;
126  typedef std::shared_ptr<hdf5Filter_Type> hdf5FilterPtr_Type;
127 #endif
128  typedef LifeV::FactorySingleton<LifeV::Factory<LifeV::FSIOperator, std::string> > FSIFactory_Type;
129  typedef LifeV::RegionMesh<LifeV::LinearTetra> mesh_Type;
130  typedef LifeV::MapEpetra map_Type;
131  typedef LifeV::FESpace<mesh_Type, map_Type> fespace_Type;
132  /*!
133  This routine sets up the problem:
134 
135  -# create the standard boundary conditions for the fluid and
136  structure problems.
137 
138  -# initialize and setup the FSIsolver
139  */
140 
141  Problem ( GetPot const& data_file ) :
142  M_Tstart (0.),
143  M_saveEvery (1)
144  {
145  using namespace LifeV;
146 
147  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "linearVenantKirchhoff", &FSIOperator::createVenantKirchhoffLinear );
148  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "exponential", &FSIOperator::createExponentialMaterialNonLinear );
149  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "neoHookean", &FSIOperator::createNeoHookeanMaterialNonLinear );
150  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "nonLinearVenantKirchhoff", &FSIOperator::createVenantKirchhoffNonLinear );
151 
152  std::cout << "register MonolithicGE : " << FSIMonolithicGE::S_register << std::endl;
153  std::cout << "register MonolithicGI : " << FSIMonolithicGI::S_register << std::endl;
154  //bool reg=FSIMonolithicGI::S_register&&FSIMonolithicGE::S_register;
155 
156  M_data = dataPtr_Type ( new data_Type() );
157  M_data->setup ( data_file );
158 
159  M_fsi = FSISolverPtr_Type ( new FSISolver( ) );
160  MPI_Barrier ( MPI_COMM_WORLD );
161 
162  M_fsi->setData ( M_data );
163  M_fsi->FSIOper()->setDataFile ( data_file ); //TO BE REMOVED!
164 
165  // Setting FESpace and DOF
166 
167  std::string fluidMeshPartitioned = data_file ( "problem/fluidMeshPartitioned", "none" );
168  std::string solidMeshPartitioned = data_file ( "problem/solidMeshPartitioned", "none" );
169 #ifdef HAVE_HDF5
170  if ( fluidMeshPartitioned.compare ( "none" ) )
171  {
172  FSIOperator::meshFilter_Type fluidMeshFilter ( data_file, fluidMeshPartitioned );
173  fluidMeshFilter.setComm ( M_fsi->FSIOper()->worldComm() );
174  FSIOperator::meshFilter_Type solidMeshFilter ( data_file, solidMeshPartitioned );
175  solidMeshFilter.setComm ( M_fsi->FSIOper( )->worldComm( ) );
176  M_fsi->FSIOper( )->partitionMeshes ( fluidMeshFilter, solidMeshFilter );
177  M_fsi->FSIOper( )->setupFEspace( );
178  M_fsi->FSIOper( )->setupDOF ( fluidMeshFilter );
179  fluidMeshFilter.closeFile( );
180  solidMeshFilter.closeFile( );
181  }
182  else
183 #endif
184  {
185  M_fsi->FSIOper( )->partitionMeshes( );
186  M_fsi->FSIOper( )->setupFEspace( );
187  M_fsi->FSIOper( )->setupDOF( );
188  }
189 
190  MPI_Barrier ( MPI_COMM_WORLD );
191 
192  M_fsi->setFluidBC ( BCh_monolithicFlux( ) );
193  M_fsi->setSolidBC ( BCh_monolithicSolid ( *M_fsi->FSIOper( ) ) );
194 
195  M_fsi->setup();
196 
197  M_fsi->setFluidBC ( BCh_monolithicFluid ( *M_fsi->FSIOper( ) ) );
198  M_fsi->setHarmonicExtensionBC ( BCh_harmonicExtension ( *M_fsi->FSIOper( ) ) );
199 
200  dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->mergeBCHandlers();
201 
202  std::string const exporterType = data_file ( "exporter/type", "ensight" );
203  std::string const fluidName = data_file ( "exporter/fluid/filename", "fluid" );
204  std::string const solidName = data_file ( "exporter/solid/filename", "solid" );
205 
206 #ifdef HAVE_HDF5
207  if (exporterType.compare ("hdf5") == 0)
208  {
209  M_exporterFluid.reset ( new hdf5Filter_Type ( data_file, fluidName) );
210  M_exporterSolid.reset ( new hdf5Filter_Type ( data_file, solidName) );
211  }
212  else
213 #endif
214  {
215  if (exporterType.compare ("none") == 0)
216  {
217  M_exporterFluid.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), fluidName, M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
218  M_exporterSolid.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->dFESpace().mesh(), solidName, M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
219  }
220  else
221  {
222  M_exporterFluid.reset ( new ensightFilter_Type ( data_file, fluidName) );
223  M_exporterSolid.reset ( new ensightFilter_Type ( data_file, solidName) );
224  }
225  }
226 
227 
228  // load using ensight/hdf5
229  M_saveEvery = data_file ("exporter/saveEvery", 1);
230 
231  // load using ensight/hdf5
232  std::string restartType (data_file ("importer/restartType", "newSimulation") );
233  std::cout << "The load state is: " << restartType << std::endl;
234 
235  if (!restartType.compare ("Stokes") )
236  {
237  initializeStokes (data_file);
238  }
239  else if (!restartType.compare ("restartFSI") )
240  {
241  restartFSI (data_file);
242  }
243  else
244  {
245  M_fsi->initialize();
246  M_Tstart = data_file ( "fluid/time_discretization/initialtime", 0.);
247 
248  M_velAndPressure.reset ( new vector_Type ( M_fsi->FSIOper()->fluid().getMap(), M_exporterFluid->mapType() ) );
249 
250  M_fluidDisp.reset ( new vector_Type ( M_fsi->FSIOper()->mmFESpace().map(), M_exporterFluid->mapType() ) );
251 
252  M_solidDisp.reset ( new vector_Type ( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ) );
253  }
254 
255 
256  M_exporterFluid->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
257  M_exporterSolid->setMeshProcId (M_fsi->FSIOper()->dFESpace().mesh(), M_fsi->FSIOper()->dFESpace().map().comm().MyPID() );
258  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "f-velocity",
259  M_fsi->FSIOper()->uFESpacePtr(), M_velAndPressure, UInt (0) );
260  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::ScalarField, "f-pressure",
261  M_fsi->FSIOper()->pFESpacePtr(), M_velAndPressure,
262  UInt (3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() ) );
263 
264  M_exporterFluid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "f-displacement",
265  M_fsi->FSIOper()->mmFESpacePtr(), M_fluidDisp, UInt (0) );
266 
267 
268 
269  // M_solidVel.reset ( new vector_Type( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ));
270  // M_WS.reset ( new vector_Type( M_fsi->FSIOper()->dFESpace().map(), M_exporterSolid->mapType() ));
271 
272  M_exporterSolid->addVariable ( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-displacement",
273  M_fsi->FSIOper()->dFESpacePtr(), M_solidDisp, UInt (0) );
274  // M_exporterSolid->addVariable( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-velocity",
275  // M_fsi->FSIOper()->dFESpacePtr(), M_solidVel, UInt(0) );
276  // M_exporterSolid->addVariable( ExporterData<FSIOperator::mesh_Type>::VectorField, "s-ws",
277  // M_fsi->FSIOper()->dFESpacePtr(), M_WS, UInt(0) );
278 
279 
280 
281  M_data->dataFluid()->dataTime()->setInitialTime ( M_Tstart );
282  M_data->dataFluid()->dataTime()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
283  M_data->dataSolid()->dataTime()->setInitialTime ( M_Tstart );
284  M_data->dataSolid()->dataTime()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
285  M_data->timeDataALE()->setInitialTime ( M_Tstart );
286  M_data->timeDataALE()->setTime ( M_data->dataFluid()->dataTime()->initialTime() );
287  }
288 
289  /*!
290  This routine runs the temporal loop
291  */
292  void
293  run()
294  {
295  boost::timer _overall_timer;
296 
297  LifeV::UInt iter = 1;
298  LifeV::UInt offset = dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->offset();
299 
300  dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->enableStressComputation (1);
301 
302  // #ifdef HAVE_HDF5
303  // if (M_exporterFluid->mapType() == LifeV::Unique)
304  // {
305  // M_exporterFluid->postProcess( M_Tstart-M_data->dataFluid()->dataTime()->getTimeStep() );//ugly way to avoid that hdf5 starts with a deformed mesh
306  // M_exporterSolid->postProcess( M_Tstart-M_data->dataSolid()->dataTime()->getTimeStep() );//ugly way to avoid that hdf5 starts with a deformed mesh
307  // }
308  // #endif
309 
310  for ( ; M_data->dataFluid()->dataTime()->canAdvance(); M_data->dataFluid()->dataTime()->updateTime(), M_data->dataSolid()->dataTime()->updateTime(), ++iter)
311  {
312  boost::timer _timer;
313 
314  if (iter % M_saveEvery == 0)
315  {
316  M_fsi->FSIOper()->exportSolidDisplacement (*M_solidDisp); // displacement(), M_offset);
317  //M_fsi->FSIOper()->exportSolidVelocity(*M_solidVel);// displacement(), M_offset);
318 
319  M_fsi->FSIOper()->exportFluidVelocityAndPressure (*M_velAndPressure);
320  M_exporterSolid->postProcess ( M_data->dataFluid()->dataTime()->time() );
321 
322  //*M_WS= *(dynamic_cast<LifeV::FSIMonolithic*>(M_fsi->FSIOper().get())->/*WS());//*/computeStress());
323 
324  *M_fluidDisp = M_fsi->FSIOper()->meshDisp();
325  M_exporterFluid->postProcess ( M_data->dataFluid()->dataTime()->time() );
326  }
327 
328  M_fsi->iterate();
329 
330 
331  M_fsi->FSIOper()->displayer().leaderPrintMax ("[fsi_run] Iteration ", iter);
332  M_fsi->FSIOper()->displayer().leaderPrintMax (" was done in : ", _timer.elapsed() );
333 
334  // std::cout << "solution norm " << iter << " : "
335  // << M_fsi->displacement().norm2() << "\n";
336 
337  }
338  if (M_data->method().compare ("monolithicGI") )
339  {
340  vectorPtr_Type solution ( new vector_Type ( (*M_fsi->FSIOper()->couplingVariableMap() ) ) );
341  M_fsi->FSIOper()->extrapolation ( *solution );
342  M_fsi->FSIOper()->iterateMesh (*solution);
343 
344  M_solidDisp->subset (*solution, offset);
345  // M_solidVel->subset(M_fsi->FSIOper()->solid().velocity(), offset);
346  // *M_solidDisp *= 1/(M_fsi->FSIOper()->solid().rescaleFactor()*M_data->dataFluid()->dataTime()->getTimeStep());
347  // *M_solidVel *= 1/(M_fsi->FSIOper()->solid().rescaleFactor()*M_data->dataFluid()->dataTime()->getTimeStep());
348 
349  *M_velAndPressure = *solution;
350  M_exporterSolid->postProcess ( M_data->dataFluid()->dataTime()->time() );
351  *M_fluidDisp = M_fsi->FSIOper()->meshMotion().disp();
352  M_exporterFluid->postProcess ( M_data->dataFluid()->dataTime()->time() );
353  }
354 
355  std::cout << "Total computation time = "
356  << _overall_timer.elapsed() << "s" << "\n";
357 
358  }
359 
360 private:
361 
362  void initializeStokes ( GetPot const& data_file);
363  void restartFSI ( GetPot const& data_file);
364 
365  FSISolverPtr_Type M_fsi;
366  dataPtr_Type M_data;
367 
368  filterPtr_Type M_exporterSolid;
369  filterPtr_Type M_exporterFluid;
370  filterPtr_Type M_importerSolid;
371  filterPtr_Type M_importerFluid;
372  vectorPtr_Type M_velAndPressure;
373  vectorPtr_Type M_fluidDisp;
374  vectorPtr_Type M_solidDisp;
375 
376  std::vector<vectorPtr_Type> M_solidStencil;
377  std::vector<vectorPtr_Type> M_fluidStencil;
378  std::vector<vectorPtr_Type> M_ALEStencil;
379 
380  //vectorPtr_Type M_solidVel;
381  LifeV::Real M_Tstart;
382  // vectorPtr_Type M_WS;
383  LifeV::UInt M_saveEvery;
384 
385  struct RESULT_CHANGED_EXCEPTION
386  {
387  public:
388  RESULT_CHANGED_EXCEPTION (LifeV::Real time)
389  {
390  std::cout << "Some modifications led to changes in the l2 norm of the solution at time" << time << std::endl;
391  }
392  };
393 };
394 
395 struct FSIChecker
396 {
397  FSIChecker ( GetPot const& _data_file ) :
398  data_file ( _data_file )
399  {}
400 
401  void operator() ()
402  {
403  std::shared_ptr<Problem> fsip;
404 
405  try
406  {
407  fsip = std::shared_ptr<Problem> ( new Problem ( data_file ) );
408 
409  fsip->run();
410  }
411  catch ( std::exception const& _ex )
412  {
413  std::cout << "caught exception : " << _ex.what() << "\n";
414  }
415 
416  }
417 
418  GetPot data_file;
419  LifeV::Vector disp;
420 };
421 
422 
423 
424 namespace LifeV
425 {
426 namespace
427 {
428 static bool regIF = (PRECFactory::instance().registerProduct ( "Ifpack", &createIfpack ) );
429 static bool regML = (PRECFactory::instance().registerProduct ( "ML", &createML ) );
430 }
431 }
432 
433 
434 int main (int argc, char** argv)
435 {
436 #ifdef HAVE_MPI
437  MPI_Init (&argc, &argv);
438 #else
439  std::cout << "% using serial Version" << std::endl;
440 #endif
441 
442  GetPot command_line (argc, argv);
443  const bool check = command_line.search (2, "-c", "--check");
444 
445  if (check)
446  {
447  GetPot data_fileGCE ("dataGCE");
448  FSIChecker _GCE_check ( data_fileGCE );
449  _GCE_check();
450 
451  GetPot data_fileCE ("dataCE");
452  FSIChecker _CE_check (data_fileCE);
453  _CE_check();
454 
455 #ifdef HAVE_MPI
456  MPI_Finalize();
457 #endif
458  return 0;
459  }
460  else
461  {
462  const std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
463  GetPot data_file (data_file_name);
464  FSIChecker _sp_check ( data_file );
465  _sp_check();
466  }
467 #ifdef HAVE_MPI
468  MPI_Finalize();
469 #endif
470 
471  return 0;
472 
473 }
474 
475 void Problem::restartFSI ( GetPot const& data_file)
476 {
477 
478  using namespace LifeV;
479 
480  typedef FSIOperator::mesh_Type mesh_Type;
481 
482  //Creating the pointer to the filter
483  std::string const importerType = data_file ( "importer/type", "ensight");
484  std::string const fluidName = data_file ( "importer/fluid/filename", "fluid");
485  std::string const solidName = data_file ( "importer/solid/filename", "solid");
486 
487  std::string const loadInitSol = data_file ( "importer/initSol", "00000");
488  std::string const loadInitSolFD = data_file ("importer/initSolFD", "-1");
489  std::string iterationString;
490 
491  M_Tstart = data_file ( "fluid/time_discretization/initialtime", 0.);
492 
493  std::cout << "The file for fluid is : " << fluidName << std::endl;
494  std::cout << "The file for solid is : " << solidName << std::endl;
495  std::cout << "The importerType is : " << importerType << std::endl;
496  std::cout << "The iteration is : " << loadInitSol << std::endl;
497  std::cout << "For the fluid disp is : " << loadInitSolFD << std::endl;
498  std::cout << "Starting time : " << M_Tstart << std::endl;
499 
500 #ifdef HAVE_HDF5
501  if (importerType.compare ("hdf5") == 0)
502  {
503  M_importerFluid.reset ( new hdf5Filter_Type ( data_file, fluidName) );
504  M_importerSolid.reset ( new hdf5Filter_Type ( data_file, solidName) );
505  }
506  else
507 #endif
508  {
509  if (importerType.compare ("none") == 0)
510  {
511  M_importerFluid.reset ( new ExporterEmpty<mesh_Type > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), "fluid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
512  M_importerSolid.reset ( new ExporterEmpty<mesh_Type > ( data_file, M_fsi->FSIOper()->dFESpace().mesh(), "solid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
513  }
514  else
515  {
516  M_importerFluid.reset ( new ensightFilter_Type ( data_file, fluidName) );
517  M_importerSolid.reset ( new ensightFilter_Type ( data_file, solidName) );
518  }
519  }
520 
521  M_importerFluid->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
522  M_importerSolid->setMeshProcId (M_fsi->FSIOper()->dFESpace().mesh(), M_fsi->FSIOper()->dFESpace().map().comm().MyPID() );
523 
524  //Each of the vectors of the stencils has the dimension of the big vector
525  //Performing a cycle of the size of the timeAdvance classes for each problem
526  //The size of TimeAdvanceClass depends on the order of the BDF (data file)
527 
528  //It should work just initializing the timeAdvance classes
529  //Three stencils are needed (Fluid-Structure-Geometric)
530 
531 
532  vectorPtr_Type fluidSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
533  vectorPtr_Type initFluid (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
534  vectorPtr_Type HarmonicSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
535  vectorPtr_Type structureSol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
536  vectorPtr_Type temporarySol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
537 
538  vectorPtr_Type vel (new vector_Type (M_fsi->FSIOper()->uFESpace().map(), M_importerFluid->mapType() ) );
539  vectorPtr_Type pressure (new vector_Type (M_fsi->FSIOper()->pFESpace().map(), M_importerFluid->mapType() ) );
540  vectorPtr_Type solidDisp (new vector_Type (M_fsi->FSIOper()->dFESpace().map(), M_importerSolid->mapType() ) );
541  vectorPtr_Type fluidDisp (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), M_importerFluid->mapType() ) );
542 
543  vectorPtr_Type firstFluidDisp (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), M_importerFluid->mapType() ) );
544 
545  UInt offset = dynamic_cast<LifeV::FSIMonolithic*> (M_fsi->FSIOper().get() )->offset();
546 
547  iterationString = loadInitSol;
548  std::cout << "Fluid size TimeAdvance:" << M_fsi->FSIOper()->fluidTimeAdvance()->size() << std::endl;
549  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->fluidTimeAdvance()->size(); iterInit++ )
550  {
551 
552  /*!
553  definition of the vector to fill with the initialization.
554  */
555  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
556  fluidSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
557  initFluid.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
558 
559 
560  LifeV::ExporterData<mesh_Type> initSolFluidVel (LifeV::ExporterData<mesh_Type>::VectorField, std::string ("f-velocity." + iterationString), M_fsi->FSIOper()->uFESpacePtr(), vel, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
561  LifeV::ExporterData<mesh_Type> initSolFluidPress (LifeV::ExporterData<mesh_Type>::ScalarField, std::string ("f-pressure." + iterationString), M_fsi->FSIOper()->pFESpacePtr(), pressure, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
562 
563  /*!load of the solutions*/
564  M_importerFluid->readVariable (initSolFluidVel);
565  M_importerFluid->readVariable (initSolFluidPress);
566 
567  fluidSol.reset ( new vector_Type (*pressure, Unique, Zero) );
568  initFluid->subset (*fluidSol, fluidSol->map(), UInt (0), (UInt) 3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() );
569  vector_Type tmpVec (*initFluid);
570  //temporarySol->subset(*fluidSol, fluidSol->map(), UInt(0), (UInt)3*M_fsi->FSIOper()->uFESpace().dof().numTotalDof());
571  fluidSol.reset ( new vector_Type (*vel, Unique, Zero) );
572  tmpVec = *fluidSol;
573  *initFluid += tmpVec;
574 
575  // std::string firstFl="firstFluid";
576  // initFluid->spy(firstFl);
577 
578  std::cout << "Norm of first Fluid sol: " << initFluid->norm2() << std::endl;
579  *temporarySol = *initFluid;
580 
581  // if(!iterInit)
582  // {
583  // *un += *temporarySol;
584  // }
585 
586  M_fluidStencil.push_back (temporarySol);
587  //Updating string name
588  int iterations = std::atoi (iterationString.c_str() );
589  iterations--;
590 
591  std::ostringstream iter;
592  iter.fill ( '0' );
593  iter << std::setw (5) << ( iterations );
594  iterationString = iter.str();
595 
596  }
597 
598  iterationString = loadInitSol;
599  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->solidTimeAdvance()->size(); iterInit++ )
600  {
601  /*!
602  definition of the vector to fill with the initialization.
603  */
604  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
605  structureSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
606 
607 
608  /*!Definition of the ExporterData, used to load the solution inside the previously defined vectors*/
609  LifeV::ExporterData<mesh_Type> initSolSolidDisp (LifeV::ExporterData<mesh_Type>::VectorField, "s-displacement." + iterationString, M_fsi->FSIOper()->dFESpacePtr(), solidDisp, UInt (0) /*offset*/, LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
610 
611  /*!load of the solutions*/
612  M_importerSolid->readVariable (initSolSolidDisp);
613 
614  structureSol->subset (*solidDisp, solidDisp->map(), (UInt) 0, offset);
615 
616  *temporarySol = *structureSol / (M_fsi->FSIOper()->solid().rescaleFactor() );
617 
618  M_solidStencil.push_back (temporarySol);
619 
620  //Updating the string name for the next iteration
621  UInt iterations = std::atoi (iterationString.c_str() );
622  iterations--;
623 
624  std::ostringstream iter;
625  iter.fill ( '0' );
626  iter << std::setw (5) << ( iterations );
627  iterationString = iter.str();
628 
629  }
630 
631  Int convectiveTerm = 0; //loadInitSol;
632  if (!M_data->dataFluid()->domainVelImplicit() )
633  {
634  convectiveTerm = 1;
635  }
636  HarmonicSol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique, Zero) );
637 
638  iterationString = loadInitSolFD;
639 
640  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->ALETimeAdvance()->size() + convectiveTerm + 1; iterInit++ )
641  {
642  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV:: Unique, Zero) );
643  /*!
644  definition of the vector to fill with the initialization.
645  */
646 
647  LifeV::ExporterData<mesh_Type> initSolFluidDisp (LifeV::ExporterData<mesh_Type>::VectorField, "f-displacement." + iterationString, M_fsi->FSIOper()->mmFESpacePtr(), fluidDisp, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
648 
649  /*!load of the solutions*/
650  M_importerFluid->readVariable (initSolFluidDisp);
651 
652  std::cout << "Reloaded Harmonic sol norm: " << fluidDisp->norm2() << std::endl;
653  if (iterInit == 0)
654  {
655  HarmonicSol->subset (*fluidDisp, fluidDisp->map(), (UInt) 0, dynamic_cast<LifeV::FSIMonolithicGI*> (M_fsi->FSIOper().get() )->mapWithoutMesh().map (Unique)->NumGlobalElements() );
656  if (!convectiveTerm)
657  {
658  *firstFluidDisp = *fluidDisp;
659  }
660  }
661  else if (convectiveTerm && iterInit == 1)
662  {
663  *firstFluidDisp = *fluidDisp;
664  }
665  else
666  {
667  M_ALEStencil.push_back (/*temporarySol*/fluidDisp);
668  }
669  //Updating the iteration String name
670  int iterations = std::atoi (iterationString.c_str() );
671  iterations--;
672 
673  std::ostringstream iter;
674  iter.fill ( '0' );
675  iter << std::setw (5) << ( iterations );
676  iterationString = iter.str();
677  }
678 
679  *M_fluidStencil[0] += *M_solidStencil[0];
680  *M_fluidStencil[0] += *HarmonicSol;
681  //this is going to be the global solutions returned by the method solution()
682 
683  M_fsi->initialize (M_fluidStencil, M_solidStencil, M_ALEStencil);
684 
685  if (!M_data->dataFluid()->domainVelImplicit() )
686  {
687  //The following is needed because (and if) of the extrapolation of the fluid domain velocity is used, i.e. M_domainVelImplicit
688  M_fsi->FSIOper()->ALETimeAdvance()->updateRHSFirstDerivative ( M_data->dataSolid()->dataTime()->timeStep() );
689  M_fsi->FSIOper()->ALETimeAdvance()->shiftRight (*firstFluidDisp);
690  }
691 
692  M_velAndPressure.reset ( new vector_Type ( M_fsi->FSIOper()->fluid().getMap(), M_importerFluid->mapType() ) );
693  M_velAndPressure->subset (*pressure, pressure->map(), UInt (0), (UInt) 3 * M_fsi->FSIOper()->uFESpace().dof().numTotalDof() );
694  *M_velAndPressure += *vel;
695 
696  M_fluidDisp.reset ( new vector_Type ( *fluidDisp, M_importerFluid->mapType() ) );
697 
698  M_solidDisp.reset ( new vector_Type ( *solidDisp, M_importerSolid->mapType() ) );
699 
700  M_data->dataFluid()->dataTime()->updateTime(), M_data->dataSolid()->dataTime()->updateTime();
701 }
702 
703 
704 void Problem::initializeStokes ( GetPot const& data_file)
705 {
706 
707  using namespace LifeV;
708 
709  typedef FSIOperator::mesh_Type mesh_Type;
710 
711 
712  //Creating the pointer to the filter
713  filterPtr_Type importer;
714  std::string const importerType = data_file ( "importer/type", "hdf5");
715  std::string const filename = data_file ( "importer/fluid/filename", "fluid");
716 
717  std::string const loadInitSol = data_file ( "importer/initSol", "00000");
718  std::string const loadInitSolFD = data_file ("importer/initSolFD", "-1");
719  std::string iterationString;
720 
721  std::cout << "The filename is : " << filename << std::endl;
722  std::cout << "The importerType is: " << importerType << std::endl;
723 
724 #ifdef HAVE_HDF5
725  if (importerType.compare ("hdf5") == 0)
726  {
727  importer.reset ( new hdf5Filter_Type ( data_file, filename) );
728  }
729  else
730 #endif
731  {
732  if (importerType.compare ("none") == 0)
733  {
734  importer.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( data_file, M_fsi->FSIOper()->uFESpace().mesh(), "fluid", M_fsi->FSIOper()->uFESpace().map().comm().MyPID() ) );
735  }
736  else
737  {
738  importer.reset ( new ensightFilter_Type ( data_file, filename) );
739  }
740  }
741 
742  importer->setMeshProcId (M_fsi->FSIOper()->uFESpace().mesh(), M_fsi->FSIOper()->uFESpace().map().comm().MyPID() );
743 
744 
745  std::vector<vectorPtr_Type> fluidStencil;;
746 
747  vectorPtr_Type temporarySol (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), LifeV::Unique) );
748 
749  //UInt offset=dynamic_cast<LifeV::FSIMonolithic*>(M_fsi->FSIOper().get())->offset();
750 
751  iterationString = loadInitSol;
752  for (UInt iterInit = 0; iterInit < M_fsi->FSIOper()->fluidTimeAdvance()->size(); ++iterInit )
753  {
754 
755  /*!
756  definition of the vector to fill with the initialization.
757  */
758  temporarySol.reset (new vector_Type (*M_fsi->FSIOper()->couplingVariableMap(), Unique) );
759 
760  vectorPtr_Type vel (new vector_Type (M_fsi->FSIOper()->uFESpace().map(), M_importerFluid->mapType() ) );
761  vectorPtr_Type pressure (new vector_Type (M_fsi->FSIOper()->pFESpace().map(), M_importerFluid->mapType() ) );
762 
763  *vel *= 0.0;
764  *pressure *= 0.0;
765 
766  LifeV::ExporterData<mesh_Type> initSolFluidVel (LifeV::ExporterData<mesh_Type>::VectorField, std::string ("f-velocity." + iterationString), M_fsi->FSIOper()->uFESpacePtr(), vel, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
767  LifeV::ExporterData<mesh_Type> initSolFluidPress (LifeV::ExporterData<mesh_Type>::ScalarField, std::string ("f-pressure." + iterationString), M_fsi->FSIOper()->pFESpacePtr(), pressure, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );
768 
769  /*!load of the solutions*/
770  importer->readVariable (initSolFluidVel);
771  importer->readVariable (initSolFluidPress);
772 
773  int iterations = std::atoi (iterationString.c_str() );
774  iterations--;
775 
776  std::ostringstream iter;
777  iter.fill ( '0' );
778  iter << std::setw (5) << ( iterations );
779  iterationString = iter.str();
780 
781  *temporarySol = *vel + *pressure;
782  fluidStencil.push_back (temporarySol);
783  }
784 
785  // M_fsi->initialize(fluidStencil);
786 
787 }
main
int main(int argc, char **argv)
Definition: lifev/level_set/tutorials/1_Volume_computation/main.cpp:172
GetPot::GetPot
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:507
anonymous_namespace{main.cpp}::regIF
static bool regIF
Definition: lifev/structure/testsuite/time_advance_ii/main.cpp:86
LifeV::MeshIO::anonymous_namespace{ParserGmsh.hpp}::elm_nodes_num
static const LifeV::UInt elm_nodes_num[]
Definition: ParserGmsh.hpp:57
FSIChecker::FSIChecker
FSIChecker(GetPot const &_data_file)
Definition: lifev/fsi/examples/challenge_VPH/main.cpp:397
GetPot::GetPot
GetPot(const char *FileName, const char *CommentStart=0x0, const char *CommentEnd=0x0, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:541
GetPot::search
bool search(unsigned No, const char *P,...)
Definition: GetPot.hpp:1217
GetPot::GetPot
GetPot(const STRING_VECTOR &FileNameList)
Definition: GetPot.hpp:645
GetPot::follow
const std::string follow(const char *Default, unsigned No, const char *Option,...)
Definition: GetPot.hpp:1510
anonymous_namespace{main.cpp}::regML
static bool regML
Definition: lifev/structure/testsuite/time_advance_ii/main.cpp:87