This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/utilities/mesh/conversion/Optional/ccm26ToFoam/ccm26ToFoam.C
2019-02-04 21:58:51 +01:00

1118 lines
30 KiB
C++

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.1
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Description
Reads CCM files as written by Prostar/ccm using ccm 2.6 (not 2.4)
- does polyhedral mesh
- does not handle 'interfaces' (couples)
- does not handle cyclics. Use createPatch to recreate these
- does not do data
- does patch names only if they are in the problem description
\*---------------------------------------------------------------------------*/
#include "ListOps.H"
#include "argList.H"
#include "foamTime.H"
#include "fvMesh.H"
#include "volFields.H"
#include "emptyPolyPatch.H"
#include "symmetryPolyPatch.H"
#include "wallPolyPatch.H"
#include "SortableList.H"
#include "cellSet.H"
#include <ccmio.h>
#include <vector>
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
static char const kDefaultState[] = "default";
static int const kVertOffset = 2;
// Determine upper-triangular order for internal faces:
labelList getInternalFaceOrder
(
const cellList& cells,
const labelList& owner,
const labelList& neighbour
)
{
labelList oldToNew(owner.size(), -1);
// First unassigned face
label newFaceI = 0;
forAll(cells, cellI)
{
const labelList& cFaces = cells[cellI];
SortableList<label> nbr(cFaces.size());
forAll(cFaces, i)
{
label faceI = cFaces[i];
label nbrCellI = neighbour[faceI];
if (nbrCellI != -1)
{
// Internal face. Get cell on other side.
if (nbrCellI == cellI)
{
nbrCellI = owner[faceI];
}
if (cellI < nbrCellI)
{
// CellI is master
nbr[i] = nbrCellI;
}
else
{
// nbrCell is master. Let it handle this face.
nbr[i] = -1;
}
}
else
{
// External face. Do later.
nbr[i] = -1;
}
}
nbr.sort();
forAll(nbr, i)
{
if (nbr[i] != -1)
{
oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
}
}
}
// Keep boundary faces in same order.
for (label faceI = newFaceI; faceI < owner.size(); faceI++)
{
oldToNew[faceI] = faceI;
}
return oldToNew;
}
void storeCellInZone
(
const label cellI,
const label cellType,
Map<label>& typeToZone,
List<dynamicLabelList >& zoneCells
)
{
if (cellType >= 0)
{
Map<label>::iterator zoneFnd = typeToZone.find(cellType);
if (zoneFnd == typeToZone.end())
{
// New type. Allocate zone for it.
zoneCells.setSize(zoneCells.size() + 1);
label zoneI = zoneCells.size()-1;
Info<< "Mapping type " << cellType << " to Foam cellZone "
<< zoneI << endl;
typeToZone.insert(cellType, zoneI);
zoneCells[zoneI].append(cellI);
}
else
{
// Existing zone for type
zoneCells[zoneFnd()].append(cellI);
}
}
}
void CheckError(CCMIOError const &err, const Foam::string& str)
{
if (err != kCCMIONoErr)
{
FatalErrorIn("CheckError")
<< str << exit(FatalError);
}
}
void ReadVertices
(
CCMIOError &err,
CCMIOID &vertices,
labelList& foamPointMap,
pointField& foamPoints
)
{
// Read the vertices. This involves reading both the vertex data and
// the map, which maps the index into the data array with the ID number.
// As we process the vertices we need to be sure to scale them by the
// appropriate scaling factor. The offset is just to show you can read
// any chunk. Normally this would be in a for loop.
CCMIOSize nVertices;
CCMIOEntitySize(&err, vertices, &nVertices, nullptr);
List<int> mapData(nVertices, 0);
List<float> verts(3*nVertices, 0);
int offset = 0;
int offsetPlusSize = offset+nVertices;
int dims = 1;
float scale;
CCMIOID mapID;
CCMIOReadVerticesf(&err, vertices, &dims, &scale, &mapID, verts.begin(),
offset, offsetPlusSize);
CCMIOReadMap(&err, mapID, mapData.begin(), offset, offsetPlusSize);
//CCMIOSize size;
//CCMIOEntityDescription(&err, vertices, &size, nullptr);
//char *desc = new char[size + 1];
//CCMIOEntityDescription(&err, vertices, nullptr, desc);
//Pout<< "label: '" << desc << "'" << endl;
//delete [] desc;
// Convert to foamPoints
foamPoints.setSize(nVertices);
foamPoints = vector::zero;
foamPointMap.setSize(nVertices);
forAll(foamPointMap, i)
{
foamPointMap[i] = mapData[i];
for (direction cmpt = 0; cmpt < dims; cmpt++)
{
foamPoints[i][cmpt] = verts[dims*i + cmpt]*scale;
}
}
}
void ReadProblem
(
CCMIOError &err,
CCMIOID& problem,
const Map<label>& prostarToFoamPatch,
Map<word>& foamCellTypeNames,
wordList& foamPatchTypes,
wordList& foamPatchNames
)
{
// ... walk through each cell type and print it...
CCMIOID next;
int i = 0;
while
(
CCMIONextEntity(nullptr, problem, kCCMIOCellType, &i, &next)
== kCCMIONoErr
)
{
char *name;
int size;
int cellType;
// ... if it has a material type. (Note that we do not pass in
// an array to get the name because we do not know how long the
// string is yet. Many parameters to CCMIO functions that
// return
// data can be nullptr if that data is not needed.)
if
(
CCMIOReadOptstr(nullptr, next, "MaterialType", &size, nullptr)
== kCCMIONoErr
)
{
name = new char[size + 1];
CCMIOReadOptstr(&err, next, "MaterialType", &size, name);
CCMIOGetEntityIndex(&err, next, &cellType);
foamCellTypeNames.insert(cellType, name);
Pout<< "Celltype:" << cellType << " name:" << name << endl;
delete [] name;
}
}
// ... walk through each region description and print it...
CCMIOID boundary;
label regionI = 0;
int k = 0;
while
(
CCMIONextEntity(nullptr, problem, kCCMIOBoundaryRegion, &k, &boundary)
== kCCMIONoErr
)
{
// Index of foam patch
label foamPatchI = -1;
// Read prostar id
int prostarI = -1;
if
(
CCMIOReadOpti(nullptr, boundary, "ProstarRegionNumber", &prostarI)
== kCCMIONoErr
)
{
Pout<< "For region:" << regionI
<< " found ProstarRegionNumber:" << prostarI << endl;
}
else
{
prostarI = regionI;
Pout<< "For region:" << regionI
<< "did not find ProstarRegionNumber entry. Assuming "
<< prostarI << endl;
}
if (prostarToFoamPatch.found(prostarI))
{
foamPatchI = prostarToFoamPatch[prostarI];
// Read boundary type
int size;
if
(
CCMIOReadOptstr(nullptr, boundary, "BoundaryType", &size, nullptr)
== kCCMIONoErr
)
{
char* s = new char[size + 1];
CCMIOReadOptstr(nullptr, boundary, "BoundaryType", &size, s);
s[size] = '\0';
foamPatchTypes[foamPatchI] = string::validate<word>(string(s));
delete [] s;
}
//foamPatchMap.append(prostarI);
// Read boundary name:
// - from BoundaryName field (Prostar)
// - from 'Label' field (ccm+?)
// - make up one from prostar id.
if
(
CCMIOReadOptstr(nullptr, boundary, "BoundaryName", &size, nullptr)
== kCCMIONoErr
)
{
char* name = new char[size + 1];
CCMIOReadOptstr(nullptr, boundary, "BoundaryName", &size, name);
name[size] = '\0';
foamPatchNames[foamPatchI] =
string::validate<word>(string(name));
delete [] name;
}
else if
(
CCMIOReadOptstr(nullptr, boundary, "Label", &size, nullptr)
== kCCMIONoErr
)
{
char* name = new char[size + 1];
CCMIOReadOptstr(nullptr, boundary, "Label", &size, name);
name[size] = '\0';
foamPatchNames[foamPatchI] =
string::validate<word>(string(name));
delete [] name;
}
else
{
foamPatchNames[foamPatchI] =
foamPatchTypes[foamPatchI]
+ Foam::name(foamPatchI);
Pout<< "Made up name:" << foamPatchNames[foamPatchI]
<< endl;
}
Pout<< "Read patch:" << foamPatchI
<< " name:" << foamPatchNames[foamPatchI]
<< " foamPatchTypes:" << foamPatchTypes[foamPatchI]
<< endl;
}
regionI++;
}
}
void ReadCells
(
CCMIOError &err,
CCMIOID& topology,
labelList& foamCellMap,
labelList& foamCellType,
Map<label>& prostarToFoamPatch,
dynamicLabelList& foamPatchSizes,
dynamicLabelList& foamPatchStarts,
labelList& foamFaceMap,
labelList& foamOwner,
labelList& foamNeighbour,
faceList& foamFaces
)
{
// Read the cells.
// ~~~~~~~~~~~~~~~
// Store cell IDs so that we know what cells are in
// this post data.
CCMIOID id;
CCMIOGetEntity(&err, topology, kCCMIOCells, 0, &id);
CCMIOSize nCells;
CCMIOEntitySize(&err, id, &nCells, nullptr);
std::vector<int> mapData(nCells);
std::vector<int> cellType(nCells);
CCMIOID mapID;
CCMIOReadCells(&err, id, &mapID, &cellType[0], 0, nCells);
CCMIOReadMap(&err, mapID, &mapData[0], 0, nCells);
CheckError(err, "Error reading cells");
foamCellMap.setSize(nCells);
foamCellType.setSize(nCells);
forAll(foamCellMap, i)
{
foamCellMap[i] = mapData[i];
foamCellType[i] = cellType[i];
}
// Read the internal faces.
// ~~~~~~~~~~~~~~~~~~~~~~~~
CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
CCMIOSize nInternalFaces;
CCMIOEntitySize(&err, id, &nInternalFaces, nullptr);
Pout<< "nInternalFaces:" << nInternalFaces << endl;
// Determine patch sizes before reading internal faces
label foamNFaces = nInternalFaces;
int index = 0;
while
(
CCMIONextEntity(nullptr, topology, kCCMIOBoundaryFaces, &index, &id)
== kCCMIONoErr
)
{
CCMIOSize size;
CCMIOEntitySize(&err, id, &size, nullptr);
Pout<< "Read kCCMIOBoundaryFaces entry with " << size
<< " faces." << endl;
foamPatchStarts.append(foamNFaces);
foamPatchSizes.append(size);
foamNFaces += size;
}
foamPatchStarts.shrink();
foamPatchSizes.shrink();
Pout<< "patchSizes:" << foamPatchSizes << endl;
Pout<< "patchStarts:" << foamPatchStarts << endl;
Pout<< "nFaces:" << foamNFaces << endl;
mapData.resize(nInternalFaces);
CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
CCMIOSize size;
CCMIOReadFaces(&err, id, kCCMIOInternalFaces, nullptr, &size, nullptr,
kCCMIOStart, kCCMIOEnd);
std::vector<int> faces(size);
CCMIOReadFaces(&err, id, kCCMIOInternalFaces, &mapID, nullptr, &faces[0],
kCCMIOStart, kCCMIOEnd);
std::vector<int> faceCells(2*nInternalFaces);
CCMIOReadFaceCells(&err, id, kCCMIOInternalFaces, &faceCells[0],
kCCMIOStart, kCCMIOEnd);
CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
CheckError(err, "Error reading internal faces");
// Copy into Foam lists
foamFaceMap.setSize(foamNFaces);
foamFaces.setSize(foamNFaces);
foamOwner.setSize(foamNFaces);
foamNeighbour.setSize(foamNFaces);
unsigned int pos = 0;
for (unsigned int faceI = 0; faceI < nInternalFaces; faceI++)
{
foamFaceMap[faceI] = mapData[faceI];
foamOwner[faceI] = faceCells[2*faceI];
foamNeighbour[faceI] = faceCells[2*faceI+1];
face& f = foamFaces[faceI];
f.setSize(faces[pos++]);
forAll(f, fp)
{
f[fp] = faces[pos++];
}
}
// Read the boundary faces.
// ~~~~~~~~~~~~~~~~~~~~~~~~
index = 0;
label regionI = 0;
while
(
CCMIONextEntity(nullptr, topology, kCCMIOBoundaryFaces, &index, &id)
== kCCMIONoErr
)
{
CCMIOSize nFaces;
CCMIOEntitySize(&err, id, &nFaces, nullptr);
mapData.resize(nFaces);
faceCells.resize(nFaces);
CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, nullptr, &size, nullptr,
kCCMIOStart, kCCMIOEnd);
faces.resize(size);
CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, &mapID, nullptr, &faces[0],
kCCMIOStart, kCCMIOEnd);
CCMIOReadFaceCells(&err, id, kCCMIOBoundaryFaces, &faceCells[0],
kCCMIOStart, kCCMIOEnd);
CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);
CheckError(err, "Error reading boundary faces");
// Read prostar id
int prostarI;
if
(
CCMIOReadOpti(nullptr, id, "ProstarRegionNumber", &prostarI)
== kCCMIONoErr
)
{
Pout<< "For region:" << regionI
<< " found ProstarRegionNumber:" << prostarI << endl;
}
else
{
prostarI = regionI;
Pout<< "For region:" << regionI
<< " did not find ProstarRegionNumber entry. Assuming "
<< prostarI << endl;
}
prostarToFoamPatch.insert(prostarI, regionI);
Pout<< "region:" << regionI
<< " ProstarRegionNumber:" << prostarI
<< " foamPatchStart:"
<< foamPatchStarts[regionI]
<< " size:"
<< foamPatchSizes[regionI]
<< endl;
// Copy into Foam list.
unsigned int pos = 0;
for (unsigned int i = 0; i < nFaces; i++)
{
label foamFaceI = foamPatchStarts[regionI] + i;
foamFaceMap[foamFaceI] = mapData[i];
foamOwner[foamFaceI] = faceCells[i];
foamNeighbour[foamFaceI] = -1;
face& f = foamFaces[foamFaceI];
f.setSize(faces[pos++]);
forAll(f, fp)
{
f[fp] = faces[pos++];
}
}
Pout<< endl;
regionI++;
}
}
// Main program:
int main(int argc, char *argv[])
{
argList::noParallel();
argList::validArgs.append("ccm26 file");
# include "setRootCase.H"
# include "createTime.H"
// Foam mesh data
// ~~~~~~~~~~~~~~
// Coordinates
pointField foamPoints;
labelList foamPointMap;
// Cell type
labelList foamCellType;
labelList foamCellMap;
// Patching info
Map<label> prostarToFoamPatch;
dynamicLabelList foamPatchSizes;
dynamicLabelList foamPatchStarts;
// Face connectivity
labelList foamFaceMap;
labelList foamOwner;
labelList foamNeighbour;
faceList foamFaces;
// Celltypes (not the names of the zones; just the material type)
// and patch type names
Map<word> foamCellTypeNames;
wordList foamPatchTypes;
wordList foamPatchNames;
{
fileName ccmFile(args.additionalArgs()[0]);
if (!isFile(ccmFile))
{
FatalErrorIn(args.executable())
<< "Cannot read file " << ccmFile
<< exit(FatalError);
}
word ccmExt = ccmFile.ext();
if (ccmExt != "ccm" && ccmExt != "ccmg")
{
FatalErrorIn(args.executable())
<< "Illegal extension " << ccmExt << " for file " << ccmFile
<< nl << "Allowed extensions are '.ccm', '.ccmg'"
<< exit(FatalError);
}
// Open the file. Because we did not initialize 'err' we need to pass
// in nullptr (which always means kCCMIONoErr) and then assign the return
// value to 'err'.).
CCMIOID root;
CCMIOError err =
CCMIOOpenFile(nullptr, ccmFile.c_str(), kCCMIORead, &root);
// We are going to assume that we have a state with a known name.
// We could instead use CCMIONextEntity() to walk through all the
// states in the file and present the list to the user for selection.
CCMIOID state;
int stateI = 0;
CCMIONextEntity(&err, root, kCCMIOState, &stateI, &state);
CheckError(err, "Error opening state");
unsigned int size;
CCMIOEntityDescription(&err, state, &size, nullptr);
char *desc = new char[size + 1];
CCMIOEntityDescription(&err, state, nullptr, desc);
Pout<< "Reading state '" << kDefaultState << "' (" << desc << ")"
<< endl;
delete [] desc;
// Find the first processor (i has previously been initialized to 0)
// and read the mesh and solution information.
int i = 0;
CCMIOID processor;
CCMIONextEntity(&err, state, kCCMIOProcessor, &i, &processor);
CCMIOID solution, vertices, topology;
CCMIOReadProcessor
(
&err,
processor,
&vertices,
&topology,
nullptr,
&solution
);
if (err != kCCMIONoErr)
{
// Maybe no solution; try again
err = kCCMIONoErr;
CCMIOReadProcessor
(
&err,
processor,
&vertices,
&topology,
nullptr,
nullptr
);
if (err != kCCMIONoErr)
{
FatalErrorIn(args.executable())
<< "Could not read the file."
<< exit(FatalError);
}
}
ReadVertices(err, vertices, foamPointMap, foamPoints);
Pout<< "nPoints:" << foamPoints.size() << endl
<< "bounding box:" << boundBox(foamPoints) << endl
<< endl;
ReadCells
(
err,
topology,
foamCellMap,
foamCellType,
prostarToFoamPatch,
foamPatchSizes,
foamPatchStarts,
foamFaceMap,
foamOwner,
foamNeighbour,
foamFaces
);
Pout<< "nCells:" << foamCellMap.size() << endl
<< "nFaces:" << foamOwner.size() << endl
<< "nPatches:" << foamPatchStarts.size() << endl
<< "nInternalFaces:" << foamPatchStarts[0] << endl
<< endl;
// Create some default patch names/types. These will be overwritten
// by any problem desciption (if it is there)
foamPatchTypes.setSize(foamPatchStarts.size());
foamPatchNames.setSize(foamPatchStarts.size());
forAll(foamPatchNames, i)
{
foamPatchNames[i] = word("patch") + Foam::name(i);
foamPatchTypes[i] = "patch";
}
// Get problem description
CCMIOID problem;
int problemI = 0;
CCMIONextEntity
(
&err,
root,
kCCMIOProblemDescription,
&problemI,
&problem
);
CheckError(err, "Error stepping to first problem description");
if (CCMIOIsValidEntity(problem)) // if we have a problem description
{
ReadProblem
(
err,
problem,
prostarToFoamPatch,
foamCellTypeNames,
foamPatchTypes,
foamPatchNames
);
}
CCMIOCloseFile(&err, vertices);
CCMIOCloseFile(&err, topology);
CCMIOCloseFile(&err, solution);
CCMIOCloseFile(&err, root);
}
Pout<< "foamPatchNames:" << foamPatchNames << endl;
Pout<< "foamOwner : min:" << min(foamOwner)
<< " max:" << max(foamOwner)
<< nl
<< "foamNeighbour : min:" << min(foamNeighbour)
<< " max:" << max(foamNeighbour)
<< nl
<< "foamCellType : min:" << min(foamCellType)
<< " max:" << max(foamCellType)
<< nl << endl;
// We now have extracted all info from CCMIO:
// - coordinates (points)
// - face to point addressing (faces)
// - face to cell addressing (owner, neighbour)
// - cell based data (cellId)
// Renumber vertex labels to Foam point labels
{
label maxCCMPointI = max(foamPointMap);
labelList toFoamPoints(invert(maxCCMPointI+1, foamPointMap));
forAll(foamFaces, faceI)
{
inplaceRenumber(toFoamPoints, foamFaces[faceI]);
}
}
// Renumber cell labels
{
label maxCCMCellI = max(foamCellMap);
labelList toFoamCells(invert(maxCCMCellI+1, foamCellMap));
inplaceRenumber(toFoamCells, foamOwner);
inplaceRenumber(toFoamCells, foamNeighbour);
}
//
// Basic mesh info complete. Now convert to Foam convention:
// - owner < neighbour
// - face vertices such that normal points away from owner
// - order faces: upper-triangular for internal faces; boundary faces after
// internal faces
//
// Set owner/neighbour so owner < neighbour
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(foamNeighbour, faceI)
{
label nbr = foamNeighbour[faceI];
label own = foamOwner[faceI];
if (nbr >= foamCellType.size() || own >= foamCellType.size())
{
FatalErrorIn(args.executable())
<< "face:" << faceI
<< " nbr:" << nbr
<< " own:" << own
<< " nCells:" << foamCellType.size()
<< exit(FatalError);
}
if (nbr >= 0)
{
if (nbr < own)
{
foamOwner[faceI] = foamNeighbour[faceI];
foamNeighbour[faceI] = own;
foamFaces[faceI] = foamFaces[faceI].reverseFace();
}
}
// And check the face
const face& f = foamFaces[faceI];
forAll(f, fp)
{
if (f[fp] < 0 || f[fp] >= foamPoints.size())
{
FatalErrorIn(args.executable()) << "face:" << faceI
<< " f:" << f << abort(FatalError);
}
}
}
// Do upper-triangular ordering
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
// Create cells (inverse of face-to-cell addressing)
cellList foamCells;
primitiveMesh::calcCells
(
foamCells,
foamOwner,
foamNeighbour,
foamCellType.size()
);
// Determine face order for upper-triangular ordering
labelList oldToNew
(
getInternalFaceOrder
(
foamCells,
foamOwner,
foamNeighbour
)
);
// Reorder faces accordingly
forAll(foamCells, cellI)
{
inplaceRenumber(oldToNew, foamCells[cellI]);
}
// Reorder faces.
inplaceReorder(oldToNew, foamFaces);
inplaceReorder(oldToNew, foamOwner);
inplaceReorder(oldToNew, foamNeighbour);
}
// Construct fvMesh
// ~~~~~~~~~~~~~~~~
fvMesh mesh
(
IOobject
(
fvMesh::defaultRegion,
runTime.constant(),
runTime
),
xferMove<pointField>(foamPoints),
xferMove<faceList>(foamFaces),
xferCopy<labelList>(foamOwner),
xferMove<labelList>(foamNeighbour)
);
// Create patches. Use patch types to determine what Foam types to generate.
List<polyPatch*> newPatches(foamPatchNames.size());
label meshFaceI = foamPatchStarts[0];
forAll(newPatches, patchI)
{
const word& patchName = foamPatchNames[patchI];
const word& patchType = foamPatchTypes[patchI];
Pout<< "Patch:" << patchName << " start at:" << meshFaceI
<< " size:" << foamPatchSizes[patchI]
<< " end at:" << meshFaceI+foamPatchSizes[patchI]
<< endl;
if (patchType == "wall")
{
newPatches[patchI] =
new wallPolyPatch
(
patchName,
foamPatchSizes[patchI],
meshFaceI,
patchI,
mesh.boundaryMesh()
);
}
else if (patchType == "symmetryplane")
{
newPatches[patchI] =
new symmetryPolyPatch
(
patchName,
foamPatchSizes[patchI],
meshFaceI,
patchI,
mesh.boundaryMesh()
);
}
else if (patchType == "empty")
{
// Note: not ccm name, introduced by us above.
newPatches[patchI] =
new emptyPolyPatch
(
patchName,
foamPatchSizes[patchI],
meshFaceI,
patchI,
mesh.boundaryMesh()
);
}
else
{
// All other ccm types become straight polyPatch:
// 'inlet', 'outlet', 'pressured'.
newPatches[patchI] =
new polyPatch
(
patchName,
foamPatchSizes[patchI],
meshFaceI,
patchI,
mesh.boundaryMesh()
);
}
meshFaceI += foamPatchSizes[patchI];
}
if (meshFaceI != foamOwner.size())
{
FatalErrorIn(args.executable())
<< "meshFaceI:" << meshFaceI
<< " nFaces:" << foamOwner.size()
<< abort(FatalError);
}
mesh.addFvPatches(newPatches);
// Construct sets and zones from types
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label maxType = max(foamCellType);
label minType = min(foamCellType);
if (maxType > minType)
{
// From foamCellType physical region to Foam cellZone
Map<label> typeToZone;
// Storage for cell zones.
List<dynamicLabelList > zoneCells(0);
forAll(foamCellType, cellI)
{
storeCellInZone
(
cellI,
foamCellType[cellI],
typeToZone,
zoneCells
);
}
mesh.cellZones().clear();
mesh.cellZones().setSize(typeToZone.size());
label nValidCellZones = 0;
forAllConstIter(Map<label>, typeToZone, iter)
{
label type = iter.key();
label zoneI = iter();
zoneCells[zoneI].shrink();
word zoneName = "cellZone_" + name(type);
Info<< "Writing zone " << type
<< " size " << zoneCells[zoneI].size()
<< " to cellZone "
<< zoneName << " and cellSet " << zoneName
<< endl;
cellSet cset(mesh, zoneName, labelHashSet(zoneCells[zoneI]));
cset.write();
mesh.cellZones().set
(
nValidCellZones,
new cellZone
(
zoneName,
zoneCells[zoneI],
nValidCellZones,
mesh.cellZones()
)
);
nValidCellZones++;
}
mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
}
Info<< "Writing mesh to " << mesh.objectRegistry::objectPath()
<< "..." << nl << endl;
// Construct field with calculated bc to hold Star cell Id.
volScalarField cellIdField
(
IOobject
(
"cellId",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("cellId", dimless, 0.0)
);
forAll(foamCellMap, cellI)
{
cellIdField[cellI] = foamCellMap[cellI];
}
// Construct field with calculated bc to hold cell type.
volScalarField cellTypeField
(
IOobject
(
"cellType",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("cellType", dimless, 0.0)
);
forAll(foamCellType, cellI)
{
cellTypeField[cellI] = foamCellType[cellI];
}
Info<< "Writing cellIds as volScalarField to " << cellIdField.objectPath()
<< "..." << nl << endl;
mesh.write();
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //