C++ labelListList::size方法代码示例

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

    CompactListList<label> cellCells(globalCellCells);

    // Decompose using default weights
    labelList decomp;
    decompose(cellCells.m(), cellCells.offsets(), cellWeights, decomp);

    return decomp;
}

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

//    // For running sequential ...
//    if (Pstream::nProcs() <= 1)
//    {
//        return scotchDecomp(decompositionDict_, mesh)
//            .decompose(globalCellCells, cellCentres, cWeights);
//    }


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

    CompactListList<label> cellCells(globalCellCells);

    // Decompose using weights
    List<int> finalDecomp;
    decomposeZeroDomains
    (
        "ptscotch",
        cellCells.m(),
        cellCells.offsets(),
        cWeights,
        finalDecomp
    );

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

void referredCell::locallyMapFaceList
(
    const labelList& points,
    const labelListList& sourceCellFaces
)
{
    faces_.setSize(sourceCellFaces.size());

    forAll(sourceCellFaces, sCF)
    {
        const labelList& sourceCellFace(sourceCellFaces[sCF]);

        labelList& localFace(faces_[sCF]);

        localFace.setSize(sourceCellFace.size());

        forAll(sourceCellFace, p)
        {
            localFace[p] = findIndex(points, sourceCellFace[p]);

            if (localFace[p] == -1)
            {
                FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
                    << "edgeList and points labelList for "
                    << "referred cell do not match: "
                    << nl << "points: " << points
                    << nl << "faces: " << sourceCellFaces
                    << abort(FatalError);
            }
        }
    }

Foam::Field<Type>::Field
(
    const UList<Type>& mapF,
    const labelListList& mapAddressing,
    const scalarListList& mapWeights
)
:
    List<Type>(mapAddressing.size())
{
    map(mapF, mapAddressing, mapWeights);
}

Foam::labelList Foam::scotchDecomp::decompose
(
    const labelListList& globalCellCells,
    const pointField& cc,
    const scalarField& cWeights
)
{
    if (cc.size() != globalCellCells.size())
    {
        FatalErrorIn
        (
            "scotchDecomp::decompose"
            "(const labelListList&, const pointField&, const scalarField&)"
        )   << "Inconsistent number of cells (" << globalCellCells.size()
            << ") and number of cell centres (" << cc.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;
    calcCSR(globalCellCells, adjncy, xadj);

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

    // Copy back to labelList
    labelList decomp(finalDecomp.size());

    forAll(decomp, i)
    {
        decomp[i] = finalDecomp[i];
    }

    return decomp;
}

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

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

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