629 lines
18 KiB
C
629 lines
18 KiB
C
|
/*---------------------------------------------------------------------------*\
|
||
|
========= |
|
||
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||
|
\\ / O peration |
|
||
|
\\ / A nd | Copyright held by original author
|
||
|
\\/ 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
|
||
|
|
||
|
Description
|
||
|
Create polyMesh from cell and patch shapes
|
||
|
|
||
|
\*---------------------------------------------------------------------------*/
|
||
|
|
||
|
#include "polyMesh.H"
|
||
|
#include "Time.H"
|
||
|
#include "primitiveMesh.H"
|
||
|
#include "DynamicList.H"
|
||
|
|
||
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||
|
|
||
|
Foam::labelListList Foam::polyMesh::cellShapePointCells
|
||
|
(
|
||
|
const cellShapeList& c
|
||
|
) const
|
||
|
{
|
||
|
List<DynamicList<label, primitiveMesh::cellsPerPoint_> >
|
||
|
pc(points().size());
|
||
|
|
||
|
// For each cell
|
||
|
forAll(c, i)
|
||
|
{
|
||
|
// For each vertex
|
||
|
const labelList& labels = c[i];
|
||
|
|
||
|
forAll(labels, j)
|
||
|
{
|
||
|
// Set working point label
|
||
|
label curPoint = labels[j];
|
||
|
DynamicList<label, primitiveMesh::cellsPerPoint_>& curPointCells =
|
||
|
pc[curPoint];
|
||
|
|
||
|
// Enter the cell label in the point's cell list
|
||
|
curPointCells.append(i);
|
||
|
}
|
||
|
}
|
||
|
|
||
|
labelListList pointCellAddr(pc.size());
|
||
|
|
||
|
forAll (pc, pointI)
|
||
|
{
|
||
|
pointCellAddr[pointI].transfer(pc[pointI].shrink());
|
||
|
}
|
||
|
|
||
|
return pointCellAddr;
|
||
|
}
|
||
|
|
||
|
|
||
|
Foam::labelList Foam::polyMesh::facePatchFaceCells
|
||
|
(
|
||
|
const faceList& patchFaces,
|
||
|
const labelListList& pointCells,
|
||
|
const faceListList& cellsFaceShapes,
|
||
|
const label patchID
|
||
|
) const
|
||
|
{
|
||
|
register bool found;
|
||
|
|
||
|
labelList FaceCells(patchFaces.size());
|
||
|
|
||
|
forAll(patchFaces, fI)
|
||
|
{
|
||
|
found = false;
|
||
|
|
||
|
const face& curFace = patchFaces[fI];
|
||
|
const labelList& facePoints = patchFaces[fI];
|
||
|
|
||
|
forAll(facePoints, pointI)
|
||
|
{
|
||
|
const labelList& facePointCells = pointCells[facePoints[pointI]];
|
||
|
|
||
|
forAll(facePointCells, cellI)
|
||
|
{
|
||
|
faceList cellFaces = cellsFaceShapes[facePointCells[cellI]];
|
||
|
|
||
|
forAll(cellFaces, cellFace)
|
||
|
{
|
||
|
if (cellFaces[cellFace] == curFace)
|
||
|
{
|
||
|
// Found the cell corresponding to this face
|
||
|
FaceCells[fI] = facePointCells[cellI];
|
||
|
|
||
|
found = true;
|
||
|
}
|
||
|
if (found) break;
|
||
|
}
|
||
|
if (found) break;
|
||
|
}
|
||
|
if (found) break;
|
||
|
}
|
||
|
|
||
|
if (!found)
|
||
|
{
|
||
|
FatalErrorIn
|
||
|
(
|
||
|
"polyMesh::facePatchFaceCells(const faceList& patchFaces,"
|
||
|
"const labelListList& pointCells,"
|
||
|
"const faceListList& cellsFaceShapes,"
|
||
|
"const label patchID)"
|
||
|
) << "face " << fI << " in patch " << patchID
|
||
|
<< " does not have neighbour cell"
|
||
|
<< " face: " << patchFaces[fI]
|
||
|
<< abort(FatalError);
|
||
|
}
|
||
|
}
|
||
|
|
||
|
return FaceCells;
|
||
|
}
|
||
|
|
||
|
|
||
|
Foam::polyMesh::polyMesh
|
||
|
(
|
||
|
const IOobject& io,
|
||
|
const pointField& points,
|
||
|
const cellShapeList& cellsAsShapes,
|
||
|
const faceListList& boundaryFaces,
|
||
|
const wordList& boundaryPatchNames,
|
||
|
const wordList& boundaryPatchTypes,
|
||
|
const word& defaultBoundaryPatchName,
|
||
|
const word& defaultBoundaryPatchType,
|
||
|
const wordList& boundaryPatchPhysicalTypes,
|
||
|
const bool syncPar
|
||
|
)
|
||
|
:
|
||
|
objectRegistry(io),
|
||
|
primitiveMesh(),
|
||
|
allPoints_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"points",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::AUTO_WRITE
|
||
|
),
|
||
|
points
|
||
|
),
|
||
|
// To be re-sliced later. HJ, 19/oct/2008
|
||
|
points_(allPoints_, allPoints_.size()),
|
||
|
allFaces_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"faces",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::AUTO_WRITE
|
||
|
),
|
||
|
0
|
||
|
),
|
||
|
faces_(allFaces_, allFaces_.size()),
|
||
|
owner_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"owner",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::AUTO_WRITE
|
||
|
),
|
||
|
0
|
||
|
),
|
||
|
neighbour_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"neighbour",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::AUTO_WRITE
|
||
|
),
|
||
|
0
|
||
|
),
|
||
|
clearedPrimitives_(false),
|
||
|
boundary_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"boundary",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::AUTO_WRITE
|
||
|
),
|
||
|
*this,
|
||
|
boundaryFaces.size() + 1 // add room for a default patch
|
||
|
),
|
||
|
bounds_(allPoints_, syncPar),
|
||
|
directions_(Vector<label>::zero),
|
||
|
pointZones_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"pointZones",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::NO_WRITE
|
||
|
),
|
||
|
*this,
|
||
|
0
|
||
|
),
|
||
|
faceZones_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"faceZones",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::NO_WRITE
|
||
|
),
|
||
|
*this,
|
||
|
0
|
||
|
),
|
||
|
cellZones_
|
||
|
(
|
||
|
IOobject
|
||
|
(
|
||
|
"cellZones",
|
||
|
instance(),
|
||
|
meshSubDir,
|
||
|
*this,
|
||
|
IOobject::NO_READ,
|
||
|
IOobject::NO_WRITE
|
||
|
),
|
||
|
*this,
|
||
|
0
|
||
|
),
|
||
|
globalMeshDataPtr_(NULL),
|
||
|
moving_(false),
|
||
|
changing_(false),
|
||
|
curMotionTimeIndex_(time().timeIndex()),
|
||
|
oldAllPointsPtr_(NULL),
|
||
|
oldPointsPtr_(NULL)
|
||
|
{
|
||
|
if (debug)
|
||
|
{
|
||
|
Info<<"Constructing polyMesh from cell and boundary shapes." << endl;
|
||
|
}
|
||
|
|
||
|
// Remove all of the old mesh files if they exist
|
||
|
removeFiles(instance());
|
||
|
|
||
|
// Calculate the faces of all cells
|
||
|
// Initialise maximum possible numer of mesh faces to 0
|
||
|
label maxFaces = 0;
|
||
|
|
||
|
// Set up a list of face shapes for each cell
|
||
|
faceListList cellsFaceShapes(cellsAsShapes.size());
|
||
|
cellList cells(cellsAsShapes.size());
|
||
|
|
||
|
forAll(cellsFaceShapes, cellI)
|
||
|
{
|
||
|
cellsFaceShapes[cellI] = cellsAsShapes[cellI].faces();
|
||
|
|
||
|
cells[cellI].setSize(cellsFaceShapes[cellI].size());
|
||
|
|
||
|
// Initialise cells to -1 to flag undefined faces
|
||
|
static_cast<labelList&>(cells[cellI]) = -1;
|
||
|
|
||
|
// Count maximum possible numer of mesh faces
|
||
|
maxFaces += cellsFaceShapes[cellI].size();
|
||
|
}
|
||
|
|
||
|
// Set size of faces array to maximum possible number of mesh faces
|
||
|
allFaces_.setSize(maxFaces);
|
||
|
|
||
|
// Initialise number of faces to 0
|
||
|
label nFaces = 0;
|
||
|
|
||
|
// set reference to point-cell addressing
|
||
|
labelListList PointCells = cellShapePointCells(cellsAsShapes);
|
||
|
|
||
|
bool found = false;
|
||
|
|
||
|
forAll(cells, cellI)
|
||
|
{
|
||
|
// Note:
|
||
|
// Insertion cannot be done in one go as the faces need to be
|
||
|
// added into the list in the increasing order of neighbour
|
||
|
// cells. Therefore, all neighbours will be detected first
|
||
|
// and then added in the correct order.
|
||
|
|
||
|
const faceList& curFaces = cellsFaceShapes[cellI];
|
||
|
|
||
|
// Record the neighbour cell
|
||
|
labelList neiCells(curFaces.size(), -1);
|
||
|
|
||
|
// Record the face of neighbour cell
|
||
|
labelList faceOfNeiCell(curFaces.size(), -1);
|
||
|
|
||
|
label nNeighbours = 0;
|
||
|
|
||
|
// For all faces ...
|
||
|
forAll(curFaces, faceI)
|
||
|
{
|
||
|
// Skip faces that have already been matched
|
||
|
if (cells[cellI][faceI] >= 0) continue;
|
||
|
|
||
|
found = false;
|
||
|
|
||
|
const face& curFace = curFaces[faceI];
|
||
|
|
||
|
// Get the list of labels
|
||
|
const labelList& curPoints = curFace;
|
||
|
|
||
|
// For all points
|
||
|
forAll(curPoints, pointI)
|
||
|
{
|
||
|
// Gget the list of cells sharing this point
|
||
|
const labelList& curNeighbours =
|
||
|
PointCells[curPoints[pointI]];
|
||
|
|
||
|
// For all neighbours
|
||
|
forAll(curNeighbours, neiI)
|
||
|
{
|
||
|
label curNei = curNeighbours[neiI];
|
||
|
|
||
|
// Reject neighbours with the lower label
|
||
|
if (curNei > cellI)
|
||
|
{
|
||
|
// Get the list of search faces
|
||
|
const faceList& searchFaces = cellsFaceShapes[curNei];
|
||
|
|
||
|
forAll(searchFaces, neiFaceI)
|
||
|
{
|
||
|
if (searchFaces[neiFaceI] == curFace)
|
||
|
{
|
||
|
// Match!!
|
||
|
found = true;
|
||
|
|
||
|
// Record the neighbour cell and face
|
||
|
neiCells[faceI] = curNei;
|
||
|
faceOfNeiCell[faceI] = neiFaceI;
|
||
|
nNeighbours++;
|
||
|
|
||
|
break;
|
||
|
}
|
||
|
}
|
||
|
if (found) break;
|
||
|
}
|
||
|
// if (found) break; HJ, not needed. Never true
|
||
|
}
|
||
|
if (found) break;
|
||
|
} // End of current points
|
||
|
} // End of current faces
|
||
|
|
||
|
// Add the faces in the increasing order of neighbours
|
||
|
for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
|
||
|
{
|
||
|
// Find the lowest neighbour which is still valid
|
||
|
label nextNei = -1;
|
||
|
label minNei = cells.size();
|
||
|
|
||
|
forAll (neiCells, ncI)
|
||
|
{
|
||
|
if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
|
||
|
{
|
||
|
nextNei = ncI;
|
||
|
minNei = neiCells[ncI];
|
||
|
}
|
||
|
}
|
||
|
|
||
|
if (nextNei > -1)
|
||
|
{
|
||
|
// Add the face to the list of faces
|
||
|
allFaces_[nFaces] = curFaces[nextNei];
|
||
|
|
||
|
// Set cell-face and cell-neighbour-face to current face label
|
||
|
cells[cellI][nextNei] = nFaces;
|
||
|
cells[neiCells[nextNei]][faceOfNeiCell[nextNei]] = nFaces;
|
||
|
|
||
|
// Stop the neighbour from being used again
|
||
|
neiCells[nextNei] = -1;
|
||
|
|
||
|
// Increment number of faces counter
|
||
|
nFaces++;
|
||
|
}
|
||
|
else
|
||
|
{
|
||
|
FatalErrorIn
|
||
|
(
|
||
|
"polyMesh::polyMesh\n"
|
||
|
"(\n"
|
||
|
" const IOobject& io,\n"
|
||
|
" const pointField& points,\n"
|
||
|
" const cellShapeList& cellsAsShapes,\n"
|
||
|
" const faceListList& boundaryFaces,\n"
|
||
|
" const wordList& boundaryPatchTypes,\n"
|
||
|
" const wordList& boundaryPatchNames,\n"
|
||
|
" const word& defaultBoundaryPatchType\n"
|
||
|
")"
|
||
|
) << "Error in internal face insertion"
|
||
|
<< abort(FatalError);
|
||
|
}
|
||
|
}
|
||
|
}
|
||
|
|
||
|
// Do boundary faces
|
||
|
|
||
|
labelList patchSizes(boundaryFaces.size(), -1);
|
||
|
labelList patchStarts(boundaryFaces.size(), -1);
|
||
|
|
||
|
forAll (boundaryFaces, patchI)
|
||
|
{
|
||
|
const faceList& patchFaces = boundaryFaces[patchI];
|
||
|
|
||
|
labelList curPatchFaceCells =
|
||
|
facePatchFaceCells
|
||
|
(
|
||
|
patchFaces,
|
||
|
PointCells,
|
||
|
cellsFaceShapes,
|
||
|
patchI
|
||
|
);
|
||
|
|
||
|
// Grab the start label
|
||
|
label curPatchStart = nFaces;
|
||
|
|
||
|
forAll (patchFaces, faceI)
|
||
|
{
|
||
|
const face& curFace = patchFaces[faceI];
|
||
|
|
||
|
const label cellInside = curPatchFaceCells[faceI];
|
||
|
|
||
|
allFaces_[nFaces] = curFace;
|
||
|
|
||
|
// get faces of the cell inside
|
||
|
const faceList& facesOfCellInside = cellsFaceShapes[cellInside];
|
||
|
|
||
|
bool found = false;
|
||
|
|
||
|
forAll (facesOfCellInside, cellFaceI)
|
||
|
{
|
||
|
if (facesOfCellInside[cellFaceI] == curFace)
|
||
|
{
|
||
|
if (cells[cellInside][cellFaceI] >= 0)
|
||
|
{
|
||
|
FatalErrorIn
|
||
|
(
|
||
|
"polyMesh::polyMesh\n"
|
||
|
"(\n"
|
||
|
" const IOobject& io,\n"
|
||
|
" const pointField& points,\n"
|
||
|
" const cellShapeList& cellsAsShapes,\n"
|
||
|
" const faceListList& boundaryFaces,\n"
|
||
|
" const wordList& boundaryPatchTypes,\n"
|
||
|
" const wordList& boundaryPatchNames,\n"
|
||
|
" const word& defaultBoundaryPatchType\n"
|
||
|
")"
|
||
|
) << "Trying to specify a boundary face " << curFace
|
||
|
<< " on the face on cell " << cellInside
|
||
|
<< " which is either an internal face or already "
|
||
|
<< "belongs to some other patch. This is face "
|
||
|
<< faceI << " of patch "
|
||
|
<< patchI << " named "
|
||
|
<< boundaryPatchNames[patchI] << "."
|
||
|
<< abort(FatalError);
|
||
|
}
|
||
|
|
||
|
found = true;
|
||
|
|
||
|
cells[cellInside][cellFaceI] = nFaces;
|
||
|
|
||
|
break;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
if (!found)
|
||
|
{
|
||
|
FatalErrorIn("polyMesh::polyMesh(... construct from shapes...)")
|
||
|
<< "face " << faceI << " of patch " << patchI
|
||
|
<< " does not seem to belong to cell " << cellInside
|
||
|
<< " which, according to the addressing, "
|
||
|
<< "should be next to it."
|
||
|
<< abort(FatalError);
|
||
|
}
|
||
|
|
||
|
// increment the counter of faces
|
||
|
nFaces++;
|
||
|
}
|
||
|
|
||
|
patchSizes[patchI] = nFaces - curPatchStart;
|
||
|
patchStarts[patchI] = curPatchStart;
|
||
|
}
|
||
|
|
||
|
// Grab "non-existing" faces and put them into a default patch
|
||
|
|
||
|
label defaultPatchStart = nFaces;
|
||
|
|
||
|
forAll(cells, cellI)
|
||
|
{
|
||
|
labelList& curCellFaces = cells[cellI];
|
||
|
|
||
|
forAll(curCellFaces, faceI)
|
||
|
{
|
||
|
if (curCellFaces[faceI] == -1) // "non-existent" face
|
||
|
{
|
||
|
curCellFaces[faceI] = nFaces;
|
||
|
allFaces_[nFaces] = cellsFaceShapes[cellI][faceI];
|
||
|
|
||
|
nFaces++;
|
||
|
}
|
||
|
}
|
||
|
}
|
||
|
|
||
|
// Reset the size of the face list
|
||
|
allFaces_.setSize(nFaces);
|
||
|
faces_.reset(allFaces_, nFaces);
|
||
|
|
||
|
// Warning: Patches can only be added once the face list is
|
||
|
// completed, as they hold a subList of the face list
|
||
|
forAll (boundaryFaces, patchI)
|
||
|
{
|
||
|
// add the patch to the list
|
||
|
boundary_.set
|
||
|
(
|
||
|
patchI,
|
||
|
polyPatch::New
|
||
|
(
|
||
|
boundaryPatchTypes[patchI],
|
||
|
boundaryPatchNames[patchI],
|
||
|
patchSizes[patchI],
|
||
|
patchStarts[patchI],
|
||
|
patchI,
|
||
|
boundary_
|
||
|
)
|
||
|
);
|
||
|
|
||
|
if
|
||
|
(
|
||
|
boundaryPatchPhysicalTypes.size()
|
||
|
&& boundaryPatchPhysicalTypes[patchI].size()
|
||
|
)
|
||
|
{
|
||
|
boundary_[patchI].physicalType() =
|
||
|
boundaryPatchPhysicalTypes[patchI];
|
||
|
}
|
||
|
}
|
||
|
|
||
|
label nAllPatches = boundaryFaces.size();
|
||
|
|
||
|
if (nFaces > defaultPatchStart)
|
||
|
{
|
||
|
WarningIn("polyMesh::polyMesh(... construct from shapes...)")
|
||
|
<< "Found " << nFaces - defaultPatchStart
|
||
|
<< " undefined faces in mesh; adding to default patch." << endl;
|
||
|
|
||
|
boundary_.set
|
||
|
(
|
||
|
nAllPatches,
|
||
|
polyPatch::New
|
||
|
(
|
||
|
defaultBoundaryPatchType,
|
||
|
defaultBoundaryPatchName,
|
||
|
nFaces - defaultPatchStart,
|
||
|
defaultPatchStart,
|
||
|
boundary_.size() - 1,
|
||
|
boundary_
|
||
|
)
|
||
|
);
|
||
|
|
||
|
nAllPatches++;
|
||
|
}
|
||
|
|
||
|
// Reset the size of the boundary
|
||
|
boundary_.setSize(nAllPatches);
|
||
|
|
||
|
// Set the primitive mesh
|
||
|
initMesh(cells);
|
||
|
|
||
|
if (syncPar)
|
||
|
{
|
||
|
// Calculate topology for the patches (processor-processor comms etc.)
|
||
|
boundary_.updateMesh();
|
||
|
|
||
|
// Calculate the geometry for the patches (transformation tensors etc.)
|
||
|
boundary_.calcGeometry();
|
||
|
}
|
||
|
|
||
|
if (debug)
|
||
|
{
|
||
|
if (checkMesh())
|
||
|
{
|
||
|
Info << "Mesh OK" << endl;
|
||
|
}
|
||
|
}
|
||
|
}
|
||
|
|
||
|
|
||
|
// ************************************************************************* //
|