Commit 9891d371 authored by mattijs's avatar mattijs
Browse files

parallel postChannel

parent 0bbc8fd2
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
EXE_LIBS = \
-lmeshTools \
-lfiniteVolume \
-lsampling
......@@ -25,156 +25,269 @@ License
\*---------------------------------------------------------------------------*/
#include "channelIndex.H"
#include "boolList.H"
#include "syncTools.H"
#include "OFstream.H"
#include "meshTools.H"
#include "Time.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * * * //
channelIndex::channelIndex(const fvMesh& m)
:
indexingDict_
(
IOobject
(
"postChannelDict",
m.time().constant(),
m,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
nx_(readLabel(indexingDict_.lookup("Nx"))),
ny_(indexingDict_.lookup("Ny")),
nz_(readLabel(indexingDict_.lookup("Nz"))),
symmetric_
(
readBool(indexingDict_.lookup("symmetric"))
),
cumNy_(ny_.size()),
nLayers_(ny_[0])
template<>
const char* Foam::NamedEnum<Foam::vector::components, 3>::names[] =
{
// initialise the layers
cumNy_[0] = ny_[0];
for (label j=1; j<ny_.size(); j++)
{
nLayers_ += ny_[j];
cumNy_[j] = ny_[j]+cumNy_[j-1];
}
}
"x",
"y",
"z"
};
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
const Foam::NamedEnum<Foam::vector::components, 3>
Foam::channelIndex::vectorComponentsNames_;
channelIndex::~channelIndex()
{}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
scalarField channelIndex::collapse
// Determines face blocking
void Foam::channelIndex::walkOppositeFaces
(
const volScalarField& vsf,
const bool asymmetric
) const
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
)
{
scalarField cs(nLayers(), 0.0);
const cellList& cells = mesh.cells();
const faceList& faces = mesh.faces();
label nBnd = mesh.nFaces() - mesh.nInternalFaces();
DynamicList<label> frontFaces(startFaces);
forAll(frontFaces, i)
{
label faceI = frontFaces[i];
blockedFace[faceI] = true;
}
forAll(cs, j)
while (returnReduce(frontFaces.size(), sumOp<label>()) > 0)
{
// sweep over all cells in this layer
for (label i=0; i<nx(); i++)
// Transfer across.
boolList isFrontBndFace(nBnd, false);
forAll(frontFaces, i)
{
for (label k=0; k<nz(); k++)
label faceI = frontFaces[i];
if (!mesh.isInternalFace(faceI))
{
cs[j] += vsf[operator()(i,j,k)];
isFrontBndFace[faceI-mesh.nInternalFaces()] = true;
}
}
syncTools::swapBoundaryFaceList(mesh, isFrontBndFace, false);
// and divide by the number of cells in the layer
cs[j] /= scalar(nx()*nz());
}
// Add
forAll(isFrontBndFace, i)
{
label faceI = mesh.nInternalFaces()+i;
if (isFrontBndFace[i] && !blockedFace[faceI])
{
blockedFace[faceI] = true;
frontFaces.append(faceI);
}
}
if (symmetric_)
{
label nlb2 = nLayers()/2;
// Transfer across cells
DynamicList<label> newFrontFaces(frontFaces.size());
if (asymmetric)
forAll(frontFaces, i)
{
for (label j=0; j<nlb2; j++)
label faceI = frontFaces[i];
{
cs[j] = 0.5*(cs[j] - cs[nLayers() - j - 1]);
const cell& ownCell = cells[mesh.faceOwner()[faceI]];
label oppositeFaceI = ownCell.opposingFaceLabel(faceI, faces);
if (oppositeFaceI == -1)
{
FatalErrorIn("channelIndex::walkOppositeFaces(..)")
<< "Face:" << faceI << " owner cell:" << ownCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFaceI])
{
blockedFace[oppositeFaceI] = true;
newFrontFaces.append(oppositeFaceI);
}
}
}
if (mesh.isInternalFace(faceI))
{
const cell& neiCell = mesh.cells()[mesh.faceNeighbour()[faceI]];
label oppositeFaceI = neiCell.opposingFaceLabel(faceI, faces);
if (oppositeFaceI == -1)
{
FatalErrorIn("channelIndex::walkOppositeFaces(..)")
<< "Face:" << faceI << " neighbour cell:" << neiCell
<< " is not a hex?" << abort(FatalError);
}
else
{
if (!blockedFace[oppositeFaceI])
{
blockedFace[oppositeFaceI] = true;
newFrontFaces.append(oppositeFaceI);
}
}
}
}
else
frontFaces.transfer(newFrontFaces);
}
}
// Calculate regions.
void Foam::channelIndex::calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
)
{
boolList blockedFace(mesh.nFaces(), false);
walkOppositeFaces
(
mesh,
startFaces,
blockedFace
);
if (false)
{
OFstream str(mesh.time().path()/"blockedFaces.obj");
label vertI = 0;
forAll(blockedFace, faceI)
{
for (label j=0; j<nlb2; j++)
if (blockedFace[faceI])
{
cs[j] = 0.5*(cs[j] + cs[nLayers() - j - 1]);
const face& f = mesh.faces()[faceI];
forAll(f, fp)
{
meshTools::writeOBJ(str, mesh.points()[f[fp]]);
}
str<< 'f';
forAll(f, fp)
{
str << ' ' << vertI+fp+1;
}
str << nl;
vertI += f.size();
}
}
cs.setSize(nlb2);
}
return cs;
// Do analysis for connected regions
cellRegion_.reset(new regionSplit(mesh, blockedFace));
Info<< "Detected " << cellRegion_().nRegions() << " layers." << nl << endl;
// Sum number of entries per region
regionCount_ = regionSum(scalarField(mesh.nCells(), 1.0));
// Average cell centres to determine ordering.
pointField regionCc
(
regionSum(mesh.cellCentres())
/ regionCount_
);
SortableList<scalar> sortComponent(regionCc.component(dir_));
sortMap_ = sortComponent.indices();
y_ = sortComponent;
if (symmetric_)
{
y_.setSize(cellRegion_().nRegions()/2);
}
}
scalarField channelIndex::y
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::channelIndex::channelIndex
(
const volVectorField& cellCentres
) const
const polyMesh& mesh,
const dictionary& dict
)
:
symmetric_(readBool(dict.lookup("symmetric"))),
dir_(vectorComponentsNames_.read(dict.lookup("component")))
{
if (symmetric_)
const polyBoundaryMesh& patches = mesh.boundaryMesh();
const wordList patchNames(dict.lookup("patches"));
label nFaces = 0;
forAll(patchNames, i)
{
scalarField Y(nLayers()/2);
label patchI = patches.findPatchID(patchNames[i]);
for (label j=0; j<nLayers()/2; j++)
if (patchI == -1)
{
Y[j] = cellCentres[operator()(0, j, 0)].y();
FatalErrorIn("channelIndex::channelIndex(const polyMesh&)")
<< "Illegal patch " << patchNames[i]
<< ". Valid patches are " << patches.name()
<< exit(FatalError);
}
return Y;
nFaces += patches[patchI].size();
}
else
labelList startFaces(nFaces);
nFaces = 0;
forAll(patchNames, i)
{
scalarField Y(nLayers());
const polyPatch& pp = patches[patches.findPatchID(patchNames[i])];
for (label j=0; j<nLayers(); j++)
forAll(pp, j)
{
Y[j] = cellCentres[operator()(0, j, 0)].y();
startFaces[nFaces++] = pp.start()+j;
}
return Y;
}
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
label channelIndex::operator()
Foam::channelIndex::channelIndex
(
const label Jx,
const label Jy,
const label Jz
) const
const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
)
:
symmetric_(symmetric),
dir_(dir)
{
label index(0);
// Calculate regions.
calcLayeredRegions(mesh, startFaces);
}
// count up `full' layers in the mesh
label j(0);
label tmpJy(Jy);
while(Jy >= cumNy_[j])
{
index += nx_*ny_[j]*nz_;
tmpJy -= ny_[j];
j++;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
index += Jx + nx_*tmpJy + nx_*ny_[j]*Jz;
return index;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// ************************************************************************* //
......@@ -26,19 +26,25 @@ Class
Foam::channelIndex
Description
does indexing for a standard (Christer-style) channel.
Assumes that the blocks are aranged in the y direction.
Does averaging of fields over layers of cells. Assumes layered mesh.
SourceFiles
channelIndex.C
channelIndexIO.C
\*---------------------------------------------------------------------------*/
#ifndef channelIndex_H
#define channelIndex_H
#include "fvCFD.H"
#include "regionSplit.H"
#include "direction.H"
#include "scalarField.H"
#include "polyMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
......@@ -47,21 +53,47 @@ SourceFiles
class channelIndex
{
// Private data
IOdictionary indexingDict_;
static const NamedEnum<vector::components, 3> vectorComponentsNames_;
const label nx_;
const labelList ny_;
const label nz_;
//- Is mesh symmetric
const bool symmetric_;
labelList cumNy_;
label nLayers_;
//- direction to sort
const direction dir_;
//- Per cell the global region
autoPtr<regionSplit> cellRegion_;
//- Per global region the number of cells (scalarField so we can use
// field algebra)
scalarField regionCount_;
//- From sorted region back to unsorted global region
labelList sortMap_;
//- Sorted component of cell centres
scalarField y_;
// Private Member Functions
void walkOppositeFaces
(
const polyMesh& mesh,
const labelList& startFaces,
boolList& blockedFace
);
void calcLayeredRegions
(
const polyMesh& mesh,
const labelList& startFaces
);
//- Disallow default bitwise copy construct and assignment
channelIndex(const channelIndex&);
void operator=(const channelIndex&);
......@@ -71,55 +103,53 @@ public:
// Constructors
channelIndex(const fvMesh& m);
//- Construct from dictionary
channelIndex(const polyMesh&, const dictionary&);
// Destructor
~channelIndex();
//- Construct from supplied starting faces
channelIndex
(
const polyMesh& mesh,
const labelList& startFaces,
const bool symmetric,
const direction dir
);
// Member Functions
// Access
//- return number of layers
label nLayers() const
{
return nLayers_;
}
//- number of cells in X direction
label nx() const
{
return nx_;
}
//- number of cells in Z direction
label nz() const
{
return nz_;
}
//- Sum field per region
template<class T>
Field<T> regionSum(const Field<T>& cellField) const;
//- collapse a field to a line
scalarField collapse
template<class T>
Field<T> collapse
(
const volScalarField& vsf,
const Field<T>& vsf,
const bool asymmetric=false
) const;
//- return the field of Y locations from the cell centres
scalarField y
(
const volVectorField& cellCentres
) const;
const scalarField& y() const
{
return y_;
}
};
// Member Operators
label operator()(const label, const label, const label) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "channelIndexTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "channelIndex.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class T>
Foam::Field<T> Foam::channelIndex::regionSum(const Field<T>& cellField) const
{
Field<T> regionField(cellRegion_().nRegions(), pTraits<T>::zero);
forAll(cellRegion_(), cellI)