LifeV
MultiscaleModelFSI3D.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 Multiscale Model FSI3D
30  *
31  * @date 19-04-2010
32  * @author Paolo Crosetto <paolo.crosetto@epfl.ch>
33  * @author Cristiano Malossi <cristiano.malossi@epfl.ch>
34  *
35  * @maintainer Cristiano Malossi <cristiano.malossi@epfl.ch>
36  */
37 
38 #include <lifev/multiscale/models/MultiscaleModelFSI3D.hpp>
39 
40 namespace LifeV
41 {
42 namespace Multiscale
43 {
44 
45 // ===================================================
46 // Constructors & Destructor
47 // ===================================================
48 MultiscaleModelFSI3D::MultiscaleModelFSI3D() :
49  multiscaleModel_Type (),
50  MultiscaleInterface (),
51  M_FSIoperator (),
52  M_data ( new data_Type() ),
53  M_exporterFluid (),
54  M_exporterSolid (),
55  M_importerFluid (),
56  M_importerSolid (),
57 #ifndef FSI_WITH_EXTERNALPRESSURE
58  M_boundaryStressFunctions (),
59  M_externalPressureScalar (),
60 #endif
61 #ifndef FSI_WITHOUT_VELOCITYPROFILE
62  M_boundaryFlowRateFunctions (),
63  M_boundaryFlowRateType (),
64 #endif
65 #ifdef FSI_WITH_BOUNDARYAREA
66  M_boundaryAreaFunctions (),
67  M_boundaryFlagsArea (),
68 #endif
69 #ifdef FSI_WITH_STRESSOUTPUT
70  M_wallTensionEstimator ( new wallTensionEstimator_Type() ),
71 #endif
72  M_fluidVelocity (),
73  M_fluidPressure (),
74  M_fluidDisplacement (),
75  M_solidVelocity (),
76  M_solidDisplacement (),
77 #ifdef FSI_WITH_STRESSOUTPUT
78  M_solidStressXX (),
79  M_solidStressXY (),
80  M_solidStressXZ (),
81  M_solidStressYY (),
82  M_solidStressYZ (),
83  M_solidStressZZ (),
84  M_solidStressVonMises (),
85 #endif
86  M_stateVariable (),
87  M_nonLinearRichardsonIteration ( 0 ),
88  M_fluidBC ( new bcInterface_Type() ),
89  M_solidBC ( new bcInterface_Type() ),
90  M_harmonicExtensionBC ( new bcInterface_Type() ),
91  M_linearFluidBC (),
92  M_linearSolidBC (),
93  M_linearRHS (),
94  M_linearSolution ()
95 {
96 
97 #ifdef HAVE_LIFEV_DEBUG
98  debugStream ( 8140 ) << "MultiscaleModelFSI3D::MultiscaleModelFSI3D() \n";
99 #endif
100 
101  M_type = FSI3D;
102 
103  FSIOperator_Type::factory_Type::instance().registerProduct ( "monolithicGE", &createFSIMonolithicGE );
104  FSIOperator_Type::factory_Type::instance().registerProduct ( "monolithicGI", &createFSIMonolithicGI );
105 
106  BlockPrecFactory::instance().registerProduct ("AdditiveSchwarz", &MonolithicBlockMatrix::createAdditiveSchwarz) ;
107  BlockPrecFactory::instance().registerProduct ("AdditiveSchwarzRN", &MonolithicBlockMatrixRN::createAdditiveSchwarzRN ) ;
108  BlockPrecFactory::instance().registerProduct ("ComposedDN", &MonolithicBlockComposedDN::createComposedDN ) ;
109  BlockPrecFactory::instance().registerProduct ("ComposedDN2", &MonolithicBlockComposedDN::createComposedDN2 );
110  BlockPrecFactory::instance().registerProduct ("ComposedDNND", &MonolithicBlockComposedDNND::createComposedDNND);
111 
112  MonolithicBlockMatrix::Factory_Type::instance().registerProduct ("AdditiveSchwarz", &MonolithicBlockMatrix::createAdditiveSchwarz ) ;
113  MonolithicBlockMatrix::Factory_Type::instance().registerProduct ("AdditiveSchwarzRN", &MonolithicBlockMatrixRN::createAdditiveSchwarzRN ) ;
114 
115  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "linearVenantKirchhoff", &FSIOperator::createVenantKirchhoffLinear );
116  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "exponential", &FSIOperator::createExponentialMaterialNonLinear );
117  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "neoHookean", &FSIOperator::createNeoHookeanMaterialNonLinear );
118  FSIOperator::solid_Type::material_Type::isotropicLaw_Type::StructureIsotropicMaterialFactory::instance().registerProduct ( "nonLinearVenantKirchhoff", &FSIOperator::createVenantKirchhoffNonLinear );
119 
120 }
121 
122 // ===================================================
123 // MultiscaleModel Methods
124 // ===================================================
125 void
126 MultiscaleModelFSI3D::setupData ( const std::string& fileName )
127 {
128 
129 #ifdef HAVE_LIFEV_DEBUG
130  debugStream ( 8140 ) << "MultiscaleModelFSI3D::setupData( fileName ) \n";
131 #endif
132 
133  multiscaleModel_Type::setupData ( fileName );
134 
135  GetPot dataFile ( fileName );
136 
137  // Load data
138  M_data->setup ( dataFile );
139  if ( M_globalData.get() )
140  {
141  setupGlobalData ( fileName );
142  }
143 
144  // Create FSI
145  M_FSIoperator = FSIOperatorPtr_Type ( FSIOperator_Type::factory_Type::instance().createObject ( M_data->method() ) );
146 
147  // Setup Communicator
148  setupCommunicator();
149 
150  // Set data
151  M_FSIoperator->setData ( M_data );
152  M_FSIoperator->setDataFile ( dataFile );
153 
154  // Setup Boundary Conditions from file
155  setupBC ( fileName );
156 
157  // Setup exporter
158  if ( M_FSIoperator->isFluid() )
159  {
160  setupExporter ( M_exporterFluid, dataFile );
161  setupImporter ( M_importerFluid, dataFile );
162  }
163  if ( M_FSIoperator->isSolid() )
164  {
165  setupExporter ( M_exporterSolid, dataFile, "_Solid" );
166  setupImporter ( M_importerSolid, dataFile, "_Solid" );
167  }
168 
169 #ifndef FSI_WITHOUT_VELOCITYPROFILE
170  std::map< std::string, FSI3DBoundaryFlowRate_Type > boundaryFlowRateMap;
171  boundaryFlowRateMap["Weak"] = Weak;
172  boundaryFlowRateMap["Semiweak"] = Semiweak;
173  boundaryFlowRateMap["Strong"] = Strong;
174 
175  M_boundaryFlowRateType.reserve ( M_boundaryFlags.size() );
176  for ( UInt j ( 0 ); j < M_boundaryFlags.size(); ++j )
177  {
178  M_boundaryFlowRateType[j] = boundaryFlowRateMap[dataFile ( "Multiscale/flowRateCouplingType", "Weak", j )];
179  }
180 #endif
181 
182 #ifdef FSI_WITH_BOUNDARYAREA
183  M_boundaryFlagsArea.reserve ( M_boundaryFlags.size() );
184  M_boundaryFlagsAreaPerturbed.reserve ( M_boundaryFlags.size() );
185  for ( UInt j ( 0 ); j < M_boundaryFlags.size(); ++j )
186  {
187  M_boundaryFlagsArea.push_back ( dataFile ( "Multiscale/couplingAreaFlags", 0, j ) );
188  M_boundaryFlagsAreaPerturbed.push_back ( false );
189  }
190 #endif
191 }
192 
193 void
194 MultiscaleModelFSI3D::setupModel()
195 {
196 
197 #ifdef HAVE_LIFEV_DEBUG
198  debugStream ( 8140 ) << "MultiscaleModelFSI3D::setupProblem() \n";
199 #endif
200 
201  // Mesh transformation (before partitioning, ideally should be done after for scalability)
202  //M_FSIoperator->fluidMesh().transformMesh( M_geometryScale, M_geometryRotate, M_geometryTranslate );
203  M_FSIoperator->solidMesh().meshTransformer().transformMesh ( M_geometryScale, M_geometryRotate, M_geometryTranslate );
204 
205  // Mesh partitioning
206  M_FSIoperator->partitionMeshes();
207 
208  // Mesh transformation (after partitioning - not working for solid)
209  M_FSIoperator->fluidLocalMesh().meshTransformer().transformMesh ( M_geometryScale, M_geometryRotate, M_geometryTranslate );
210  //M_FSIoperator->solidLocalMesh().meshTransformer()->transformMesh( M_geometryScale, M_geometryRotate, M_geometryTranslate );
211 
212  // Setup FEspace & DOF
213  M_FSIoperator->setupFEspace();
214  M_FSIoperator->setupDOF();
215 
216  // Setup FSI Interface Boundary Conditions (by giving the operator to BCInterface)
217  M_fluidBC->setPhysicalSolver ( M_FSIoperator );
218  M_solidBC->setPhysicalSolver ( M_FSIoperator );
219  M_harmonicExtensionBC->setPhysicalSolver ( M_FSIoperator );
220 
221  // Setup Fluid & Solid solver
222  M_FSIoperator->setupFluidSolid ( M_FSIoperator->imposedFluxes() );
223  M_FSIoperator->setupSystem();
224 
225 #ifdef FSI_WITH_BOUNDARYAREA
226  // Setup the boundary area functions
227  for ( boundaryAreaFunctionsContainerIterator_Type i = M_boundaryAreaFunctions.begin(); i < M_boundaryAreaFunctions.end(); ++i )
228  {
229  ( *i )->setup();
230  }
231 #endif
232 
233 #ifdef FSI_WITH_STRESSOUTPUT
234  M_wallTensionEstimator->setup ( M_FSIoperator->solid().data(), M_FSIoperator->dFESpacePtr(), M_FSIoperator->dFESpaceETPtr(), M_comm, static_cast<UInt> ( 0 ) );
235 #endif
236 
237  // Setup Exporters
238  if ( M_FSIoperator->isFluid() )
239  {
240  setExporterFluid ( M_exporterFluid );
241  }
242 
243  if ( M_FSIoperator->isSolid() )
244  {
245  setExporterSolid ( M_exporterSolid );
246  }
247 
248  //Setup solution
249  initializeSolution();
250 
251  //Setup linear model
252  setupLinearModel();
253 }
254 
255 void
256 MultiscaleModelFSI3D::buildModel()
257 {
258 
259 #ifdef HAVE_LIFEV_DEBUG
260  debugStream ( 8140 ) << "MultiscaleModelFSI3D::buildModel() \n";
261 #endif
262 
263  // Display data
264  // if ( M_comm->MyPID() == 0 )
265  // M_data->showMe();
266 
267  M_FSIoperator->buildSystem();
268  M_FSIoperator->updateSystem();
269 
270  // Update BCInterface solver variables
271  updateBC();
272 }
273 
274 void
275 MultiscaleModelFSI3D::updateModel()
276 {
277 
278 #ifdef HAVE_LIFEV_DEBUG
279  debugStream ( 8140 ) << "MultiscaleModelFSI3D::updateModel() \n";
280 #endif
281 
282  // Update System
283  M_FSIoperator->updateSystem();
284 
285  // Update BCInterface solver variables
286  updateBC();
287 
288  // TODO This is a workaround of Paolo Crosetto to make it possible to perform subiterations
289  // in the multiscale when using 3D FSI models. In the future this should be replaced with
290  // a more proper implementation.
291  M_FSIoperator->precPtrView()->setRecompute ( 1, true );
292  M_nonLinearRichardsonIteration = 0;
293 }
294 
295 void
296 MultiscaleModelFSI3D::solveModel()
297 {
298 
299 #ifdef HAVE_LIFEV_DEBUG
300  debugStream ( 8140 ) << "MultiscaleModelFSI3D::solveModel() \n";
301 #endif
302 
303  displayModelStatus ( "Solve" );
304 
305  // TODO This is a workaround of Paolo Crosetto to make it possible to perform subiterations
306  // in the multiscale when using 3D FSI models. In the future this should be replaced with
307  // a more proper implementation.
308  if ( M_nonLinearRichardsonIteration != 0 )
309  {
310  M_FSIoperator->resetRHS();
311  M_FSIoperator->updateRHS();
312  M_FSIoperator->applyBoundaryConditions();
313  }
314 
315  // Non-linear Richardson solver
316  UInt maxSubIterationNumber = M_data->maxSubIterationNumber();
317  M_FSIoperator->extrapolation ( *M_stateVariable );
318 
319  NonLinearRichardson ( *M_stateVariable, *M_FSIoperator,
320  M_data->absoluteTolerance(), M_data->relativeTolerance(),
321  maxSubIterationNumber, M_data->errorTolerance(),
322  M_data->NonLinearLineSearch(),
323  M_nonLinearRichardsonIteration,
324  1
325  );
326 
327  // TODO This is a workaround of Paolo Crosetto to make it possible to perform subiterations
328  // in the multiscale when using 3D FSI models. In the future this should be replaced with
329  // a more proper implementation.
330  M_FSIoperator->precPtrView()->setRecompute ( 1, false );
331  M_nonLinearRichardsonIteration = 1;
332 
333 #ifdef FSI_WITH_BOUNDARYAREA
334  for ( UInt j ( 0 ); j < M_boundaryFlagsAreaPerturbed.size(); ++j )
335  {
336  M_boundaryFlagsAreaPerturbed[j] = false;
337  }
338 #endif
339 }
340 
341 void
342 MultiscaleModelFSI3D::updateSolution()
343 {
344 
345 #ifdef HAVE_LIFEV_DEBUG
346  debugStream ( 8140 ) << "MultiscaleModelFSI3D::updateSolution() \n";
347 #endif
348 
349  M_FSIoperator->updateSolution ( *M_stateVariable );
350 }
351 
352 void
353 MultiscaleModelFSI3D::saveSolution()
354 {
355 
356 #ifdef HAVE_LIFEV_DEBUG
357  debugStream ( 8140 ) << "MultiscaleModelFSI3D::saveSolution() \n";
358 #endif
359 
360  exportFluidSolution();
361  if ( M_FSIoperator->isFluid() )
362  {
363  M_exporterFluid->postProcess ( M_data->dataFluid()->dataTime()->time() );
364  }
365 
366  exportSolidSolution();
367  if ( M_FSIoperator->isSolid() )
368  {
369  M_exporterSolid->postProcess ( M_data->dataSolid()->dataTime()->time() );
370  }
371 
372 #ifdef HAVE_HDF5
373  if ( M_data->dataFluid()->dataTime()->isLastTimeStep() )
374  {
375  if ( M_FSIoperator->isFluid() )
376  {
377  M_exporterFluid->closeFile();
378  }
379  if ( M_FSIoperator->isSolid() )
380  {
381  M_exporterSolid->closeFile();
382  }
383  }
384 #endif
385 
386 }
387 
388 void
389 MultiscaleModelFSI3D::showMe()
390 {
391  if ( M_comm->MyPID() == 0 )
392  {
393  multiscaleModel_Type::showMe();
394 
395  std::cout << "FSI method = " << M_data->method() << std::endl << std::endl;
396 
397  std::cout << "Velocity FE order = " << M_FSIoperator->dataFluid()->uOrder() << std::endl
398  << "Pressure FE order = " << M_FSIoperator->dataFluid()->pOrder() << std::endl
399  << "Structure FE order = " << M_FSIoperator->dataSolid()->order() << std::endl << std::endl;
400 
401  std::cout << "Velocity DOF = " << 3 * M_FSIoperator->uFESpace().dof().numTotalDof() << std::endl
402  << "Pressure DOF = " << M_FSIoperator->pFESpace().dof().numTotalDof() << std::endl
403  << "Harmonic ext. DOF = " << M_FSIoperator->mmFESpace().dof().numTotalDof() << std::endl
404  << "Structure DOF = " << M_FSIoperator->dFESpace().dof().numTotalDof() << std::endl << std::endl;
405 
406  std::cout << "Fluid mesh maxH = " << MeshUtility::MeshStatistics::computeSize ( * ( M_FSIoperator->uFESpace().mesh() ) ).maxH << std::endl
407  << "Fluid mesh meanH = " << MeshUtility::MeshStatistics::computeSize ( * ( M_FSIoperator->uFESpace().mesh() ) ).meanH << std::endl
408  << "Solid mesh maxH = " << MeshUtility::MeshStatistics::computeSize ( * ( M_FSIoperator->dFESpace().mesh() ) ).maxH << std::endl
409  << "Solid mesh meanH = " << MeshUtility::MeshStatistics::computeSize ( * ( M_FSIoperator->dFESpace().mesh() ) ).meanH << std::endl << std::endl;
410  }
411 }
412 
413 Real
414 MultiscaleModelFSI3D::checkSolution() const
415 {
416  return M_stateVariable->norm2();
417 }
418 
419 // ===================================================
420 // MultiscaleInterface Methods
421 // ===================================================
422 void
423 MultiscaleModelFSI3D::imposeBoundaryFlowRate ( const multiscaleID_Type& boundaryID, const function_Type& function )
424 {
425  BCFunctionBase base;
426 
427 #ifndef FSI_WITHOUT_VELOCITYPROFILE
428  switch ( M_boundaryFlowRateType[boundaryID] )
429  {
430  case Strong:
431  {
432  boundaryFlowRateFunctionPtr_Type boundaryFlowRateFunction ( new boundaryFlowRateFunction_Type() );
433  boundaryFlowRateFunction->setModel ( this );
434  boundaryFlowRateFunction->setFluidFlag ( boundaryFlag ( boundaryID ) );
435  boundaryFlowRateFunction->setFunction ( function);
436  boundaryFlowRateFunction->setBoundaryFlowRateType ( Strong );
437 
438  M_boundaryFlowRateFunctions.push_back ( boundaryFlowRateFunction );
439 
440  base.setFunction ( std::bind ( &FSI3DBoundaryFlowRateFunction::function, M_boundaryFlowRateFunctions.back(), std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
441  M_fluidBC->handler()->addBC ( "CouplingFlowRate_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), boundaryFlag ( boundaryID ), EssentialVertices, Full, base, 3 );
442 
443  break;
444  }
445  case Semiweak:
446  {
447  // TODO: implementation is still required for the switching of the boundary condition.
448  // Possible idea: we can use the same function and change just the type of BC inside the BCHandler from Flux to Essential.
449  boundaryFlowRateFunctionPtr_Type boundaryFlowRateFunction ( new boundaryFlowRateFunction_Type() );
450  boundaryFlowRateFunction->setModel ( this );
451  boundaryFlowRateFunction->setFluidFlag ( boundaryFlag ( boundaryID ) );
452  boundaryFlowRateFunction->setFunction ( function);
453  boundaryFlowRateFunction->setBoundaryFlowRateType ( Semiweak );
454 
455  M_boundaryFlowRateFunctions.push_back ( boundaryFlowRateFunction );
456 
457  base.setFunction ( std::bind ( &FSI3DBoundaryFlowRateFunction::function, M_boundaryFlowRateFunctions.back(), std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
458  M_fluidBC->handler()->addBC ( "CouplingFlowRate_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), boundaryFlag ( boundaryID ), Flux, Full, base, 3 );
459 
460  break;
461  }
462  case Weak:
463  default:
464 
465  base.setFunction ( function );
466  M_fluidBC->handler()->addBC ( "CouplingFlowRate_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), boundaryFlag ( boundaryID ), Flux, Full, base, 3 );
467 
468  break;
469  }
470 #else
471  base.setFunction ( function );
472  M_fluidBC->handler()->addBC ( "CouplingFlowRate_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), boundaryFlag ( boundaryID ), Flux, Full, base, 3 );
473 #endif
474 }
475 
476 void
477 MultiscaleModelFSI3D::imposeBoundaryMeanNormalStress ( const multiscaleID_Type& boundaryID, const function_Type& function )
478 {
479  BCFunctionBase base;
480 #ifdef FSI_WITH_EXTERNALPRESSURE
481  base.setFunction ( function );
482 #else
483  boundaryStressFunctionPtr_Type boundaryStressFunction ( new boundaryStressFunction_Type() );
484  boundaryStressFunction->setDelta ( M_data->dataSolid()->externalPressure() );
485  boundaryStressFunction->setFunction ( function);
486 
487  M_boundaryStressFunctions.push_back ( boundaryStressFunction );
488 
489  base.setFunction ( std::bind ( &FSI3DBoundaryStressFunction::function, M_boundaryStressFunctions.back(), std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
490 #endif
491  M_fluidBC->handler()->addBC ( "BoundaryStress_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), boundaryFlag ( boundaryID ), Natural, Normal, base );
492 }
493 
494 void
495 MultiscaleModelFSI3D::imposeBoundaryArea ( const multiscaleID_Type& boundaryID, const function_Type& function )
496 {
497 #ifdef FSI_WITH_BOUNDARYAREA
498  boundaryAreaFunctionPtr_Type boundaryAreaFunction ( new boundaryAreaFunction_Type() );
499  boundaryAreaFunction->setModel ( this );
500  boundaryAreaFunction->setFluidFlag ( boundaryFlag ( boundaryID ) );
501  boundaryAreaFunction->setFunction ( function);
502 
503  M_boundaryAreaFunctions.push_back ( boundaryAreaFunction );
504 
505  BCFunctionBase base;
506  base.setFunction ( std::bind ( &FSI3DBoundaryAreaFunction::function, M_boundaryAreaFunctions.back(), std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
507  M_solidBC->handler()->addBC ( "BoundaryArea_Model_" + number2string ( M_ID ) + "_BoundaryID_" + number2string ( boundaryID ), M_boundaryFlagsArea[boundaryID], EssentialEdges, Full, base, 3 );
508 #endif
509 }
510 
511 Real
512 MultiscaleModelFSI3D::boundaryMeanNormalStress ( const multiscaleID_Type& boundaryID ) const
513 {
514 #ifdef FSI_WITH_EXTERNALPRESSURE
515  return M_FSIoperator->fluid().meanNormalStress ( boundaryFlag ( boundaryID ), *M_fluidBC->handler(), *M_stateVariable );
516 #else
517  return M_FSIoperator->fluid().meanNormalStress ( boundaryFlag ( boundaryID ), *M_fluidBC->handler(), *M_stateVariable ) - M_externalPressureScalar;
518 #endif
519 }
520 
521 Real
522 MultiscaleModelFSI3D::boundaryMeanTotalNormalStress ( const multiscaleID_Type& boundaryID ) const
523 {
524 #ifdef FSI_WITH_EXTERNALPRESSURE
525  return M_FSIoperator->fluid().meanTotalNormalStress ( boundaryFlag ( boundaryID ), *M_fluidBC->handler(), *M_stateVariable );
526 #else
527  return M_FSIoperator->fluid().meanTotalNormalStress ( boundaryFlag ( boundaryID ), *M_fluidBC->handler(), *M_stateVariable ) - M_externalPressureScalar;
528 #endif
529 }
530 
531 Real
532 MultiscaleModelFSI3D::boundaryDeltaFlowRate ( const multiscaleID_Type& boundaryID, bool& solveLinearSystem )
533 {
534  solveLinearModel ( solveLinearSystem );
535 
536  return M_FSIoperator->fluid().flux ( boundaryFlag ( boundaryID ), *M_linearSolution );
537 }
538 
539 Real
540 MultiscaleModelFSI3D::boundaryDeltaMeanNormalStress ( const multiscaleID_Type& boundaryID, bool& solveLinearSystem )
541 {
542  solveLinearModel ( solveLinearSystem );
543 
544  return M_FSIoperator->fluid().linearMeanNormalStress ( boundaryFlag ( boundaryID ), *M_linearFluidBC, *M_linearSolution );
545 }
546 
547 Real
548 MultiscaleModelFSI3D::boundaryDeltaMeanTotalNormalStress ( const multiscaleID_Type& boundaryID, bool& solveLinearSystem )
549 {
550  solveLinearModel ( solveLinearSystem );
551 
552  return M_FSIoperator->fluid().linearMeanTotalNormalStress ( boundaryFlag ( boundaryID ), *M_linearFluidBC, *M_stateVariable, *M_linearSolution );
553 }
554 
555 // ===================================================
556 // Get Methods
557 // ===================================================
558 Real
559 MultiscaleModelFSI3D::boundaryPressure ( const multiscaleID_Type& boundaryID ) const
560 {
561 #ifdef FSI_WITH_EXTERNALPRESSURE
562  return M_FSIoperator->fluid().pressure ( boundaryFlag ( boundaryID ), *M_stateVariable );
563 #else
564  return M_FSIoperator->fluid().pressure ( boundaryFlag ( boundaryID ), *M_stateVariable ) + M_externalPressureScalar;
565 #endif
566 }
567 
568 Real
569 MultiscaleModelFSI3D::boundaryTotalPressure ( const multiscaleID_Type& boundaryID ) const
570 {
571 #ifdef FSI_WITH_EXTERNALPRESSURE
572  return M_FSIoperator->fluid().pressure ( boundaryFlag ( boundaryID ), *M_stateVariable )
573  + M_FSIoperator->fluid().kineticNormalStress ( boundaryFlag ( boundaryID ), *M_stateVariable );
574 #else
575  return M_FSIoperator->fluid().pressure ( boundaryFlag ( boundaryID ), *M_stateVariable )
576  + M_FSIoperator->fluid().kineticNormalStress ( boundaryFlag ( boundaryID ), *M_stateVariable ) + M_externalPressureScalar;
577 #endif
578 }
579 
580 Real
581 MultiscaleModelFSI3D::externalPressure() const
582 {
583 #ifndef FSI_WITH_EXTERNALPRESSURE
584  return M_externalPressureScalar;
585 #else
586  return M_data->dataSolid()->externalPressure();
587 #endif
588 }
589 
590 // ===================================================
591 // Private Methods
592 // ===================================================
593 void
594 MultiscaleModelFSI3D::setupGlobalData ( const std::string& fileName )
595 {
596  GetPot dataFile ( fileName );
597 
598  // Global data time
599  M_data->dataFluid()->setTimeData ( M_globalData->dataTime() );
600  M_data->dataSolid()->setTimeData ( M_globalData->dataTime() );
601  M_data->setTimeDataALE ( M_globalData->dataTime() );
602 
603  // Fluid global physical quantities
604  if ( !dataFile.checkVariable ( "fluid/physics/density" ) )
605  {
606  M_data->dataFluid()->setDensity ( M_globalData->fluidDensity() );
607  }
608  if ( !dataFile.checkVariable ( "fluid/physics/viscosity" ) )
609  {
610  M_data->dataFluid()->setViscosity ( M_globalData->fluidViscosity() );
611  }
612 
613  // Solid global physical quantities
614  if ( !dataFile.checkVariable ( "solid/physics/density" ) )
615  {
616  M_data->dataSolid()->setDensity ( M_globalData->solidDensity() );
617  }
618  if ( !dataFile.checkVariable ( "solid/physics/externalPressure" ) )
619  {
620  M_data->dataSolid()->setExternalPressure ( M_globalData->solidExternalPressure() );
621  }
622 
623  std::vector< UInt > materialFlags;
624  if ( !dataFile.checkVariable ( "solid/physics/material_flag" ) )
625  {
626  materialFlags.push_back ( 1 );
627  }
628  else
629  for ( UInt i ( 0 ) ; i < dataFile.vector_variable_size ( "solid/physics/material_flag" ) ; ++i )
630  {
631  materialFlags.push_back ( dataFile ( "solid/physics/material_flag", 1, i ) );
632  }
633 
634  for ( std::vector< UInt >::const_iterator i = materialFlags.begin(); i != materialFlags.end() ; ++i )
635  {
636  if ( !dataFile.checkVariable ( "solid/physics/poisson" ) )
637  {
638  M_data->dataSolid()->setPoisson ( M_globalData->solidPoissonCoefficient(), *i );
639  }
640  if ( !dataFile.checkVariable ( "solid/physics/young" ) )
641  {
642  M_data->dataSolid()->setYoung ( M_globalData->solidYoungModulus(), *i );
643  }
644  }
645 }
646 
647 void
648 MultiscaleModelFSI3D::initializeSolution()
649 {
650 
651 #ifdef HAVE_LIFEV_DEBUG
652  debugStream ( 8140 ) << "MultiscaleModelFSI3D::initializeSolution() \n";
653 #endif
654 
655 #ifndef FSI_WITH_EXTERNALPRESSURE
656  // Initialize the external pressure scalar
657  M_externalPressureScalar = M_data->dataSolid()->externalPressure();
658  M_data->dataSolid()->setExternalPressure ( 0.0 );
659 #endif
660 
661  if ( multiscaleProblemStep > 0 )
662  {
663  M_importerFluid->setMeshProcId ( M_FSIoperator->uFESpace().mesh(), M_FSIoperator->uFESpace().map().comm().MyPID() );
664  M_importerSolid->setMeshProcId ( M_FSIoperator->dFESpace().mesh(), M_FSIoperator->dFESpace().map().comm().MyPID() );
665 
666  M_importerFluid->addVariable ( IOData_Type::VectorField, "Velocity (fluid)", M_FSIoperator->uFESpacePtr(), M_fluidVelocity, static_cast <UInt> (0) );
667  M_importerFluid->addVariable ( IOData_Type::ScalarField, "Pressure (fluid)", M_FSIoperator->pFESpacePtr(), M_fluidPressure, static_cast <UInt> (0) );
668  M_importerFluid->addVariable ( IOData_Type::VectorField, "Displacement (fluid)", M_FSIoperator->mmFESpacePtr(), M_fluidDisplacement, static_cast <UInt> (0) );
669  M_importerSolid->addVariable ( IOData_Type::VectorField, "Displacement (solid)", M_FSIoperator->dFESpacePtr(), M_solidDisplacement, static_cast <UInt> (0) );
670 
671  // Read fluid velocity and pressure + solid displacement
672  for ( UInt i (0); i < M_FSIoperator->fluidTimeAdvance()->size() ; ++i )
673  {
674  UInt iterationImported (0);
675  iterationImported = M_importerFluid->importFromTime ( M_data->dataFluid()->dataTime()->initialTime() - i * M_data->dataFluid()->dataTime()->timeStep() );
676  iterationImported = M_importerSolid->importFromTime ( M_data->dataSolid()->dataTime()->initialTime() - i * M_data->dataSolid()->dataTime()->timeStep() );
677  if ( i == 0 )
678  {
679  M_exporterFluid->setTimeIndex ( iterationImported + 1 );
680  M_exporterSolid->setTimeIndex ( iterationImported + 1 );
681  }
682 
683 #ifndef FSI_WITH_EXTERNALPRESSURE
684  // Remove external pressure
685  *M_fluidPressure -= M_externalPressureScalar;
686 #endif
687 
688  M_FSIoperator->setFluidVectorInStencil ( M_fluidVelocity, M_fluidPressure, i );
689  M_FSIoperator->setSolidVectorInStencil ( M_solidDisplacement, i );
690  }
691 
692  // Read last solid displacement
693  M_importerSolid->importFromTime ( M_data->dataSolid()->dataTime()->initialTime() - M_FSIoperator->fluidTimeAdvance()->size() * M_data->dataSolid()->dataTime()->timeStep() );
694  M_FSIoperator->setSolidVectorInStencil ( M_solidDisplacement, M_FSIoperator->fluidTimeAdvance()->size() );
695 
696  // Read fluid displacement
697  for ( UInt i (0); i < M_FSIoperator->ALETimeAdvance()->size() ; ++i )
698  {
699  M_importerFluid->importFromTime ( M_data->dataFluid()->dataTime()->initialTime() - (i + 1) * M_data->dataFluid()->dataTime()->timeStep() );
700 
701  M_FSIoperator->setALEVectorInStencil ( M_fluidDisplacement, i, false );
702  }
703 
704  // Read first fluid displacement
705  M_importerFluid->importFromTime ( M_data->dataSolid()->dataTime()->initialTime() );
706 
707  //This is ugly but it's the only way I have figured out at the moment
708  if ( M_data->method().compare ("monolithicGI") == 0 )
709  {
710  //Don't be scared by the ten. The goal of 10 is just to make the first if fail
711  M_FSIoperator->setALEVectorInStencil ( M_fluidDisplacement, 10, true );
712  }
713 
714  //Setting the vector in the stencil
715  M_FSIoperator->ALETimeAdvance()->shiftRight ( *M_fluidDisplacement );
716 
717  // Finalize restart
718  M_FSIoperator->finalizeRestart();
719 
720  // Re-initialize fluid and velocity vectors
721  exportFluidSolution();
722  exportSolidSolution();
723 
724 #ifdef HAVE_HDF5
725  if ( M_FSIoperator->isFluid() )
726  {
727  M_importerFluid->closeFile();
728  }
729  if ( M_FSIoperator->isSolid() )
730  {
731  M_importerSolid->closeFile();
732  }
733 #endif
734 
735  }
736  else
737  {
738  // Initialize containers to zero
739  *M_fluidVelocity *= 0.0;
740  *M_fluidPressure = M_data->dataSolid()->externalPressure();
741  *M_solidDisplacement *= 0.0;
742  *M_fluidDisplacement *= 0.0;
743 
744  for ( UInt i (0); i < M_FSIoperator->fluidTimeAdvance()->size() ; ++i )
745  {
746  M_FSIoperator->setFluidVectorInStencil ( M_fluidVelocity, M_fluidPressure, i );
747  M_FSIoperator->setSolidVectorInStencil ( M_solidDisplacement, i );
748  }
749 
750  M_FSIoperator->setSolidVectorInStencil ( M_solidDisplacement, M_FSIoperator->fluidTimeAdvance()->size() );
751 
752  for ( UInt i (0); i < M_FSIoperator->ALETimeAdvance()->size() ; ++i )
753  {
754  M_FSIoperator->setALEVectorInStencil ( M_fluidDisplacement, i, false );
755  }
756 
757  //This is ugly but it's the only way I have figured out at the moment
758  if ( M_data->method().compare ("monolithicGI") == 0 )
759  {
760  //Don't be scared by the ten. The goal of 10 is just to make the first if fail
761  M_FSIoperator->setALEVectorInStencil ( M_fluidDisplacement, 10, true );
762  }
763 
764  //Setting the vector in the stencil
765  M_FSIoperator->ALETimeAdvance()->shiftRight ( *M_fluidDisplacement );
766 
767  // Finalize restart
768  M_FSIoperator->finalizeRestart();
769  }
770 
771  // Initialize state variables
772  M_stateVariable.reset ( new vector_Type ( M_FSIoperator->fluidTimeAdvance()->singleElement (0) ) );
773 }
774 
775 void
776 MultiscaleModelFSI3D::setupCommunicator()
777 {
778  M_FSIoperator->setFluid ( true );
779  M_FSIoperator->setSolid ( true );
780 
781  M_FSIoperator->setFluidLeader ( 0 );
782  M_FSIoperator->setSolidLeader ( 0 );
783 
784  M_FSIoperator->setComm ( M_comm, M_comm );
785 }
786 
787 void
788 MultiscaleModelFSI3D::setupBC ( const std::string& fileName )
789 {
790  if ( M_FSIoperator->isFluid() )
791  {
792  M_fluidBC->createHandler();
793  M_fluidBC->fillHandler ( fileName, "fluid" );
794 
795  M_FSIoperator->setFluidBC ( M_fluidBC->handler() );
796 
797  M_harmonicExtensionBC->createHandler();
798  M_harmonicExtensionBC->fillHandler ( fileName, "mesh_motion" );
799 
800  M_FSIoperator->setHarmonicExtensionBC ( M_harmonicExtensionBC->handler() );
801  }
802 
803  if ( M_FSIoperator->isSolid() )
804  {
805  M_solidBC->createHandler();
806  M_solidBC->fillHandler ( fileName, "solid" );
807 
808  M_FSIoperator->setSolidBC ( M_solidBC->handler() );
809  }
810 }
811 
812 void
813 MultiscaleModelFSI3D::updateBC()
814 {
815 #ifndef FSI_WITHOUT_VELOCITYPROFILE
816  for ( boundaryFlowRateFunctionsContainerIterator_Type i = M_boundaryFlowRateFunctions.begin() ; i != M_boundaryFlowRateFunctions.end() ; ++i )
817  {
818  ( *i )->updateParameters();
819  }
820 #endif
821 
822  if ( M_FSIoperator->isFluid() )
823  {
824  M_fluidBC->updatePhysicalSolverVariables();
825  M_harmonicExtensionBC->updatePhysicalSolverVariables();
826  }
827  if ( M_FSIoperator->isSolid() )
828  {
829  M_solidBC->updatePhysicalSolverVariables();
830  }
831 }
832 
833 void
834 MultiscaleModelFSI3D::setupExporter ( IOFilePtr_Type& exporter, const GetPot& dataFile, const std::string& label )
835 {
836  const std::string exporterType = dataFile ( "exporter/type", "ensight" );
837 #ifdef HAVE_HDF5
838  if ( exporterType.compare ( "hdf5" ) == 0 )
839  {
840  exporter.reset ( new hdf5IOFile_Type() );
841  }
842  else
843 #endif
844  exporter.reset ( new ensightIOFile_Type() );
845 
846  exporter->setDataFromGetPot ( dataFile );
847  exporter->setPrefix ( multiscaleProblemPrefix + "_Model_" + number2string ( M_ID ) + label + "_" + number2string ( multiscaleProblemStep ) );
848  exporter->setPostDir ( multiscaleProblemFolder );
849 }
850 
851 void
852 MultiscaleModelFSI3D::setupImporter ( IOFilePtr_Type& importer, const GetPot& dataFile, const std::string& label )
853 {
854  const std::string importerType = dataFile ( "exporter/type", "ensight" );
855 #ifdef HAVE_HDF5
856  if ( importerType.compare ( "hdf5" ) == 0 )
857  {
858  importer.reset ( new hdf5IOFile_Type() );
859  }
860  else
861 #endif
862  importer.reset ( new ensightIOFile_Type() );
863 
864  importer->setDataFromGetPot ( dataFile );
865  importer->setPrefix ( multiscaleProblemPrefix + "_Model_" + number2string ( M_ID ) + label + "_" + number2string ( multiscaleProblemStep - 1 ) );
866  importer->setPostDir ( multiscaleProblemFolder );
867 }
868 
869 void
870 MultiscaleModelFSI3D::setExporterFluid ( const IOFilePtr_Type& exporter )
871 {
872  M_fluidVelocity.reset ( new vector_Type ( M_FSIoperator->uFESpace().map(), M_exporterFluid->mapType() ) );
873  M_fluidPressure.reset ( new vector_Type ( M_FSIoperator->pFESpace().map(), M_exporterFluid->mapType() ) );
874  M_fluidDisplacement.reset ( new vector_Type ( M_FSIoperator->mmFESpace().map(), M_exporterFluid->mapType() ) );
875 
876  exporter->setMeshProcId ( M_FSIoperator->uFESpace().mesh(), M_FSIoperator->uFESpace().map().comm().MyPID() );
877 
878  exporter->addVariable ( IOData_Type::VectorField, "Velocity (fluid)", M_FSIoperator->uFESpacePtr(), M_fluidVelocity, static_cast<UInt> ( 0 ) );
879  exporter->addVariable ( IOData_Type::ScalarField, "Pressure (fluid)", M_FSIoperator->pFESpacePtr(), M_fluidPressure, static_cast<UInt> ( 0 ) );
880  exporter->addVariable ( IOData_Type::VectorField, "Displacement (fluid)", M_FSIoperator->mmFESpacePtr(), M_fluidDisplacement, static_cast<UInt> ( 0 ) );
881 }
882 
883 void
884 MultiscaleModelFSI3D::setExporterSolid ( const IOFilePtr_Type& exporter )
885 {
886  M_solidDisplacement.reset ( new vector_Type ( M_FSIoperator->dFESpace().map(), M_exporterSolid->mapType() ) );
887  M_solidVelocity.reset ( new vector_Type ( M_FSIoperator->dFESpace().map(), M_exporterSolid->mapType() ) );
888 
889 #ifdef FSI_WITH_STRESSOUTPUT
890  M_solidStressXX.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
891  M_solidStressXY.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
892  M_solidStressXZ.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
893  M_solidStressYY.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
894  M_solidStressYZ.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
895  M_solidStressZZ.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
896  M_solidStressVonMises.reset ( new vector_Type ( M_wallTensionEstimator->feSpaceScalar().map(), M_exporterSolid->mapType() ) );
897 #endif
898 
899  exporter->setMeshProcId ( M_FSIoperator->dFESpace().mesh(), M_FSIoperator->dFESpace().map().comm().MyPID() );
900 
901  exporter->addVariable ( IOData_Type::VectorField, "Displacement (solid)", M_FSIoperator->dFESpacePtr(), M_solidDisplacement, static_cast<UInt> ( 0 ) );
902  exporter->addVariable ( IOData_Type::VectorField, "Velocity (solid)", M_FSIoperator->dFESpacePtr(), M_solidVelocity, static_cast<UInt> ( 0 ) );
903 
904 #ifdef FSI_WITH_STRESSOUTPUT
905  exporter->addVariable ( IOData_Type::ScalarField, "Stress XX (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressXX, static_cast<UInt> ( 0 ) );
906  exporter->addVariable ( IOData_Type::ScalarField, "Stress XY (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressXY, static_cast<UInt> ( 0 ) );
907  exporter->addVariable ( IOData_Type::ScalarField, "Stress XZ (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressXZ, static_cast<UInt> ( 0 ) );
908  exporter->addVariable ( IOData_Type::ScalarField, "Stress YY (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressYY, static_cast<UInt> ( 0 ) );
909  exporter->addVariable ( IOData_Type::ScalarField, "Stress YZ (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressYZ, static_cast<UInt> ( 0 ) );
910  exporter->addVariable ( IOData_Type::ScalarField, "Stress ZZ (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressZZ, static_cast<UInt> ( 0 ) );
911  exporter->addVariable ( IOData_Type::ScalarField, "Stress Von Mises (solid)", M_wallTensionEstimator->feSpaceScalarPtr(), M_solidStressVonMises, static_cast<UInt> ( 0 ) );
912 #endif
913 }
914 
915 void
916 MultiscaleModelFSI3D::exportFluidSolution()
917 {
918  if ( M_FSIoperator->isFluid() )
919  {
920  M_FSIoperator->exportFluidVelocity ( *M_fluidVelocity );
921  M_FSIoperator->exportFluidPressure ( *M_fluidPressure );
922  M_FSIoperator->exportFluidDisplacement ( *M_fluidDisplacement );
923 
924 #ifndef FSI_WITH_EXTERNALPRESSURE
925  // Add the external pressure
926  *M_fluidPressure += M_externalPressureScalar;
927 #endif
928  }
929 }
930 
931 void
932 MultiscaleModelFSI3D::exportSolidSolution()
933 {
934  if ( M_FSIoperator->isSolid() )
935  {
936  M_FSIoperator->exportSolidDisplacement ( *M_solidDisplacement );
937  M_FSIoperator->exportSolidVelocity ( *M_solidVelocity );
938 
939 #ifdef FSI_WITH_STRESSOUTPUT
940  M_wallTensionEstimator->setDisplacement ( *M_solidDisplacement );
941  M_wallTensionEstimator->analyzeTensionsRecoveryVonMisesStress ( );
942 
943  M_wallTensionEstimator->exportStressXX ( *M_solidStressXX );
944  M_wallTensionEstimator->exportStressXY ( *M_solidStressXY );
945  M_wallTensionEstimator->exportStressXZ ( *M_solidStressXZ );
946  M_wallTensionEstimator->exportStressYY ( *M_solidStressYY );
947  M_wallTensionEstimator->exportStressYZ ( *M_solidStressYZ );
948  M_wallTensionEstimator->exportStressZZ ( *M_solidStressZZ );
949  M_wallTensionEstimator->exportStressVonMises ( *M_solidStressVonMises );
950 #endif
951  }
952 }
953 
954 void
955 MultiscaleModelFSI3D::setupLinearModel()
956 {
957 
958 #ifdef HAVE_LIFEV_DEBUG
959  debugStream ( 8140 ) << "MultiscaleModelFSI3D::setupLinearModel() \n";
960 #endif
961 
962  // The linear BCHandler is a copy of the original BCHandler with all BCFunctions giving zero
963  M_linearFluidBC.reset ( new bc_Type ( *M_fluidBC->handler() ) );
964  M_linearSolidBC.reset ( new bc_Type ( *M_solidBC->handler() ) );
965 
966  // Set all the BCFunctions to zero
967  BCFunctionBase bcBaseDeltaZero;
968  bcBaseDeltaZero.setFunction ( std::bind ( &MultiscaleModelFSI3D::bcFunctionDeltaZero, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
969 
970  for ( bc_Type::bcBaseIterator_Type i = M_linearFluidBC->begin() ; i != M_linearFluidBC->end() ; ++i )
971  {
972  i->setBCFunction ( bcBaseDeltaZero );
973  }
974 
975  for ( bc_Type::bcBaseIterator_Type i = M_linearSolidBC->begin() ; i != M_linearSolidBC->end() ; ++i )
976  {
977  i->setBCFunction ( bcBaseDeltaZero );
978  }
979 
980  // Setup linear solution & the RHS
981  M_linearSolution.reset ( new vector_Type ( *M_FSIoperator->couplingVariableMap() ) );
982  M_linearRHS.reset ( new vector_Type ( *M_FSIoperator->couplingVariableMap() ) );
983 }
984 
985 void
986 MultiscaleModelFSI3D::updateLinearModel()
987 {
988 
989 #ifdef HAVE_LIFEV_DEBUG
990  debugStream ( 8140 ) << "MultiscaleModelFSI3D::updateLinearModel() \n";
991 #endif
992 
993  //Create the RHS
994  *M_linearRHS *= 0;
995  M_FSIoperator->bcManageVectorRHS ( M_linearFluidBC, M_linearSolidBC, *M_linearRHS );
996 }
997 
998 void
999 MultiscaleModelFSI3D::solveLinearModel ( bool& solveLinearSystem )
1000 {
1001 
1002 #ifdef HAVE_LIFEV_DEBUG
1003  debugStream ( 8140 ) << "MultiscaleModelFSI3D::solveLinearModel() \n";
1004 #endif
1005 
1006  if ( !solveLinearSystem )
1007  {
1008  return;
1009  }
1010 
1011  imposePerturbation();
1012 
1013  updateLinearModel();
1014 
1015  //Solve the linear problem
1016  displayModelStatus ( "Solve linear" );
1017  M_FSIoperator->solveJac ( *M_linearSolution, *M_linearRHS, 0. );
1018 
1019  resetPerturbation();
1020 
1021  //This flag avoid recomputation of the same system
1022  solveLinearSystem = false;
1023 }
1024 
1025 void
1026 MultiscaleModelFSI3D::imposePerturbation()
1027 {
1028 
1029 #ifdef HAVE_LIFEV_DEBUG
1030  debugStream ( 8140 ) << "MultiscaleModelFSI3D::imposePerturbation() \n";
1031 #endif
1032 
1033  for ( multiscaleCouplingsContainerConstIterator_Type i = M_couplings.begin(); i < M_couplings.end(); ++i )
1034  if ( ( *i )->isPerturbed() )
1035  {
1036  BCFunctionBase bcBaseDeltaOne;
1037 
1038  multiscaleID_Type boundaryID ( ( *i )->boundaryID ( ( *i )->modelGlobalToLocalID ( M_ID ) ) );
1039  if ( M_boundaryFlagsAreaPerturbed[boundaryID] == true )
1040  {
1041  for ( boundaryAreaFunctionsContainerIterator_Type j = M_boundaryAreaFunctions.begin(); j < M_boundaryAreaFunctions.end(); ++j )
1042  if ( ( *j )->fluidFlag() == boundaryFlag ( boundaryID ) )
1043  {
1044  bcBaseDeltaOne.setFunction ( std::bind ( &FSI3DBoundaryAreaFunction::functionLinear, *j, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
1045  M_linearSolidBC->findBCWithFlag ( M_boundaryFlagsArea[boundaryID] ).setBCFunction ( bcBaseDeltaOne );
1046  }
1047  }
1048  else
1049  {
1050  bcBaseDeltaOne.setFunction ( std::bind ( &MultiscaleModelFSI3D::bcFunctionDeltaOne, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
1051  M_linearFluidBC->findBCWithFlag ( boundaryFlag ( boundaryID ) ).setBCFunction ( bcBaseDeltaOne );
1052  }
1053 
1054  break;
1055  }
1056 }
1057 
1058 void
1059 MultiscaleModelFSI3D::resetPerturbation()
1060 {
1061 
1062 #ifdef HAVE_LIFEV_DEBUG
1063  debugStream ( 8140 ) << "MultiscaleModelFSI3D::resetPerturbation() \n";
1064 #endif
1065 
1066  for ( multiscaleCouplingsContainerConstIterator_Type i = M_couplings.begin(); i < M_couplings.end(); ++i )
1067  if ( ( *i )->isPerturbed() )
1068  {
1069  BCFunctionBase bcBaseDeltaZero;
1070  bcBaseDeltaZero.setFunction ( std::bind ( &MultiscaleModelFSI3D::bcFunctionDeltaZero, this, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4, std::placeholders::_5 ) );
1071 
1072  multiscaleID_Type boundaryID ( ( *i )->boundaryID ( ( *i )->modelGlobalToLocalID ( M_ID ) ) );
1073  if ( M_boundaryFlagsAreaPerturbed[boundaryID] == true )
1074  {
1075  M_linearSolidBC->findBCWithFlag ( M_boundaryFlagsArea[boundaryID] ).setBCFunction ( bcBaseDeltaZero );
1076  }
1077  else
1078  {
1079  M_linearFluidBC->findBCWithFlag ( boundaryFlag ( boundaryID ) ).setBCFunction ( bcBaseDeltaZero );
1080  M_boundaryFlagsAreaPerturbed[boundaryID] = true;
1081  }
1082 
1083  break;
1084  }
1085 }
1086 
1087 } // Namespace multiscale
1088 } // Namespace LifeV
FSI_WITH_BOUNDARYAREA
#define FSI_WITH_BOUNDARYAREA
Definition: MultiscaleModelFSI3D.hpp:54
FSI_WITH_STRESSOUTPUT
#define FSI_WITH_STRESSOUTPUT
Definition: MultiscaleModelFSI3D.hpp:57