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

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 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;
    }

}

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
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 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;
}

//.........这里部分代码省略.........
    // Absorbing
    bc_Type::bcFunctionSolverDefinedPtr_Type absorbing ( new OneDFSIFunctionSolverDefinedAbsorbing ( OneDFSI::right, OneDFSI::W2 ) );
    bcFunction_Type absorbingFunction ( boost::bind ( &OneDFSIFunctionSolverDefinedAbsorbing::operator(),
                                                      dynamic_cast<OneDFSIFunctionSolverDefinedAbsorbing*> ( & ( *absorbing ) ), _1, _2 ) );

    // BC to test A_from_P conversion
    //Constant constantArea( 1.00 );
    //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;

void Heart::computeRhs ( vector_Type& rhs,
                         HeartBidomainSolver< mesh_Type >& electricModel,
                         boost::shared_ptr< HeartIonicSolver< mesh_Type > > ionicModel,
                         HeartBidomainData& data )
{
    bool verbose = (M_heart_fct->M_comm->MyPID() == 0);
    if (verbose)
    {
        std::cout << "  f-  Computing Rhs ...        " << "\n" << std::flush;
    }
    LifeChrono chrono;
    chrono.start();

    //! u, w with repeated map
    vector_Type uVecRep (electricModel.solutionTransmembranePotential(), Repeated);
    ionicModel->updateRepeated();

    VectorElemental elvec_Iapp ( electricModel.potentialFESpace().fe().nbFEDof(), 2 ),
                    elvec_u ( electricModel.potentialFESpace().fe().nbFEDof(), 1 ),
                    elvec_Iion ( electricModel.potentialFESpace().fe().nbFEDof(), 1 );
    for (UInt iVol = 0; iVol < electricModel.potentialFESpace().mesh()->numVolumes(); ++iVol)
    {
        electricModel.potentialFESpace().fe().updateJacQuadPt ( electricModel.potentialFESpace().mesh()->volumeList ( iVol ) );
        elvec_u.zero();
        elvec_Iion.zero();
        elvec_Iapp.zero();

        UInt eleIDu = electricModel.potentialFESpace().fe().currentLocalId();
        UInt nbNode = ( UInt ) electricModel.potentialFESpace().fe().nbFEDof();
        for ( UInt iNode = 0 ; iNode < nbNode ; iNode++ )
        {
            Int ig = electricModel.potentialFESpace().dof().localToGlobalMap ( eleIDu, iNode );
            elvec_u.vec() [ iNode ] = uVecRep[ig];
        }

        UInt eleID = electricModel.potentialFESpace().fe().currentLocalId();
        ionicModel->updateElementSolution (eleID);
        ionicModel->computeIonicCurrent (data.membraneCapacitance(), elvec_Iion, elvec_u, electricModel.potentialFESpace() );

        //! Computing Iapp
        AssemblyElemental::source (M_heart_fct->stimulus(),
                                   elvec_Iapp,
                                   electricModel.potentialFESpace().fe(),
                                   data.time(), 0);
        AssemblyElemental::source (M_heart_fct->stimulus(),
                                   elvec_Iapp,
                                   electricModel.potentialFESpace().fe(),
                                   data.time(),
                                   1);
        UInt totalUDof  = electricModel.potentialFESpace().map().map (Unique)->NumGlobalElements();

        for ( UInt iNode = 0 ; iNode < nbNode ; iNode++ )
        {
            Int ig = electricModel.potentialFESpace().dof().localToGlobalMap ( eleIDu, iNode );
            rhs.sumIntoGlobalValues (ig, elvec_Iapp.vec() [iNode] +
                                     data.volumeSurfaceRatio() * elvec_Iion.vec() [iNode] );
            rhs.sumIntoGlobalValues (ig + totalUDof,
                                     -elvec_Iapp.vec() [iNode + nbNode] -
                                     data.volumeSurfaceRatio() * elvec_Iion.vec() [iNode] );
        }
    }
    rhs.globalAssemble();

    rhs += electricModel.matrMass() * data.volumeSurfaceRatio() *
           data.membraneCapacitance() * electricModel.BDFIntraExtraPotential().time_der (data.timeStep() );

    MPI_Barrier (MPI_COMM_WORLD);

    chrono.stop();
    if (verbose)
    {
        std::cout << "done in " << chrono.diff() << " s." << std::endl;
    }
}

Real
darcy_nonlinear::run()
{
    using namespace dataProblem;

    // Life chonos
    LifeChrono chronoTotal;
    LifeChrono chronoReadAndPartitionMesh;
    LifeChrono chronoBoundaryCondition;
    LifeChrono chronoFiniteElementSpace;
    LifeChrono chronoProblem;
    LifeChrono chronoProcess;
    LifeChrono chronoError;

    // Displayer to print on screen
    boost::shared_ptr < Displayer > displayer ( new Displayer ( Members->comm ) );

    // Start chronoTotal for measure the total time for the computation
    chronoTotal.start();

    // Reading from data file
    GetPot dataFile ( Members->data_file_name.c_str() );

    //
    // The Darcy Solver
    //

    displayer->leaderPrint ( "The Darcy solver\n" );

    // Start chronoReadAndPartitionMesh for measure the total time for the creation of the local meshes
    chronoReadAndPartitionMesh.start();

    // Create the data file
    darcyDataPtr_Type darcyData ( new darcyData_Type );

    // Set up the data
    darcyData->setup ( dataFile );

    // Create a Teuchos parameter list for the linear algebra.
    darcyData_Type::paramListPtr_Type linearAlgebraList = Teuchos::rcp ( new darcyData_Type::paramList_Type );
    linearAlgebraList = Teuchos::getParametersFromXmlFile ( Members->xml_file_name );

    // Set the parameter list into the data darcy.
    darcyData->setLinearAlgebraList ( linearAlgebraList );

    // Create the mesh file handler
    MeshData meshData;

    // Set up the mesh file
    meshData.setup ( dataFile,  Members->discretization_section + "/space_discretization");

    // Create the the mesh
    regionMeshPtr_Type fullMeshPtr ( new regionMesh_Type ( Members->comm ) );

    // Select if the mesh is structured or not
    if ( meshData.meshType() != "structured" )
    {
        // Set up the mesh
        readMesh ( *fullMeshPtr, meshData );
    }
    else
    {
        // Section of the structured mesh
        const std::string structuredSection = Members->discretization_section + "/space_discretization/";

        // Set up the structured mesh
        regularMesh2D ( *fullMeshPtr, 0,
                        dataFile ( ( structuredSection + "nx" ).data(), 4 ),
                        dataFile ( ( structuredSection + "ny" ).data(), 4 ),
                        dataFile ( ( structuredSection + "verbose" ).data(), false ),
                        dataFile ( ( structuredSection + "lx" ).data(), 1. ),
                        dataFile ( ( structuredSection + "ly" ).data(), 1. ) );
    }

    // Create the local mesh ( local scope used to delete the meshPart object )
    regionMeshPtr_Type meshPtr;
    {
        // Create the partitioner
        MeshPartitioner< regionMesh_Type >  meshPart;

        // Partition the mesh using ParMetis
        meshPart.doPartition ( fullMeshPtr, Members->comm );

        // Get the local mesh
        meshPtr = meshPart.meshPartition();
    }

    // Stop chronoReadAndPartitionMesh
    chronoReadAndPartitionMesh.stop();

    // The leader process print chronoReadAndPartitionMesh
    displayer->leaderPrint ( "Time for read and partition the mesh ",
                             chronoReadAndPartitionMesh.diff(), "\n" );

    // Create the boundary conditions

    // Start chronoBoundaryCondition for measure the total time for create the boundary conditions
    chronoBoundaryCondition.start();

    bcHandlerPtr_Type bcDarcy ( new bcHandler_Type );
//.........这里部分代码省略.........

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

    LifeChrono chrono;

    this->M_displayer->leaderPrint (" \n*********************************\n  ");
    this->M_displayer->leaderPrint ("   Performing the analysis recovering the tensions..., ", this->M_dataMaterial->solidType() );
    this->M_displayer->leaderPrint (" \n*********************************\n  ");

    solutionVect_Type patchArea (* (this->M_displacement), Unique, Add);
    patchArea *= 0.0;

    super::constructPatchAreaVector ( patchArea );

    //Before assembling the reconstruction process is done
    solutionVect_Type patchAreaR (patchArea, Repeated);

    QuadratureRule fakeQuadratureRule;

    Real refElemArea (0); //area of reference element
    //compute the area of reference element
    for (UInt iq = 0; iq < this->M_FESpace->qr().nbQuadPt(); iq++)
    {
        refElemArea += this->M_FESpace->qr().weight (iq);
    }

    Real wQuad (refElemArea / this->M_FESpace->refFE().nbDof() );

    //Setting the quadrature Points = DOFs of the element and weight = 1
    std::vector<GeoVector> coords = this->M_FESpace->refFE().refCoor();
    std::vector<Real> weights (this->M_FESpace->fe().nbFEDof(), wQuad);
    fakeQuadratureRule.setDimensionShape ( shapeDimension (this->M_FESpace->refFE().shape() ), this->M_FESpace->refFE().shape() );
    fakeQuadratureRule.setPoints (coords, weights);

    //Set the new quadrature rule
    this->M_FESpace->setQuadRule (fakeQuadratureRule);

    UInt totalDof = this->M_FESpace->dof().numTotalDof();
    VectorElemental dk_loc (this->M_FESpace->fe().nbFEDof(), this->M_FESpace->fieldDim() );

    //Vectors for the deformation tensor
    std::vector<matrix_Type> vectorDeformationF (this->M_FESpace->fe().nbFEDof(), * (this->M_deformationF) );
    //Copying the displacement field into a vector with repeated map for parallel computations
    solutionVect_Type dRep (* (this->M_displacement), Repeated);

    VectorElemental elVecTens (this->M_FESpace->fe().nbFEDof(), this->M_FESpace->fieldDim() );

    chrono.start();

    //Loop on each volume
    for ( UInt i = 0; i < this->M_FESpace->mesh()->numVolumes(); ++i )
    {
        this->M_FESpace->fe().updateFirstDerivQuadPt ( this->M_FESpace->mesh()->volumeList ( i ) );
        elVecTens.zero();

        this->M_marker = this->M_FESpace->mesh()->volumeList ( i ).markerID();

        UInt eleID = this->M_FESpace->fe().currentLocalId();

        //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;
                dk_loc[iloc + iComp * this->M_FESpace->fe().nbFEDof()] = dRep[ig];
            }
        }

        //Compute the element tensor F
        AssemblyElementalStructure::computeLocalDeformationGradientWithoutIdentity ( dk_loc, vectorDeformationF, this->M_FESpace->fe() );

        //Compute the local vector of the principal stresses
        for ( UInt nDOF = 0; nDOF < ( UInt ) this->M_FESpace->fe().nbFEDof(); nDOF++ )
        {
            UInt  iloc = this->M_FESpace->fe().patternFirst ( nDOF );
            vector_Type localDisplacement (this->M_FESpace->fieldDim(), 0.0);

            for ( UInt coor = 0; coor < this->M_FESpace->fieldDim(); coor++ )
            {
                localDisplacement[coor] = iloc + coor * this->M_FESpace->fe().nbFEDof();
            }

            this->M_sigma->Scale (0.0);
            this->M_firstPiola->Scale (0.0);
            this->M_cofactorF->Scale (0.0);
            M_deformationCylindricalF->Scale (0.0);

            moveToCylindricalCoordinates (vectorDeformationF[nDOF], 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
//.........这里部分代码省略.........

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

    LifeChrono chrono;

    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);

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

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;
        }
    }
}