LifeV
ETRobinMembraneSolver.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
30 
31  @author Claudia Colciago <claudia.colciago@epfl.ch>
32  @date 08-11-2012
33  */
34 
35 // Tell the compiler to ignore specific kind of warnings:
36 #pragma GCC diagnostic ignored "-Wunused-variable"
37 #pragma GCC diagnostic ignored "-Wunused-parameter"
38 
39 #include <Epetra_ConfigDefs.h>
40 #ifdef EPETRA_MPI
41 #include <mpi.h>
42 #include <Epetra_MpiComm.h>
43 #else
44 #include <Epetra_SerialComm.h>
45 #endif
46 
47 //Tell the compiler to restore the warning previously silented
48 #pragma GCC diagnostic warning "-Wunused-variable"
49 #pragma GCC diagnostic warning "-Wunused-parameter"
50 
51 
52 #include <lifev/core/array/MatrixEpetra.hpp>
53 #include <lifev/core/array/MapEpetra.hpp>
54 #include <lifev/core/mesh/MeshData.hpp>
55 #include <lifev/core/mesh/MeshPartitioner.hpp>
56 #include <lifev/navier_stokes/fem/TimeAdvanceBDFNavierStokes.hpp>
57 
58 #include <lifev/core/fem/BCManage.hpp>
59 #ifdef HAVE_HDF5
60 #include <lifev/core/filter/ExporterHDF5.hpp>
61 #endif
62 #include <lifev/core/filter/ExporterEnsight.hpp>
63 
64 #include <lifev/eta/expression/Integrate.hpp>
65 #include <lifev/navier_stokes/solver/OseenData.hpp>
66 #include <lifev/navier_stokes/solver/StabilizationIP.hpp>
67 
68 
69 #include "ETRobinMembraneSolver.hpp"
70 #include "ud_functions.hpp"
71 
72 #include <iostream>
73 
74 //#define FLUX 1
75 
76 using namespace LifeV;
77 
78 
79 const Real PI = 3.141592653589793;
80 
81 
82 struct ETRobinMembraneSolver::Private
83 {
84  typedef std::function<Real ( Real const&, Real const&, Real const&, Real const&, ID const& ) > fct_type;
85 
86  std::string data_file_name;
87 
88  // Data for the fluid
89  Real Re;
90  Real nu; /**< viscosity (in m^2/s) */
91  Real H; /**< height and width of the domain (in m) */
92  Real D; /**< diameter of the cylinder (in m) */
93  Real density;
94  Real R; //radius
95 
96 
97  //Data for the membrnae
98  Real rhos;
99  Real Hs;
100  Real ni;
101  Real E;
102 
103  //Discretization choices
104  UInt transpirationOrder;
105  std::string stabilization;
106  std::string initial_sol;
107  Real numLM;
108  bool useFlowRate;
109 
110  std::shared_ptr<Epetra_Comm> comm;
111 
112 };
113 
114 ETRobinMembraneSolver::ETRobinMembraneSolver ( int argc, char** argv )
115  :
116  M_d ( new Private )
117 {
118  GetPot command_line (argc, argv);
119  std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
120  GetPot dataFile ( data_file_name );
121 
122  M_d->data_file_name = data_file_name;
123 
124  M_d->Re = dataFile ( "fluid/physics/Re", 1. );
125  M_d->nu = dataFile ( "fluid/physics/viscosity", 1. );
126  M_d->H = dataFile ( "fluid/physics/H", 20. );
127  M_d->D = dataFile ( "fluid/physics/D", 1. );
128  M_d->R = M_d->D / 2;
129  M_d->density = dataFile ( "fluid/physics/density", 1. );
130 
131  M_d->rhos = dataFile ( "membrane/physics/density_sol", 1. );
132  M_d->E = dataFile ( "membrane/physics/young", 1. );
133  M_d->ni = dataFile ( "membrane/physics/poisson", 1. );
134  M_d->Hs = dataFile ( "membrane/physics/wall_thickness", 1. );
135 
136 
137  M_d->initial_sol = (std::string) dataFile ( "fluid/problem/initial_sol", "none");
138  M_d->stabilization = (std::string) dataFile ( "fluid/problem/stabilization", "none");
139  M_d->numLM = dataFile ( "fluid/problem/numLM" , 0 );
140  M_d->useFlowRate = 0;
141  if ( M_d->numLM > 0 )
142  {
143  M_d->useFlowRate = 1 ;
144  }
145  M_d->transpirationOrder = dataFile ( "fluid/problem/transpiration_order" , 0 );
146 
147 #ifdef EPETRA_MPI
148  std::cout << "mpi initialization ... " << std::endl;
149 
150  int ntasks = 0;
151  M_d->comm.reset ( new Epetra_MpiComm ( MPI_COMM_WORLD ) );
152  if (!M_d->comm->MyPID() )
153  {
154  std::cout << "My PID = " << M_d->comm->MyPID() << " out of " << ntasks << " running." << std::endl;
155  std::cout << "fluid density = " << M_d->density << std::endl
156  << "nu = " << M_d->nu << std::endl
157  << "poisson = " << M_d->ni << std::endl
158  << "initial solution = " << M_d->initial_sol << std::endl
159  << "Young modulus = " << M_d->E << std::endl
160  << "Wall thickness = " << M_d->Hs << std::endl;
161 
162  }
163 
164 #else
165  M_d->comm.reset ( new Epetra_SerialComm() );
166 #endif
167 
168 }
169 
170 void
171 ETRobinMembraneSolver::run()
172 
173 {
174  bool verbose = (M_d->comm->MyPID() == 0);
175 
176  //------------Creating file to store data------------------------------
177 
178  GetPot dataFile ( M_d->data_file_name );
179 
180  // NS
181  Real NSSupg (dataFile ("fluid/space_discretization/supg_coeff", 1e-3) ); // 1e-3
182  Real NSPspg (dataFile ("fluid/space_discretization/pspg_coeff", 1e-3) ); // 1e-3
183  Real NSdivdiv (dataFile ("fluid/space_discretization/dividiv_coeff", 5e-2) );
184  Real OUTLET (dataFile ("fluid/problem/boundary_flags/outlet", 3) );
185  Real INLET (dataFile ("fluid/problem/boundary_flags/inlet", 2) );
186  Real WALL (dataFile ("fluid/problem/boundary_flags/wall", 100) );
187  Real RING (dataFile ("fluid/problem/boundary_flags/ring_in", 20) );
188  Real RING2 (dataFile ("fluid/problem/boundary_flags/ring_out", 30) );
189 
190  std::shared_ptr<OseenData> oseenData (new OseenData() );
191  oseenData->setup ( dataFile );
192 
193  MeshData meshData;
194  meshData.setup (dataFile, "fluid/space_discretization");
195 
196  Real LameI = (/*M_d->Hs * */M_d->E * M_d->ni ) / ( ( 1 - M_d->ni * M_d->ni ) );
197  Real LameII = /*M_d->Hs * */M_d->E / ( 2 * ( 1 + M_d->ni ) );
198 
199  if (verbose)
200  {
201  std::cout << " Outlet : " << OUTLET << std::endl;
202  std::cout << " Wall : " << WALL << std::endl;
203  }
204  if (verbose)
205  {
206  std::cout << " LameI : " << LameI << std::endl;
207  std::cout << " LameII : " << LameII << std::endl;
208  }
209 
210  //-----------------------------------------------------------------------
211 
212  //----------Creating Mesh and FESpaces-----------------------------------
213 
214  std::shared_ptr< mesh_type > fullMeshPtr (new mesh_type);
215  readMesh (*fullMeshPtr, meshData);
216 
217  MeshPartitioner< mesh_type > meshPart (fullMeshPtr, M_d->comm);
218 
219  std::string uOrder = dataFile ( "fluid/space_discretization/vel_order", "P1");
220  std::string pOrder = dataFile ( "fluid/space_discretization/press_order", "P1");
221 
222  if (verbose)
223  {
224  std::cout << "Building the FE spaces ... " << std::flush;
225  }
226 
227  M_uFESpace.reset ( new FESpace< mesh_type, MapEpetra > ( meshPart, uOrder, 3, M_d->comm ) );
228  M_uCompFESpace.reset ( new FESpace< mesh_type, MapEpetra > ( meshPart, uOrder, 1, M_d->comm ) );
229  M_pFESpace.reset ( new FESpace< mesh_type, MapEpetra > ( meshPart, pOrder, 1, M_d->comm ) );
230 
231  M_ETuFESpace.reset ( new ETFESpace< mesh_type, MapEpetra, 3, 3 > ( meshPart, & (M_uFESpace->refFE() ), M_d->comm ) );
232  M_ETpFESpace.reset ( new ETFESpace< mesh_type, MapEpetra, 3, 1 > ( meshPart, & (M_pFESpace->refFE() ), M_d->comm ) );
233 
234  if (verbose)
235  {
236  std::cout << "ok." << std::endl;
237  }
238 
239  UInt totalVelDof = M_uFESpace->map().map ( Unique )->NumGlobalElements();
240  UInt totalPressDof = M_pFESpace->map().map ( Unique )->NumGlobalElements();
241 
242  if (verbose)
243  {
244  std::cout << "Total Velocity DOF = " << totalVelDof << std::endl;
245  }
246  if (verbose)
247  {
248  std::cout << "Total Pressure DOF = " << totalPressDof << std::endl;
249  }
250  if (verbose)
251  {
252  std::cout << "Total FLux DOF = " << M_d->numLM << std::endl;
253  }
254 
255  MapEpetra fluidMap ( M_uFESpace->map() );
256  MapEpetra fluxMap ( M_d->numLM, M_d->comm );
257  MapEpetra fullMap ( M_uFESpace->map() + M_pFESpace->map() + fluxMap );
258 
259  DOFInterface3Dto3D interfaceDOF ( M_uFESpace->refFE() , M_uFESpace->dof() );
260  interfaceDOF.update ( *M_uFESpace->mesh() , WALL , *M_uFESpace->mesh() , WALL , 0);
261  createInterfaceMap ( interfaceDOF.localDofMap() , M_uFESpace->dof() );
262 
263  //------------------FLUID PROBLEM--------------------------------
264 
265  if (verbose)
266  {
267  std::cout << "Calling the solver constructors ... ";
268  }
269 
270  SolverAztecOO NSSolver;
271  NSSolver.setCommunicator (M_d->comm);
272  NSSolver.setDataFromGetPot (dataFile, "solver");
273  NSSolver.setupPreconditioner (dataFile, "prec");
274 
275  if (verbose)
276  {
277  std::cout << "done." << std::endl;
278  }
279 
280  //---------------------------------------------------------
281 
282  //--------------------Creating BDF objects---------------------------
283 
284  if (verbose)
285  {
286  std::cout << "Calling the TimeAdvance constructors ... ";
287  }
288 
289  Real dt = oseenData->dataTime()->timeStep();
290  Real t0 = oseenData->dataTime()->initialTime();
291  Real tFinal = oseenData->dataTime()->endTime();
292 
293  TimeAdvanceBDFNavierStokes<vector_type> fluidTimeAdvance;
294  fluidTimeAdvance.setup (oseenData->dataTimeAdvance()->orderBDF() );
295 
296  TimeAdvanceBDF<vector_type> dispTimeAdvance;
297  dispTimeAdvance.setup (oseenData->dataTimeAdvance()->orderBDF() );
298 
299  //TimeAdvance for the fluid
300  fluidTimeAdvance.bdfVelocity().setTimeStep (oseenData->dataTime()->timeStep() );
301  if (verbose)
302  {
303  fluidTimeAdvance.showMe();
304  }
305 
306  //TimeAdvance for the displacement
307  dispTimeAdvance.setTimeStep (oseenData->dataTime()->timeStep() );
308  if (verbose)
309  {
310  dispTimeAdvance.showMe();
311  }
312 
313  Real alpha = fluidTimeAdvance.bdfVelocity().coefficientFirstDerivative (t0);
314 
315  if (verbose)
316  {
317  std::cout << " done. \n";
318  }
319 
320  //---------------------------------------------------------
321 
322  //--------------Boundary Conditions-------------------------
323 
324  if (verbose)
325  {
326  std::cout << "Calling the BCHandler constructor ... ";
327  }
328 
329  //creating BCHandler objects
330  BCHandler bcHFluid;
331  BCFunctionBase vel_in (flatNormalVelInlet /*linearVelInletCylinder*/);
332  BCFunctionBase flux_in (linearInletCylinder);
333  BCFunctionBase uZero ( fZero );
334  BCFunctionBase press_out ( p0 );
335 
336  if (M_d->useFlowRate == true )
337  {
338  bcHFluid.addBC ("InFlow" , INLET, Flux, Normal, flux_in );
339  bcHFluid.addBC ("InFlow_2" , OUTLET, Natural, Full, uZero, 3 );
340 
341  //In order to have a well defined problem you need to set conditions on boundaries of boundaries
342  bcHFluid.addBC ("Ring" , RING, EssentialEdges, Full, uZero, 3 );
343  bcHFluid.addBC ("Ring3" , RING2, EssentialEdges, Full, uZero, 3 );
344 
345  bcHFluid.setOffset ("InFlow", totalVelDof + totalPressDof);
346 
347  }
348  else
349  {
350 
351  bcHFluid.addBC ("InFlow" , INLET, Essential, Full, vel_in, 3 );
352  bcHFluid.addBC ("InFlow_2" , OUTLET, Natural, Normal, press_out );
353 
354  //In order to have a well defined problem you need to set conditions on boundaries of boundaries
355  bcHFluid.addBC ("Ring" , RING, EssentialEdges, Full, uZero, 3 );
356  bcHFluid.addBC ("Ring3" , RING2, EssentialEdges, Full, uZero, 3 );
357 
358  }
359 
360  if (verbose)
361  {
362  std::cout << " BC done \n";
363  }
364 
365  //----------------------------------------------------------------------------
366 
367  vector_block_type NSSolution (M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap, Unique);
368  vector_type velocitySolution (M_ETuFESpace->map(), Repeated);
369  vector_type dispSolution (M_ETuFESpace->map(), Repeated);
370 
371  NSSolution *= 0.0;
372 
373  //-------------Creating Exporter Object------------------------------------------------
374 
375  if (verbose)
376  {
377  std::cout << "Calling the Exporter constructor ... ";
378  }
379 
380  std::string const exporterFileName = dataFile ( "exporter/filename", "cylinder");
381  UInt exportEach = dataFile ("exporter/each", 1);
382 
383  ExporterHDF5< mesh_type > exporter ( dataFile, meshPart.meshPartition(), exporterFileName, M_d->comm->MyPID() );
384 
385  vectorPtr_type velAndPressureExporter ( new vector_type (NSSolution, exporter.mapType() ) );
386  vectorPtr_type dispExporter (new vector_type (dispSolution, Repeated ) );
387 
388  exporter.addVariable ( ExporterData< mesh_type >::VectorField, "f-velocity",
389  M_uFESpace, velAndPressureExporter, UInt (0) );
390 
391  exporter.addVariable ( ExporterData< mesh_type >::ScalarField, "f-pressure",
392  M_pFESpace, velAndPressureExporter, UInt (3 * M_uFESpace->dof().numTotalDof() ) );
393 
394  exporter.addVariable ( ExporterData< mesh_type >::VectorField, "s-displacement",
395  M_uFESpace, dispExporter, UInt (0) ) ;
396  if (verbose)
397  {
398  std::cout << " done. \n";
399  }
400 
401  //-----------------------------------------------------------------------------
402 
403  //--------------------------Initialization--------------------------------------
404 
405 
406 
407  if (M_d->initial_sol == "restart")
408  {
409  ASSERT (oseenData->dataTimeAdvance()->orderBDF() == 1, "Restart works only for BDF1");
410 
411  if (verbose)
412  {
413  std::cout << std::endl;
414  }
415  if (verbose)
416  {
417  std::cout << "Restoring the previous solution ... " << std::endl;
418  }
419 
420  std::string filename = dataFile ("importer/filename", "cylinder");
421 
422  LifeV::ExporterHDF5<mesh_type> importer ( dataFile, filename );
423  importer.setMeshProcId ( M_uFESpace->mesh(), M_d->comm->MyPID() );
424 
425  vectorPtr_type velAndPressureImporter ( new vector_type (NSSolution, importer.mapType() ) );
426  vectorPtr_type dispImporter (new vector_type (dispSolution, importer.mapType() ) );
427 
428  importer.addVariable ( ExporterData<mesh_type>::VectorField,
429  "f-velocity",
430  M_uFESpace,
431  velAndPressureImporter,
432  UInt ( 0 ) );
433 
434  importer.addVariable ( ExporterData<mesh_type>::ScalarField,
435  "f-pressure",
436  M_pFESpace,
437  velAndPressureImporter,
438  3 * M_uFESpace->dof().numTotalDof() );
439 
440  importer.addVariable ( ExporterData<mesh_type>::VectorField,
441  "s-displacement",
442  M_uFESpace,
443  dispImporter,
444  UInt ( 0 ) );
445 
446  exporter.setTimeIndex ( importer.importFromTime (t0) );
447 
448  Real norm = velAndPressureImporter->norm2();
449  *velAndPressureExporter = *velAndPressureImporter;
450  *dispExporter = *dispImporter;
451  if (verbose)
452  {
453  std::cout << " f- restart solution norm = " << norm << std::endl;
454  }
455 
456 
457  }
458 
459  velocitySolution.subset ( *velAndPressureExporter );
460  dispSolution = *dispExporter;
461 
462  fluidTimeAdvance.bdfVelocity().setInitialCondition ( * (velAndPressureExporter) );
463  dispTimeAdvance.setInitialCondition ( dispSolution );
464 
465  fluidTimeAdvance.bdfVelocity().updateRHSContribution ( oseenData->dataTime()->timeStep() );
466  dispTimeAdvance.updateRHSContribution (oseenData->dataTime()->timeStep() );
467 
468  if (verbose)
469  {
470  std::cout << " done \n";
471  }
472 
473 
474  //-------------------------------------------------------------------------------------
475 
476  std::shared_ptr<matrix_block_type> NSMatrixConstant (new matrix_block_type ( M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap ) );
477  *NSMatrixConstant *= 0.0;
478  std::shared_ptr<matrix_block_type> NSMatrixSteadyStokes (new matrix_block_type ( M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap ) );
479  *NSMatrixSteadyStokes *= 0.0;
480 
481  //----------------------Temporal Loop----------------------------------
482 
483  if (verbose)
484  {
485  std::cout << std::endl;
486  }
487  if (verbose)
488  {
489  std::cout << " ### Simulation times ### " << std::endl;
490  }
491  if (verbose)
492  {
493  std::cout << " From " << t0 << " to " << tFinal << std::endl;
494  }
495  if (verbose)
496  {
497  std::cout << " Time step: " << dt << std::endl;
498  }
499  if (verbose)
500  {
501  std::cout << std::endl;
502  }
503 
504  Real currentTime (t0);
505  UInt niter (0);
506 
507  *NSMatrixConstant *= 0.0;
508 
509  QuadratureBoundary myBDQR (buildTetraBDQR (quadRuleTria4pt) );
510  vector_type robinExt ( M_pFESpace->map(), Unique);
511  M_pFESpace->interpolate ( uZero , robinExt, 0);
512  vector_type hSolid (M_pFESpace->map(), Unique);
513  hSolid += M_d->Hs; // * robinExt;
514  //robinExt *= 0.0;
515 
516  {
517 
518  using namespace ExpressionAssembly;
519 
520  integrate (
521  elements (M_ETuFESpace->mesh() ), // Mesh
522 
523  M_uFESpace->qr(), // QR
524 
525  M_ETuFESpace,
526  M_ETuFESpace,
527 
528  // Viscous Term
529  M_d->nu * dot (grad (phi_i) , grad (phi_j) )
530 
531  )
532  >> NSMatrixSteadyStokes->block (0, 0);
533 
534 
535 
536  integrate (
537  elements (M_ETuFESpace->mesh() ), // Mesh
538 
539  M_uFESpace->qr(), // QR
540 
541  M_ETuFESpace,
542  M_ETpFESpace,
543 
544  value (-1.0) *phi_j * div (phi_i)
545 
546 
547  )
548  >> NSMatrixSteadyStokes->block (0, 1);
549 
550  integrate (
551  elements (M_ETuFESpace->mesh() ), // Mesh
552 
553  M_uFESpace->qr(), // QR
554 
555  M_ETpFESpace,
556  M_ETuFESpace,
557 
558  value (1.0) *phi_i * div (phi_j)
559 
560  )
561  >> NSMatrixSteadyStokes->block (1, 0);
562 
563  NSMatrixSteadyStokes->globalAssemble();
564 
565  integrate (
566  elements (M_ETuFESpace->mesh() ), // Mesh
567 
568  M_uFESpace->qr(), // QR
569 
570  M_ETuFESpace,
571  M_ETuFESpace,
572 
573  // Intertial Term
574  M_d->density
575  * ( value (alpha / dt) * dot (phi_i, phi_j) )
576 
577  )
578  >> NSMatrixConstant->block (0, 0);
579 
580 
581  integrate ( boundary (M_ETuFESpace->mesh(), WALL),
582  myBDQR,
583 
584  M_ETuFESpace,
585  M_ETuFESpace,
586 
587  //Boundary Mass
588  ( value (M_d->rhos /*M_d->Hs*/ * alpha / dt) * value ( M_ETpFESpace, hSolid ) + value ( M_ETpFESpace, robinExt ) * ( 0.1 + dt ) ) * dot ( phi_j, phi_i )
589 
590  )
591  >> NSMatrixConstant->block (0, 0);
592 
593  }
594 
595 
596  MatrixSmall<3, 3> Eye;
597  Eye *= 0.0;
598  Eye[0][0] = 1;
599  Eye[1][1] = 1;
600  Eye[2][2] = 1;
601 
602  {
603  using namespace ::LifeV::ExpressionAssembly;
604 
605 
606  integrate ( boundary (M_ETuFESpace->mesh(), WALL),
607  myBDQR,
608 
609  M_ETuFESpace,
610  M_ETuFESpace,
611 
612  //Boundary Stiffness
613  ( dt / alpha ) *
614  2 * LameII * value ( M_ETpFESpace, hSolid ) *
615  0.5 * dot ( ( grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) )
616  + transpose (grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) ),
617  ( grad (phi_i) + ( (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) ) ) +
618 
619  ( dt / alpha ) * value ( M_ETpFESpace, hSolid ) *
620  LameI * dot ( value ( Eye ) , ( grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) ) )
621  * dot ( value ( Eye ) , ( grad (phi_i) + (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) )
622 
623  )
624  >> NSMatrixConstant->block (0, 0);
625 
626  }
627 
628  *NSMatrixConstant += *NSMatrixSteadyStokes;
629  NSMatrixConstant->globalAssemble();
630 
631  details::StabilizationIP<mesh_type, DOF> M_ipStabilization;
632  M_ipStabilization.setFeSpaceVelocity ( *M_uFESpace );
633  M_ipStabilization.setViscosity ( M_d->nu );
634 
635  // Parameters from J. Michalik
636  M_ipStabilization.setGammaBeta ( dataFile ( "ipstab/gammaBeta", 1.0) );
637  M_ipStabilization.setGammaDiv ( dataFile ( "ipstab/gammaDiv", 0.2) );
638  M_ipStabilization.setGammaPress ( dataFile ( "ipstab/gammaPress", 0.5) );
639 
640  vector_type steadyResidual ( fullMap, Unique );
641 
642  // if (M_d->initial_sol == "steady")
643  // {
644 
645  // #ifdef FLUX
646  // std::shared_ptr<matrix_block_type> NSMatrix (new matrix_block_type ( M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap ) );
647  // *NSMatrix *= 0.0;
648 
649  // #else
650 
651  // std::shared_ptr<matrix_block_type> NSMatrix (new matrix_block_type ( M_ETuFESpace->map() | M_ETpFESpace->map() ) );
652  // *NSMatrix *= 0.0;
653  // #endif
654 
655  // vector_type NSRhsUnique ( fullMap, Unique );
656 
657  // std::shared_ptr<matrix_type> stabMatrix (new matrix_type ( fullMap ));
658  // M_ipStabilization.apply ( *stabMatrix, velocitySolution, false );
659  // stabMatrix->globalAssemble();
660  // *NSMatrix += *NSMatrixSteadyStokes;
661  // *NSMatrix += *stabMatrix;
662 
663  // NSMatrix->globalAssemble();
664 
665  // bcHFluid.bcUpdate ( *meshPart.meshPartition(), M_uFESpace->feBd(), M_uFESpace->dof() );
666  // // bcManage (*NSMatrix, NSRhsUnique,
667  // // *M_uFESpace->mesh(), M_uFESpace->dof(),
668  // // bcHFluid, M_uFESpace->feBd(), 1.0, currentTime);
669 
670  // std::shared_ptr<matrix_type> NSMatrixDiri( new matrix_type( fullMap ) );
671  // *NSMatrixDiri += *NSMatrix;
672  // NSMatrixDiri->globalAssemble();
673 
674  // BCHandler bcHInit;
675  // bcHInit.addBC ("Wall" , WALL, Essential, Full, uZero, 3 );
676  // bcHInit.bcUpdate ( *meshPart.meshPartition(), M_uFESpace->feBd(), M_uFESpace->dof() );
677  // // bcManage (*NSMatrixDiri, NSRhsUnique,
678  // // *M_uFESpace->mesh(), M_uFESpace->dof(),
679  // // bcHInit, M_uFESpace->feBd(), 1.0, currentTime);
680 
681  // NSSolver.setMatrix (*NSMatrixDiri);
682  // NSSolver.solveSystem (NSRhsUnique, NSSolution, NSMatrixDiri);
683 
684  // *velAndPressureExporter = NSSolution;
685 
686  // steadyResidual = (-1) * (*NSMatrix * NSSolution);
687 
688  // vector_type dispRhsInit( M_ETuFESpace->map(), Unique );
689  // dispRhsInit.subset( steadyResidual );
690 
691  // vector_type dispSolutionInit( M_ETuFESpace->map(), Unique );
692  // std::shared_ptr<matrix_type> DispMatrixInit (new matrix_type ( M_ETuFESpace->map() ) );
693  // *DispMatrixInit *= 0.0;
694  // DispMatrixInit->insertOneDiagonal();
695  // DispMatrixInit->insertValueDiagonal( -1 , (*M_interfaceMap) );
696 
697  // {
698  // using namespace ::LifeV::ExpressionAssembly;
699 
700 
701  // integrate ( boundary (M_ETuFESpace->mesh(), WALL),
702  // myBDQR,
703 
704  // M_ETuFESpace,
705  // M_ETuFESpace,
706 
707  // //Boundary Stiffness
708  // //Boundary Mass
709  // value( M_ETpFESpace, robinExt ) * dot ( phi_j, phi_i )
710 
711  // + 2 * LameII *
712  // 0.5 * dot ( ( grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) )
713  // + transpose (grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) ),
714  // ( grad (phi_i) + ( (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) ) ) +
715 
716  // LameI * dot ( value ( Eye ) , ( grad (phi_j) + (-1) * grad (phi_j) * outerProduct ( Nface, Nface ) ) )
717  // * dot ( value ( Eye ) , ( grad (phi_i) + (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) )
718 
719  // )
720  // >> *DispMatrixInit;
721 
722  // }
723 
724  // DispMatrixInit->globalAssemble();
725 
726  // DispMatrixInit->spy("dispMatrix");
727 
728  // BCHandler bcHDiri;
729  // bcHDiri.addBC ("Inlet" , INLET, Essential, Full, uZero, 3 );
730  // bcHDiri.addBC ("Outlet" , OUTLET, Essential, Full, uZero, 3 );
731  // bcHDiri.addBC ("Ring" , RING2, EssentialEdges, Full, uZero, 3 );
732  // bcHDiri.addBC ("Ring" , RING, EssentialEdges, Full, uZero, 3 );
733  // bcHDiri.bcUpdate ( *meshPart.meshPartition(), M_uFESpace->feBd(), M_uFESpace->dof() );
734  // // bcManage (*DispMatrixInit, dispRhsInit,
735  // // *M_uFESpace->mesh(), M_uFESpace->dof(),
736  // // bcHDiri, M_uFESpace->feBd(), 1.0, currentTime);
737 
738  // NSSolver.setMatrix (*DispMatrixInit);
739  // NSSolver.solveSystem (dispRhsInit, dispSolutionInit , DispMatrixInit);
740 
741  // *dispExporter = dispSolutionInit;
742 
743  // dispSolution = dispSolutionInit;
744 
745  // dispTimeAdvance.shiftRight ( dispSolution );
746 
747  // dispTimeAdvance.updateRHSContribution (oseenData->dataTime()->timeStep() );
748 
749  // velocitySolution.subset (NSSolution);
750 
751  // fluidTimeAdvance.bdfVelocity().shiftRight ( NSSolution );
752 
753  // fluidTimeAdvance.bdfVelocity().updateRHSContribution ( oseenData->dataTime()->timeStep() );
754  // }
755 
756  exporter.postProcess (t0);
757 
758 
759  while ( currentTime < tFinal)
760  {
761  LifeChrono ChronoIteration;
762  ChronoIteration.start();
763 
764  currentTime += dt;
765  niter += 1;
766 
767  if (verbose)
768  {
769  std::cout << std::endl;
770  }
771  if (verbose)
772  {
773  std::cout << "----------------------------" << std::endl;
774  }
775  if (verbose)
776  {
777  std::cout << " Time : " << currentTime << std::endl;
778  }
779  if (verbose)
780  {
781  std::cout << " Iter : " << niter << std::endl;
782  }
783  if (verbose)
784  {
785  std::cout << "----------------------------" << std::endl;
786  }
787  if (verbose)
788  {
789  std::cout << std::endl;
790  }
791 
792  vector_type velocityExtrapolated (velocitySolution, Repeated);
793  fluidTimeAdvance.bdfVelocity().extrapolation ( velocityExtrapolated );
794  vector_type velocityBdfRHS (velocitySolution, Repeated);
795  velocityBdfRHS = fluidTimeAdvance.bdfVelocity().rhsContributionFirstDerivative();
796 
797  vector_type dispExtrapolated (dispSolution, Repeated);
798  dispTimeAdvance.extrapolation ( dispExtrapolated );
799  vector_type dispBdfRHS (fluidMap, Repeated);
800  dispBdfRHS = dispTimeAdvance.rhsContributionFirstDerivative();
801  vector_type dtDisp (fluidMap, Repeated);
802  dtDisp = ( alpha / dt ) * dispSolution - dispTimeAdvance.rhsContributionFirstDerivative();
803 
804  std::shared_ptr<matrix_block_type> NSMatrix (new matrix_block_type ( M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap ) );
805  *NSMatrix *= 0.0;
806 
807 #define DIVDIV_TEST value(NSdivdiv) * h_K * div(phi_i)
808 
809 #define SUPG_TEST value(NSSupg) * h_K * (grad(phi_i)*eval(normalize,value(M_ETuFESpace,velocitySolution)))
810 
811 #define PSPG_TEST value(NSPspg) * h_K * grad(phi_i)
812 
813 #define TRANSP_GRADGRAD grad(phi_j) * ( grad(M_ETuFESpace, dispExtrapolated) ) * ( Eye + (-1) * outerProduct(Nface, Nface) )
814 
815  if (verbose)
816  {
817  std::cout << "[Navier-Stokes] Assembling the matrix ... " << std::flush;
818  }
819 
820  LifeChrono ChronoItem;
821  ChronoItem.start();
822 
823  {
824  std::shared_ptr<NormalizeFct> normalize (new NormalizeFct);
825  using namespace ::LifeV::ExpressionAssembly;
826 
827  integrate (
828  elements (M_ETuFESpace->mesh() ), // Mesh
829  M_uFESpace->qr(), // QR
830 
831  M_ETuFESpace,
832  M_ETuFESpace,
833 
834  // Advection Term
835  M_d->density
836  * ( dot (grad (phi_j) * value (M_ETuFESpace, velocityExtrapolated) , phi_i) )
837 
838  )
839  >> NSMatrix->block (0, 0);
840 
841 
842  if (M_d->stabilization == "sup")
843  {
844 
845  integrate (
846  elements (M_ETuFESpace->mesh() ), // Mesh
847  M_uFESpace->qr(), // QR
848 
849  M_ETuFESpace,
850  M_ETuFESpace,
851 
852  // SUPG
853  dot (
854  grad (phi_j) *value (M_ETuFESpace, velocityExtrapolated) * M_d->density
855  + value (alpha / dt) * M_d->density * phi_j
856  , SUPG_TEST)
857 
858  // Div div
859  + DIVDIV_TEST
860  *div (phi_j)
861 
862 
863  )
864  >> NSMatrix->block (0, 0);
865 
866 
867  integrate (
868  elements (M_ETuFESpace->mesh() ), // Mesh
869  M_uFESpace->qr(), // QR
870 
871  M_ETuFESpace,
872  M_ETpFESpace,
873 
874  // SUPG
875  dot ( grad (phi_j) , SUPG_TEST)
876 
877  )
878  >> NSMatrix->block (0, 1);
879 
880 
881  integrate (
882  elements (M_ETuFESpace->mesh() ), // Mesh
883  M_uFESpace->qr(), // QR
884 
885  M_ETpFESpace,
886  M_ETuFESpace,
887 
888  // PSPG
889  dot (
890  grad (phi_j) *value (M_ETuFESpace, velocityExtrapolated) * value ( M_d->density ) + value (1.0 / dt) * value ( M_d->density ) * phi_j
891  , PSPG_TEST)
892  )
893  >> NSMatrix->block (1, 0);
894 
895 
896  integrate (
897  elements (M_ETuFESpace->mesh() ), // Mesh
898 
899  M_uFESpace->qr(), // QR
900 
901  M_ETpFESpace,
902  M_ETpFESpace,
903 
904  // PSPG
905  dot (
906  grad (phi_j)
907  , PSPG_TEST)
908 
909  )
910  >> NSMatrix->block (1, 1);
911  }
912 
913  }
914 
915  if ( M_d->transpirationOrder == 1 )
916  {
917 
918  {
919  using namespace ::LifeV::ExpressionAssembly;
920 
921 
922  integrate ( boundary (M_ETuFESpace->mesh(), WALL),
923  myBDQR,
924 
925  M_ETuFESpace,
926  M_ETuFESpace,
927 
928  //Boundary Stiffness
929  ( dt / alpha ) *
930  2 * LameII * value ( M_ETpFESpace, hSolid ) *
931  0.5 * dot ( TRANSP_GRADGRAD
932  + transpose ( TRANSP_GRADGRAD ),
933  ( grad (phi_i) + ( (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) ) ) +
934 
935  ( dt / alpha ) *
936  LameI * value ( M_ETpFESpace, hSolid ) * dot ( value ( Eye ) , TRANSP_GRADGRAD )
937  * dot ( value ( Eye ) , ( grad (phi_i) + (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) )
938 
939  )
940  >> NSMatrix->block (0, 0);
941 
942 
943  integrate ( boundary (M_ETuFESpace->mesh(), WALL),
944  myBDQR,
945 
946  M_ETuFESpace,
947  M_ETuFESpace,
948 
949  //Boundary Mass
950  value (M_d->rhos* /*M_d->Hs*/ alpha / dt) * value ( M_ETpFESpace, hSolid) * dot ( grad (phi_j) * value ( M_ETuFESpace, dispExtrapolated) , phi_i )
951  + dot ( grad (phi_j) * value ( M_ETuFESpace, dtDisp ) , phi_i )
952  + value ( M_ETpFESpace, robinExt ) * ( 0.1 + dt / alpha ) * dot ( grad (phi_j) * value ( M_ETuFESpace, dispExtrapolated) , phi_i )
953 
954  )
955  >> NSMatrix->block (0, 0);
956 
957 
958 
959  }
960 
961  }
962 
963 
964  if (M_d->stabilization == "ip")
965  {
966 
967  std::shared_ptr<matrix_type> stabMatrix (new matrix_type ( fullMap ) );
968  M_ipStabilization.apply ( *stabMatrix, velocityExtrapolated, false );
969  stabMatrix->globalAssemble();
970  *NSMatrix += *stabMatrix;
971  }
972 
973 
974 
975  ChronoItem.stop();
976  if (verbose)
977  {
978  std::cout << ChronoItem.diff() << " s" << std::endl;
979  }
980 
981 
982  if (verbose)
983  {
984  std::cout << "[Navier-Stokes] Adding constant parts ... " << std::flush;
985  }
986  ChronoItem.start();
987 
988  *NSMatrix += *NSMatrixConstant;
989 
990  ChronoItem.stop();
991  if (verbose)
992  {
993  std::cout << ChronoItem.diff() << " s" << std::endl;
994  }
995 
996 
997  if (verbose)
998  {
999  std::cout << "[Navier-Stokes] Assembling the rhs ... " << std::flush;
1000  }
1001  ChronoItem.start();
1002 
1003  vector_block_type NSRhs ( M_ETuFESpace->map() | M_ETpFESpace->map() | fluxMap, Repeated );
1004  NSRhs *= 0.0;
1005 
1006  {
1007  std::shared_ptr<NormalizeFct> normalize (new NormalizeFct);
1008 
1009  using namespace ExpressionAssembly;
1010 
1011  integrate (
1012  elements (M_ETuFESpace->mesh() ), // Mesh
1013 
1014  M_uFESpace->qr(), // QR
1015 
1016  M_ETuFESpace,
1017 
1018  //Inertial Term
1019  M_d->density * dot (value (M_ETuFESpace, velocityBdfRHS), phi_i )
1020 
1021  )
1022  >> NSRhs.block (0);
1023 
1024  if (M_d->stabilization == "sup")
1025  {
1026 
1027  integrate (
1028  elements (M_ETuFESpace->mesh() ), // Mesh
1029 
1030  M_uFESpace->qr(), // QR
1031 
1032  M_ETuFESpace,
1033 
1034  // SUPG
1035  dot (
1036  M_d->density * value (M_ETuFESpace, velocityBdfRHS)
1037  , SUPG_TEST)
1038 
1039  )
1040  >> NSRhs.block (0);
1041 
1042  integrate (
1043  elements (M_ETuFESpace->mesh() ), // Mesh
1044 
1045  M_uFESpace->qr(), // QR
1046 
1047  M_ETpFESpace,
1048 
1049 
1050  // PSPG
1051  dot (
1052  M_d->density * value (M_ETuFESpace, velocityBdfRHS)
1053  , PSPG_TEST)
1054 
1055 
1056  )
1057  >> NSRhs.block (1);
1058  }
1059 
1060  }
1061 
1062  ChronoItem.stop();
1063  if (verbose)
1064  {
1065  std::cout << ChronoItem.diff() << " s" << std::endl;
1066  }
1067 
1068  if (verbose)
1069  {
1070  std::cout << "[Navier-Stokes] Boundary Intergrals in the rhs ... " << std::flush;
1071  }
1072  ChronoItem.start();
1073 
1074  {
1075  using namespace ExpressionAssembly;
1076 
1077  integrate ( boundary (M_ETuFESpace->mesh(), WALL),
1078  myBDQR,
1079 
1080  M_ETuFESpace,
1081 
1082  //BOUNDARY STIFFNESS
1083  value (-1) * ( dt / alpha ) *
1084  2 * LameII * value ( M_ETpFESpace, hSolid ) *
1085  0.5 * dot ( ( grad (M_ETuFESpace, dispBdfRHS ) + (-1) * grad (M_ETuFESpace, dispBdfRHS ) * outerProduct ( Nface, Nface ) )
1086  + transpose (grad (M_ETuFESpace, dispBdfRHS ) + (-1) * grad (M_ETuFESpace, dispBdfRHS ) * outerProduct ( Nface, Nface ) ),
1087  ( grad (phi_i) + ( (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) ) ) +
1088 
1089  value (-1) * ( dt / alpha ) *
1090  LameI * value ( M_ETpFESpace, hSolid ) * dot ( value ( Eye ) , ( grad (M_ETuFESpace, dispBdfRHS ) + (-1) * grad (M_ETuFESpace, dispBdfRHS ) * outerProduct ( Nface, Nface ) ) )
1091  * dot ( value ( Eye ) , ( grad (phi_i) + (-1) * grad (phi_i) * outerProduct ( Nface, Nface ) ) )
1092  ) >> NSRhs.block (0);
1093 
1094  integrate (
1095  boundary (M_ETuFESpace->mesh(), WALL), // Mesh
1096 
1097  myBDQR, // QR
1098 
1099  M_ETuFESpace,
1100 
1101  //BOUNDARY MASS
1102  ( /*M_d->Hs*/value ( M_ETpFESpace, hSolid ) * M_d->rhos * alpha) * dot (value (M_ETuFESpace, velocityBdfRHS), phi_i )
1103  - dt * dot ( value ( M_ETpFESpace, robinExt) * value (M_ETuFESpace, dispBdfRHS ) , phi_i )
1104 
1105  )
1106  >> NSRhs.block (0);
1107 
1108  }
1109 
1110  if ( M_d->transpirationOrder == 1 )
1111  {
1112  {
1113  using namespace ExpressionAssembly;
1114 
1115  integrate (
1116  boundary (M_ETuFESpace->mesh(), WALL), // Mesh
1117 
1118  myBDQR, // QR
1119 
1120  M_ETuFESpace,
1121 
1122  //BOUNDARY MASS
1123  ( /*M_d->Hs*/value ( M_ETpFESpace, hSolid) * M_d->rhos * alpha) * dot (grad (M_ETuFESpace, velocityBdfRHS) * value ( M_ETuFESpace, dispExtrapolated ), phi_i )
1124 
1125 
1126  )
1127  >> NSRhs.block (0);
1128  }
1129  }
1130 
1131  ChronoItem.stop();
1132  if (verbose)
1133  {
1134  std::cout << ChronoItem.diff() << " s" << std::endl;
1135  }
1136 
1137  if (verbose)
1138  {
1139  std::cout << "[Navier-Stokes] Closing the matrix and the rhs ... " << std::flush;
1140  }
1141 
1142  ChronoItem.start();
1143 
1144  NSMatrix->globalAssemble();
1145 
1146  NSRhs.globalAssemble();
1147  vector_block_type NSRhsUnique ( NSRhs, Unique );
1148  //NSRhsUnique += (-1) * steadyResidual;
1149  //steadyResidual *= 0.0;
1150 
1151  ChronoItem.stop();
1152  if (verbose)
1153  {
1154  std::cout << ChronoItem.diff() << " s" << std::endl;
1155  }
1156 
1157 
1158  if (verbose)
1159  {
1160  std::cout << "[Navier-Stokes] Applying boundary conditions ... " << std::flush;
1161  }
1162 
1163  bcHFluid.bcUpdate ( *meshPart.meshPartition(), M_uFESpace->feBd(), M_uFESpace->dof() );
1164 
1165  std::shared_ptr<matrix_type> NSMatrixNoBlock (new matrix_type ( fullMap ) );
1166  *NSMatrixNoBlock += *NSMatrix;
1167  NSMatrixNoBlock->globalAssemble();
1168  bcManage (*NSMatrixNoBlock, NSRhsUnique,
1169  *M_uFESpace->mesh(), M_uFESpace->dof(),
1170  bcHFluid, M_uFESpace->feBd(), 1.0, currentTime);
1171 
1172  ChronoItem.stop();
1173  if (verbose)
1174  {
1175  std::cout << ChronoItem.diff() << " s" << std::endl;
1176  }
1177 
1178  if (verbose)
1179  {
1180  std::cout << "[Navier-Stokes] Solving the system " << std::endl;
1181  }
1182 
1183  NSSolver.setMatrix (*NSMatrixNoBlock);
1184 
1185  NSSolver.solveSystem (NSRhsUnique, NSSolution, NSMatrixNoBlock);
1186 
1187  fluidTimeAdvance.bdfVelocity().shiftRight ( NSSolution );
1188 
1189  velocitySolution.subset (NSSolution);
1190 
1191  vector_type dispInterface (M_interfaceMap, Repeated);
1192 
1193  dispInterface = dt / alpha * ( velocitySolution + dispBdfRHS );
1194 
1195  dispSolution *= 0.0;
1196 
1197  dispSolution.subset (dispInterface, *M_interfaceMap, 0, 0);
1198 
1199  dispTimeAdvance.shiftRight ( dispSolution );
1200 
1201  fluidTimeAdvance.bdfVelocity().updateRHSContribution ( oseenData->dataTime()->timeStep() );
1202  dispTimeAdvance.updateRHSContribution (oseenData->dataTime()->timeStep() );
1203 
1204  *velAndPressureExporter = NSSolution;
1205  *dispExporter = dispSolution;
1206 
1207  //dispExporter->spy ("dispExporter");
1208 
1209  if (niter % exportEach == 0)
1210  {
1211  exporter.postProcess (currentTime);
1212  }
1213 
1214  ChronoIteration.stop();
1215  if (verbose)
1216  {
1217  std::cout << std::endl << " Total iteration time : " << ChronoIteration.diff() << " s" << std::endl;
1218  }
1219 
1220  } // end time loop
1221 
1222 
1223  exporter.closeFile();
1224 
1225 }
1226 
1227 
1228 void ETRobinMembraneSolver::createInterfaceMap ( std::map<ID, ID> const& locDofMap, const DOF& dof )
1229 {
1230  Displayer disp (M_d->comm);
1231  disp.leaderPrint ("Building the Interface Map ... ");
1232 
1233  std::vector<int> dofInterfaceFluid;
1234 
1235  typedef std::map<ID, ID>::const_iterator iterator_Type;
1236 
1237  //std::map<ID, ID> const& locDofMap = M_dofStructureToHarmonicExtension->locDofMap();
1238  dofInterfaceFluid.reserve ( locDofMap.size() );
1239 
1240  for (UInt dim = 0; dim < nDimensions; ++dim)
1241  for ( iterator_Type i = locDofMap.begin(); i != locDofMap.end(); ++i )
1242  {
1243  dofInterfaceFluid.push_back (i->second + dim * dof.numTotalDof() ); // in solid numerotation
1244  }
1245 
1246  int* pointerToDofs (0);
1247  if (dofInterfaceFluid.size() > 0)
1248  {
1249  pointerToDofs = &dofInterfaceFluid[0];
1250  }
1251 
1252  M_interfaceMap.reset ( new MapEpetra ( -1,
1253  static_cast<int> (dofInterfaceFluid.size() ),
1254  pointerToDofs,
1255  M_d->comm ) );
1256  disp.leaderPrint ("done\n");
1257  M_d->comm->Barrier();
1258 
1259 }
TRANSP_GRADGRAD
#define TRANSP_GRADGRAD
LifeV::operator*
VectorEpetra operator*(const VectorEpetra::data_type &scalar, const VectorEpetra &vector)
Definition: VectorEpetra.cpp:1135
GetPot::operator()
double operator()(const char *VarName, const double &Default) const
Definition: GetPot.hpp:2033
LifeV::VectorEpetra::VectorEpetra
VectorEpetra(const VectorEpetra &vector, const MapEpetraType &mapType)
Copy constructor.
Definition: VectorEpetra.cpp:94
LifeV::ETRobinMembraneSolver::mesh_type
RegionMesh< LinearTetra > mesh_type
Definition: ETRobinMembraneSolver.hpp:89
LifeV::quadRuleTria4pt
const QuadratureRule quadRuleTria4pt(pt_tria_4pt, 3, "Quadrature rule 4 points on a triangle", TRIANGLE, 4, 3)
LifeV::ETRobinMembraneSolver::vectorPtr_type
std::shared_ptr< vector_type > vectorPtr_type
Definition: ETRobinMembraneSolver.hpp:98
GetPot::GetPot
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:507
DIVDIV_TEST
#define DIVDIV_TEST
LifeV::ETRobinMembraneSolver::Private::D
Real D
diameter of the cylinder (in m)
Definition: ETRobinMembraneSolver.cpp:92
SUPG_TEST
#define SUPG_TEST
LifeV::VectorEpetra::VectorEpetra
VectorEpetra(const MapEpetra &map, const MapEpetraType &mapType=Unique, const combineMode_Type combineMode=Add)
Constructor - Using Maps.
Definition: VectorEpetra.cpp:61
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::MapEpetra::importer
Epetra_Import const & importer()
Getter for the Epetra_Import.
Definition: MapEpetra.cpp:394
ASSERT
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
LifeV::VectorEpetra::operator+
const VectorEpetra operator+(const VectorEpetra &vector) const
Addition operator.
Definition: VectorEpetra.cpp:365
LifeV::MeshPartitioner::M_isOnProc
std::shared_ptr< std::vector< Int > > M_isOnProc
Definition: MeshPartitioner.hpp:388
LifeV::TimeAdvanceBDF
class TimeAdvanceBDF - Backward differencing formula time discretization for the first and the second...
Definition: TimeAdvanceBDF.hpp:140
LifeV::MeshData
MeshData - class for handling spatial discretization.
Definition: MeshData.hpp:72
PSPG_TEST
#define PSPG_TEST
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::ETRobinMembraneSolver::ETRobinMembraneSolver
ETRobinMembraneSolver(int argc, char **argv)
Definition: ETRobinMembraneSolver.cpp:114
GetPot::operator()
int operator()(const char *VarName, int Default) const
Definition: GetPot.hpp:2021
LifeV::ETRobinMembraneSolver::vector_block_type
VectorEpetraStructured vector_block_type
Definition: ETRobinMembraneSolver.hpp:93
LifeV::ETRobinMembraneSolver::Private::fct_type
std::function< Real(Real const &, Real const &, Real const &, Real const &, ID const &) > fct_type
Definition: ETRobinMembraneSolver.cpp:84
LifeV::VectorEpetra::operator*=
VectorEpetra & operator*=(const data_type &scalar)
Multiplication operator.
Definition: VectorEpetra.cpp:426
PI
const Real PI
Definition: ETRobinMembraneSolver.cpp:79
LifeV::VectorEpetra::operator=
VectorEpetra & operator=(const VectorEpetra &vector)
Affectation operator.
Definition: VectorEpetra.cpp:216
LifeV::ETRobinMembraneSolver::Private::H
Real H
height and width of the domain (in m)
Definition: ETRobinMembraneSolver.cpp:91
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191