LifeV
lifev/structure/testsuite/anisotropicLaw/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  \file main.cpp
28 
29  This test is the case of traction of a cube. It does not use the symmetry BCs
30  This test uses the FESpace which is the standard in LifeV and the ETFESpace
31  The FESpace is used for the BCs of Neumann type since in the ET branch there
32  is not the integration over the boundary faces.
33 
34  \author Paolo Tricerri <paolo.tricerri@epfl.ch>
35  \date 1861-03-17
36  */
37 
38 #ifdef TWODIM
39 #error test_structure cannot be compiled in 2D
40 #endif
41 
42 // Tell the compiler to ignore specific kind of warnings:
43 #pragma GCC diagnostic ignored "-Wunused-variable"
44 #pragma GCC diagnostic ignored "-Wunused-parameter"
45 
46 #include <Epetra_ConfigDefs.h>
47 #ifdef EPETRA_MPI
48 #include <mpi.h>
49 #include <Epetra_MpiComm.h>
50 #else
51 #include <Epetra_SerialComm.h>
52 #endif
53 
54 //Tell the compiler to restore the warning previously silented
55 #pragma GCC diagnostic warning "-Wunused-variable"
56 #pragma GCC diagnostic warning "-Wunused-parameter"
57 
58 #include <lifev/core/LifeV.hpp>
59 #include <lifev/core/algorithm/PreconditionerIfpack.hpp>
60 #include <lifev/core/algorithm/PreconditionerML.hpp>
61 
62 #include <lifev/core/array/MapEpetra.hpp>
63 
64 #include <lifev/core/fem/TimeAdvance.hpp>
65 #include <lifev/core/fem/TimeAdvanceBDF.hpp>
66 
67 #include <lifev/core/mesh/MeshData.hpp>
68 #include <lifev/core/mesh/MeshPartitioner.hpp>
69 
70 #include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>
71 
72 #include <lifev/structure/solver/StructuralConstitutiveLaw.hpp>
73 #include <lifev/structure/solver/StructuralOperator.hpp>
74 #include <lifev/structure/solver/isotropic/VenantKirchhoffMaterialLinear.hpp>
75 #include <lifev/structure/solver/isotropic/VenantKirchhoffMaterialNonLinear.hpp>
76 #include <lifev/structure/solver/isotropic/ExponentialMaterialNonLinear.hpp>
77 #include <lifev/structure/solver/isotropic/VenantKirchhoffMaterialNonLinearPenalized.hpp>
78 #include <lifev/structure/solver/isotropic/SecondOrderExponentialMaterialNonLinear.hpp>
79 #include <lifev/structure/solver/isotropic/NeoHookeanMaterialNonLinear.hpp>
80 
81 
82 #include <lifev/structure/solver/anisotropic/StructuralAnisotropicConstitutiveLaw.hpp>
83 #include <lifev/structure/solver/anisotropic/HolzapfelMaterialNonLinear.hpp>
84 #include <lifev/structure/solver/anisotropic/DistributedHolzapfelMaterialNonLinear.hpp>
85 
86 
87 #include <lifev/core/filter/ExporterEnsight.hpp>
88 #ifdef HAVE_HDF5
89 #include <lifev/core/filter/ExporterHDF5.hpp>
90 #endif
91 #include <lifev/core/filter/ExporterEmpty.hpp>
92 
93 //Includes for the Expression Template
94 #include <lifev/eta/fem/ETFESpace.hpp>
95 
96 #include <iostream>
97 #include <fstream>
98 
99 #include "ud_functions.hpp"
100 
101 
102 using namespace LifeV;
103 
104 int returnValue = EXIT_FAILURE;
106 
107 namespace
108 {
109 static bool regIF = (PRECFactory::instance().registerProduct ( "Ifpack", &createIfpack ) );
110 static bool regML = (PRECFactory::instance().registerProduct ( "ML", &createML ) );
111 }
112 
113 
115 {
116  std::string stringList = list;
117  std::set<UInt> setList;
118  if ( list == "" )
119  {
120  return setList;
121  }
122  size_t commaPos = 0;
123  while ( commaPos != std::string::npos )
124  {
125  commaPos = stringList.find ( "," );
126  setList.insert ( atoi ( stringList.substr ( 0, commaPos ).c_str() ) );
127  stringList = stringList.substr ( commaPos + 1 );
128  }
129  setList.insert ( atoi ( stringList.c_str() ) );
130  return setList;
131 }
132 
133 
134 class Structure
135 {
136 public:
137 
138  // Public typedefs
144 
147 
148  // typedefs for fibers
149  // std function for fiber direction
150  typedef std::function<Real ( Real const&, Real const&, Real const&, Real const&, ID const& ) > fiberFunction_Type;
152 
155 
157 
158  /** @name Constructors, destructor
159  */
160  //@{
161  Structure ( int argc,
162  char** argv,
163  std::shared_ptr<Epetra_Comm> structComm );
164 
166  {}
167  //@}
168 
169  //@{
170  void run()
171  {
172  run3d();
173  }
174  void CheckResultHolzapfelModel (const Real& dispNorm, const Real& time);
175  void CheckResultDistributedModel (const Real& dispNorm, const Real& time);
176  void resultChanged (Real time);
177  //@}
178 
179 protected:
180 
181 private:
182 
183  /**
184  * run the 2D driven cylinder simulation
185  */
186  void run2d();
187 
188  /**
189  * run the 3D driven cylinder simulation
190  */
191  void run3d();
192 
193 private:
194  struct Private;
196 };
197 
198 
199 
200 struct Structure::Private
201 {
203  rho (1), poisson (1), young (1), bulk (1), alpha (1), gamma (1)
204  {}
205  double rho, poisson, young, bulk, alpha, gamma;
206 
208 
210 };
211 
212 
213 
214 Structure::Structure ( int argc,
215  char** argv,
216  std::shared_ptr<Epetra_Comm> structComm) :
217  parameters ( new Private() )
218 {
219  GetPot command_line (argc, argv);
220  std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
221  GetPot dataFile ( data_file_name );
222  parameters->data_file_name = data_file_name;
223 
224  parameters->comm = structComm;
225  int ntasks = parameters->comm->NumProc();
226 
227  if (!parameters->comm->MyPID() )
228  {
229  std::cout << "My PID = " << parameters->comm->MyPID() << " out of " << ntasks << " running." << std::endl;
230  }
231 }
232 
233 
234 
235 void
236 Structure::run2d()
237 {
238  std::cout << "2D cylinder test case is not available yet\n";
239 }
240 
241 
242 
243 void
244 Structure::run3d()
245 {
246  bool verbose = (parameters->comm->MyPID() == 0);
247 
248  //! Number of boundary conditions for the velocity and mesh motion
249  std::shared_ptr<BCHandler> BCh ( new BCHandler() );
250 
251  //! dataElasticStructure
252  GetPot dataFile ( parameters->data_file_name.c_str() );
253 
254  std::shared_ptr<StructuralConstitutiveLawData> dataStructure (new StructuralConstitutiveLawData( ) );
255  dataStructure->setup (dataFile);
256 
257  dataStructure->showMe();
258 
259  MeshData meshData;
260  meshData.setup (dataFile, "solid/space_discretization");
261 
262  std::shared_ptr<RegionMesh<LinearTetra> > fullMeshPtr (new RegionMesh<LinearTetra> ( parameters->comm ) );
263  std::shared_ptr<RegionMesh<LinearTetra> > localMeshPtr (new RegionMesh<LinearTetra> ( parameters->comm ) );
264  readMesh (*fullMeshPtr, meshData);
265 
266  MeshPartitioner< RegionMesh<LinearTetra> > meshPart ( fullMeshPtr, parameters->comm );
267  localMeshPtr = meshPart.meshPartition();
268 
269  std::string dOrder = dataFile ( "solid/space_discretization/order", "P1");
270 
271  //Mainly used for BCs assembling (Neumann type)
272  solidFESpacePtr_Type dFESpace ( new solidFESpace_Type (localMeshPtr, dOrder, 3, parameters->comm) );
273  solidETFESpacePtr_Type dETFESpace ( new solidETFESpace_Type (meshPart, & (dFESpace->refFE() ), & (dFESpace->fe().geoMap() ), parameters->comm) );
274 
275  // Setting the fibers
276  vectorFiberFunctionPtr_Type pointerToVectorOfFamilies( new vectorFiberFunction_Type( ) );
277  (*pointerToVectorOfFamilies).resize( dataStructure->numberFibersFamilies( ) );
278 
279  listOfFiberDirections_Type fiberDirections;
280  fiberDirections.resize( dataStructure->numberFibersFamilies( ) );
281 
282  std::cout << "Size of the number of families: " << (*pointerToVectorOfFamilies).size() << std::endl;
283 
284  fibersDirectionList setOfFiberFunctions;
285  setOfFiberFunctions.setupFiberDefinitions( dataStructure->numberFibersFamilies( ) );
286 
287  if (verbose)
288  {
289  std::cout << std::endl;
290  }
291 
292  std::string timeAdvanceMethod = dataFile ( "solid/time_discretization/method", "BDF");
293 
294  timeAdvance_Type timeAdvance ( TimeAdvanceFactory::instance().createObject ( timeAdvanceMethod ) );
295 
296  UInt OrderDev = 2;
297 
298  timeAdvance->setup (dataStructure->dataTimeAdvance()->orderBDF() , OrderDev);
299 
300  timeAdvance->setTimeStep (dataStructure->dataTime()->timeStep() );
301  //timeAdvance->showMe();
302 
303  //! #################################################################################
304  //! BOUNDARY CONDITIONS
305  //! #################################################################################
306  std::vector <ID> compx (1), compy (1), compz (1), compxy (2), compxz (2), compyz (2);
307  compx[0] = 0;
308  compy[0] = 1, compz[0] = 2;
309  compxy[0] = 0;
310  compxy[1] = 1;
311  compxz[0] = 0;
312  compxz[1] = 2;
313  compyz[0] = 1;
314  compyz[1] = 2;
315 
316  BCFunctionBase zero (bcZero);
317  BCFunctionBase nonZero;
318 
319  nonZero.setFunction (bcNonZero);
320 
321  //! =================================================================================
322  //! BC for StructuredCube4_test_structuralsolver.mesh
323  //! =================================================================================
324  BCh->addBC ("EdgesIn", 20, Natural, Component, nonZero, compy);
325  BCh->addBC ("EdgesIn", 40, Essential, Component, zero, compy);
326 
327  //! Symmetry BC
328  BCh->addBC ("EdgesIn", 50, EssentialVertices, Component, zero, compxy);
329  BCh->addBC ("EdgesIn", 30, EssentialVertices, Component, zero, compyz);
330  BCh->addBC ("EdgesIn", 80, EssentialVertices, Component, zero, compxz);
331  BCh->addBC ("EdgesIn", 100, EssentialVertices, Full, zero, 3);
332 
333  BCh->addBC ("EdgesIn", 7, Essential, Component , zero, compx);
334  BCh->addBC ("EdgesIn", 3, Essential, Component , zero, compz);
335  //! =================================================================================
336 
337  //! 1. Constructor of the structuralSolver
338  StructuralOperator< RegionMesh<LinearTetra> > solid;
339 
340  //! 2. Setup of the structuralSolver
341  solid.setup (dataStructure,
342  dFESpace,
343  dETFESpace,
344  BCh,
345  parameters->comm);
346 
347  //! 3. Setting data from getPot
348  solid.setDataFromGetPot (dataFile);
349 
350  // Setting the vector of fibers functions
351  for( UInt k(1); k <= pointerToVectorOfFamilies->size( ); k++ )
352  {
353  // Setting up the name of the function to define the family
354  std::string family="Family";
355  // adding the number of the family
356  std::string familyNumber;
357  std::ostringstream number;
358  number << ( k );
359  familyNumber = number.str();
360 
361  // Name of the function to create
362  std::string creationString = family + familyNumber;
363  (*pointerToVectorOfFamilies)[ k-1 ].reset( new fiberFunction_Type() );
364  (*pointerToVectorOfFamilies)[ k-1 ] = setOfFiberFunctions.fiberDefinition( creationString );
365 
366  fiberDirections[ k-1 ].reset( new vector_Type(solid.displacement(), Unique) );
367  }
368 
369  //! 3.b Setting the fibers in the abstract class of Anisotropic materials
370  solid.material()->anisotropicLaw()->setupFiberDirections( pointerToVectorOfFamilies );
371 
372  //! 4. Building system using TimeAdvance class
373  double timeAdvanceCoefficient = timeAdvance->coefficientSecondDerivative ( 0 ) / (dataStructure->dataTime()->timeStep() * dataStructure->dataTime()->timeStep() );
374  solid.buildSystem (timeAdvanceCoefficient);
375 
376 
377  //dataStructure->showMe();
378  //! =================================================================================
379  //! Temporal data and initial conditions
380  //! =================================================================================
381 
382  //! 5. Initial data
383  Real dt = dataStructure->dataTime()->timeStep();
384 
385  vectorPtr_Type rhs (new vector_Type (solid.displacement(), Unique) );
386  vectorPtr_Type disp (new vector_Type (solid.displacement(), Unique) );
387  vectorPtr_Type vel (new vector_Type (solid.displacement(), Unique) );
388  vectorPtr_Type acc (new vector_Type (solid.displacement(), Unique) );
389 
390  if (verbose)
391  {
392  std::cout << "S- initialization ... ";
393  }
394 
395  std::vector<vectorPtr_Type> uv0;
396 
397  if (timeAdvanceMethod == "Newmark")
398  {
399  uv0.push_back (disp);
400  uv0.push_back (vel);
401  uv0.push_back (acc);
402  }
403 
404  vectorPtr_Type initialDisplacement (new vector_Type (solid.displacement(), Unique) );
405 
406  if ( !dataStructure->solidTypeIsotropic().compare ("secondOrderExponential") )
407  {
408  dFESpace->interpolate ( static_cast<solidFESpace_Type::function_Type> ( d0 ), *initialDisplacement, 0.0 );
409  }
410 
411  if (timeAdvanceMethod == "BDF")
412  {
413  Real tZero = dataStructure->dataTime()->initialTime();
414 
415  for ( UInt previousPass = 0; previousPass < timeAdvance->size() ; previousPass++)
416  {
417  Real previousTimeStep = tZero - previousPass * dt;
418  std::cout << "BDF " << previousTimeStep << "\n";
419  if ( !dataStructure->solidTypeIsotropic().compare ("secondOrderExponential") )
420  {
421  uv0.push_back (initialDisplacement);
422  }
423  else
424  {
425  uv0.push_back (disp);
426  }
427  }
428  }
429 
430  timeAdvance->setInitialCondition (uv0);
431 
432  timeAdvance->setTimeStep ( dt );
433 
434  timeAdvance->updateRHSContribution ( dt );
435 
436  if ( !dataStructure->solidTypeIsotropic().compare ("secondOrderExponential") )
437  {
438  solid.initialize ( initialDisplacement );
439  }
440  else
441  {
442  solid.initialize ( disp );
443  }
444 
445  MPI_Barrier (MPI_COMM_WORLD);
446 
447  if (verbose )
448  {
449  std::cout << "ok." << std::endl;
450  }
451 
452  std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporter;
453 
454  std::string const exporterType = dataFile ( "exporter/type", "ensight");
455 #ifdef HAVE_HDF5
456  if (exporterType.compare ("hdf5") == 0)
457  {
458  exporter.reset ( new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, "structure" ) );
459  }
460  else
461 #endif
462  {
463  if (exporterType.compare ("none") == 0)
464  {
465  exporter.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(), "structure", parameters->comm->MyPID() ) );
466  }
467 
468  else
469  {
470  exporter.reset ( new ExporterEnsight<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(), "structure", parameters->comm->MyPID() ) );
471  }
472  }
473 
474  exporter->setPostDir ( "./" );
475  exporter->setMeshProcId ( meshPart.meshPartition(), parameters->comm->MyPID() );
476 
477  vectorPtr_Type solidDisp ( new vector_Type (solid.displacement(), exporter->mapType() ) );
478  vectorPtr_Type solidVel ( new vector_Type (solid.displacement(), exporter->mapType() ) );
479  vectorPtr_Type solidAcc ( new vector_Type (solid.displacement(), exporter->mapType() ) );
480  //vectorPtr_Type rhsVector ( new vector_Type (solid.displacement(), exporter->mapType() ) );
481 
482  exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "displacement", dFESpace, solidDisp, UInt (0) );
483  exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "velocity", dFESpace, solidVel, UInt (0) );
484  exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "acceleration", dFESpace, solidAcc, UInt (0) );
485  //exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "forcingTerm", dFESpace, rhsVector, UInt (0) );
486 
487  // Adding the fibers vectors
488  // Setting the vector of fibers functions
489  for( UInt k(1); k <= pointerToVectorOfFamilies->size( ); k++ )
490  {
491  // Setting up the name of the function to define the family
492  std::string family="Family-";
493  // adding the number of the family
494  std::string familyNumber;
495  std::ostringstream number;
496  number << ( k );
497  familyNumber = number.str();
498 
499  // Name of the function to create
500  std::string creationString = family + familyNumber;
501  exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, creationString, dFESpace, fiberDirections[ k-1 ], UInt (0) );
502 
503  // Extracting the fibers vectors
504  *(fiberDirections[ k-1 ]) = solid.material()->anisotropicLaw()->ithFiberVector( k );
505  }
506 
507 
508  exporter->postProcess ( 0 );
509  std::cout.precision(16);
510 
511  // int IDPointX = 618;
512  // int IDPointY = 331;
513 
514  /*
515  //!--------------------------------------------------------------------------------------------
516  //! MATLAB FILE WITH DISPLACEMENT OF A CHOSEN POINT
517  //!--------------------------------------------------------------------------------------------
518 
519  ofstream file_comp( "Displacement_components_NL.m" );
520  if ( !file_comp )
521  {
522  std::cout <<" Unable to open file! You need to specify the output folder in the data file " << std::endl;
523  }
524 
525 
526  //int IDPoint = 401; // StructuredCube8
527  //int IDPoint = 2593; // StructuredCube16
528 
529  //int IDPoint = 74;// cube4x4.mesh
530  //int IDPoint = 315;// cube8x8.mesh
531  //int IDPoint = 1526;// cube16x16.mesh
532 
533  file_comp << " % TEST NONLINEAR ELASTICITY" << endl;
534  file_comp << " % Displacement components of ID node " << IDPoint << " :" << endl;
535  file_comp << " % Each row is a time step" << endl;
536  file_comp << " % First column = comp x, second = comp y and third = comp z. " << endl;
537  file_comp << endl;
538  file_comp << " SolidDisp_NL = [ " ;
539 
540  for ( UInt k = IDPoint - 1; k <= solid.displacement().size() - 1; k = k + solid.displacement().size()/nDimensions )
541  {
542  file_comp<< solid.displacement()[ k ] << " ";
543  }
544 
545  file_comp<< endl;
546  */
547  //!--------------------------------------------------------------------------------------------
548  //!The update of the RHS is done by the TimeAdvance class
549  //solid.updateSystem();
550  //! =================================================================================
551 
552  Real normVect;
553  normVect = solid.displacement().norm2();
554 
555  if( verbose )
556  {
557  std::cout << "The norm 2 of the displacement field is: " << normVect << std::endl;
558  }
559 
560  //! =============================================================================
561  //! Temporal loop
562  //! =============================================================================
563  // for (Real time = dt; time <= T; time += dt)
564  for (dataStructure->dataTime()->setTime ( dt ) ; dataStructure->dataTime()->canAdvance( ); dataStructure->dataTime()->updateTime( ) )
565  {
566  returnValue = EXIT_FAILURE;
567 
568  if (verbose)
569  {
570  std::cout << std::endl;
571  std::cout << "S- Now we are at time " << dataStructure->dataTime()->time() << " s." << std::endl;
572  }
573 
574  //! 6. Updating right-hand side
575  *rhs *= 0;
576  timeAdvance->updateRHSContribution ( dt );
577  *rhs += *solid.massMatrix() * timeAdvance->rhsContributionSecondDerivative() / timeAdvanceCoefficient;
578  solid.setRightHandSide ( *rhs );
579 
580  //! 7. Iterate --> Calling Newton
581  solid.iterate ( BCh );
582 
583  timeAdvance->shiftRight ( solid.displacement() );
584 
585  *solidDisp = solid.displacement();
586  *solidVel = timeAdvance->firstDerivative();
587  *solidAcc = timeAdvance->secondDerivative();
588 
589  // This vector is to export the forcing term
590  //*rhsVector = solid.rhsCopy();
591 
592  exporter->postProcess ( dataStructure->dataTime()->time() );
593 
594  /* This part lets to save the displacement at one point of the mesh and to check the result
595  w.r.t. manufactured solution.
596  //!--------------------------------------------------------------------------------------------------
597  //! MATLAB FILE WITH DISPLACEMENT OF A CHOOSEN POINT
598  //!--------------------------------------------------------------------------------------------------
599  cout <<"******* DISPLACEMENT COMPONENTS of ID node "<< IDPoint << " *******"<< std::endl;
600  for ( UInt k = IDPoint - 1 ; k <= solid.displacement().size() - 1; k = k + solid.displacement().size()/nDimensions )
601  {
602  file_comp<< solid.displacement()[ k ] << " ";
603  cout.precision(16);
604  cout <<"*********************************************************"<< std::endl;
605  cout <<" solid.disp()[ "<< k <<" ] = "<< solid.displacement()[ k ] << std::endl;
606  cout <<"*********************************************************"<< std::endl;
607  }
608  file_comp<< endl;
609 
610  cout <<"*********************************************************"<< std::endl;
611  cout <<" solid.disp()[ "<< IDPointX - 1 <<" ] = "<< solid.displacement()[ IDPointX - 1 ] << std::endl;
612  cout <<" solid.disp()[ "<< IDPointX - 1 + dFESpace->dof().numTotalDof() <<" ] = "<< solid.displacement()[ IDPointX - 1 + dFESpace->dof().numTotalDof() ] << std::endl;
613  cout <<"*********************************************************"<< std::endl;
614  */
615 
616  Real normVect;
617  normVect = solid.displacement().norm2();
618  std::cout << "The norm 2 of the displacement field is: " << normVect << std::endl;
619 
620  // Check results
621  if ( !dataStructure->solidTypeAnisotropic().compare ("holzapfel") )
622  CheckResultHolzapfelModel (normVect, dataStructure->dataTime()->time() );
623  else
624  CheckResultDistributedModel (normVect, dataStructure->dataTime()->time() );
625 
626  //!--------------------------------------------------------------------------------------------------
627 
628  MPI_Barrier (MPI_COMM_WORLD);
629  }
630 
631 }
632 
633 void Structure::CheckResultHolzapfelModel (const Real& dispNorm, const Real& time)
634 {
635  if ( time == 0.05 && std::fabs (dispNorm - 0.84960668) <= 1e-7 )
636  {
637  this->resultChanged (time);
638  }
639  if ( time == 0.1 && std::fabs (dispNorm - 0.84981715) <= 1e-7 )
640  {
641  this->resultChanged (time);
642  }
643 }
644 
645 void Structure::CheckResultDistributedModel (const Real& dispNorm, const Real& time)
646 {
647  if ( time == 0.05 && std::fabs (dispNorm - 1.7747561302) <= 1e-7 )
648  {
649  this->resultChanged (time);
650  }
651  if ( time == 0.1 && std::fabs (dispNorm - 1.7757263996) <= 1e-7 )
652  {
653  this->resultChanged (time);
654  }
655 }
656 
657 void Structure::resultChanged (Real time)
658 {
659  std::cout << "Correct value at time: " << time << std::endl;
660  returnValue = EXIT_SUCCESS;
661 }
662 
663 
664 
665 int
666 main ( int argc, char** argv )
667 {
668 
669 #ifdef HAVE_MPI
670  MPI_Init (&argc, &argv);
671  std::shared_ptr<Epetra_MpiComm> Comm (new Epetra_MpiComm ( MPI_COMM_WORLD ) );
672  if ( Comm->MyPID() == 0 )
673  {
674  std::cout << "% using MPI" << std::endl;
675  }
676 #else
677  std::shared_ptr<Epetra_SerialComm> Comm ( new Epetra_SerialComm() );
678  std::cout << "% using serial Version" << std::endl;
679 #endif
680 
681  Structure structure ( argc, argv, Comm );
682  structure.run();
683 
684 #ifdef HAVE_MPI
685  MPI_Finalize();
686 #endif
687  return returnValue ;
688 }
void assignFunction(bcBase_Type &base)
Assign the function to the base of the BCHandler.
Structure(int argc, char **argv, std::shared_ptr< Epetra_Comm > structComm)
int main(int argc, char **argv)
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
Definition: GetPot.hpp:507
std::shared_ptr< fiberFunction_Type > fiberFunctionPtr_Type
std::shared_ptr< hdf5Filter_Type > hdf5FilterPtr_Type
FESpace< RegionMesh< LinearTetra >, MapEpetra > solidFESpace_Type
void run3d()
run the 3D driven cylinder simulation
ETFESpace< RegionMesh< LinearTetra >, MapEpetra, 3, 3 > solidETFESpace_Type
void run2d()
run the 2D driven cylinder simulation
MeshData - class for handling spatial discretization.
Definition: MeshData.hpp:72
void CheckResultDistributedModel(const Real &dispNorm, const Real &time)
std::shared_ptr< vector_Type > vectorPtr_Type
double Real
Generic real data.
Definition: LifeV.hpp:175
std::set< UInt > parseList(const std::string &list)
std::shared_ptr< solidETFESpace_Type > solidETFESpacePtr_Type
std::shared_ptr< TimeAdvance< vector_Type > > timeAdvance_Type
std::vector< vectorPtr_Type > listOfFiberDirections_Type
void resultChanged(Real time)
StructuralOperator< RegionMesh< LinearTetra > >::vector_Type vector_Type
void CheckResultHolzapfelModel(const Real &dispNorm, const Real &time)
std::shared_ptr< solidFESpace_Type > solidFESpacePtr_Type
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191
std::shared_ptr< vectorFiberFunction_Type > vectorFiberFunctionPtr_Type
std::function< Real(Real const &, Real const &, Real const &, Real const &, ID const &) > fiberFunction_Type
std::vector< fiberFunctionPtr_Type > vectorFiberFunction_Type