39 #error test_structure cannot be compiled in 2D
43 #pragma GCC diagnostic ignored "-Wunused-variable" 44 #pragma GCC diagnostic ignored "-Wunused-parameter" 46 #include <Epetra_ConfigDefs.h> 49 #include <Epetra_MpiComm.h> 51 #include <Epetra_SerialComm.h> 55 #pragma GCC diagnostic warning "-Wunused-variable" 56 #pragma GCC diagnostic warning "-Wunused-parameter" 58 #include <lifev/core/LifeV.hpp> 59 #include <lifev/core/algorithm/PreconditionerIfpack.hpp> 60 #include <lifev/core/algorithm/PreconditionerML.hpp> 62 #include <lifev/core/array/MapEpetra.hpp> 64 #include <lifev/core/fem/TimeAdvance.hpp> 65 #include <lifev/core/fem/TimeAdvanceBDF.hpp> 67 #include <lifev/core/mesh/MeshData.hpp> 68 #include <lifev/core/mesh/MeshPartitioner.hpp> 70 #include <lifev/structure/solver/StructuralConstitutiveLawData.hpp> 72 #include <lifev/structure/solver/StructuralConstitutiveLaw.hpp> 73 #include <lifev/structure/solver/StructuralOperator.hpp> 74 #include <lifev/structure/solver/isotropic/NeoHookeanMaterialNonLinear.hpp> 75 #include <lifev/structure/solver/anisotropic/StructuralAnisotropicConstitutiveLaw.hpp> 76 #include <lifev/structure/solver/anisotropic/HolzapfelMaterialNonLinear.hpp> 79 #include <lifev/core/filter/ExporterEnsight.hpp> 81 #include <lifev/core/filter/ExporterHDF5.hpp> 83 #include <lifev/core/filter/ExporterEmpty.hpp> 86 #include <lifev/eta/fem/ETFESpace.hpp> 93 using namespace LifeV;
107 std::string stringList = list;
108 std::set<UInt> setList;
114 while ( commaPos != std::string::npos )
116 commaPos = stringList.find (
"," );
117 setList.insert ( atoi ( stringList.substr ( 0, commaPos ).c_str() ) );
118 stringList = stringList.substr ( commaPos + 1 );
120 setList.insert ( atoi ( stringList.c_str() ) );
158 std::shared_ptr<Epetra_Comm> structComm );
208 std::shared_ptr<Epetra_Comm> structComm) :
212 std::string data_file_name = command_line.follow (
"data", 2,
"-f",
"--file");
213 GetPot dataFile ( data_file_name );
214 parameters->data_file_name = data_file_name;
216 parameters->comm = structComm;
217 int ntasks = parameters->comm->NumProc();
219 if (!parameters->comm->MyPID() )
221 std::cout <<
"My PID = " << parameters->comm->MyPID() <<
" out of " << ntasks <<
" running." << std::endl;
230 std::cout <<
"2D cylinder test case is not available yet\n";
238 bool verbose = (parameters->comm->MyPID() == 0);
241 std::shared_ptr<BCHandler> BCh (
new BCHandler() );
244 GetPot dataFile ( parameters->data_file_name.c_str() );
246 std::shared_ptr<StructuralConstitutiveLawData> dataStructure (
new StructuralConstitutiveLawData( ) );
247 dataStructure->setup (dataFile);
249 dataStructure->showMe();
252 meshData.setup (dataFile,
"solid/space_discretization");
254 std::shared_ptr<RegionMesh<LinearTetra> > fullMeshPtr (
new RegionMesh<LinearTetra> ( parameters->comm ) );
255 readMesh (*fullMeshPtr, meshData);
257 MeshPartitioner< RegionMesh<LinearTetra> > meshPart ( fullMeshPtr, parameters->comm );
258 std::shared_ptr<RegionMesh<LinearTetra> > localMeshPtr (
new RegionMesh<LinearTetra> ( parameters->comm ) );
259 localMeshPtr = meshPart.meshPartition();
261 std::string dOrder = dataFile (
"solid/space_discretization/order",
"P1");
264 solidFESpacePtr_Type dFESpace (
new solidFESpace_Type (localMeshPtr, dOrder, 3, parameters->comm) );
265 solidETFESpacePtr_Type dETFESpace (
new solidETFESpace_Type (meshPart, & (dFESpace->refFE() ), & (dFESpace->fe().geoMap() ), parameters->comm) );
271 if( !dataStructure->constitutiveLaw().compare(
"anisotropic") )
274 pointerToVectorOfFamilies.reset(
new vectorFiberFunction_Type( ) );
275 (*pointerToVectorOfFamilies).resize( dataStructure->numberFibersFamilies( ) );
277 fiberDirections.resize( dataStructure->numberFibersFamilies( ) );
284 setOfFiberFunctions.setupFiberDefinitions( dataStructure->numberFibersFamilies( ) );
287 std::string timeAdvanceMethod = dataFile (
"solid/time_discretization/method",
"BDF");
289 timeAdvance_Type timeAdvance ( TimeAdvanceFactory::instance().createObject ( timeAdvanceMethod ) );
293 timeAdvance->setup (dataStructure->dataTimeAdvance()->orderBDF() , OrderDev);
295 timeAdvance->setTimeStep (dataStructure->dataTime()->timeStep() );
301 std::vector <ID> compx (1), compy (1), compz (1), compxy (2), compxz (2), compyz (2);
303 compy[0] = 1, compz[0] = 2;
311 BCFunctionBase zero (bcZero);
312 BCFunctionBase nonZero;
313 BCFunctionBase displacementImposed;
315 nonZero.setFunction (bcNonZero);
316 displacementImposed.setFunction( analyticalDisplacement );
321 BCh->addBC (
"EdgesIn", 20, Natural, Component, nonZero, compy);
322 BCh->addBC (
"EdgesIn", 40, Essential, Component, zero, compy);
325 BCh->addBC (
"EdgesIn", 50, EssentialVertices, Component, zero, compxy);
326 BCh->addBC (
"EdgesIn", 30, EssentialVertices, Component, zero, compyz);
327 BCh->addBC (
"EdgesIn", 80, EssentialVertices, Component, zero, compxz);
328 BCh->addBC (
"EdgesIn", 100, EssentialVertices, Full, zero, 3);
330 BCh->addBC (
"EdgesIn", 7, Essential, Component, zero, compx);
331 BCh->addBC (
"EdgesIn", 3, Essential, Component, zero, compz);
336 StructuralOperator< RegionMesh<LinearTetra> > solid;
339 solid.setup (dataStructure,
346 solid.setDataFromGetPot (dataFile);
349 bool bodyForce = dataFile
( "solid/physics/bodyForce" , false );
352 solid.setHavingSourceTerm( bodyForce );
355 solid.setSourceTerm( bodyTerm );
358 if( !dataStructure->constitutiveLaw().compare(
"anisotropic") )
361 for( UInt k(1); k <= pointerToVectorOfFamilies->size( ); k++ )
364 std::string family=
"Family";
366 std::string familyNumber;
367 std::ostringstream number;
369 familyNumber = number.str();
372 std::string creationString = family + familyNumber;
373 (*pointerToVectorOfFamilies)[ k-1 ].reset(
new analyticalFunction_Type() );
374 (*pointerToVectorOfFamilies)[ k-1 ] = setOfFiberFunctions.fiberDefinition( creationString );
376 fiberDirections[ k-1 ].reset(
new vector_Type(solid.displacement(), Unique) );
380 solid.material()->anisotropicLaw()->setupFiberDirections( pointerToVectorOfFamilies );
384 double timeAdvanceCoefficient = timeAdvance->coefficientSecondDerivative ( 0 ) / (dataStructure->dataTime()->timeStep() * dataStructure->dataTime()->timeStep() );
385 solid.buildSystem (timeAdvanceCoefficient);
387 vectorPtr_Type analyticDispl (
new vector_Type ( dFESpace->map() ) );
389 dFESpace->interpolate (
static_cast<solidFESpace_Type::function_Type> ( analyticalDisplacement ), *analyticDispl, 0.0 );
397 Real dt = dataStructure->dataTime()->timeStep();
404 std::vector<vectorPtr_Type> uv0;
406 vectorPtr_Type initialDisplacement (
new vector_Type (solid.displacement(), Unique) );
408 if (timeAdvanceMethod ==
"BDF")
410 Real tZero = dataStructure->dataTime()->initialTime();
412 for ( UInt previousPass = 0; previousPass < timeAdvance->size() ; previousPass++)
414 Real previousTimeStep = tZero - previousPass * dt;
415 std::cout <<
"BDF " << previousTimeStep <<
"\n";
416 uv0.push_back (disp);
420 timeAdvance->setInitialCondition (uv0);
422 timeAdvance->setTimeStep ( dt );
424 timeAdvance->updateRHSContribution ( dt );
426 solid.initialize ( disp );
428 MPI_Barrier (MPI_COMM_WORLD);
430 std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporter;
431 std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporterCheck;
433 std::string
const exporterType = dataFile (
"exporter/type",
"ensight");
434 std::string
const exporterNameFile = dataFile (
"exporter/nameFile",
"structure");
435 std::string
const exporterCheckName = dataFile (
"exporter/nameCheck",
"verifyVectors");
438 if (exporterType.compare (
"hdf5") == 0)
440 exporter.reset (
new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, exporterNameFile ) );
441 exporterCheck.reset (
new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, exporterCheckName ) );
446 if (exporterType.compare (
"none") == 0)
448 exporter.reset (
new ExporterEmpty<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(),
"structure", parameters->comm->MyPID() ) );
453 exporter.reset (
new ExporterEnsight<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(),
"structure", parameters->comm->MyPID() ) );
457 exporter->setPostDir (
"./" );
458 exporter->setMeshProcId ( meshPart.meshPartition(), parameters->comm->MyPID() );
459 exporterCheck->setPostDir (
"./" );
460 exporterCheck->setMeshProcId ( meshPart.meshPartition(), parameters->comm->MyPID() );
462 vectorPtr_Type imposedSolution (
new vector_Type ( dFESpace->map() ) );
469 exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"displacement", dFESpace, solidDisp, UInt (0) );
470 exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"velocity", dFESpace, solidVel, UInt (0) );
471 exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"acceleration", dFESpace, solidAcc, UInt (0) );
474 exporterCheck->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField,
"imposedSolution", dFESpace, imposedSolution, UInt (0) );
476 *imposedSolution = *analyticDispl;
478 if( !dataStructure->constitutiveLaw().compare(
"anisotropic") )
482 for( UInt k(1); k <= pointerToVectorOfFamilies->size( ); k++ )
485 std::string family=
"Family-";
487 std::string familyNumber;
488 std::ostringstream number;
490 familyNumber = number.str();
493 std::string creationString = family + familyNumber;
494 exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, creationString, dFESpace, fiberDirections[ k-1 ], UInt (0) );
497 *(fiberDirections[ k-1 ]) = solid.material()->anisotropicLaw()->ithFiberVector( k );
501 exporter->postProcess ( 0 );
502 exporterCheck->postProcess ( 0 );
503 std::cout.precision(16);
506 exporterCheck->postProcess ( 1.0 );
523 for (dataStructure->dataTime()->setTime ( dt ) ; dataStructure->dataTime()->canAdvance( ); dataStructure->dataTime()->updateTime( ) )
528 timeAdvance->updateRHSContribution ( dt );
529 *rhs += *solid.massMatrix() * timeAdvance->rhsContributionSecondDerivative() / timeAdvanceCoefficient;
534 if( !solid.havingSourceTerm() )
536 solid.setRightHandSide ( *rhs );
540 solid.updateRightHandSideWithBodyForce( dataStructure->dataTime()->time(), *rhs );
547 solid.iterate ( BCh );
549 timeAdvance->shiftRight ( solid.displacement() );
551 *solidDisp = solid.displacement();
552 *solidVel = timeAdvance->firstDerivative();
553 *solidAcc = timeAdvance->secondDerivative();
555 exporter->postProcess ( dataStructure->dataTime()->time() );
558 MPI_Barrier (MPI_COMM_WORLD);
560 Real errorAbs, error;
563 errorAbs = dFESpace->l2Error(
static_cast<solidFESpace_Type::function_Type> ( analyticalDisplacement ),
564 solid.displacement(), 0.0, &error);
566 std::cout <<
"L2 relative error: " << error << std::endl;
567 std::cout <<
"L2 abs error: " << errorAbs << std::endl;
575 MPI_Init (&argc, &argv);
576 std::shared_ptr<Epetra_MpiComm> Comm (
new Epetra_MpiComm ( MPI_COMM_WORLD ) );
577 if ( Comm->MyPID() == 0 )
579 std::cout <<
"% using MPI" << std::endl;
582 std::shared_ptr<Epetra_SerialComm> Comm (
new Epetra_SerialComm() );
583 std::cout <<
"% using serial Version" << std::endl;
586 Structure structure ( argc, argv, Comm );
std::shared_ptr< Epetra_Comm > comm
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)
std::function< VectorSmall< 3 > Real const &, Real const &, Real const &, Real const &) > bodyFunction_Type
int main(int argc, char **argv)
GetPot(const int argc_, char **argv_, const char *FieldSeparator=0x0)
static const LifeV::UInt elm_nodes_num[]
std::shared_ptr< hdf5Filter_Type > hdf5FilterPtr_Type
FESpace< RegionMesh< LinearTetra >, MapEpetra > solidFESpace_Type
std::shared_ptr< analyticalFunction_Type > analyticalFunctionPtr_Type
void run3d()
run the 3D driven cylinder simulation
ETFESpace< RegionMesh< LinearTetra >, MapEpetra, 3, 3 > solidETFESpace_Type
std::string data_file_name
void run2d()
run the 2D driven cylinder simulation
MeshData - class for handling spatial discretization.
std::shared_ptr< bodyFunction_Type > bodyFunctionPtr_Type
std::shared_ptr< Private > parameters
std::shared_ptr< vector_Type > vectorPtr_Type
double Real
Generic real data.
std::set< UInt > parseList(const std::string &list)
std::shared_ptr< solidETFESpace_Type > solidETFESpacePtr_Type
std::function< Real(Real const &, Real const &, Real const &, Real const &, ID const &) > analyticalFunction_Type
std::shared_ptr< TimeAdvance< vector_Type > > timeAdvance_Type
std::vector< vectorPtr_Type > listOfFiberDirections_Type
bool operator()(const char *VarName, bool Default) const
StructuralOperator< RegionMesh< LinearTetra > >::vector_Type vector_Type
std::shared_ptr< solidFESpace_Type > solidFESpacePtr_Type
uint32_type UInt
generic unsigned integer (used mainly for addressing)
std::shared_ptr< vectorFiberFunction_Type > vectorFiberFunctionPtr_Type
std::vector< fiberFunctionPtr_Type > vectorFiberFunction_Type