Foam::labelList Foam::metisDecomp::decompose
(
    const labelListList& globalCellCells,
    const pointField& cellCentres,
    const scalarField& cellWeights
)
{
    if (cellCentres.size() != globalCellCells.size())
    {
        FatalErrorIn
        (
            "metisDecomp::decompose"
            "(const pointField&, const labelListList&, const scalarField&)"
        )   << "Inconsistent number of cells (" << globalCellCells.size()
            << ") and number of cell centres (" << cellCentres.size()
            << ")." << exit(FatalError);
    }


    // Make Metis CSR (Compressed Storage Format) storage
    //   adjncy      : contains neighbours (= edges in graph)
    //   xadj(celli) : start of information in adjncy for celli

    List<int> adjncy;
    List<int> xadj;
    scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);


    // Decompose using default weights
    List<int> finalDecomp;
    decompose(adjncy, xadj, cellWeights, finalDecomp);

    // Copy back to labelList
    labelList decomp(finalDecomp.size());
    forAll(decomp, i)
    {
        decomp[i] = finalDecomp[i];
    }
    return decomp;
}

autoPtr<mapDistribute> buildMap
(
    const T& mesh,
    labelListList& pointPoints
)
{
    pointPoints.setSize(mesh.vertexCount());

    globalIndex globalIndexing(mesh.vertexCount());

    for
    (
        typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
        vit != mesh.finite_vertices_end();
        ++vit
    )
    {
        if (!vit->real())
        {
            continue;
        }

        std::list<typename T::Vertex_handle> adjVerts;
        mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));

        DynamicList<label> indices(adjVerts.size());

        for
        (
            typename std::list<typename T::Vertex_handle>::const_iterator
                adjVertI = adjVerts.begin();
            adjVertI != adjVerts.end();
            ++adjVertI
        )
        {
            typename T::Vertex_handle vh = *adjVertI;

            if (!vh->farPoint())
            {
                indices.append
                (
                    globalIndexing.toGlobal(vh->procIndex(), vh->index())
                );
            }
        }

        pointPoints[vit->index()].transfer(indices);
    }

    List<Map<label>> compactMap;

    return autoPtr<mapDistribute>
    (
        new mapDistribute
        (
            globalIndexing,
            pointPoints,
            compactMap
        )
    );
}

// Given a subset of cells determine the new global indices. The problem
// is in the cells from neighbouring processors which need to be renumbered.
void Foam::multiLevelDecomp::subsetGlobalCellCells
(
    const label nDomains,
    const label domainI,
    const labelList& dist,

    const labelListList& cellCells,
    const labelList& set,
    labelListList& subCellCells,
    labelList& cutConnections
) const
{
    // Determine new index for cells by inverting subset
    labelList oldToNew(invert(cellCells.size(), set));

    globalIndex globalCells(cellCells.size());

    // Subset locally the elements for which I have data
    subCellCells = UIndirectList<labelList>(cellCells, set);

    // Get new indices for neighbouring processors
    List<Map<label>> compactMap;
    mapDistribute map(globalCells, subCellCells, compactMap);
    map.distribute(oldToNew);
    labelList allDist(dist);
    map.distribute(allDist);

    // Now we have:
    // oldToNew : the locally-compact numbering of all our cellCells. -1 if
    //            cellCell is not in set.
    // allDist  : destination domain for all our cellCells
    // subCellCells : indexes into oldToNew and allDist

    // Globally compact numbering for cells in set.
    globalIndex globalSubCells(set.size());

    // Now subCellCells contains indices into oldToNew which are the
    // new locations of the neighbouring cells.

    cutConnections.setSize(nDomains);
    cutConnections = 0;

    forAll(subCellCells, subCelli)
    {
        labelList& cCells = subCellCells[subCelli];

        // Keep the connections to valid mapped cells
        label newI = 0;
        forAll(cCells, i)
        {
            // Get locally-compact cell index of neighbouring cell
            label nbrCelli = oldToNew[cCells[i]];
            if (nbrCelli == -1)
            {
                cutConnections[allDist[cCells[i]]]++;
            }
            else
            {
                // Reconvert local cell index into global one

                // Get original neighbour
                label celli = set[subCelli];
                label oldNbrCelli = cellCells[celli][i];
                // Get processor from original neighbour
                label proci = globalCells.whichProcID(oldNbrCelli);
                // Convert into global compact numbering
                cCells[newI++] = globalSubCells.toGlobal(proci, nbrCelli);
            }
        }

void Pstream::exchange
(
    const List<Container>& sendBufs,
    List<Container>& recvBufs,
    labelListList& sizes,
    const int tag,
    const bool block
)
{
    if (!contiguous<T>())
    {
        FatalErrorIn
        (
            "Pstream::exchange(..)"
        )   << "Continuous data only." << Foam::abort(FatalError);
    }

    if (sendBufs.size() != UPstream::nProcs())
    {
        FatalErrorIn
        (
            "Pstream::exchange(..)"
        )   << "Size of list:" << sendBufs.size()
            << " does not equal the number of processors:"
            << UPstream::nProcs()
            << Foam::abort(FatalError);
    }

    sizes.setSize(UPstream::nProcs());
    labelList& nsTransPs = sizes[UPstream::myProcNo()];
    nsTransPs.setSize(UPstream::nProcs());

    forAll(sendBufs, procI)
    {
        nsTransPs[procI] = sendBufs[procI].size();
    }

    // Send sizes across. Note: blocks.
    combineReduce(sizes, UPstream::listEq(), tag);

    if (Pstream::parRun())
    {
        label startOfRequests = Pstream::nRequests();

        // Set up receives
        // ~~~~~~~~~~~~~~~

        recvBufs.setSize(sendBufs.size());
        forAll(sizes, procI)
        {
            label nRecv = sizes[procI][UPstream::myProcNo()];

            if (procI != Pstream::myProcNo() && nRecv > 0)
            {
                recvBufs[procI].setSize(nRecv);
                UIPstream::read
                (
                    UPstream::nonBlocking,
                    procI,
                    reinterpret_cast<char*>(recvBufs[procI].begin()),
                    nRecv*sizeof(T),
                    tag
                );
            }
        }


        // Set up sends
        // ~~~~~~~~~~~~

        forAll(sendBufs, procI)
        {
            if (procI != Pstream::myProcNo() && sendBufs[procI].size() > 0)
            {
                if
                (
                   !UOPstream::write
                    (
                        UPstream::nonBlocking,
                        procI,
                        reinterpret_cast<const char*>(sendBufs[procI].begin()),
                        sendBufs[procI].size()*sizeof(T),
                        tag
                    )
                )
                {
                    FatalErrorIn("Pstream::exchange(..)")
                        << "Cannot send outgoing message. "
                        << "to:" << procI << " nBytes:"
                        << label(sendBufs[procI].size()*sizeof(T))
                        << Foam::abort(FatalError);
                }
            }
        }


        // Wait for all to finish
        // ~~~~~~~~~~~~~~~~~~~~~~

        if (block)
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
#   include "setRootCase.H"
#   include "createTime.H"
#   include "createMesh.H"

    //----------------------------------------------------------------------//
    //- open abaqus input file
    //----------------------------------------------------------------------//
    Info << "Opening Abaqus input file" << endl << endl;
    OFstream abaqusInp("abaqusMesh.inp");

    //- input header
    Info << "Writing input file header" << endl << endl;
    abaqusInp << "*Heading\n"
        << "** Job name: OpenFoamMesh Model name: OpenFoamMesh\n"
        << "** Generated by: Abaqus/CAE 6.9-2\n"
        << "*Preprint, echo=NO, model=NO, history=NO, contact=NO\n"
        << "**\n"
        << "** PARTS\n"
        << "**\n"
        << "*Part, name=OpenFoamMeshPart"
        << endl;


    //----------------------------------------------------------------------//
    //- write nodes
    //----------------------------------------------------------------------//
    Info << "Writing Nodes" << endl << endl;
    abaqusInp << "*Node" << endl;
    const pointField& points = mesh.points();
    forAll(points, pointi)
    {
        abaqusInp << "\t" << (pointi+1)
            << ",\t" << (points[pointi]).x()
            << ",\t" << (points[pointi]).y()
            << ",\t" << (points[pointi]).z()
            << endl;
    }


    //----------------------------------------------------------------------//
    //- determine abaqusCellPoints
    //----------------------------------------------------------------------//
    //- for hex 1st order elements, abaqus orders the points
    //- where nodes 1 2 3 4 are clockwise from the outside of the cell
    //- and then 5 6 7 8 are apposite 1 2 3 4 respectively
    Info << "Determining Abaqus element node ordering" << endl << endl;

    const cellList& cells = mesh.cells();
    const faceList& faces = mesh.faces();
    const labelListList pointPoints = mesh.pointPoints();
    const labelListList cellPoints = mesh.cellPoints();
    const labelList faceOwner = mesh.faceOwner();
    const vectorField faceNormals = mesh.Sf()/mesh.magSf();

    labelListList abaqusCellPoints(cellPoints.size(), List<label>(8, 1));

    forAll(cells, celli)
    {
        //- face0 will be our reference face
        //- face0 seems to be always owned by the cell
        //- so the face normal points out of the cell
        label refFace = cells[celli][0];

        //- insert first four abaqusCellPoints
        abaqusCellPoints[celli][0] = (faces[refFace][3] + 1);
        abaqusCellPoints[celli][1] = (faces[refFace][2] + 1);
        abaqusCellPoints[celli][2] = (faces[refFace][1] + 1);
        abaqusCellPoints[celli][3] = (faces[refFace][0] + 1);

        //- now find the opposite face in the cell
        //Info << "Finding oppFace" << endl << endl;
        labelList refFacePoints = faces[refFace];
        //- compare each faces points to the refFace, the opposite
        //- face will share points with the refFace
        label oppFace = -1;

        //- for all the cell faces
        forAll(cells[celli], facei)
        {
            bool faceHasNoCommonPoints = true;
            label globalFaceLabel = cells[celli][facei];

            //- for all the face points
            forAll(faces[globalFaceLabel], pointi)
            {
                label globalPointLabel = faces[globalFaceLabel][pointi];

                //- compare each point with all the refFace points
                forAll(faces[refFace], refFacePointi)
                {
                    label refFaceGlobalPointLabel = faces[refFace][refFacePointi];
                    if(globalPointLabel == refFaceGlobalPointLabel)
                      {
                        faceHasNoCommonPoints = false;
                      }
                }

	void CalculateDragForce
	(
		cfdemCloud& sm,
		const volScalarField& alpf_,
		const volVectorField& Uf_,
		const volScalarField& rho_,
		const bool& verbose_,
		vectorField& DragForce_,
		const labelListList& particleList_
	)
	{
		
		// get viscosity field
		#ifdef comp
		    const volScalarField nufField = sm.turbulence().mu()/rho_;
		#else
		    const volScalarField& nufField = sm.turbulence().nu();
		#endif

		// Local variables	
		label  cellI(-1);
		vector drag(0,0,0);
		vector Ufluid(0,0,0);
		
		vector position(0,0,0);
		scalar voidfraction(1);
		
		vector Up(0,0,0);
		vector Ur(0,0,0);
		scalar ds(0);
		
		scalar nuf(0);
		scalar rhof(0);
		
		vector WenYuDrag(0,0,0);
		
		interpolationCellPoint<scalar> voidfractionInterpolator_(alpf_);
		interpolationCellPoint<vector> UInterpolator_(Uf_);	
				
		// 
		//_AO_Parallel
		DragForce_.resize(particleList_.size());
					
	    	for(int ii =0; ii < particleList_.size(); ii++)
		{
			int index = particleList_[ii][0];
			cellI = sm.cellIDs()[index][0];
			drag = vector(0,0,0);
			Ufluid = vector(0,0,0);
			WenYuDrag = vector(0,0,0);
			DragForce_[ii] = vector(0,0,0);
			    
			if (cellI > -1) // particle Found
			{
				position = sm.position(index);
			
				if ( alpf_[cellI] > 1. ) Pout << " voidfraction > 1 " << alpf_[cellI] << endl;
			
				voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
				Ufluid = UInterpolator_.interpolate(position,cellI);
			
				if ( voidfraction > 1. ) 
				{
						Pout << " Int. voidfraction > 1 " << " value= " << voidfraction;
						voidfraction = alpf_[cellI];
						Pout << " mod. value = " << voidfraction << endl;
				}		
				Up = sm.velocity(index);
			
				Ur = Ufluid-Up;				
				ds = 2*sm.radius(index);		
				rhof = rho_[cellI];
				nuf = nufField[cellI];
			
				// Drag force
				WenYuDragForce(Ur,ds,rhof,nuf,voidfraction,WenYuDrag);
						
				if(verbose_ && index <= 1)
				{
					Info << "" << endl;
					Pout << " index = " << index << endl;
					Pout << " position = " << position << endl; 
					Pout << " Up = " << Up << endl;
					Pout << " Ur = " << Ur << endl;
					Pout << " dp = " << ds << endl;
					Pout << " rho = " << rhof << endl;
					Pout << " nuf = " << nuf << endl;
					Pout << " voidfraction = " << voidfraction << endl;
					Pout << " drag = " << WenYuDrag << endl;
					Info << " " << endl;
				}
			}	
			for(int j=0;j<3;j++) DragForce_[ii][j] = WenYuDrag[j];
		}	
	}

	void EulerianParticleVelocityForce
	(
		cfdemCloud& sm,			
		const fvMesh& mesh,
		volVectorField& Uf_,
		volVectorField&	Up_,
		volScalarField& rho_,
		volScalarField& alpf_,
		volScalarField& Pg_,
		volVectorField& MappedDragForce_,
		const labelListList& particleList_,
		const bool& weighting_
	)
	{		
		// Neighbouring cells
		CPCCellToCellStencil neighbourCells(mesh);
				
		// get viscosity field
		#ifdef comp
		    const volScalarField nufField = sm.turbulence().mu()/rho_;
		#else
		    const volScalarField& nufField = sm.turbulence().nu();
		#endif

		// Gas pressure gradient
		volVectorField gradPg_ = fvc::grad(Pg_);
		interpolationCellPoint<vector> gradPgInterpolator_(gradPg_);

		// Local variables	
		label  cellID(-1);
		vector drag(0,0,0);
		vector Ufluid(0,0,0);
		
		vector position(0,0,0);
		scalar voidfraction(1);
		
		vector Up(0,0,0);
		vector Ur(0,0,0);
		scalar ds(0);
		
		scalar nuf(0);
		scalar rhof(0);
		
		vector WenYuDrag(0,0,0);
		
		interpolationCellPoint<scalar> voidfractionInterpolator_(alpf_);
		interpolationCellPoint<vector> UInterpolator_(Uf_);	
		
		scalar dist_s(0);
		scalar sumWeights(0);
		
		scalarField               weightScalar(27,scalar(0.0));
		Field <Field <scalar> >   particleWeights(particleList_.size(),weightScalar);
		
		//Info << " particle size " << particleList_.size() << endl;
		
		// Number of particle in a cell
		scalarField np(mesh.cells().size(),scalar(0));
		
		// Particle volume
		scalar Volp(0);
		vector gradPg_int(0,0,0);
		
		for(int ii = 0; ii < particleList_.size(); ii++)
		{
			int index = particleList_[ii][0];
			
			cellID = sm.cellIDs()[index][0];
			position = sm.position(index);			    

                        Ufluid = UInterpolator_.interpolate(position,cellID); 
			Up = sm.velocity(index);
                        Ur = Ufluid-Up;

                        ds = 2*sm.radius(index);

                        // Calculate WenYu Drag 
                        voidfraction = voidfractionInterpolator_.interpolate(position,cellID);
                        nuf = nufField[cellID];
                        rhof = rho_[cellID];	
                        WenYuDragForce(Ur,ds,rhof,nuf,voidfraction,WenYuDrag);	
    
        		Volp = ds*ds*ds*M_PI/6;
			gradPg_int = gradPgInterpolator_.interpolate(position,cellID);
			
			//if (cellID > -1)  // particle centre is in domain
            		//{
				if(weighting_)
				{
					labelList& cellsNeigh = neighbourCells[cellID];
					sumWeights = 0;
					dist_s = 0;

					//Info << " index = " << index << " ii = " << ii << " cellID = " << cellID << endl;

					forAll(cellsNeigh,jj)
					{
						// Find distances between particle and neighbouring cells					
						dist_s = mag(sm.mesh().C()[cellsNeigh[jj]]-position)/pow(sm.mesh().V()[cellsNeigh[jj]],1./3.);

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