C++ LifeChrono::diff方法代码示例

void ExporterEnsight<MeshType>::import (const Real& time)
{
    this->M_timeSteps.push_back (time);
    ++this->M_steps;

    // typedef std::list< exporterData_Type >::iterator Iterator;

    this->computePostfix();

    assert ( this->M_postfix.find ( "*" ) == std::string::npos );

    if (!this->M_procId)
    {
        std::cout << "  X-  ExporterEnsight importing ..." << std::endl;
    }

    LifeChrono chrono;
    chrono.start();
    for (typename super::dataVectorIterator_Type i = this->M_dataVector.begin(); i != this->M_dataVector.end(); ++i)
    {
        this->readVariable (*i);
    }
    chrono.stop();
    if (!this->M_procId)
    {
        std::cout << "      done in " << chrono.diff() << " s." << std::endl;
    }

}

int
PreconditionerComposed::push_back (operatorPtr_Type& oper,
                                   const bool useInverse,
                                   const bool useTranspose
                                  )
{
    if (!M_prec.get() )
    {
        M_prec.reset (new prec_Type (M_displayer.comm() ) );
    }
    M_operVector.push_back (oper);
    LifeChrono chrono;
    epetraPrecPtr_Type prec;

    this->M_displayer.leaderPrint ( std::string ("ICP-  Preconditioner type:                     ") + M_prec->Operator() [M_operVector.size() - 1]->preconditionerType() + std::string ("\n") );
    this->M_displayer.leaderPrint ( "ICP-  Computing preconditioner ...             " );
    chrono.start();
    createPrec (oper, M_prec->OperatorView() [M_operVector.size() - 1]);
    chrono.stop();
    this->M_displayer.leaderPrintMax ("done in ", chrono.diff() );
    M_prec->replace (prec, useInverse, useTranspose); // \TODO to reset as push_back
    if ( M_prec->Operator().size() == M_operVector.size() )
    {
        this->M_preconditionerCreated = true;
    }
    return EXIT_SUCCESS;
}

Int SolverAztecOO::solveSystem ( const vector_type& rhsFull,
                                 vector_type&       solution,
                                 matrix_ptrtype&    baseMatrixForPreconditioner )

{

    bool retry ( true );

    LifeChrono chrono;

    M_displayer->leaderPrint ( "SLV-  Setting up the solver ...                \n" );

    if ( baseMatrixForPreconditioner.get() == 0 )
    {
        M_displayer->leaderPrint ( "SLV-  Warning: baseMatrixForPreconditioner is empty     \n" );
    }

    if ( !isPreconditionerSet() || !M_reusePreconditioner  )
    {
        buildPreconditioner ( baseMatrixForPreconditioner );
        // do not retry if I am recomputing the preconditioner
        retry = false;
    }
    else
    {
        M_displayer->leaderPrint ( "SLV-  Reusing precond ...                 \n" );
    }

    Int numIter = solveSystem ( rhsFull, solution, M_preconditioner );

    // If we do not want to retry, return now.
    // otherwise rebuild the preconditioner and solve again:
    if ( numIter < 0  && retry )
    {
        chrono.start();

        M_displayer->leaderPrint ( "SLV-  Iterative solver failed, numiter = " , - numIter );
        M_displayer->leaderPrint ( "SLV-  maxIterSolver = " , M_maxIter );
        M_displayer->leaderPrint ( "SLV-  retrying:          " );

        buildPreconditioner ( baseMatrixForPreconditioner );

        chrono.stop();
        M_displayer->leaderPrintMax ( "done in " , chrono.diff() );
        // Solving again, but only once (retry = false)
        numIter = solveSystem ( rhsFull, solution, M_preconditioner );

        if ( numIter < 0 )
        {
            M_displayer->leaderPrint ( " ERROR: Iterative solver failed again.\n" );
        }
    }

    if ( std::abs (numIter) > M_maxIterForReuse )
    {
        resetPreconditioner();
    }

    return numIter;
}

void
WallTensionEstimatorCylindricalCoordinates<Mesh >::analyzeTensionsRecoveryCauchyStressesCylindrical ( void )
{

    LifeChrono chrono;

    chrono.start();
    UInt dim = this->M_FESpace->dim();

    constructGlobalStressVector();


    for ( UInt iDOF = 0; iDOF < ( UInt ) this->M_FESpace->dof().numTotalDof(); iDOF++ )
    {

        if ( this->M_displacement->blockMap().LID ( static_cast<EpetraInt_Type> (iDOF) ) != -1 ) // The Global ID is on the calling processors
        {

            (* (this->M_sigma) ).Scale (0.0);

            //Extracting the gradient of U on the current DOF
            for ( UInt iComp = 0; iComp < this->M_FESpace->fieldDim(); ++iComp )
            {
                Int LIDid = this->M_displacement->blockMap().LID ( static_cast<EpetraInt_Type> (iDOF + iComp * dim + this->M_offset ) );
                Int GIDid = this->M_displacement->blockMap().GID (static_cast<EpetraInt_Type> (LIDid) );
                (* (this->M_sigma) ) (iComp, 0) = (*this->M_sigmaX) (GIDid); // (d_xX,d_yX,d_zX)
                (* (this->M_sigma) ) (iComp, 1) = (*this->M_sigmaY) (GIDid); // (d_xY,d_yY,d_zY)
                (* (this->M_sigma) ) (iComp, 2) = (*this->M_sigmaZ) (GIDid); // (d_xZ,d_yZ,d_zZ)
            }

            //Compute the eigenvalue
            AssemblyElementalStructure::computeEigenvalues (* (this->M_sigma), this->M_eigenvaluesR, this->M_eigenvaluesI);

            //The Cauchy tensor is symmetric and therefore, the eigenvalues are real
            //Check on the imaginary part of eigen values given by the Lapack method
            Real sum (0);
            for ( int i = 0; i < this->M_eigenvaluesI.size(); i++ )
            {
                sum += std::abs (this->M_eigenvaluesI[i]);
            }
            ASSERT_PRE ( sum < 1e-6 , "The eigenvalues of the Cauchy stress tensors have to be real!" );

            std::sort ( this->M_eigenvaluesR.begin(), this->M_eigenvaluesR.end() );

            //Save the eigenvalues in the global vector
            for ( UInt icoor = 0; icoor < this->M_FESpace->fieldDim(); ++icoor )
            {
                Int LIDid = this->M_displacement->blockMap().LID ( static_cast<EpetraInt_Type> (iDOF + icoor * dim + this->M_offset) );
                Int GIDid = this->M_displacement->blockMap().GID (static_cast<EpetraInt_Type> (LIDid) );
                (* (this->M_globalEigenvalues) ) (GIDid) = this->M_eigenvaluesR[icoor];
            }

        }
    }

    chrono.stop();
    this->M_displayer->leaderPrint ("Analysis done in: ", chrono.diff() );

}

void ExporterEnsight<MeshType>::postProcess (const Real& time)
{
    // writing the geo file and the list of global IDs, but only upon the first instance
    if ( M_firstTimeStep )
    {
        if (!this->M_multimesh)
        {
            writeAsciiGeometry ( this->M_postDir + this->M_prefix + this->M_me + ".geo" );
        }

        writeGlobalIDs ( this->M_postDir + super::M_prefix + "_globalIDs" +
                         this->M_me + ".scl" );

        M_firstTimeStep = false;
    }
    // prepare the file postfix
    this->computePostfix();

    // the postfix will be full of stars, if this time step is not going to generate a snapshot
    std::size_t found ( this->M_postfix.find ( "*" ) );
    if ( found == std::string::npos )
    {
        if (!this->M_procId)
        {
            std::cout << "  X-  ExporterEnsight post-processing ...        " << std::flush;
        }
        LifeChrono chrono;
        chrono.start();
        for (typename super::dataVectorIterator_Type i = this->M_dataVector.begin();
                i != this->M_dataVector.end(); ++i)
        {
            // the "regime" attribute needs to be valid
            if ( i->regime() != exporterData_Type::NullRegime )
            {
                writeAscii (*i);
            }
            // if the solution is steady, we do not need to export it at each time step
            if (i->regime() == exporterData_Type::SteadyRegime)
            {
                i->setRegime ( exporterData_Type::NullRegime );
            }
        }
        // write an updated case file
        writeCase (time);

        // write an updated geo file, if needed
        if (this->M_multimesh)
        {
            writeAsciiGeometry ( this->M_postDir + this->M_prefix + this->M_postfix + this->M_me + ".geo" );
        }
        chrono.stop();
        if (!this->M_procId)
        {
            std::cout << "      done in " << chrono.diff() << " s." << std::endl;
        }
    }

}

void
LinearSolver::buildPreconditioner()
{
    LifeChrono chrono;
    Real condest ( -1 );

    if ( M_preconditioner )
    {
        if ( M_matrix.get() == 0 )
        {
            M_displayer->leaderPrint ( "SLV-  ERROR: LinearSolver requires a matrix to build the preconditioner!\n" );
            exit ( 1 );
        }
        else
        {
            chrono.start();
            if ( !M_silent )
            {
                M_displayer->leaderPrint ( "SLV-  Computing the preconditioner...\n" );
            }
            if ( M_baseMatrixForPreconditioner.get() == 0 )
            {
                if ( !M_silent )
                {
                    M_displayer->leaderPrint ( "SLV-  Build the preconditioner using the problem matrix\n" );
                }
                M_preconditioner->buildPreconditioner ( M_matrix );
            }
            else
            {
                if ( !M_silent )
                {
                    M_displayer->leaderPrint ( "SLV-  Build the preconditioner using the base matrix provided\n" );
                }
                M_preconditioner->buildPreconditioner ( M_baseMatrixForPreconditioner );
            }
            condest = M_preconditioner->condest();
            chrono.stop();
            if ( !M_silent )
            {
                M_displayer->leaderPrintMax ( "SLV-  Preconditioner computed in " , chrono.diff(), " s." );
            }
            if ( !M_silent )
            {
                M_displayer->leaderPrint ( "SLV-  Estimated condition number               " , condest, "\n" );
            }
        }
    }
}

void SolverAztecOO::buildPreconditioner ( matrix_ptrtype& preconditioner )
{
    LifeChrono chrono;
    Real condest (-1);

    chrono.start();

    M_displayer->leaderPrint ( "SLV-  Computing the precond ...                " );

    M_preconditioner->buildPreconditioner ( preconditioner );

    condest = M_preconditioner->condest();
    chrono.stop();

    M_displayer->leaderPrintMax ( "done in " , chrono.diff() );
    M_displayer->leaderPrint ( "SLV-  Estimated condition number               " , condest, "\n" );
}

void ExporterVTK<MeshType>::import(const Real& /*time*/)
{
    this->computePostfix();

    assert( this->M_postfix != "*****" );

    if (!this->M_procId) std::cout << "  X-  ExporterVTK importing ..."<< std::endl;

    LifeChrono chrono;
    chrono.start();
    for (typename super::dataVectorIterator_Type iData=this->M_dataVector.begin();
                    iData != this->M_dataVector.end(); ++iData)
    {
        this->readVTUFiles(*iData);
    }
    chrono.stop();
    if (!this->M_procId) std::cout << "      done in " << chrono.diff() << " s." << std::endl;

}

int
PreconditionerComposed::replace(operatorPtr_Type& oper,
                            const UInt index,
                            const bool useInverse,
                            const bool useTranspose)
{
    ASSERT(index <= M_operVector.size(), "PreconditionerComposed::replace: index too large");

    M_operVector[index] = oper;
    LifeChrono chrono;
    //ifpack_prec_type prec;
    this->M_displayer.leaderPrint( std::string("ICP-  Preconditioner type:                     ") + M_prec->Operator()[index]->preconditionerType() + std::string("\n") );
    this->M_displayer.leaderPrint( "ICP-  Computing preconditioner ...             " );
    chrono.start();
    createPrec(oper, M_prec->OperatorView()[index]);
    chrono.stop();
    this->M_displayer.leaderPrintMax("done in ", chrono.diff());

    M_prec->replace(M_prec->Operator()[index], index, useInverse, useTranspose);

    return EXIT_SUCCESS;
}

//.........这里部分代码省略.........
        std::cout << "ok." << std::endl;
    }

    fluid.setUp (dataFile);

    // Initialization
    Real dt     = oseenData->dataTime()->timeStep();
    Real t0     = oseenData->dataTime()->initialTime();
    Real tFinal = oseenData->dataTime()->endTime();

    boost::shared_ptr< Exporter<mesh_Type > > exporter;
    boost::shared_ptr< Exporter<mesh_Type > > importer;

    std::string const exporterType =  dataFile ( "exporter/type", "hdf5");
    std::string const exporterName =  dataFile ( "exporter/filename", "ethiersteinman");
    std::string const exportDir    =  dataFile ( "exporter/exportDir", "./");

    std::string const importerType =  dataFile ( "importer/type", "ensight");
    std::string const importerName =  dataFile ( "importer/filename", "ethiersteinman");
    std::string const importDir    =  dataFile ( "importer/importDir", "importDir/");

#ifdef HAVE_HDF5
    if (exporterType.compare ("hdf5") == 0)
    {
        exporter.reset ( new ExporterHDF5<mesh_Type > ( dataFile, exporterName ) );
    }
    else
#endif
        exporter.reset ( new ExporterEnsight<mesh_Type > ( dataFile, exporterName ) );

    exporter->setPostDir ( exportDir ); // This is a test to see if M_post_dir is working
    exporter->setMeshProcId ( meshPtr, d->comm->MyPID() );

#ifdef HAVE_HDF5
    if (importerType.compare ("hdf5") == 0)
    {
        importer.reset ( new ExporterHDF5<mesh_Type > ( dataFile, importerName ) );
    }
    else
#endif
        importer.reset ( new ExporterEnsight<mesh_Type > ( dataFile, importerName ) );

    // todo this will not work with the ExporterEnsight filter (it uses M_importDir, a private variable)
    importer->setPostDir ( importDir ); // This is a test to see if M_post_dir is working
    importer->setMeshProcId ( meshPtr, d->comm->MyPID() );

    vectorPtr_Type velAndPressureExport ( new vector_Type (*fluid.solution(), exporter->mapType() ) );
    vectorPtr_Type velAndPressureImport ( new vector_Type (*fluid.solution(), importer->mapType() ) );

    if ( exporter->mapType() == Unique )
    {
        velAndPressureExport->setCombineMode (Zero);
    }

    importer->addVariable ( ExporterData<mesh_Type>::VectorField, "velocity", uFESpacePtr,
                            velAndPressureImport, UInt (0) );
    importer->addVariable ( ExporterData<mesh_Type>::ScalarField, "pressure", pFESpacePtr,
                            velAndPressureImport, UInt (3 * uFESpacePtr->dof().numTotalDof() ) );
    importer->import ( t0 );

    *velAndPressureExport = *velAndPressureImport;


    vectorPtr_Type P0pres ( new vector_Type (p0FESpacePtr->map() ) );
    MPI_Barrier (MPI_COMM_WORLD);
    computeP0pressure (pFESpacePtr, p0FESpacePtr, uFESpacePtr, *velAndPressureImport, *P0pres, t0);

    exporter->addVariable ( ExporterData<mesh_Type>::VectorField, "velocity", uFESpacePtr,
                            velAndPressureExport, UInt (0) );

    exporter->addVariable ( ExporterData<mesh_Type>::ScalarField, "pressure", pFESpacePtr,
                            velAndPressureExport, UInt (3 * uFESpacePtr->dof().numTotalDof() ) );

    exporter->addVariable (ExporterData<mesh_Type>::ScalarField, "P0pressure", p0FESpacePtr,
                           P0pres, UInt (0),
                           ExporterData<mesh_Type>::SteadyRegime, ExporterData<mesh_Type>::Cell );
    exporter->postProcess ( t0 );

    // Temporal loop
    LifeChrono chrono;
    int iter = 1;

    for ( Real time = t0 + dt ; time <= tFinal + dt / 2.; time += dt, iter++)
    {
        chrono.stop();
        importer->import ( time );

        *velAndPressureExport = *velAndPressureImport;
        MPI_Barrier (MPI_COMM_WORLD);
        computeP0pressure (pFESpacePtr, p0FESpacePtr, uFESpacePtr, *velAndPressureImport, *P0pres, time);

        exporter->postProcess ( time );

        chrono.stop();
        if (verbose)
        {
            std::cout << "Total iteration time " << chrono.diff() << " s." << std::endl;
        }
    }
}

//.........这里部分代码省略.........

    std::vector<vectorPtr_Type> uv0;

    if (TimeAdvanceMethod =="Newmark")
    {
      uv0.push_back(U);
      uv0.push_back(V);
    }
    if (TimeAdvanceMethod =="BDF")
    {
        for ( UInt previousPass=0; previousPass < dataProblem->dataTimeAdvance()->orderBDF() ; previousPass++)
        {
            Real previousTimeStep = -previousPass*dt;
// <<<<<<< HEAD
            feSpace->interpolate(uexact, *U, previousTimeStep );
            uv0.push_back(U);
// =======
            // feSpace->interpolate( static_cast<FESpace_type::function_Type>( uexact ), *U, previousTimeStep );
            // uv0.push_back(*U);
        //>>>>>>> master
        }
    }

    //the uv0[0] should be the displacement
    //the uv0[1] should be the velocity

    //timeAdvance->setInitialCondition(uv0[0],uv0[1]);
    timeAdvance->setInitialCondition(uv0);

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

    timeAdvance->showMe();

    feSpace->interpolate( static_cast<FESpace_type::function_Type>( uexact ), *Exact , 0);
    feSpace->interpolate( static_cast<FESpace_type::function_Type>( v0 ), *vExact , 0);

    *U = timeAdvance->solution();
    *V = timeAdvance->firstDerivative();

    for (Real time = dt; time <= T; time += dt)
    {
        dataProblem->dataTime()->setTime(time);

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

        //evaluate rhs
        rhs *=0;
        timeAdvance->updateRHSContribution(dt);
        rhsV = timeAdvance->rhsContributionFirstDerivative();
        feSpace->l2ScalarProduct(source_in, rhs, time);
        rhs += problem.matrMass() *rhsV;

        //updateRHS
        problem.updateRHS(rhs);

        //solver system
        problem.iterate( bcH );    // Computes the matrices and solves the system

        //update unknowns of timeAdvance
        timeAdvance->shiftRight(*problem.solution());

        //evaluate uexact solution
        feSpace->interpolate( static_cast<FESpace_type::function_Type>( uexact ), *Exact , time);
        feSpace->interpolate( static_cast<FESpace_type::function_Type>( v0 ), *vExact , time);
        *U =  timeAdvance->solution();
        *V  =timeAdvance->firstDerivative();

        //postProcess
        exporter->postProcess(time);

        // Error L2 and H1 Norms
        AnalyticalSol uExact;
        vector_Type u (*problem.solution(), Repeated);

        Real H1_Error,  H1_RelError,  L2_Error, L2_RelError;

        L2_Error = feSpace->l2Error(uexact, u, time ,&L2_RelError);
        H1_Error = feSpace->h1Error(uExact, u, time ,&H1_RelError);

        //save the norm
        if (verbose)
        {
            out_norm.open("norm.txt", std::ios::app);
            out_norm << time             << "   "
            << L2_Error       << "   "
            << H1_Error      << "   "
            << L2_RelError << "   "
            << H1_RelError << "\n";
            out_norm.close();
        }

        MPI_Barrier(MPI_COMM_WORLD);
    }
    chrono.stop();
    if (verbose) std::cout << "Total iteration time " << chrono.diff() << " s." << std::endl;
}

Int main ( Int argc, char** argv )
{
    //! Initializing Epetra communicator
    MPI_Init (&argc, &argv);
    boost::shared_ptr<Epetra_Comm>  Comm ( new Epetra_MpiComm (MPI_COMM_WORLD) );
    if ( Comm->MyPID() == 0 )
    {
        std::cout << "% using MPI" << std::endl;
    }

    //********************************************//
    // Import parameters from an xml list. Use    //
    // Teuchos to create a list from a given file //
    // in the execution directory.                //
    //********************************************//

    std::cout << "Importing parameters list...";
    Teuchos::ParameterList FHNParameterList = * ( Teuchos::getParametersFromXmlFile ( "FitzHughNagumoParameters.xml" ) );
    std::cout << " Done!" << std::endl;


    //********************************************//
    // Creates a new model object representing the//
    // model from Fitz-Hugh Nagumo. The           //
    // model input are the parameters. Pass  the  //
    // parameter list in the constructor          //
    //********************************************//
    std::cout << "Building Constructor for Fitz-Hugh Nagumo Model with parameters ... ";
    IonicFitzHughNagumo  model ( FHNParameterList );
    std::cout << " Done!" << std::endl;


    //********************************************//
    // Show the parameters of the model as well as//
    // other informations  about the object.      //
    //********************************************//
    model.showMe();

    //********************************************//
    // Simulation starts on t=0 and ends on t=TF. //
    // The timestep is given by dt                //
    //********************************************//

    Real TF ( FHNParameterList.get ("TF", 300.0) );
    Real dt ( FHNParameterList.get ("dt", 0.01) );

    std::cout << "Time parameters : " << std::endl;
    std::cout << "TF = " << TF << std::endl;
    std::cout << "dt = " << dt << std::endl;

    std::string filename = "test_expeuler.txt";


    std::ofstream output (filename.c_str() );

    //********************************************//
    // Time loop starts.                          //
    //********************************************//
    std::cout << "Time loop starts...\n\n\n";

    LifeChrono chrono;
    chrono.start();

    Real valueToTest;

    valueToTest = EulerExplicit (dt, TF, model, FHNParameterList.get ("Iapp", 2000.0), output);


    chrono.stop();
    std::cout << "\n...Time loop ends.\n";
    std::cout << "\nElapsed time : " << chrono.diff() << std::endl;
    std::cout << "Solution written on file: " << filename << "\n";
    //********************************************//
    // Close exported file.                       //
    //********************************************//
    output.close();

    //! Finalizing Epetra communicator
    MPI_Finalize();

    Real returnValue;
    Real err = std::abs (valueToTest - SolutionTestNorm) / std::abs (SolutionTestNorm);
    std::cout << std::setprecision (20) << "\nError: " << err << "\nSolution norm: " << valueToTest << "\n";
    if ( err > 1e-12 )
    {
        returnValue = EXIT_FAILURE; // Norm of solution did not match
    }
    else
    {
        returnValue = EXIT_SUCCESS;
    }
    return ( returnValue );
}

void ExporterVTK<MeshType>::postProcess (const Real& time)
{
    // typedef std::list< ExporterData >::iterator Iterator;

    this->computePostfix();

    std::size_t found ( this->M_postfix.find ( "*" ) );
    if ( found == string::npos )
    {
        if (this->M_procId == 0)
        {
            std::cout << "  X-  ExporterVTK post-processing" << std::flush;
        }

        LifeChrono chrono;
        chrono.start();

        this->M_timeSteps.push_back (time);

        for (typename super::dataVectorIterator_Type iData = this->M_dataVector.begin();
                iData != this->M_dataVector.end(); ++iData)
        {
            std::ofstream vtkFile;
            std::stringstream buffer ("");

            // a unique PVTU file + a time collection is produced by the leader process
            if (this->M_procId == 0)
            {
                composePVTUStream (*iData, buffer);

                std::string vtkPFileName ( this->M_prefix + "_" + iData->variableName() +
                                           this->M_postfix + ".pvtu" );
                std::string vtkPFileNameWithDir (this->M_postDir + vtkPFileName);

                std::ofstream vtkPFile;
                vtkPFile.open ( vtkPFileNameWithDir.c_str() );
                ASSERT (vtkPFile.is_open(), "There is an error while opening " + vtkPFileName );
                ASSERT (vtkPFile.good(), "There is an error while writing to " + vtkPFileName );
                vtkPFile << buffer.str();
                vtkPFile.close();

                buffer.str ("");

                this->M_pvtuFiles[iData->variableName()].push_back (vtkPFileName);

            }


            // redundant. should be done just the first time
            std::map<UInt, UInt> ltgPointsMap, gtlPointsMap;
            std::vector<Vector> pointsCoords;
            createPointsMaps ( iData->feSpacePtr(), gtlPointsMap, ltgPointsMap, pointsCoords );

            composeVTUHeaderStream ( gtlPointsMap.size(), buffer );
            composeVTUGeoStream ( iData->feSpacePtr(), gtlPointsMap, pointsCoords, buffer );

            composeTypeDataHeaderStream (iData->where(), buffer);
            composeDataArrayStream (*iData, ltgPointsMap, buffer);
            composeTypeDataFooterStream (iData->where(), buffer);

            composeVTUFooterStream ( buffer );

            // each process writes its own file
            std::string filename ( this->M_postDir + this->M_prefix + "_" + iData->variableName() +
                                   this->M_postfix + "." + this->M_procId + ".vtu" );
            vtkFile.open ( filename.c_str() );
            ASSERT (vtkFile.is_open(), "There is an error while opening " + filename );
            ASSERT (vtkFile.good(), "There is an error while writing to " + filename );
            vtkFile << buffer.str();
            vtkFile.close();

            buffer.str ("");

        }
        chrono.stop();
        if (!this->M_procId)
        {
            std::cout << "...done in " << chrono.diff() << " s." << std::endl;
        }
    }
}

//.........这里部分代码省略.........
    //bcFunction_Type constantAreaFunction( boost::bind( &Constant::operator(), &constantArea, _1 ) );

    //Constant constantPressure( 24695.0765959599 );
    //bcFunction_Type constantPressureFunction( boost::bind( &Constant::operator(), &constantPressure, _1 ) );

    oneDModel.bc().setBC ( OneDFSI::left,  OneDFSI::first, OneDFSI::Q,  sinusoidalFunction  );
    oneDModel.bc().setBC ( OneDFSI::right, OneDFSI::first, OneDFSI::W2, absorbingFunction );

    //oneDModel.bc().setBC( OneDFSI::right, OneDFSI::first, OneDFSI::A, constantAreaFunction );
    //oneDModel.bc().setBC( OneDFSI::right, OneDFSI::first, OneDFSI::P,   constantPressureFunction );

    oneDModel.setupModel();

    absorbing->setSolution ( oneDModel.solution() );
    absorbing->setFluxSource ( oneDModel.flux(), oneDModel.source() );

    // *********************************
    // Tempolar loop
    // *********************************
    std::cout << "\nTemporal loop:" << std::endl;

    LifeChrono chronoTotal;
    LifeChrono chronoSystem;
    LifeChrono chronoIteration;

    Int count = 0;
    chronoTotal.start();
    for ( ; oneDModel.data().dataTime()->canAdvance() ; oneDModel.data().dataTime()->updateTime(), ++count )
    {
        std::cout << std::endl << "--------- Iteration " << count << " time = " << oneDModel.data().dataTime()->time() << std::endl;

        chronoIteration.start();

        if ( oneDModel.data().dataTime()->isFirstTimeStep() )
        {
            oneDModel.buildModel();
        }
        else
        {
            oneDModel.updateModel();
        }

        chronoSystem.start();

        oneDModel.solveModel();

        chronoSystem.stop();

        oneDModel.updateSolution();

        //Save solution
        if ( count % 50 == 0 || oneDModel.data().dataTime()->isLastTimeStep() )
        {
            oneDModel.saveSolution();
        }

        chronoIteration.stop();

        std::cout << " System solved in " << chronoSystem.diff() << " s, (total time " << chronoIteration.diff() << " s)." << std::endl;
    }

    chronoTotal.stop();
    std::cout << std::endl << " Simulation ended successfully in " << chronoTotal.diff()  << " s" << std::endl;

#ifdef HAVE_MPI
    std::cout << "MPI Finalization" << std::endl;
    MPI_Finalize();
#endif

    if ( check )
    {
        bool ok = true;
        Int rightNodeID = oneDModel.solver()->boundaryDOF ( OneDFSI::right );

        ok = ok && checkValue ( 0.999998  , (*oneDModel.solution ("A") ) [rightNodeID]);
        ok = ok && checkValue (-0.00138076, (*oneDModel.solution ("Q") ) [rightNodeID]);
        ok = ok && checkValue (-0.00276153, (*oneDModel.solution ("W1") ) [rightNodeID]);
        ok = ok && checkValue ( 0.00000000, (*oneDModel.solution ("W2") ) [rightNodeID]);

        ok = ok && checkValue ( 0.999999  , (*oneDModel.solution ("A") ) [rightNodeID - 1]);
        ok = ok && checkValue (-0.00040393, (*oneDModel.solution ("Q") ) [rightNodeID - 1]);
        ok = ok && checkValue (-0.00080833, (*oneDModel.solution ("W1") ) [rightNodeID - 1]);
        ok = ok && checkValue ( 0.00000045, (*oneDModel.solution ("W2") ) [rightNodeID - 1]);

        if (ok)
        {
            std::cout << "Test successful" << std::endl;
            return 0;
        }
        else
        {
            std::cout << "Test unsuccessful" << std::endl;
            return -1;
        }
    }
    else
    {
        return 0;
    }
}

//.........这里部分代码省略.........
    this->M_displayer->leaderPrint (" \n*********************************\n  ");
    this->M_displayer->leaderPrint ("   Performing the analysis recovering the displacement..., ", this->M_dataMaterial->solidType() );
    this->M_displayer->leaderPrint (" \n*********************************\n  ");

    solutionVectPtr_Type grDisplX ( new solutionVect_Type (* (this->M_FESpace->mapPtr() ) ) );
    solutionVectPtr_Type grDisplY ( new solutionVect_Type (* (this->M_FESpace->mapPtr() ) ) );
    solutionVectPtr_Type grDisplZ ( new solutionVect_Type (* (this->M_FESpace->mapPtr() ) ) );

    //Compute the deformation gradient tensor F of the displ field
    super::computeDisplacementGradient ( grDisplX, grDisplY, grDisplZ);

    solutionVect_Type grXRep (*grDisplX, Repeated);
    solutionVect_Type grYRep (*grDisplY, Repeated);
    solutionVect_Type grZRep (*grDisplZ, Repeated);
    solutionVect_Type dRep (* (this->M_displacement), Repeated);

    //For each of the DOF, the Cauchy tensor is computed.
    //Therefore the tensor C,P, \sigma are computed for each DOF

    chrono.start();

    for ( UInt i (0); i < this->M_FESpace->mesh()->numVolumes(); ++i )
    {

        //Setup quantities
        this->M_FESpace->fe().updateFirstDerivQuadPt ( this->M_FESpace->mesh()->volumeList ( i ) );
        UInt eleID = this->M_FESpace->fe().currentLocalId();
        this->M_marker = this->M_FESpace->mesh()->volumeList ( i ).markerID();

        //store the local \grad(u)
        matrix_Type gradientDispl ( this->M_FESpace->fieldDim(), this->M_FESpace->fieldDim() );
        gradientDispl.Scale ( 0.0 );

        //Extracting the local displacement
        for ( UInt iNode = 0; iNode < ( UInt ) this->M_FESpace->fe().nbFEDof(); iNode++ )
        {
            UInt  iloc = this->M_FESpace->fe().patternFirst ( iNode );

            for ( UInt iComp = 0; iComp < this->M_FESpace->fieldDim(); ++iComp )
            {
                UInt ig = this->M_FESpace->dof().localToGlobalMap ( eleID, iloc ) + iComp * this->M_FESpace->dim() + this->M_offset;

                gradientDispl (iComp, 0) = grXRep[ig];
                gradientDispl (iComp, 1) = grYRep[ig];
                gradientDispl (iComp, 2) = grZRep[ig];
            }

            //Reinitialization of matrices and arrays
            ( * (this->M_cofactorF) ).Scale (0.0);
            ( * (this->M_firstPiola) ).Scale (0.0);
            ( * (this->M_sigma) ).Scale (0.0);

            //Moving to cylindrical coordinates
            moveToCylindricalCoordinates ( gradientDispl, iloc, *M_deformationCylindricalF );

            //Compute the rightCauchyC tensor
            AssemblyElementalStructure::computeInvariantsRightCauchyGreenTensor (this->M_invariants, *M_deformationCylindricalF, * (this->M_cofactorF) );

            //Compute the first Piola-Kirchhoff tensor
            this->M_material->computeLocalFirstPiolaKirchhoffTensor (* (this->M_firstPiola), *M_deformationCylindricalF, * (this->M_cofactorF), this->M_invariants, this->M_marker);

            //Compute the Cauchy tensor
            AssemblyElementalStructure::computeCauchyStressTensor (* (this->M_sigma), * (this->M_firstPiola), this->M_invariants[3], *M_deformationCylindricalF);

            //Compute the eigenvalue
            AssemblyElementalStructure::computeEigenvalues (* (this->M_sigma), this->M_eigenvaluesR, this->M_eigenvaluesI);

            //The Cauchy tensor is symmetric and therefore, the eigenvalues are real
            //Check on the imaginary part of eigen values given by the Lapack method
            Real sum (0);
            for ( UInt j (0); j < this->M_eigenvaluesI.size(); ++j )
            {
                sum += std::abs ( this->M_eigenvaluesI[j] );
            }

            ASSERT ( sum < 1e-6 , "The eigenvalues of the Cauchy stress tensors have to be real!" );

            std::sort ( this->M_eigenvaluesR.begin(), this->M_eigenvaluesR.end() );

            std::cout << "Saving eigenvalues" << i << std::endl;

            //Save the eigenvalues in the global vector
            for ( UInt icoor = 0; icoor < this->M_FESpace->fieldDim(); ++icoor )
            {
                UInt ig = this->M_FESpace->dof().localToGlobalMap ( eleID, iloc ) + icoor * this->M_FESpace->dim() + this->M_offset;
                (* (this->M_globalEigenvalues) ) (ig) = this->M_eigenvaluesR[icoor];
                // Int LIDid = this->M_displacement->blockMap().LID(ig);
                // if( this->M_globalEigenvalues->blockMap().LID(ig) != -1  )
                // {
                //     Int GIDid = this->M_globalEigenvalues->blockMap().GID(LIDid);

                // }
            }
        }
    }

    chrono.stop();
    this->M_displayer->leaderPrint ("Analysis done in: ", chrono.diff() );

}