lifev/structure/examples/example_tractionWithSymmetry/main.cpp

Go to the documentation of this file.

//@HEADER
/*
*******************************************************************************

    Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
    Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University

    This file is part of LifeV.

    LifeV is free software; you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    LifeV is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with LifeV.  If not, see <http://www.gnu.org/licenses/>.

*******************************************************************************
*/
//@HEADER
/**
   \file main.cpp
   \author Christophe Prud'homme <christophe.prudhomme@epfl.ch>
   \date 2005-04-16
 */

#ifdef TWODIM
#error test_structure cannot be compiled in 2D
#endif

// Tell the compiler to ignore specific kind of warnings:
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"

#include <Epetra_ConfigDefs.h>
#ifdef EPETRA_MPI
#include <mpi.h>
#include <Epetra_MpiComm.h>
#else
#include <Epetra_SerialComm.h>
#endif

//Tell the compiler to restore the warning previously silented
#pragma GCC diagnostic warning "-Wunused-variable"
#pragma GCC diagnostic warning "-Wunused-parameter"

#include <lifev/core/LifeV.hpp>
#include <lifev/core/algorithm/PreconditionerIfpack.hpp>
#include <lifev/core/algorithm/PreconditionerML.hpp>


//Include fils which were in the structure.cpp file
#include <lifev/core/array/MapEpetra.hpp>

#include <lifev/core/fem/TimeAdvance.hpp>
#include <lifev/core/fem/TimeAdvanceNewmark.hpp>
#include <lifev/core/fem/TimeAdvanceBDF.hpp>

#include <lifev/core/mesh/MeshData.hpp>
#include <lifev/core/mesh/MeshPartitioner.hpp>

#include <lifev/structure/solver/StructuralConstitutiveLawData.hpp>

#include <lifev/structure/solver/StructuralConstitutiveLaw.hpp>
#include <lifev/structure/solver/StructuralOperator.hpp>
#include <lifev/structure/solver/VenantKirchhoffMaterialLinear.hpp>
#include <lifev/structure/solver/VenantKirchhoffMaterialNonLinear.hpp>
#include <lifev/structure/solver/VenantKirchhoffMaterialNonLinearPenalized.hpp>
#include <lifev/structure/solver/NeoHookeanMaterialNonLinear.hpp>
#include <lifev/structure/solver/ExponentialMaterialNonLinear.hpp>
#include <lifev/structure/solver/SecondOrderExponentialMaterialNonLinear.hpp>

#include <lifev/core/filter/ExporterEnsight.hpp>
#ifdef HAVE_HDF5
#include <lifev/core/filter/ExporterHDF5.hpp>
#endif
#include <lifev/core/filter/ExporterEmpty.hpp>

#include <iostream>


using namespace LifeV;

int returnValue = EXIT_FAILURE;
enum TimeScheme { BDF_ORDER_ONE = 1, BDF_ORDER_TWO = 2, BDF_ORDER_THREE = 3 };

namespace
{
static bool regIF = (PRECFactory::instance().registerProduct ( "Ifpack", &createIfpack ) );
static bool regML = (PRECFactory::instance().registerProduct ( "ML", &createML ) );
}


std::set<UInt> parseList ( const std::string& list )
{
    std::string stringList = list;
    std::set<UInt> setList;
    if ( list == "" )
    {
        return setList;
    }
    size_t commaPos = 0;
    while ( commaPos != std::string::npos )
    {
        commaPos = stringList.find ( "," );
        setList.insert ( atoi ( stringList.substr ( 0, commaPos ).c_str() ) );
        stringList = stringList.substr ( commaPos + 1 );
    }
    setList.insert ( atoi ( stringList.c_str() ) );
    return setList;
}


class Structure
{
public:

    /** @name Constructors, destructor
     */
    //@{
    Structure ( int                                   argc,
                char**                                argv,
                std::shared_ptr<Epetra_Comm>        structComm );

    ~Structure()
    {}
    //@}

    //@{
    void run()
    {
        run3d();
    }
    void CheckResultLE (const Real& dispNorm, const Real& time);
    void CheckResultSVK (const Real& dispNorm, const Real& time);
    void CheckResultEXP (const Real& dispNorm, const Real& time);
    void CheckResultNH (const Real& dispNorm, const Real& time);
    void resultChanged (Real time);
    //@}

protected:

private:

    /**
     * run the 2D driven cylinder simulation
     */
    void run2d();

    /**
     * run the 3D driven cylinder simulation
     */
    void run3d();

private:
    struct Private;
    std::shared_ptr<Private> parameters;
};



struct Structure::Private
{
    Private() :
        rho (1), poisson (1), young (1), bulk (1), alpha (1), gamma (1)
    {}
    typedef std::function<Real ( Real const&, Real const&, Real const&, Real const&, ID const& ) > fct_type;
    double rho, poisson, young, bulk, alpha, gamma;

    std::string data_file_name;

    std::shared_ptr<Epetra_Comm>     comm;

    static Real bcZero (const Real& /*t*/, const Real&  /*X*/, const Real& /*Y*/, const Real& /*Z*/, const ID& /*i*/)
    {
        return  0.;
    }

    static Real bcNonZero (const Real& /*t*/, const Real&  /*X*/, const Real& /*Y*/, const Real& /*Z*/, const ID& /*i*/)
    {
        return  70000; //136400.0
    }

    static Real d0 (const Real& /*t*/, const Real& x, const Real& y, const Real& z, const ID& i)
    {
        switch (i)
        {
            case 0:
                return - 0.058549 * ( x - 0.5 );
                break;
            case 1:
                return  ( 0.256942 / 2 ) * y;
                break;
            case 2:
                return - 0.058549 * ( z + 0.5);
                break;
            default:
                ERROR_MSG ("This entry is not allowed: ud_functions.hpp");
                return 0.;
                break;
        }
    }


};



Structure::Structure ( int                                   argc,
                       char**                                argv,
                       std::shared_ptr<Epetra_Comm>        structComm) :
    parameters ( new Private() )
{
    GetPot command_line (argc, argv);
    std::string data_file_name = command_line.follow ("data", 2, "-f", "--file");
    GetPot dataFile ( data_file_name );
    parameters->data_file_name = data_file_name;

    parameters->rho     = dataFile ( "solid/physics/density", 1. );
    parameters->young   = dataFile ( "solid/physics/young",   1. );
    parameters->poisson = dataFile ( "solid/physics/poisson", 1. );
    parameters->bulk    = dataFile ( "solid/physics/bulk",    1. );
    parameters->alpha   = dataFile ( "solid/physics/alpha",   1. );
    parameters->gamma   = dataFile ( "solid/physics/gamma",   1. );

    std::cout << "density = " << parameters->rho     << std::endl
              << "young   = " << parameters->young   << std::endl
              << "poisson = " << parameters->poisson << std::endl
              << "bulk    = " << parameters->bulk    << std::endl
              << "alpha   = " << parameters->alpha   << std::endl
              << "gamma   = " << parameters->gamma   << std::endl;

    parameters->comm = structComm;
    int ntasks = parameters->comm->NumProc();

    if (!parameters->comm->MyPID() )
    {
        std::cout << "My PID = " << parameters->comm->MyPID() << " out of " << ntasks << " running." << std::endl;
    }
}



void
Structure::run2d()
{
    std::cout << "2D cylinder test case is not available yet\n";
}



void
Structure::run3d()
{
    typedef StructuralOperator< RegionMesh<LinearTetra> >::vector_Type  vector_Type;
    typedef std::shared_ptr<vector_Type> vectorPtr_Type;
    typedef std::shared_ptr< TimeAdvance< vector_Type > >       timeAdvance_type;
    typedef RegionMesh<LifeV::LinearTetra >                       mesh_Type;
    typedef LifeV::ExporterEmpty<mesh_Type >                      emptyExporter_Type;
    typedef std::shared_ptr<emptyExporter_Type>                 emptyExporterPtr_Type;

    typedef LifeV::ExporterEnsight<mesh_Type >                    ensightFilter_Type;
    typedef std::shared_ptr<ensightFilter_Type>                 ensightFilterPtr_Type;

#ifdef HAVE_HDF5
    typedef LifeV::ExporterHDF5<mesh_Type >                       hdf5Filter_Type;
    typedef std::shared_ptr<hdf5Filter_Type>                    hdf5FilterPtr_Type;
#endif

    bool verbose = (parameters->comm->MyPID() == 0);

    //! Number of boundary conditions for the velocity and mesh motion
    std::shared_ptr<BCHandler> BCh ( new BCHandler() );

    //! dataElasticStructure
    GetPot dataFile ( parameters->data_file_name.c_str() );

    std::shared_ptr<StructuralConstitutiveLawData> dataStructure (new StructuralConstitutiveLawData( ) );
    dataStructure->setup (dataFile);

    MeshData             meshData;
    meshData.setup (dataFile, "solid/space_discretization");

    std::shared_ptr<RegionMesh<LinearTetra> > fullMeshPtr (new RegionMesh<LinearTetra> (  parameters->comm  ) );
    readMesh (*fullMeshPtr, meshData);

    MeshPartitioner< RegionMesh<LinearTetra> > meshPart ( fullMeshPtr, parameters->comm );

    std::string dOrder =  dataFile ( "solid/space_discretization/order", "P1");

    typedef FESpace< RegionMesh<LinearTetra>, MapEpetra > solidFESpace_type;
    typedef std::shared_ptr<solidFESpace_type> solidFESpace_ptrtype;

    typedef ETFESpace< RegionMesh<LinearTetra>, MapEpetra, 3, 3 >       solidETFESpace_Type;
    typedef std::shared_ptr<solidETFESpace_Type>                      solidETFESpacePtr_Type;


    solidFESpace_ptrtype dFESpace ( new solidFESpace_type (meshPart, dOrder, 3, parameters->comm) );
    solidETFESpacePtr_Type dETFESpace ( new solidETFESpace_Type (meshPart, & (dFESpace->refFE() ), & (dFESpace->fe().geoMap() ), parameters->comm) );
    if (verbose)
    {
        std::cout << std::endl;
    }

    std::string timeAdvanceMethod =  dataFile ( "solid/time_discretization/method", "Newmark");

    std::cout << "Method: " << timeAdvanceMethod << std::endl;

    timeAdvance_type  timeAdvance ( TimeAdvanceFactory::instance().createObject ( timeAdvanceMethod ) );

    UInt OrderDev = 2;

    //! initialization of parameters of time Advance method:
    if (timeAdvanceMethod == "Newmark")
    {
        timeAdvance->setup ( dataStructure->dataTimeAdvance()->coefficientsNewmark() , OrderDev);
    }

    if (timeAdvanceMethod == "BDF")
    {
        timeAdvance->setup (dataStructure->dataTimeAdvance()->orderBDF() , OrderDev);
    }

    timeAdvance->setTimeStep (dataStructure->dataTime()->timeStep() );

    //! #################################################################################
    //! BOUNDARY CONDITIONS
    //! #################################################################################
    std::vector <ID> compx (1), compy (1), compz (1), compxy (2), compxz (2), compyz (2);
    compx[0] = 0;
    compy[0] = 1, compz[0] = 2;
    compxy[0] = 0;
    compxy[1] = 1;
    compxz[0] = 0;
    compxz[1] = 2;
    compyz[0] = 1;
    compyz[1] = 2;

    BCFunctionBase zero (Private::bcZero);
    BCFunctionBase nonZero (Private::bcNonZero);

    //! =================================================================================
    //! BC for StructuredCube4_test_structuralsolver.mesh
    //! =================================================================================
    BCh->addBC ("EdgesIn",      20,  Natural,   Component, nonZero, compy);
    BCh->addBC ("EdgesIn",      40,  Essential, Component, zero,    compy);

    BCh->addBC ("SymmetryX",    5,  Essential, Component, zero,    compx);
    BCh->addBC ("SymmetryY",    3,  Essential, Component, zero,    compz);
    BCh->addBC ("edgeone",      100,  EssentialVertices, Full, zero, 3);
    BCh->addBC ("edgetwo",      50,  EssentialVertices, Component, zero,    compxy);
    BCh->addBC ("edgetwo",      30,  EssentialVertices, Component, zero,    compyz);
    BCh->addBC ("edgetwo",      80,  EssentialVertices, Component, zero,    compxz);


    //! 1. Constructor of the structuralSolver
    StructuralOperator< RegionMesh<LinearTetra> > solid;

    //! 2. Setup of the structuralSolver
    solid.setup (dataStructure,
                 dFESpace,
                 dETFESpace,
                 BCh,
                 parameters->comm);

    //! 3. Setting data from getPot
    solid.setDataFromGetPot (dataFile);

    //! 4. Building system using TimeAdvance class
    double timeAdvanceCoefficient = timeAdvance->coefficientSecondDerivative ( 0 ) / (dataStructure->dataTime()->timeStep() * dataStructure->dataTime()->timeStep() );
    solid.buildSystem (timeAdvanceCoefficient);


    dataStructure->showMe();
    //! =================================================================================
    //! Temporal data and initial conditions
    //! =================================================================================

    //! 5. Initial data
    Real dt = dataStructure->dataTime()->timeStep();
    //Real T  = dataStructure->dataTime()->endTime();

    vectorPtr_Type rhs (new vector_Type (solid.displacement(), Unique) );
    vectorPtr_Type disp (new vector_Type (solid.displacement(), Unique) );
    vectorPtr_Type vel (new vector_Type (solid.displacement(), Unique) );
    vectorPtr_Type acc (new vector_Type (solid.displacement(), Unique) );

    //Initialization of TimeAdvance
    std::string const restart =  dataFile ( "importer/restart", "none");
    std::vector<vectorPtr_Type> solutionStencil;
    solutionStencil.resize ( timeAdvance->size() );

    std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > importerSolid;
    if ( restart.compare ( "none" ) )
    {
        //Reading fileNames - setting data for reading
        std::string const importerType =  dataFile ( "importer/type", "ensight");
        std::string const fileName     =  dataFile ( "importer/filename", "structure");
        std::string const initialLoaded     =  dataFile ( "importer/initialSol", "NO_DEFAULT_VALUE");
        //LifeV::Real initialTime        =  dataFile ( "importer/initialTime", 0.0);

        //Creating the importer
#ifdef HAVE_HDF5
        if ( !importerType.compare ("hdf5") )
        {
            importerSolid.reset ( new  hdf5Filter_Type ( dataFile, fileName) );
        }
        else
#endif
        {
            if ( !importerType.compare ("none") )
            {
                importerSolid.reset ( new emptyExporter_Type ( dataFile, dFESpace->mesh(), "solid", dFESpace->map().comm().MyPID() ) );
            }
            else
            {
                importerSolid.reset ( new  ensightFilter_Type ( dataFile, fileName) );
            }
        }

        importerSolid->setMeshProcId (dFESpace->mesh(), dFESpace->map().comm().MyPID() );

        //Creation of Exporter to check the loaded solution (working only for HDF5)
        // std::string expVerFile = "verificationDisplExporter";
        // LifeV::ExporterHDF5<RegionMesh<LinearTetra> > exporter( dataFile, meshPart.meshPartition(), expVerFile, parameters->comm->MyPID());
        // vectorPtr_Type vectVer ( new vector_Type(solid.displacement(),  LifeV::Unique ) );

        // exporter.addVariable( ExporterData<mesh_Type >::VectorField, "displVer", dFESpace, vectVer, UInt(0) );

        // exporter.postProcess(0.0);

        //Reading the displacement field and the timesteps need by the TimeAdvance class
        vectorPtr_Type solidDisp (new vector_Type (dFESpace->map(), importerSolid->mapType() ) );

        std::string iterationString;

        std::cout << "size TimeAdvance:" << timeAdvance->size() << std::endl;

        //Loading the stencil
        iterationString = initialLoaded;
        for (UInt iterInit = 0; iterInit < timeAdvance->size(); iterInit++ )
        {
            solidDisp.reset (new vector_Type (dFESpace->map(), importerSolid->mapType() ) );
            *solidDisp *= 0.0;

            LifeV::ExporterData<mesh_Type> solidDataReader (LifeV::ExporterData<mesh_Type>::VectorField, std::string ("displacement." + iterationString), dFESpace, solidDisp, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );

            importerSolid->readVariable (solidDataReader);

            std::cout << "Norm of the " << iterInit + 1 << "-th solution : " << solidDisp->norm2() << std::endl;

            //Exporting the just loaded solution (debug purposes)
            // Real currentLoading(iterInit + 1.0);
            // *vectVer = *solidDisp;
            // exporter.postProcess( currentLoading );

            //solutionStencil.push_back( solidDisp );
            solutionStencil[ iterInit ] = solidDisp;
            //initializing the displacement field in the StructuralSolver class with the first solution
            if ( !iterInit )
            {
                solid.initialize ( solidDisp );
            }

            //Updating string name
            int iterations = std::atoi (iterationString.c_str() );
            iterations--;

            std::ostringstream iter;
            iter.fill ( '0' );
            iter << std::setw (5) << ( iterations );
            iterationString = iter.str();

        }

        importerSolid->closeFile();

        //Putting the vector in the TimeAdvance Stencil
        timeAdvance->setInitialCondition (solutionStencil);
    }
    else
    {
        if (verbose)
        {
            std::cout << "S- initialization ... ";
        }

        std::vector<vectorPtr_Type> uv0;

        if (timeAdvanceMethod == "Newmark")
        {
            uv0.push_back (disp);
            uv0.push_back (vel);
            uv0.push_back (acc);
        }

        vectorPtr_Type initialDisplacement (new vector_Type (solid.displacement(), Unique) );
        dFESpace->interpolate ( static_cast<solidFESpace_type::function_Type> ( Private::d0 ), *initialDisplacement, 0.0 );

        if (timeAdvanceMethod == "BDF")
        {
            Real tZero = dataStructure->dataTime()->initialTime();

            for ( UInt previousPass = 0; previousPass < timeAdvance->size() ; previousPass++)
            {
                Real previousTimeStep = tZero - previousPass * dt;
                std::cout << "BDF " << previousTimeStep << "\n";
                //uv0.push_back(disp);
                uv0.push_back (initialDisplacement);
            }
        }

        timeAdvance->setInitialCondition (uv0);

        timeAdvance->setTimeStep ( dt );

        timeAdvance->updateRHSContribution ( dt );


        // debug purposes
        // vector_Type RHSfirstTerm( timeAdvance->rhsContributionFirstDerivative() );
        // vector_Type RHSsecondTerm( timeAdvance->rhsContributionSecondDerivative() );

        // vector_Type velInitialize( timeAdvance->firstDerivative() );
        // vector_Type accInitialize( timeAdvance->secondDerivative() );


        // std::cout << "Norm RHS 1: " << RHSfirstTerm.norm2() << std::endl;
        // std::cout << "Norm RHS 2: " << RHSsecondTerm.norm2() << std::endl;


        // std::cout << "Velocity 1: " << velInitialize.norm2() << std::endl;
        // std::cout << "Acceleration 2: " << accInitialize.norm2() << std::endl;

        //In the case of non-zero displacement
        //solid.initialize( disp );
        solid.initialize ( initialDisplacement );
        //Let's verify that the set displacement is the one I expect
        //Creation of Exporter to check the loaded solution (working only for HDF5)
        std::string expVerFile = "verificationDisplExporter";
        LifeV::ExporterHDF5<RegionMesh<LinearTetra> > exporter ( dataFile, meshPart.meshPartition(), expVerFile, parameters->comm->MyPID() );
        vectorPtr_Type vectVer ( new vector_Type (solid.displacement(),  LifeV::Unique ) );

        exporter.addVariable ( ExporterData<mesh_Type >::VectorField, "displVer", dFESpace, vectVer, UInt (0) );

        //Let's get the initial displacement and velocity
        exporter.postProcess (0.0);

        *vectVer = *initialDisplacement; //*disp;
        exporter.postProcess ( 1.0 );

    }
    MPI_Barrier (MPI_COMM_WORLD);

    if (verbose )
    {
        std::cout << "ok." << std::endl;
    }

    std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporter;
    //    std::shared_ptr< Exporter<RegionMesh<LinearTetra> > > exporterCheck;

    std::string const exporterType =  dataFile ( "exporter/type", "ensight");
    std::string const exporterName =  dataFile ( "exporter/name", "structure");
    //    std::string const exporterCheckName =  dataFile ( "exporter/nameCheck", "verifyVectors");
#ifdef HAVE_HDF5
    if (exporterType.compare ("hdf5") == 0)
    {
        exporter.reset ( new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, exporterName ) );
        //        exporterCheck.reset ( new ExporterHDF5<RegionMesh<LinearTetra> > ( dataFile, exporterCheckName ) );
    }
    else
#endif
    {
        if (exporterType.compare ("none") == 0)
        {
            exporter.reset ( new ExporterEmpty<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(), exporterName, parameters->comm->MyPID() ) );
        }

        else
        {
            exporter.reset ( new ExporterEnsight<RegionMesh<LinearTetra> > ( dataFile, meshPart.meshPartition(), exporterName, parameters->comm->MyPID() ) );
        }
    }

    exporter->setPostDir ( "./" );
    exporter->setMeshProcId ( meshPart.meshPartition(), parameters->comm->MyPID() );

    // exporterCheck->setPostDir ( "./" );
    // exporterCheck->setMeshProcId ( meshPart.meshPartition(), parameters->comm->MyPID() );

    vectorPtr_Type solidDisp ( new vector_Type (solid.displacement(),  exporter->mapType() ) );
    vectorPtr_Type solidVel  ( new vector_Type (solid.displacement(),  exporter->mapType() ) );
    vectorPtr_Type solidAcc  ( new vector_Type (solid.displacement(),  exporter->mapType() ) );

    // vectorPtr_Type rhsCopy ( new vector_Type (solid.displacement(),  exporterCheck->mapType() ) );
    // vectorPtr_Type residualCopy ( new vector_Type (solid.displacement(),  exporterCheck->mapType() ) );

    exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "displacement", dFESpace, solidDisp, UInt (0) );
    exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "velocity",     dFESpace, solidVel,  UInt (0) );
    exporter->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "acceleration", dFESpace, solidAcc,  UInt (0) );

    // exporterCheck->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "rhs", dFESpace, rhsCopy,  UInt (0) );
    // exporterCheck->addVariable ( ExporterData<RegionMesh<LinearTetra> >::VectorField, "residual", dFESpace, residualCopy,  UInt (0) );

    //Let's get the initial quantities
    *solidDisp = solid.displacement();
    *solidVel = timeAdvance->firstDerivative();
    *solidAcc = timeAdvance->secondDerivative();

    exporter->postProcess ( 0 );
    //exporterCheck->postProcess ( 0 );


    //!--------------------------------------------------------------------------------------------
    //! MATLAB FILE WITH DISPLACEMENT OF A CHOSEN POINT
    //!--------------------------------------------------------------------------------------------
    std::cout.precision (16);

    // Used to debug
    int IDPoint = 158;
    // This one depends on the direction of the deformation
    int indexPoint = IDPoint - 1 + dFESpace->dof().numTotalDof();
    //!--------------------------------------------------------------------------------------------
    //!The update of the RHS is done by the TimeAdvance class
    //solid.updateSystem();
    //! =================================================================================



    //! =============================================================================
    //! Temporal loop
    //! =============================================================================
    //    for (Real time = dt; time <= T; time += dt)
    for (dataStructure->dataTime()->setTime ( dt ) ; dataStructure->dataTime()->canAdvance( ); dataStructure->dataTime()->updateTime( ) )
    {
        //      dataStructure->dataTime()->setTime(time);

        if (verbose)
        {
            std::cout << std::endl;
            std::cout << "S- Now we are at time " << dataStructure->dataTime()->time() << " s." << std::endl;
        }

        //! 6. Updating right-hand side
        *rhs *= 0;
        timeAdvance->updateRHSContribution ( dt );
        *rhs += *solid.massMatrix() * timeAdvance->rhsContributionSecondDerivative() / timeAdvanceCoefficient;
        solid.setRightHandSide ( *rhs );

        //! 7. Iterate --> Calling Newton
        solid.iterate ( BCh );

        timeAdvance->shiftRight ( solid.displacement() );

        *solidDisp = solid.displacement();
        *solidVel  = timeAdvance->firstDerivative();
        *solidAcc  = timeAdvance->secondDerivative();

        // *rhsCopy = solid.rhsCopy();
        // *residualCopy = solid.residualCopy();

        exporter->postProcess ( dataStructure->dataTime()->time() );
        //exporterCheck->postProcess ( dataStructure->dataTime()->time() );

        cout << "*********************************************************" << std::endl;
        cout << " solid.disp()[ " << indexPoint  << " ] = " <<  solid.displacement() [ indexPoint ]  << std::endl;
        cout << "*********************************************************" << std::endl;


        Real normVect;
        normVect =  solid.displacement().norm2();
        std::cout << "The norm 2 of the displacement field is: " << normVect << std::endl;

        //!--------------------------------------------------------------------------------------------------

        MPI_Barrier (MPI_COMM_WORLD);
    }
}



void Structure::CheckResultLE (const Real& dispNorm, const Real& time)
{
    if ( time == 0.1  && std::fabs (dispNorm - 0.276527) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.2  && std::fabs (dispNorm - 0.276536) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.3  && std::fabs (dispNorm - 0.276529) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.4  && std::fabs (dispNorm - 0.276531) <= 1e-5 )
    {
        this->resultChanged (time);
    }
}

void Structure::CheckResultSVK (const Real& dispNorm, const Real& time)
{
    if ( time == 0.1  && std::fabs (dispNorm - 0.263348) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.2  && std::fabs (dispNorm - 0.263350) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.3  && std::fabs (dispNorm - 0.263350) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.4  && std::fabs (dispNorm - 0.263351) <= 1e-5 )
    {
        this->resultChanged (time);
    }
}
void Structure::CheckResultEXP (const Real& dispNorm, const Real& time)
{
    if ( time == 0.1  && std::fabs (dispNorm - 0.284844) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.2  && std::fabs (dispNorm - 0.284853) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.3  && std::fabs (dispNorm - 0.284846) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.4  && std::fabs (dispNorm - 0.284848) <= 1e-5 )
    {
        this->resultChanged (time);
    }
}

void Structure::CheckResultNH (const Real& dispNorm, const Real& time)
{
    if ( time == 0.1  && std::fabs (dispNorm - 0.286120) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.2  && std::fabs (dispNorm - 0.286129) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.3  && std::fabs (dispNorm - 0.286122) <= 1e-5 )
    {
        this->resultChanged (time);
    }
    if ( time == 0.4  && std::fabs (dispNorm - 0.286123) <= 1e-5 )
    {
        this->resultChanged (time);
    }
}



void Structure::resultChanged (Real time)
{
    std::cout << "Correct value at time: " << time << std::endl;
    returnValue = EXIT_SUCCESS;
}



int
main ( int argc, char** argv )
{

#ifdef HAVE_MPI
    MPI_Init (&argc, &argv);
    std::shared_ptr<Epetra_MpiComm> Comm (new Epetra_MpiComm ( MPI_COMM_WORLD ) );
    if ( Comm->MyPID() == 0 )
    {
        std::cout << "% using MPI" << std::endl;
    }
#else
    std::shared_ptr<Epetra_SerialComm> Comm ( new Epetra_SerialComm() );
    std::cout << "% using serial Version" << std::endl;
#endif

    Structure structure ( argc, argv, Comm );
    structure.run();

#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return returnValue ;
}