// ===================================================
// Methods
// ===================================================
void
BCInterfaceFunctionSolverDefined< BCHandler, FSIOperator >::exportData ( dataPtr_Type& data )
{

#ifdef HAVE_LIFEV_DEBUG
    debugStream ( 5025 ) << "BCInterfaceFunctionSolverDefined::exportData" << "\n";
#endif

    data->setName ( M_name );
    data->setFlag ( M_flag );
    data->setType ( M_type );
    data->setMode ( M_mode );
    data->setComponentsVector ( M_componentsVector );
}

void
BCInterfaceFunctionParserFileSolver< BcHandlerType, PhysicalSolverType >::setData ( const dataPtr_Type& data )
{

#ifdef HAVE_LIFEV_DEBUG
    debugStream ( 5024 ) << "BCInterfaceFunctionFileSolver::setData" << "\n";
#endif

    functionParserFile_Type::setData ( data );

    functionParserSolver_Type::M_boundaryID = data->boundaryID();

    functionParserSolver_Type::createAccessList ( data );
}

void 
Combinations< _T >::
getNextCombination( const dataPtr_Type & dataPtr )
{
    assert( dataPtr->size() == this->M_K );
    
    CombinationsID::tuple_Type tupleIDs;
    
    this->getNextCombinationID( tupleIDs );
    
    for ( UInt iK(0); iK < this->M_K; ++iK )
    {
	(*dataPtr)[ iK ] = (*this->M_dataPtr)[ tupleIDs[ iK ] ];
    }
    
    return;
}

// ===================================================
// Set Methods
// ===================================================
void
BCInterfaceFunctionSolverDefined< BCHandler, StructuralOperator<RegionMesh <LinearTetra> > >::setData ( const dataPtr_Type& data )
{

#ifdef HAVE_LIFEV_DEBUG
    debugStream ( 5025 ) << "BCInterfaceFunctionSolverDefined::setData" << "\n";
#endif

    //Set mapFunction
    std::map< std::string, Solid3DFunction > mapFunction;
    mapFunction["RobinWall"] = RobinWall;

    // Retrieving the strings
    M_solid3DFunction = mapFunction[ data->baseString() ];

    M_name = data->name();
    M_flag = data->flag();
    M_type = data->type();
    M_mode = data->mode();
    M_componentsVector = data->componentsVector();

    if ( M_solid3DFunction == RobinWall )
    {
        factory_Type factory;
        M_vectorFunctionRobin.reserve (2);
        dataPtr_Type temporaryData ( new data_Type ( *data ) );

        // Create the mass term function
        temporaryData->setRobinBaseAlpha();
        M_vectorFunctionRobin.push_back ( factory.createFunctionParser ( temporaryData ) );

        // Create the RHS
        temporaryData->setRobinBaseBeta();
        M_vectorFunctionRobin.push_back ( factory.createFunctionParser ( temporaryData ) );
    }
}

void NeoHookeanMaterialNonLinear<Mesh>::computeStiffness ( const vector_Type&       sol,
                                                           Real                     /*factor*/,
                                                           const dataPtr_Type&      dataMaterial,
                                                           const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
                                                           const displayerPtr_Type& displayer )
{
    this->M_stiff.reset (new vector_Type (*this->M_localMap) );

    displayer->leaderPrint (" \n******************************************************************\n  ");
    displayer->leaderPrint (" Non-Linear S-  Computing the Neo-Hookean nonlinear stiffness vector"     );
    displayer->leaderPrint (" \n******************************************************************\n  ");

    UInt totalDof   = this->M_FESpace->dof().numTotalDof();
    UInt dim        = this->M_FESpace->dim();

    VectorElemental dk_loc ( this->M_FESpace->fe().nbFEDof(), nDimensions );
    //vector_Type disp(sol);

    vector_Type dRep (sol, Repeated);

    mapIterator_Type it;

    for ( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++ )
    {

        //Given the marker pointed by the iterator, let's extract the material parameters
        UInt marker = it->first;

        Real mu     = dataMaterial->mu (marker);
        Real bulk   = dataMaterial->bulk (marker);

        for ( UInt j (0); j < it->second.size(); j++ )
        {
            this->M_FESpace->fe().updateFirstDerivQuadPt ( * (it->second[j]) );

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

            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 < nDimensions; ++iComp )
                {
                    UInt ig = this->M_FESpace->dof().localToGlobalMap ( eleID, iloc ) + iComp * dim + this->M_offset;
                    dk_loc[ iloc + iComp * this->M_FESpace->fe().nbFEDof() ] = dRep[ig];
                }
            }

            this->M_elvecK->zero();

            computeKinematicsVariables ( dk_loc );

            //! Stiffness for non-linear terms of the Neo-Hookean model
            /*!
              The results of the integrals are stored at each step into elvecK, until to build K matrix of the bilinear form
            */
            //! Volumetric part
            /*!
              Source term Pvol: int { bulk /2* (J1^2 - J1  + log(J1) ) * 1/J1 * (CofF1 : \nabla v) }
            */
            AssemblyElementalStructure::source_Pvol ( 0.5 * bulk, (*M_CofFk), (*M_Jack),
                                                      *this->M_elvecK,  this->M_FESpace->fe() );

            //! Isochoric part
            /*!
              Source term P1iso_NH
            */
            AssemblyElementalStructure::source_P1iso_NH ( mu, (*M_CofFk) , (*M_Fk),  (*M_Jack),  (*M_trCisok) ,
                                                          *this->M_elvecK,  this->M_FESpace->fe() );

            for ( UInt ic = 0; ic < nDimensions; ++ic )
            {
                /*!
                  M_elvecK is assemble into *vec_stiff vector that is recall
                  from updateSystem(matrix_ptrtype& mat_stiff, vector_ptr_type& vec_stiff)
                */
                assembleVector ( *this->M_stiff,
                                 *this->M_elvecK,
                                 this->M_FESpace->fe(),
                                 this->M_FESpace->dof(),
                                 ic, this->M_offset +  ic * totalDof );
            }
        }
    }

    this->M_stiff->globalAssemble();
}

void VenantKirchhoffMaterialLinear<Mesh>::computeLinearStiff(dataPtr_Type& dataMaterial,
                                                             const mapMarkerVolumesPtr_Type mapsMarkerVolumes)
{
    //  std::cout<<"compute LinearStiff Matrix start\n";

    UInt totalDof = this->M_FESpace->dof().numTotalDof();
    // Number of displacement components
    UInt nc = nDimensions;

    //Compute the linear part of the Stiffness Matrix.
    //In the case of Linear Material it is the Stiffness Matrix.
    //In the case of NonLinear Materials it must be added of the non linear part.

    mapIterator_Type it;

    for( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++ )
    {

        //Given the marker pointed by the iterator, let's extract the material parameters
        UInt marker = it->first;

        Real mu = dataMaterial->mu(marker);
        Real lambda = dataMaterial->lambda(marker);

        //Given the parameters I loop over the volumes with that marker
        for ( UInt j(0); j < it->second.size(); j++ )
        {
            this->M_FESpace->fe().updateFirstDerivQuadPt( *(it->second[j]) );

            this->M_elmatK->zero();

            //These methods are implemented in AssemblyElemental.cpp
            //They have been kept in AssemblyElemental in order to avoid repetitions
            stiff_strain( 2*mu, *this->M_elmatK, this->M_FESpace->fe() );// here in the previous version was 1. (instead of 2.)
            stiff_div   ( lambda, *this->M_elmatK, this->M_FESpace->fe() );// here in the previous version was 0.5 (instead of 1.)

            //this->M_elmatK->showMe();

            // assembling
            for ( UInt ic = 0; ic < nc; ic++ )
            {
                for ( UInt jc = 0; jc < nc; jc++ )
                {
                    assembleMatrix( *this->M_linearStiff,
                                    *this->M_elmatK,
                                    this->M_FESpace->fe(),
                                    this->M_FESpace->fe(),
                                    this->M_FESpace->dof(),
                                    this->M_FESpace->dof(),
                                    ic,  jc,
                                    this->M_offset +ic*totalDof, this->M_offset + jc*totalDof );

                }
            }


        }

    }

    this->M_linearStiff->globalAssemble();

    //Initialization of the pointer M_stiff to what is pointed by M_linearStiff
    this->M_stiff = this->M_linearStiff;
    //   std::cout<<"compute LinearStiff Matrix end\n";
    this->M_jacobian = this->M_linearStiff;
}

// ===================================================
// Methods
// ===================================================
void
FSISolver::setData( const dataPtr_Type& data )
{
    M_data = data;

    int rank, numtasks;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &numtasks);

    bool fluid = false;
    bool solid = false;

    int  fluidLeader(0);
    int  solidLeader(0);

    if ( ( data->method().compare("monolithicGE") && data->method().compare("monolithicGI") ) )
    {
        MPI_Group  originGroup, newGroup;
        MPI_Comm_group(MPI_COMM_WORLD, &originGroup);

        if ( numtasks == 1 )
        {
            std::cout << "Serial Fluid/Structure computation" << std::endl;
            fluid = true;
            solid = true;
            solidLeader = 0;
            fluidLeader = solidLeader;

            M_epetraWorldComm.reset( new Epetra_MpiComm(MPI_COMM_WORLD));
            M_epetraComm = M_epetraWorldComm;
        }
        else
        {
            std::vector<int> members(numtasks);

            solidLeader = 0;
            fluidLeader = 1-solidLeader;

            if (rank == solidLeader)
            {
                members[0] = solidLeader;
                /* int ierr = */
                MPI_Group_incl(originGroup, 1, &members[0], &newGroup);
                solid = true;
            }
            else
            {
                for (Int ii = 0; ii <= numtasks; ++ii)
                {
                    if ( ii < solidLeader)
                        members[ii] = ii;
                    else if ( ii > solidLeader)
                        members[ii - 1] = ii;
                }

                /* int ierr = */ MPI_Group_incl(originGroup, numtasks - 1, &members[0], &newGroup);
                fluid = true;
            }

            MPI_Comm* localComm = new MPI_Comm;
            MPI_Comm_create(MPI_COMM_WORLD, newGroup, localComm);
            M_localComm.reset(localComm);

            M_epetraComm.reset(new Epetra_MpiComm(*M_localComm.get()));
            M_epetraWorldComm.reset(new Epetra_MpiComm(MPI_COMM_WORLD));
        }
    }
    else // Monolithic or FullMonolithic
    {
        fluid = true;
        solid = true;
        solidLeader = 0;
        fluidLeader = solidLeader;

        M_epetraWorldComm.reset( new Epetra_MpiComm(MPI_COMM_WORLD));
        M_epetraComm = M_epetraWorldComm;
    }

#ifdef DEBUG
    if ( fluid )
    {
        debugStream(6220) << M_epetraComm->MyPID()
        << " ( " << rank << " ) "
        << " out of " << M_epetraComm->NumProc()
        << " ( " << numtasks << " ) "
        << " is fluid." << std::endl;
    }
    if ( solid )
    {
        debugStream(6220) << M_epetraComm->MyPID()
        << " ( " << rank << " ) "
        << " out of " << M_epetraComm->NumProc()
        << " ( " << numtasks << " ) "
        << " is solid." << std::endl;
    }
#endif

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

inline void
BCInterfaceFunctionParserFile< BcHandlerType, PhysicalSolverType >::setData ( const dataPtr_Type& data )
{

#ifdef HAVE_LIFEV_DEBUG
    debugStream ( 5022 ) << "BCInterfaceFunctionFile::loadData            fileName: " << data->baseString() << "\n";
#endif

    std::vector< std::string > stringsVector;
    boost::split ( stringsVector, data->baseString(), boost::is_any_of ( "[" ) );

    //Load data from file
    GetPot dataFile ( stringsVector[0] );

    //Set variables
    UInt variablesNumber = dataFile.vector_variable_size ( "variables" );

    M_variables.clear();
    M_variables.reserve ( variablesNumber );

    std::vector< Real > scale;
    scale.reserve ( variablesNumber );

    for ( UInt j ( 0 ); j < variablesNumber; ++j )
    {
        M_variables.push_back ( dataFile ( "variables", "unknown", j ) );
        scale.push_back ( dataFile ( "scale", 1.0, j ) );
    }

#ifdef HAVE_LIFEV_DEBUG
    std::stringstream output;
    output << "BCInterfaceFunctionFile::loadData           variables: ";
    for ( UInt j (0); j < variablesNumber; ++j )
    {
        output << M_variables[j] << "  ";
    }

    output << "\n                                                           scale: ";
    for ( UInt j (0); j < variablesNumber; ++j )
    {
        output << scale[j] << "  ";
    }

    debugStream ( 5022 ) << output.str() << "\n";
#endif

    //Load loop flag
    M_loop = dataFile ( "loop", false );

    //Load data
    UInt dataLines = dataFile.vector_variable_size ( "data" ) / variablesNumber;

    M_data.clear();
    for ( UInt j ( 0 ); j < variablesNumber; ++j )
    {
        M_data[M_variables[j]].reserve ( dataLines );
    }

    for ( UInt i ( 0 ); i < dataLines; ++i )
        for ( UInt j ( 0 ); j < variablesNumber; ++j )
        {
            M_data[M_variables[j]].push_back ( scale[j] * dataFile ( "data", 0.0, i * variablesNumber + j ) );
        }

#ifdef HAVE_LIFEV_DEBUG
    output.str ("");
    output << "                                                 loop: " << M_loop << "\n";
    output << "                                                 data:";
    for ( UInt i (0); i < dataLines; ++i )
    {
        if (i > 0)
        {
            output << "                                                                 ";
        }

        for ( UInt j (0); j < variablesNumber; ++j )
        {
            output << " " << M_data[ M_variables[j] ][i];
        }
        output << "\n";
    }
    debugStream ( 5022 ) << output.str();
#endif

    //Initialize iterator
    M_dataIterator = M_data[M_variables[0]].begin();

    //Update the data container (IT IS A COPY!) with the correct base string for the BCInterfaceFunctionParser
    if ( stringsVector.size() < 2 )
    {
        data->setBaseString ( dataFile ( "function", "Undefined" ) );
    }
    else
    {
        boost::replace_all ( stringsVector[1], "]", "" );
        data->setBaseString ( dataFile ( ( "function" + stringsVector[1] ).c_str(), "Undefined" ) );
    }

    // Now data contains the real base string
    functionParser_Type::setData ( data );
//.........这里部分代码省略.........

// ===================================================
// Inline conversion methods
// ===================================================
inline Real
OneDFSIPhysics::fromPToA( const Real& P, const Real& timeStep, const UInt& iNode, const bool& elasticExternalNodes ) const
{
    if ( !M_dataPtr->viscoelasticWall() || ( ( iNode == 0 || iNode == M_dataPtr->numberOfNodes() - 1 ) && elasticExternalNodes ) )
        return ( M_dataPtr->area0( iNode ) * OneDFSI::pow20( ( P - externalPressure() ) / M_dataPtr->beta0( iNode ) + 1, 1 / M_dataPtr->beta1( iNode ) )  );
    else
    {
        // Newton method to solve the non linear equation
        Real tolerance(1e-6);
        Real maxIT(100);
        UInt i(0);

        Real A( M_dataPtr->area0( iNode ) );
        Real newtonUpdate(0);
        for ( ; i < maxIT ; ++i )
        {
            if ( std::abs( pressure( A, timeStep, iNode, elasticExternalNodes ) - P ) < tolerance )
                break;

            newtonUpdate = ( pressure( A, timeStep, iNode, elasticExternalNodes ) - P ) / dPdA( A, timeStep, iNode, elasticExternalNodes );
            if ( A - newtonUpdate <= 0 )
                A /= 2.0; // Bisection
            else
                A -= newtonUpdate; // Newton
        }
        if ( i == maxIT )
        {
            std::cout << "!!! Warning: conversion fromPToA below tolerance !!! " << std::endl;
            std::cout << "Tolerance: " << tolerance << "; Residual: " << std::abs( pressure( A, timeStep, iNode, elasticExternalNodes ) - P ) << std::endl;
        }

        return A;
    }
}

inline Real
OneDFSIPhysics::viscoelasticPressure( const Real& A, const Real& timeStep, const UInt& iNode, const bool& elasticExternalNodes ) const
{
    if ( !M_dataPtr->viscoelasticWall() || ( ( iNode == 0 || iNode == M_dataPtr->numberOfNodes() - 1 ) && elasticExternalNodes ) )
        return 0;
    else
        return M_dataPtr->viscoelasticCoefficient( iNode ) / ( A * std::sqrt( A ) ) * dAdt( A, timeStep, iNode );
}

inline Real
OneDFSIPhysics::dPTdU( const Real& A, const Real& Q, const Real& timeStep, const ID& id, const UInt& iNode ) const
{
    if ( id == 0 ) // dPt/dA
        return dPdA( A, timeStep, iNode ) - M_dataPtr->densityRho() * Q * Q / ( A * A * A );

    if ( id == 1 ) // dPt/dQ
        return M_dataPtr->densityRho() * Q / ( A * A );

    ERROR_MSG("Total pressure's differential function has only 2 components.");
    return -1.;
}

void Problem::readLastVectorALETimeAdvance ( vectorPtr_Type fluidDisp,
                                             const std::string loadInitSol)
{
    using namespace LifeV;

    typedef FSIOperator::mesh_Type        mesh_Type;

    //We still need to load the last vector for ALE
    std::string iterationString = loadInitSol;
    fluidDisp.reset (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), LifeV::Unique) );

    //Setting the exporterData to read: ALE problem
    LifeV::ExporterData<mesh_Type> initSolFluidDisp  (LifeV::ExporterData<mesh_Type>::VectorField, "f-displacement." + iterationString, M_fsi->FSIOper()->mmFESpacePtr(), fluidDisp, UInt (0), LifeV::ExporterData<mesh_Type>::UnsteadyRegime );

    //Initializing
    *fluidDisp *= 0.0;

    //Reading
    M_importerFluid->readVariable (initSolFluidDisp); //Fluid df

    //Output
    std::cout << "Norm of the df " << fluidDisp->norm2() << std::endl;

    //This is ugly but it's the only way I have figured out at the moment
    if ( M_data->method().compare ("monolithicGI") == 0 )
    {
        //Don't be scared by the ten. The goal of 10 is just to make the first if fail
        M_fsi->FSIOper()->setALEVectorInStencil ( fluidDisp, 10, true );
    }

    //Setting the vector in the stencil
    M_fsi->FSIOper()->ALETimeAdvance()->shiftRight ( *fluidDisp );
}

void Problem::initializeWithVectors ( void )
{

    using namespace LifeV;
    // vectors to store the solutions we want.
    vectorPtr_Type vel;
    vectorPtr_Type pressure;
    vectorPtr_Type solidDisp;
    vectorPtr_Type fluidDisp;

    vel.reset (new vector_Type (M_fsi->FSIOper()->uFESpace().map(), LifeV::Unique) );
    pressure.reset (new vector_Type (M_fsi->FSIOper()->pFESpace().map(), LifeV::Unique) );
    solidDisp.reset (new vector_Type (M_fsi->FSIOper()->dFESpace().map(), LifeV::Unique) );
    fluidDisp.reset (new vector_Type (M_fsi->FSIOper()->mmFESpace().map(), LifeV::Unique) );

    // In this case we want to initialize only the pressure
    M_fsi->FSIOper()->pFESpacePtr()->interpolate ( static_cast<FESpace<RegionMesh<LinearTetra>, MapEpetra> ::function_Type> ( pressureInitial ), *pressure, 0.0 );

    *vel *= 0.0;
    *solidDisp *= 0.0;
    *fluidDisp *= 0.0;

    UInt iterInit;

    // Filling the stencils
    for (iterInit = 0; iterInit < M_fsi->FSIOper()->fluidTimeAdvance()->size(); iterInit++ )
    {
        //We send the vectors to the FSIMonolithic class using the interface of FSIOper
        M_fsi->FSIOper()->setVectorInStencils (vel, pressure, solidDisp, iterInit );
    }

    // This was in readLastVectorSolidStencil
    M_fsi->FSIOper()->setSolidVectorInStencil ( solidDisp, iterInit );

    // Ale part
    for (iterInit = 0; iterInit < M_fsi->FSIOper()->ALETimeAdvance()->size(); iterInit++ )
    {
        //Setting the vector in the stencil
        M_fsi->FSIOper()->setALEVectorInStencil ( fluidDisp, iterInit, false );
    }

    //Initializing the vector for the RHS terms of the formulas
    M_fsi->FSIOper()->finalizeRestart();

    // This was read the last vector from ALE
    //This is ugly but it's the only way I have figured out at the moment
    if ( M_data->method().compare ("monolithicGI") == 0 )
    {
        //Don't be scared by the ten. The goal of 10 is just to make the first if fail
        M_fsi->FSIOper()->setALEVectorInStencil ( fluidDisp, 10, true );
    }

    //Setting the vector in the stencil
    M_fsi->FSIOper()->ALETimeAdvance()->shiftRight ( *fluidDisp );

}

inline Real
OneDFSIPhysics::dAdP( const Real& P, const Real& timeStep, const UInt& iNode, const bool& elasticExternalNodes ) const
{
    if ( !M_dataPtr->viscoelasticWall() || ( ( iNode == 0 || iNode == M_dataPtr->numberOfNodes() - 1 ) && elasticExternalNodes ) )
    {
        return M_dataPtr->area0( iNode ) / ( M_dataPtr->beta0( iNode ) * M_dataPtr->beta1( iNode ) )
                                      * OneDFSI::pow10( 1 + ( P - externalPressure() )
                                      / M_dataPtr->beta0( iNode ), 1 / M_dataPtr->beta1( iNode ) - 1 );
    }
    else
    {
        // Finite difference approach
        return ( fromPToA( P + M_dataPtr->jacobianPerturbationStress(), timeStep, iNode, elasticExternalNodes ) - fromPToA( P, timeStep, iNode, elasticExternalNodes ) )
               / M_dataPtr->jacobianPerturbationStress();
    }
}

inline Real
OneDFSIPhysics::dPdAelastic( const Real& A, const UInt& iNode ) const
{
    return M_dataPtr->beta0( iNode ) * M_dataPtr->beta1( iNode ) * OneDFSI::pow05( A / M_dataPtr->area0( iNode ), M_dataPtr->beta1( iNode ) ) / A;
}

// ===================================================
// Inline methods
// ===================================================
inline Real
OneDFSIPhysics::celerity0( const UInt& iNode ) const
{
    return std::sqrt( M_dataPtr->beta0( iNode ) * M_dataPtr->beta1( iNode ) / M_dataPtr->densityRho() );
}

inline Real
OneDFSIPhysics::totalPressure( const Real& A, const Real& Q, const UInt& iNode ) const
{
    return elasticPressure( A, iNode ) + M_dataPtr->densityRho() / 2 * Q * Q / ( A * A );
}