/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2015-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
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 3 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, see _sidePatch
- extrude edges inbetween different faceZones as a
(nonuniformTransform)cyclic \_\
- extrudes into master direction (i.e. away from the owner cell
if flipMap is false)
\verbatim
Internal face extrusion
-----------------------
+-------------+
| |
| |
+---AAAAAAA---+
| |
| |
+-------------+
AAA=faceZone to extrude.
For the case of no flipMap the extrusion starts at owner and extrudes
into the space of the neighbour:
+CCCCCCC+
| | <= extruded mesh
+BBBBBBB+
+-------------+
| |
| (neighbour) |
|___CCCCCCC___| <= original mesh (with 'baffles' added)
| BBBBBBB |
|(owner side) |
| |
+-------------+
BBB=mapped between owner on original mesh and new extrusion.
(zero offset)
CCC=mapped between neighbour on original mesh and new extrusion
(offset due to the thickness of the extruded mesh)
For the case of flipMap the extrusion is the other way around: from the
neighbour side into the owner side.
Boundary face extrusion
-----------------------
+--AAAAAAA--+
| |
| |
+-----------+
AAA=faceZone to extrude. E.g. slave side is owner side (no flipmap)
becomes
+CCCCCCC+
| | <= extruded mesh
+BBBBBBB+
+--BBBBBBB--+
| | <= original mesh
| |
+-----------+
BBB=mapped between original mesh and new extrusion
CCC=polypatch
Notes:
- when extruding cyclics with only one cell inbetween it does not
detect this as a cyclic since the face is the same face. It will
only work if the coupled edge extrudes a different face so if there
are more than 1 cell inbetween.
\endverbatim
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "fvMesh.H"
#include "polyTopoChange.H"
#include "OFstream.H"
#include "meshTools.H"
#include "mappedWallPolyPatch.H"
#include "createShellMesh.H"
#include "syncTools.H"
#include "cyclicPolyPatch.H"
#include "wedgePolyPatch.H"
#include "nonuniformTransformCyclicPolyPatch.H"
#include "extrudeModel.H"
#include "globalIndex.H"
#include "faceSet.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "pointFields.H"
//#include "ReadFields.H"
#include "fvMeshTools.H"
#include "OBJstream.H"
#include "PatchTools.H"
#include "processorMeshes.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
label findPatchID(const UList& newPatches, const word& name)
{
forAll(newPatches, i)
{
if (newPatches[i]->name() == name)
{
return i;
}
}
return -1;
}
#ifdef FULLDEBUG
void printPatches(Ostream& os, const UList& newPatches)
{
for (const polyPatch* ppPtr : newPatches)
{
const polyPatch& pp = *ppPtr;
os << " name:" << pp.name() << " index:" << pp.index()
<< " start:" << pp.start() << " size:" << pp.size()
<< " type:" << pp.type() << nl;
}
}
#endif
template
label addPatch
(
const polyBoundaryMesh& patches,
const word& patchName,
DynamicList& newPatches
)
{
label patchi = findPatchID(newPatches, patchName);
if (patchi != -1)
{
if (isA(*newPatches[patchi]))
{
// Already there
return patchi;
}
else
{
FatalErrorInFunction
<< "Already have patch " << patchName
<< " but of type " << newPatches[patchi]->type()
<< exit(FatalError);
}
}
patchi = newPatches.size();
label startFacei = 0;
if (patchi > 0)
{
const polyPatch& pp = *newPatches.last();
startFacei = pp.start()+pp.size();
}
#ifdef FULLDEBUG
Pout<< "addPatch : starting newPatches:"
<< " patch:" << patchi << " startFace:" << startFacei << nl;
printPatches(Pout, newPatches);
Pout<< "*** end of addPatch:" << endl;
#endif
newPatches.append
(
polyPatch::New
(
PatchType::typeName,
patchName,
0, // size
startFacei, // nFaces
patchi,
patches
).ptr()
);
return patchi;
}
template
label addPatch
(
const polyBoundaryMesh& patches,
const word& patchName,
const dictionary& dict,
DynamicList& newPatches
)
{
label patchi = findPatchID(newPatches, patchName);
if (patchi != -1)
{
if (isA(*newPatches[patchi]))
{
// Already there
return patchi;
}
else
{
FatalErrorInFunction
<< "Already have patch " << patchName
<< " but of type " << newPatches[patchi]->type()
<< exit(FatalError);
}
}
patchi = newPatches.size();
label startFacei = 0;
if (patchi > 0)
{
const polyPatch& pp = *newPatches.last();
startFacei = pp.start()+pp.size();
}
#ifdef FULLDEBUG
Pout<< "addPatch : starting newPatches:"
<< " patch:" << patchi << " startFace:" << startFacei << nl;
printPatches(Pout, newPatches);
Pout<< "*** end of addPatch:" << endl;
#endif
dictionary patchDict(dict);
patchDict.set("type", PatchType::typeName);
patchDict.set("nFaces", 0);
patchDict.set("startFace", startFacei);
newPatches.append
(
polyPatch::New
(
patchName,
patchDict,
patchi,
patches
).ptr()
);
return patchi;
}
// Remove zero-sized patches
void deleteEmptyPatches(fvMesh& mesh)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
wordList masterNames;
if (Pstream::master())
{
masterNames = patches.names();
}
Pstream::broadcast(masterNames);
labelList oldToNew(patches.size(), -1);
label usedI = 0;
label notUsedI = patches.size();
// Add all the non-empty, non-processor patches
forAll(masterNames, masterI)
{
label patchi = patches.findPatchID(masterNames[masterI]);
if (patchi != -1)
{
if (isA(patches[patchi]))
{
// Similar named processor patch? Not 'possible'.
if (patches[patchi].empty())
{
Pout<< "Deleting processor patch " << patchi
<< " name:" << patches[patchi].name()
<< endl;
oldToNew[patchi] = --notUsedI;
}
else
{
oldToNew[patchi] = usedI++;
}
}
else
{
// Common patch.
if (returnReduceAnd(patches[patchi].empty()))
{
Pout<< "Deleting patch " << patchi
<< " name:" << patches[patchi].name()
<< endl;
oldToNew[patchi] = --notUsedI;
}
else
{
oldToNew[patchi] = usedI++;
}
}
}
}
// Add remaining patches at the end
forAll(patches, patchi)
{
if (oldToNew[patchi] == -1)
{
// Unique to this processor. Note: could check that these are
// only processor patches.
if (patches[patchi].empty())
{
Pout<< "Deleting processor patch " << patchi
<< " name:" << patches[patchi].name()
<< endl;
oldToNew[patchi] = --notUsedI;
}
else
{
oldToNew[patchi] = usedI++;
}
}
}
fvMeshTools::reorderPatches(mesh, oldToNew, usedI, true);
}
// Check zone either all internal or all external faces
void checkZoneInside
(
const polyMesh& mesh,
const wordList& zoneNames,
const labelList& zoneID,
const labelList& extrudeMeshFaces,
const boolList& isInternal
)
{
forAll(zoneNames, i)
{
if (isInternal[i])
{
Info<< "Zone " << zoneNames[i] << " has internal faces" << endl;
}
else
{
Info<< "Zone " << zoneNames[i] << " has boundary faces" << endl;
}
}
forAll(extrudeMeshFaces, i)
{
label facei = extrudeMeshFaces[i];
label zoneI = zoneID[i];
if (isInternal[zoneI] != mesh.isInternalFace(facei))
{
FatalErrorInFunction
<< "Zone " << zoneNames[zoneI]
<< " is not consistently all internal or all boundary faces."
<< " Face " << facei << " at " << mesh.faceCentres()[facei]
<< " is the first occurrence."
<< exit(FatalError);
}
}
}
// Calculate global pp faces per pp edge.
labelListList globalEdgeFaces
(
const polyMesh& mesh,
const globalIndex& globalFaces,
const primitiveFacePatch& pp,
const labelList& ppMeshEdges
)
{
// From mesh edge to global pp face labels.
labelListList globalEdgeFaces(ppMeshEdges.size());
const labelListList& edgeFaces = pp.edgeFaces();
forAll(edgeFaces, edgeI)
{
// Store pp face and processor as unique tag.
globalEdgeFaces[edgeI] = globalFaces.toGlobal(edgeFaces[edgeI]);
}
// Synchronise across coupled edges.
syncTools::syncEdgeList
(
mesh,
ppMeshEdges,
globalEdgeFaces,
ListOps::uniqueEqOp(),
labelList() // null value
);
return globalEdgeFaces;
}
// Find a patch face that is not extruded. Return -1 if not found.
label findUncoveredPatchFace
(
const fvMesh& mesh,
const labelUIndList& extrudeMeshFaces, // mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
extrudeFaceSet.insert(extrudeMeshFaces);
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
forAll(eFaces, i)
{
label facei = eFaces[i];
label patchi = pbm.whichPatch(facei);
if
(
patchi != -1
&& !pbm[patchi].coupled()
&& !extrudeFaceSet.found(facei)
)
{
return facei;
}
}
return -1;
}
// Same as findUncoveredPatchFace, except explicitly checks for cyclic faces
label findUncoveredCyclicPatchFace
(
const fvMesh& mesh,
const labelUIndList& extrudeMeshFaces, // mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
extrudeFaceSet.insert(extrudeMeshFaces);
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
const labelList& eFaces = mesh.edgeFaces()[meshEdgeI];
forAll(eFaces, i)
{
label facei = eFaces[i];
label patchi = pbm.whichPatch(facei);
if
(
patchi != -1
&& isA(pbm[patchi])
&& !extrudeFaceSet.found(facei)
)
{
return facei;
}
}
return -1;
}
// Calculate per edge min and max zone
void calcEdgeMinMaxZone
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const labelList& extrudeMeshEdges,
const labelList& zoneID,
const mapDistribute& extrudeEdgeFacesMap,
const labelListList& extrudeEdgeGlobalFaces,
labelList& minZoneID,
labelList& maxZoneID
)
{
// Get zoneIDs in extrudeEdgeGlobalFaces order
labelList mappedZoneID(zoneID);
extrudeEdgeFacesMap.distribute(mappedZoneID);
// Get min and max zone per edge
minZoneID.setSize(extrudeEdgeGlobalFaces.size(), labelMax);
maxZoneID.setSize(extrudeEdgeGlobalFaces.size(), labelMin);
forAll(extrudeEdgeGlobalFaces, edgeI)
{
const labelList& eFaces = extrudeEdgeGlobalFaces[edgeI];
if (eFaces.size())
{
forAll(eFaces, i)
{
label zoneI = mappedZoneID[eFaces[i]];
minZoneID[edgeI] = Foam::min(minZoneID[edgeI], zoneI);
maxZoneID[edgeI] = Foam::max(maxZoneID[edgeI], zoneI);
}
}
}
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
minZoneID,
minEqOp(),
labelMax // null value
);
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
maxZoneID,
maxEqOp(),
labelMin // null value
);
}
// Count the number of faces in patches that need to be created. Calculates:
// zoneSidePatch[zoneI] : the number of side faces to be created
// zoneZonePatch[zoneA,zoneB] : the number of faces inbetween zoneA and B
// Since this only counts we're not taking the processor patches into
// account.
void countExtrudePatches
(
const fvMesh& mesh,
const label nZones,
const primitiveFacePatch& extrudePatch,
const labelList& extrudeMeshFaces,
const labelList& extrudeMeshEdges,
const labelListList& extrudeEdgeGlobalFaces,
const labelList& minZoneID,
const labelList& maxZoneID,
labelList& zoneSidePatch,
labelList& zoneZonePatch
)
{
// Check on master edge for use of zones. Since we only want to know
// whether they are being used at all no need to accurately count on slave
// edge as well. Just add all together at the end of this routine so it
// gets detected at least.
forAll(extrudePatch.edgeFaces(), edgeI)
{
const labelList& eFaces = extrudePatch.edgeFaces()[edgeI];
if (eFaces.size() == 2)
{
// Internal edge - check if inbetween different zones.
if (minZoneID[edgeI] != maxZoneID[edgeI])
{
zoneZonePatch[minZoneID[edgeI]*nZones+maxZoneID[edgeI]]++;
}
}
else if
(
eFaces.size() == 1
&& extrudeEdgeGlobalFaces[edgeI].size() == 2
)
{
// Coupled edge - check if inbetween different zones.
if (minZoneID[edgeI] != maxZoneID[edgeI])
{
const edge& e = extrudePatch.edges()[edgeI];
const pointField& pts = extrudePatch.localPoints();
WarningInFunction
<< "Edge " << edgeI
<< "at " << pts[e[0]] << pts[e[1]]
<< " is a coupled edge and inbetween two different zones "
<< minZoneID[edgeI] << " and " << maxZoneID[edgeI] << endl
<< " This is currently not supported." << endl;
zoneZonePatch[minZoneID[edgeI]*nZones+maxZoneID[edgeI]]++;
}
}
else
{
// One or more than two edge-faces.
// Check whether we are on a mesh edge with external patches. If
// so choose any uncovered one. If none found put face in
// undetermined zone 'side' patch
label facei = findUncoveredPatchFace
(
mesh,
labelUIndList(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (facei == -1)
{
zoneSidePatch[minZoneID[edgeI]]++;
}
}
}
// Synchronise decision. Actual numbers are not important, just make
// sure that they're > 0 on all processors.
Pstream::listReduce(zoneSidePatch, sumOp());
Pstream::listReduce(zoneZonePatch, sumOp());
}
void addCouplingPatches
(
const fvMesh& mesh,
const word& regionName,
const word& shellRegionName,
const wordList& zoneNames,
const wordList& zoneShadowNames,
const boolList& isInternal,
const labelList& zoneIDs,
DynamicList& newPatches,
labelList& interRegionTopPatch,
labelList& interRegionBottomPatch
)
{
Pout<< "Adding coupling patches:" << nl << nl
<< "patchID\tpatch\ttype" << nl
<< "-------\t-----\t----"
<< endl;
interRegionTopPatch.setSize(zoneNames.size(), -1);
interRegionBottomPatch.setSize(zoneNames.size(), -1);
label nOldPatches = newPatches.size();
forAll(zoneNames, zoneI)
{
word interName
(
regionName
+"_to_"
+shellRegionName
+'_'
+zoneNames[zoneI]
);
if (isInternal[zoneI])
{
interRegionTopPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
interName + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
interName + "_bottom",
newPatches
);
Pout<< interRegionBottomPatch[zoneI]
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
<< nl;
}
else if (zoneShadowNames.size() == 0)
{
interRegionTopPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
zoneNames[zoneI] + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
interName,
newPatches
);
Pout<< interRegionBottomPatch[zoneI]
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
<< nl;
}
else //patch using shadow face zones.
{
interRegionTopPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
zoneShadowNames[zoneI] + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
interName,
newPatches
);
Pout<< interRegionBottomPatch[zoneI]
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionBottomPatch[zoneI]]->type()
<< nl;
}
}
Pout<< "Added " << newPatches.size()-nOldPatches
<< " inter-region patches." << nl
<< endl;
}
// Sets sidePatch[edgeI] to interprocessor or cyclic patch. Adds any
// coupled patches if necessary.
void addCoupledPatches
(
const fvMesh& mesh,
const primitiveFacePatch& extrudePatch,
const labelList& extrudeMeshFaces,
const labelList& extrudeMeshEdges,
const mapDistribute& extrudeEdgeFacesMap,
const labelListList& extrudeEdgeGlobalFaces,
labelList& sidePatchID,
DynamicList& newPatches
)
{
// Calculate opposite processor for coupled edges (only if shared by
// two procs). Note: could have saved original globalEdgeFaces structure.
// Get procID in extrudeEdgeGlobalFaces order
labelList procID(extrudeEdgeGlobalFaces.size(), Pstream::myProcNo());
extrudeEdgeFacesMap.distribute(procID);
labelList minProcID(extrudeEdgeGlobalFaces.size(), labelMax);
labelList maxProcID(extrudeEdgeGlobalFaces.size(), labelMin);
forAll(extrudeEdgeGlobalFaces, edgeI)
{
const labelList& eFaces = extrudeEdgeGlobalFaces[edgeI];
if (eFaces.size())
{
forAll(eFaces, i)
{
label proci = procID[eFaces[i]];
minProcID[edgeI] = Foam::min(minProcID[edgeI], proci);
maxProcID[edgeI] = Foam::max(maxProcID[edgeI], proci);
}
}
}
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
minProcID,
minEqOp(),
labelMax // null value
);
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
maxProcID,
maxEqOp(),
labelMin // null value
);
Pout<< "Adding processor or cyclic patches:" << nl << nl
<< "patchID\tpatch" << nl
<< "-------\t-----"
<< endl;
label nOldPatches = newPatches.size();
sidePatchID.setSize(extrudePatch.edgeFaces().size(), -1);
forAll(extrudePatch.edgeFaces(), edgeI)
{
const labelList& eFaces = extrudePatch.edgeFaces()[edgeI];
if
(
eFaces.size() == 1
&& extrudeEdgeGlobalFaces[edgeI].size() == 2
)
{
// coupled boundary edge. Find matching patch.
label nbrProci = minProcID[edgeI];
if (nbrProci == Pstream::myProcNo())
{
nbrProci = maxProcID[edgeI];
}
if (nbrProci == Pstream::myProcNo())
{
// Cyclic patch since both procs the same. This cyclic should
// already exist in newPatches so no adding necessary.
label facei = findUncoveredCyclicPatchFace
(
mesh,
labelUIndList(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (facei != -1)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
label newPatchi = findPatchID
(
newPatches,
patches[patches.whichPatch(facei)].name()
);
sidePatchID[edgeI] = newPatchi;
}
else
{
FatalErrorInFunction
<< "Unable to determine coupled patch addressing"
<< abort(FatalError);
}
}
else
{
// Processor patch
word name
(
processorPolyPatch::newName(Pstream::myProcNo(), nbrProci)
);
sidePatchID[edgeI] = findPatchID(newPatches, name);
if (sidePatchID[edgeI] == -1)
{
dictionary patchDict;
patchDict.add("myProcNo", Pstream::myProcNo());
patchDict.add("neighbProcNo", nbrProci);
sidePatchID[edgeI] = addPatch
(
mesh.boundaryMesh(),
name,
patchDict,
newPatches
);
Pout<< sidePatchID[edgeI] << '\t' << name
<< nl;
}
}
}
}
Pout<< "Added " << newPatches.size()-nOldPatches
<< " coupled patches." << nl
<< endl;
}
void addZoneSidePatches
(
const fvMesh& mesh,
const wordList& zoneNames,
const word& oneDPolyPatchType,
DynamicList& newPatches,
labelList& zoneSidePatch
)
{
Pout<< "Adding patches for sides on zones:" << nl << nl
<< "patchID\tpatch" << nl
<< "-------\t-----"
<< endl;
label nOldPatches = newPatches.size();
forAll(zoneSidePatch, zoneI)
{
if (!oneDPolyPatchType.empty())
{
// Reuse single empty patch.
word patchName;
if (oneDPolyPatchType == "empty")
{
patchName = "oneDEmptyPatch";
zoneSidePatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
patchName,
newPatches
);
}
else if (oneDPolyPatchType == "wedge")
{
patchName = "oneDWedgePatch";
zoneSidePatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
patchName,
newPatches
);
}
else
{
FatalErrorInFunction
<< "Type " << oneDPolyPatchType << " does not exist "
<< exit(FatalError);
}
Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
}
else if (zoneSidePatch[zoneI] > 0)
{
word patchName = zoneNames[zoneI] + "_" + "side";
zoneSidePatch[zoneI] = addPatch
(
mesh.boundaryMesh(),
patchName,
newPatches
);
Pout<< zoneSidePatch[zoneI] << '\t' << patchName << nl;
}
}
Pout<< "Added " << newPatches.size()-nOldPatches << " zone-side patches."
<< nl << endl;
}
void addInterZonePatches
(
const fvMesh& mesh,
const wordList& zoneNames,
const bool oneD,
labelList& zoneZonePatch_min,
labelList& zoneZonePatch_max,
DynamicList& newPatches
)
{
Pout<< "Adding inter-zone patches:" << nl << nl
<< "patchID\tpatch" << nl
<< "-------\t-----"
<< endl;
dictionary transformDict;
transformDict.add
(
"transform",
cyclicPolyPatch::transformTypeNames[cyclicPolyPatch::NOORDERING]
);
label nOldPatches = newPatches.size();
if (!oneD)
{
forAll(zoneZonePatch_min, minZone)
{
for (label maxZone = minZone; maxZone < zoneNames.size(); maxZone++)
{
label index = minZone*zoneNames.size()+maxZone;
if (zoneZonePatch_min[index] > 0)
{
word minToMax =
zoneNames[minZone]
+ "_to_"
+ zoneNames[maxZone];
word maxToMin =
zoneNames[maxZone]
+ "_to_"
+ zoneNames[minZone];
{
transformDict.set("neighbourPatch", maxToMin);
zoneZonePatch_min[index] =
addPatch
(
mesh.boundaryMesh(),
minToMax,
transformDict,
newPatches
);
Pout<< zoneZonePatch_min[index] << '\t' << minToMax
<< nl;
}
{
transformDict.set("neighbourPatch", minToMax);
zoneZonePatch_max[index] =
addPatch
(
mesh.boundaryMesh(),
maxToMin,
transformDict,
newPatches
);
Pout<< zoneZonePatch_max[index] << '\t' << maxToMin
<< nl;
}
}
}
}
}
Pout<< "Added " << newPatches.size()-nOldPatches << " inter-zone patches."
<< nl << endl;
}
tmp calcOffset
(
const primitiveFacePatch& extrudePatch,
const createShellMesh& extruder,
const polyPatch& pp
)
{
vectorField::subField fc = pp.faceCentres();
auto toffsets = tmp::New(fc.size());
auto& offsets = toffsets.ref();
forAll(fc, i)
{
label meshFacei = pp.start()+i;
label patchFacei = mag(extruder.faceToFaceMap()[meshFacei])-1;
point patchFc = extrudePatch[patchFacei].centre
(
extrudePatch.points()
);
offsets[i] = patchFc - fc[i];
}
return toffsets;
}
void setCouplingInfo
(
fvMesh& mesh,
const labelList& zoneToPatch,
const word& sampleRegion,
const mappedPatchBase::sampleMode mode,
const List& offsets,
const List& bbs
)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
List newPatches
(
patches.size(),
static_cast(nullptr)
);
forAll(zoneToPatch, zoneI)
{
const label patchi = zoneToPatch[zoneI];
if (patchi != -1)
{
const polyPatch& pp = patches[patchi];
if (isA(pp))
{
const vector avgOffset = gAverage(offsets[zoneI]);
const scalar mergeSqrDist =
gMax(magSqr(offsets[zoneI]-avgOffset));
// Create with uniform offset initially
auto mappedPtr = autoPtr::New
(
pp.name(),
pp.size(),
pp.start(),
patchi,
sampleRegion, // sampleRegion
mode, // sampleMode
pp.name(), // samplePatch
avgOffset, // uniform offset
patches
);
Info<< "Adding on " << mesh.name() << " coupling patch "
<< pp.name() << " with ";
// Verify uniformity of offset
// (same check as blockMesh geom merge)
if (mergeSqrDist < magSqr(10*SMALL*bbs[zoneI].span()))
{
Info<< "uniform offset " << avgOffset << endl;
}
else
{
Info<< "non-uniform offset" << endl;
(*mappedPtr).setOffset(offsets[zoneI]);
}
newPatches[patchi] = mappedPtr.release();
}
}
}
forAll(newPatches, patchi)
{
if (!newPatches[patchi])
{
newPatches[patchi] = patches[patchi].clone(patches).ptr();
//newPatches[patchi].index() = patchi;
}
}
#ifdef FULLDEBUG
Pout<< "*** setCouplingInfo addFvPAtches:" << nl;
printPatches(Pout, newPatches);
Pout<< "*** setCouplingInfo end of addFvPAtches:" << endl;
#endif
mesh.removeFvBoundary();
mesh.addFvPatches(newPatches, true);
}
// Extrude and write geometric properties
void extrudeGeometricProperties
(
const polyMesh& mesh,
const primitiveFacePatch& extrudePatch,
const createShellMesh& extruder,
const polyMesh& regionMesh,
const extrudeModel& model
)
{
const pointIOField patchFaceCentres
(
IOobject
(
"patchFaceCentres",
mesh.pointsInstance(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ
)
);
const pointIOField patchEdgeCentres
(
IOobject
(
"patchEdgeCentres",
mesh.pointsInstance(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ
)
);
//forAll(extrudePatch.edges(), edgeI)
//{
// const edge& e = extrudePatch.edges()[edgeI];
// Pout<< "Edge:" << e.centre(extrudePatch.localPoints()) << nl
// << "read:" << patchEdgeCentres[edgeI]
// << endl;
//}
// Determine edge normals on original patch
labelList patchEdges;
labelList coupledEdges;
bitSet sameEdgeOrientation;
PatchTools::matchEdges
(
extrudePatch,
mesh.globalData().coupledPatch(),
patchEdges,
coupledEdges,
sameEdgeOrientation
);
pointField patchEdgeNormals
(
PatchTools::edgeNormals
(
mesh,
extrudePatch,
patchEdges,
coupledEdges
)
);
pointIOField faceCentres
(
IOobject
(
"faceCentres",
regionMesh.pointsInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
regionMesh.nFaces()
);
// Work out layers. Guaranteed in columns so no fancy parallel bits.
forAll(extruder.faceToFaceMap(), facei)
{
if (extruder.faceToFaceMap()[facei] != 0)
{
// 'horizontal' face
label patchFacei = mag(extruder.faceToFaceMap()[facei])-1;
label celli = regionMesh.faceOwner()[facei];
if (regionMesh.isInternalFace(facei))
{
celli = Foam::max(celli, regionMesh.faceNeighbour()[facei]);
}
// Calculate layer from cell numbering (see createShellMesh)
label layerI = (celli % model.nLayers());
if
(
!regionMesh.isInternalFace(facei)
&& extruder.faceToFaceMap()[facei] > 0
)
{
// Top face
layerI++;
}
// Recalculate based on extrusion model
faceCentres[facei] = model
(
patchFaceCentres[patchFacei],
extrudePatch.faceNormals()[patchFacei],
layerI
);
}
else
{
// 'vertical face
label patchEdgeI = extruder.faceToEdgeMap()[facei];
label layerI =
(
regionMesh.faceOwner()[facei]
% model.nLayers()
);
// Extrude patch edge centre to this layer
point pt0 = model
(
patchEdgeCentres[patchEdgeI],
patchEdgeNormals[patchEdgeI],
layerI
);
// Extrude patch edge centre to next layer
point pt1 = model
(
patchEdgeCentres[patchEdgeI],
patchEdgeNormals[patchEdgeI],
layerI+1
);
// Interpolate
faceCentres[facei] = 0.5*(pt0+pt1);
}
}
pointIOField cellCentres
(
IOobject
(
"cellCentres",
regionMesh.pointsInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
regionMesh.nCells()
);
forAll(extruder.cellToFaceMap(), celli)
{
label patchFacei = extruder.cellToFaceMap()[celli];
// Calculate layer from cell numbering (see createShellMesh)
label layerI = (celli % model.nLayers());
// Recalculate based on extrusion model
point pt0 = model
(
patchFaceCentres[patchFacei],
extrudePatch.faceNormals()[patchFacei],
layerI
);
point pt1 = model
(
patchFaceCentres[patchFacei],
extrudePatch.faceNormals()[patchFacei],
layerI+1
);
// Interpolate
cellCentres[celli] = 0.5*(pt0+pt1);
}
// Bit of checking
if (false)
{
OBJstream faceStr(regionMesh.time().path()/"faceCentres.obj");
OBJstream cellStr(regionMesh.time().path()/"cellCentres.obj");
forAll(faceCentres, facei)
{
Pout<< "Model :" << faceCentres[facei] << endl
<< "regionMesh:" << regionMesh.faceCentres()[facei] << endl;
faceStr.writeLine
(
faceCentres[facei],
regionMesh.faceCentres()[facei]
);
}
forAll(cellCentres, celli)
{
Pout<< "Model :" << cellCentres[celli] << endl
<< "regionMesh:" << regionMesh.cellCentres()[celli] << endl;
cellStr.writeLine
(
cellCentres[celli],
regionMesh.cellCentres()[celli]
);
}
}
Info<< "Writing geometric properties for mesh " << regionMesh.name()
<< " to " << regionMesh.pointsInstance() << nl
<< endl;
bool ok = faceCentres.write() && cellCentres.write();
if (!ok)
{
FatalErrorInFunction
<< "Failed writing " << faceCentres.objectPath()
<< " and " << cellCentres.objectPath()
<< exit(FatalError);
}
}
int main(int argc, char *argv[])
{
argList::addNote
(
"Create region mesh by extruding a faceZone or faceSet"
);
#include "addRegionOption.H"
#include "addOverwriteOption.H"
argList::addOption
(
"dict", "file", "Alternative extrudeToRegionMeshDict"
);
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedMesh.H"
if (mesh.boundaryMesh().checkParallelSync(true))
{
List allNames(Pstream::nProcs());
allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names();
Pstream::allGatherList(allNames);
FatalErrorInFunction
<< "Patches are not synchronised on all processors."
<< " Per processor patches " << allNames
<< exit(FatalError);
}
const word oldInstance = mesh.pointsInstance();
const bool overwrite = args.found("overwrite");
const word dictName("extrudeToRegionMeshDict");
#include "setSystemMeshDictionaryIO.H"
IOdictionary dict(dictIO);
// Point generator
autoPtr model(extrudeModel::New(dict));
// Region
const word shellRegionName(dict.get("region"));
// Faces to extrude - either faceZones or faceSets (boundary faces only)
wordList zoneNames;
wordList zoneShadowNames;
const bool hasZones = dict.found("faceZones");
if (hasZones)
{
dict.readEntry("faceZones", zoneNames);
dict.readIfPresent("faceZonesShadow", zoneShadowNames);
// Check
if (dict.found("faceSets"))
{
FatalIOErrorInFunction(dict)
<< "Please supply faces to extrude either through 'faceZones'"
<< " or 'faceSets' entry. Found both."
<< exit(FatalIOError);
}
}
else
{
dict.readEntry("faceSets", zoneNames);
dict.readIfPresent("faceSetsShadow", zoneShadowNames);
}
mappedPatchBase::sampleMode sampleMode =
mappedPatchBase::sampleModeNames_.get("sampleMode", dict);
const bool oneD(dict.get("oneD"));
bool oneDNonManifoldEdges(false);
word oneDPatchType(emptyPolyPatch::typeName);
if (oneD)
{
oneDNonManifoldEdges = dict.getOrDefault("nonManifold", false);
oneDPatchType = dict.get("oneDPolyPatchType");
}
const bool adaptMesh(dict.get("adaptMesh"));
const bool addSidePatches(dict.getOrDefault("addSidePatches", true));
if (hasZones)
{
Info<< "Extruding zones " << zoneNames
<< " on mesh " << regionName
<< " into shell mesh " << shellRegionName
<< endl;
}
else
{
Info<< "Extruding faceSets " << zoneNames
<< " on mesh " << regionName
<< " into shell mesh " << shellRegionName
<< endl;
}
if (shellRegionName == regionName)
{
FatalIOErrorInFunction(dict)
<< "Cannot extrude into same region as mesh." << endl
<< "Mesh region : " << regionName << endl
<< "Shell region : " << shellRegionName
<< exit(FatalIOError);
}
if (oneD)
{
if (oneDNonManifoldEdges)
{
Info<< "Extruding as 1D columns with sides in patch type "
<< oneDPatchType
<< " and connected points (except on non-manifold areas)."
<< endl;
}
else
{
Info<< "Extruding as 1D columns with sides in patch type "
<< oneDPatchType
<< " and duplicated points (overlapping volumes)."
<< endl;
}
}
if (addSidePatches && zoneNames.size() > 1)
{
Info<< "Extruding edges on more than one faceZone into boundary faces"
<< endl;
}
else
{
Info<< "Extruding internal edges into internal faces" << endl;
}
//// Read objects in time directory
//IOobjectList objects(mesh, runTime.timeName());
//
//// Read vol fields.
//
//PtrList vsFlds;
//ReadFields(mesh, objects, vsFlds);
//
//PtrList vvFlds;
//ReadFields(mesh, objects, vvFlds);
//
//PtrList vstFlds;
//ReadFields(mesh, objects, vstFlds);
//
//PtrList vsymtFlds;
//ReadFields(mesh, objects, vsymtFlds);
//
//PtrList vtFlds;
//ReadFields(mesh, objects, vtFlds);
//
//// Read surface fields.
//
//PtrList ssFlds;
//ReadFields(mesh, objects, ssFlds);
//
//PtrList svFlds;
//ReadFields(mesh, objects, svFlds);
//
//PtrList sstFlds;
//ReadFields(mesh, objects, sstFlds);
//
//PtrList ssymtFlds;
//ReadFields(mesh, objects, ssymtFlds);
//
//PtrList stFlds;
//ReadFields(mesh, objects, stFlds);
//
//// Read point fields.
//
//PtrList psFlds;
//ReadFields(pointMesh::New(mesh), objects, psFlds);
//
//PtrList pvFlds;
//ReadFields(pointMesh::New(mesh), objects, pvFlds);
// Create dummy fv* files
fvMeshTools::createDummyFvMeshFiles(mesh, shellRegionName, true);
word meshInstance;
if (!overwrite)
{
++runTime;
meshInstance = runTime.timeName();
}
else
{
meshInstance = oldInstance;
}
Info<< "Writing meshes to " << meshInstance << nl << endl;
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Extract faces to extrude
// ~~~~~~~~~~~~~~~~~~~~~~~~
// Note: zoneID are regions of extrusion. They are not mesh.faceZones
// indices.
// From extrude zone to mesh zone (or -1 if extruding faceSets)
labelList meshZoneID;
// Per extrude zone whether contains internal or external faces
boolList isInternal(zoneNames.size(), false);
labelList extrudeMeshFaces;
faceList zoneFaces;
labelList zoneID;
boolList zoneFlipMap;
// Shadow
labelList zoneShadowIDs; // from extrude shadow zone to mesh zone
labelList extrudeMeshShadowFaces;
boolList zoneShadowFlipMap;
labelList zoneShadowID;
if (hasZones)
{
const faceZoneMesh& faceZones = mesh.faceZones();
meshZoneID.setSize(zoneNames.size());
forAll(zoneNames, i)
{
meshZoneID[i] = faceZones.findZoneID(zoneNames[i]);
if (meshZoneID[i] == -1)
{
FatalIOErrorInFunction(dict)
<< "Cannot find zone " << zoneNames[i] << endl
<< "Valid zones are " << faceZones.names()
<< exit(FatalIOError);
}
}
// Collect per face information
label nExtrudeFaces = 0;
forAll(meshZoneID, i)
{
nExtrudeFaces += faceZones[meshZoneID[i]].size();
}
extrudeMeshFaces.setSize(nExtrudeFaces);
zoneFaces.setSize(nExtrudeFaces);
zoneID.setSize(nExtrudeFaces);
zoneFlipMap.setSize(nExtrudeFaces);
nExtrudeFaces = 0;
forAll(meshZoneID, i)
{
const faceZone& fz = faceZones[meshZoneID[i]];
const primitiveFacePatch& fzp = fz();
forAll(fz, j)
{
extrudeMeshFaces[nExtrudeFaces] = fz[j];
zoneFaces[nExtrudeFaces] = fzp[j];
zoneID[nExtrudeFaces] = i;
zoneFlipMap[nExtrudeFaces] = fz.flipMap()[j];
nExtrudeFaces++;
if (mesh.isInternalFace(fz[j]))
{
isInternal[i] = true;
}
}
}
// Shadow zone
// ~~~~~~~~~~~
if (zoneShadowNames.size())
{
zoneShadowIDs.setSize(zoneShadowNames.size());
forAll(zoneShadowNames, i)
{
zoneShadowIDs[i] = faceZones.findZoneID(zoneShadowNames[i]);
if (zoneShadowIDs[i] == -1)
{
FatalIOErrorInFunction(dict)
<< "Cannot find zone " << zoneShadowNames[i] << endl
<< "Valid zones are " << faceZones.names()
<< exit(FatalIOError);
}
}
label nShadowFaces = 0;
forAll(zoneShadowIDs, i)
{
nShadowFaces += faceZones[zoneShadowIDs[i]].size();
}
extrudeMeshShadowFaces.setSize(nShadowFaces);
zoneShadowFlipMap.setSize(nShadowFaces);
zoneShadowID.setSize(nShadowFaces);
nShadowFaces = 0;
forAll(zoneShadowIDs, i)
{
const faceZone& fz = faceZones[zoneShadowIDs[i]];
forAll(fz, j)
{
extrudeMeshShadowFaces[nShadowFaces] = fz[j];
zoneShadowFlipMap[nShadowFaces] = fz.flipMap()[j];
zoneShadowID[nShadowFaces] = i;
nShadowFaces++;
}
}
}
}
else
{
meshZoneID.setSize(zoneNames.size(), -1);
// Load faceSets
PtrList zones(zoneNames.size());
forAll(zoneNames, i)
{
Info<< "Loading faceSet " << zoneNames[i] << endl;
zones.set(i, new faceSet(mesh, zoneNames[i]));
}
// Collect per face information
label nExtrudeFaces = 0;
forAll(zones, i)
{
nExtrudeFaces += zones[i].size();
}
extrudeMeshFaces.setSize(nExtrudeFaces);
zoneFaces.setSize(nExtrudeFaces);
zoneID.setSize(nExtrudeFaces);
zoneFlipMap.setSize(nExtrudeFaces);
nExtrudeFaces = 0;
forAll(zones, i)
{
const faceSet& fz = zones[i];
for (const label facei : fz)
{
if (mesh.isInternalFace(facei))
{
FatalIOErrorInFunction(dict)
<< "faceSet " << fz.name()
<< "contains internal faces."
<< " This is not permitted."
<< exit(FatalIOError);
}
extrudeMeshFaces[nExtrudeFaces] = facei;
zoneFaces[nExtrudeFaces] = mesh.faces()[facei];
zoneID[nExtrudeFaces] = i;
zoneFlipMap[nExtrudeFaces] = false;
nExtrudeFaces++;
if (mesh.isInternalFace(facei))
{
isInternal[i] = true;
}
}
}
// Shadow zone
// ~~~~~~~~~~~
PtrList shadowZones(zoneShadowNames.size());
if (zoneShadowNames.size())
{
zoneShadowIDs.setSize(zoneShadowNames.size(), -1);
forAll(zoneShadowNames, i)
{
shadowZones.set(i, new faceSet(mesh, zoneShadowNames[i]));
}
label nShadowFaces = 0;
for (const faceSet& fz : shadowZones)
{
nShadowFaces += fz.size();
}
if (nExtrudeFaces != nShadowFaces)
{
FatalIOErrorInFunction(dict)
<< "Extruded faces " << nExtrudeFaces << endl
<< "is different from shadow faces. " << nShadowFaces
<< "This is not permitted " << endl
<< exit(FatalIOError);
}
extrudeMeshShadowFaces.setSize(nShadowFaces);
zoneShadowFlipMap.setSize(nShadowFaces);
zoneShadowID.setSize(nShadowFaces);
nShadowFaces = 0;
forAll(shadowZones, i)
{
const faceSet& fz = shadowZones[i];
for (const label facei : fz)
{
if (mesh.isInternalFace(facei))
{
FatalIOErrorInFunction(dict)
<< "faceSet " << fz.name()
<< "contains internal faces."
<< " This is not permitted."
<< exit(FatalIOError);
}
extrudeMeshShadowFaces[nShadowFaces] = facei;
zoneShadowFlipMap[nShadowFaces] = false;
zoneShadowID[nShadowFaces] = i;
nShadowFaces++;
}
}
}
}
const primitiveFacePatch extrudePatch(std::move(zoneFaces), mesh.points());
Pstream::listCombineReduce(isInternal, orEqOp());
// Check zone either all internal or all external faces
checkZoneInside(mesh, zoneNames, zoneID, extrudeMeshFaces, isInternal);
const pointField& extrudePoints = extrudePatch.localPoints();
const faceList& extrudeFaces = extrudePatch.localFaces();
const labelListList& edgeFaces = extrudePatch.edgeFaces();
Info<< "extrudePatch :"
<< " faces:" << extrudePatch.size()
<< " points:" << extrudePatch.nPoints()
<< " edges:" << extrudePatch.nEdges()
<< nl
<< endl;
// Determine per-extrude-edge info
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Corresponding mesh edges
const labelList extrudeMeshEdges
(
extrudePatch.meshEdges
(
mesh.edges(),
mesh.pointEdges()
)
);
const globalIndex globalExtrudeFaces(extrudePatch.size());
// Global pp faces per pp edge.
labelListList extrudeEdgeGlobalFaces
(
globalEdgeFaces
(
mesh,
globalExtrudeFaces,
extrudePatch,
extrudeMeshEdges
)
);
List> compactMap;
const mapDistribute extrudeEdgeFacesMap
(
globalExtrudeFaces,
extrudeEdgeGlobalFaces,
compactMap
);
// Determine min and max zone per edge
labelList edgeMinZoneID;
labelList edgeMaxZoneID;
calcEdgeMinMaxZone
(
mesh,
extrudePatch,
extrudeMeshEdges,
zoneID,
extrudeEdgeFacesMap,
extrudeEdgeGlobalFaces,
edgeMinZoneID,
edgeMaxZoneID
);
if (!addSidePatches)
{
// Disabling inter-zone multiple boundary faces by setting 'left' and
// 'right' zone to be the same
edgeMaxZoneID = edgeMinZoneID;
}
DynamicList regionPatches(patches.size());
// Copy all non-local patches since these are used on boundary edges of
// the extrusion
forAll(patches, patchi)
{
if (!isA(patches[patchi]))
{
label newPatchi = regionPatches.size();
regionPatches.append
(
patches[patchi].clone
(
patches,
newPatchi,
0, // size
0 // start
).ptr()
);
}
}
// Add interface patches
// ~~~~~~~~~~~~~~~~~~~~~
// From zone to interface patch (region side)
labelList interRegionTopPatch;
labelList interRegionBottomPatch;
addCouplingPatches
(
mesh,
regionName,
shellRegionName,
zoneNames,
zoneShadowNames,
isInternal,
meshZoneID,
regionPatches,
interRegionTopPatch,
interRegionBottomPatch
);
// From zone to interface patch (mesh side)
labelList interMeshTopPatch;
labelList interMeshBottomPatch;
if (adaptMesh)
{
// Add coupling patches to mesh
// 1. Clone existing global patches
DynamicList newPatches(patches.size());
forAll(patches, patchi)
{
if (!isA(patches[patchi]))
{
autoPtr clonedPatch(patches[patchi].clone(patches));
clonedPatch->index() = newPatches.size();
newPatches.append(clonedPatch.ptr());
}
}
// 2. Add new patches
addCouplingPatches
(
mesh,
regionName,
shellRegionName,
zoneNames,
zoneShadowNames,
isInternal,
meshZoneID,
newPatches,
interMeshTopPatch,
interMeshBottomPatch
);
// 3. Clone processor patches
forAll(patches, patchi)
{
if (isA(patches[patchi]))
{
autoPtr clonedPatch(patches[patchi].clone(patches));
clonedPatch->index() = newPatches.size();
newPatches.append(clonedPatch.ptr());
}
}
#ifdef FULLDEBUG
Pout<< "*** adaptMesh : addFvPAtches:" << nl;
printPatches(Pout, newPatches);
Pout<< "*** end of adaptMesh : addFvPAtches:" << endl;
#endif
// Add to mesh
mesh.clearOut();
mesh.removeFvBoundary();
mesh.addFvPatches(newPatches, true);
//!Note: from this point on mesh patches differs from regionPatches
}
// Patch per extruded face
labelList extrudeTopPatchID(extrudePatch.size());
labelList extrudeBottomPatchID(extrudePatch.size());
forAll(zoneID, facei)
{
extrudeTopPatchID[facei] = interRegionTopPatch[zoneID[facei]];
extrudeBottomPatchID[facei] = interRegionBottomPatch[zoneID[facei]];
}
// Count how many patches on special edges of extrudePatch are necessary
// - zoneXXX_sides
// - zoneXXX_zoneYYY
labelList zoneSidePatch(zoneNames.size(), Zero);
// Patch to use for minZone
labelList zoneZonePatch_min(zoneNames.size()*zoneNames.size(), Zero);
// Patch to use for maxZone
labelList zoneZonePatch_max(zoneNames.size()*zoneNames.size(), Zero);
countExtrudePatches
(
mesh,
zoneNames.size(),
extrudePatch, // patch
extrudeMeshFaces, // mesh face per patch face
extrudeMeshEdges, // mesh edge per patch edge
extrudeEdgeGlobalFaces, // global indexing per patch edge
edgeMinZoneID, // minZone per patch edge
edgeMaxZoneID, // maxZone per patch edge
zoneSidePatch, // per zone-side num edges that extrude into it
zoneZonePatch_min // per zone-zone num edges that extrude into it
);
// Now we'll have:
// zoneSidePatch[zoneA] : number of faces needed on the side of zoneA
// zoneZonePatch_min[zoneA,zoneB] : number of faces needed inbetween A,B
// Add the zone-side patches.
addZoneSidePatches
(
mesh,
zoneNames,
(oneD ? oneDPatchType : word::null),
regionPatches,
zoneSidePatch
);
// Add the patches inbetween zones
addInterZonePatches
(
mesh,
zoneNames,
oneD,
zoneZonePatch_min,
zoneZonePatch_max,
regionPatches
);
// Sets sidePatchID[edgeI] to interprocessor patch. Adds any
// interprocessor or cyclic patches if necessary.
labelList sidePatchID;
addCoupledPatches
(
mesh,
extrudePatch,
extrudeMeshFaces,
extrudeMeshEdges,
extrudeEdgeFacesMap,
extrudeEdgeGlobalFaces,
sidePatchID,
regionPatches
);
// // Add all the newPatches to the mesh and fields
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// {
// forAll(newPatches, patchi)
// {
// Pout<< "Adding patch " << patchi
// << " name:" << newPatches[patchi]->name()
// << endl;
// }
// //label nOldPatches = mesh.boundary().size();
// mesh.clearOut();
// mesh.removeFvBoundary();
// mesh.addFvPatches(newPatches, true);
// //// Add calculated fvPatchFields for the added patches
// //for
// //(
// // label patchi = nOldPatches;
// // patchi < mesh.boundary().size();
// // patchi++
// //)
// //{
// // Pout<< "ADDing calculated to patch " << patchi
// // << endl;
// // addCalculatedPatchFields(mesh);
// //}
// //Pout<< "** Added " << mesh.boundary().size()-nOldPatches
// // << " patches." << endl;
// }
// Set patches to use for edges to be extruded into boundary faces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// In order of edgeFaces: per edge, per originating face the
// patch to use for the side face (from the extruded edge).
// If empty size create an internal face.
labelListList extrudeEdgePatches(extrudePatch.nEdges());
// Is edge a non-manifold edge
bitSet nonManifoldEdge(extrudePatch.nEdges());
// Note: logic has to be same as in countExtrudePatches.
forAll(edgeFaces, edgeI)
{
const labelList& eFaces = edgeFaces[edgeI];
labelList& ePatches = extrudeEdgePatches[edgeI];
if (oneD)
{
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
{
ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
}
if (oneDNonManifoldEdges)
{
//- Set nonManifoldEdge[edgeI] for non-manifold edges only
// The other option is to have non-manifold edges everywhere
// and generate space overlapping columns of cells.
if (eFaces.size() != 2)
{
nonManifoldEdge.set(edgeI);
}
}
else
{
nonManifoldEdge.set(edgeI);
}
}
else if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
if (addSidePatches && (zone0 != zone1)) // || (cos(angle) > blabla))
{
label minZone = Foam::min(zone0,zone1);
label maxZone = Foam::max(zone0,zone1);
label index = minZone*zoneNames.size()+maxZone;
ePatches.setSize(eFaces.size());
if (zone0 == minZone)
{
ePatches[0] = zoneZonePatch_min[index];
ePatches[1] = zoneZonePatch_max[index];
}
else
{
ePatches[0] = zoneZonePatch_max[index];
ePatches[1] = zoneZonePatch_min[index];
}
nonManifoldEdge.set(edgeI);
}
}
else if (sidePatchID[edgeI] != -1)
{
// Coupled extrusion
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
{
ePatches[i] = sidePatchID[edgeI];
}
}
else
{
label facei = findUncoveredPatchFace
(
mesh,
labelUIndList(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (facei != -1)
{
label newPatchi = findPatchID
(
regionPatches,
patches[patches.whichPatch(facei)].name()
);
ePatches.setSize(eFaces.size(), newPatchi);
}
else
{
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
{
ePatches[i] = zoneSidePatch[zoneID[eFaces[i]]];
}
}
nonManifoldEdge.set(edgeI);
}
}
// Assign point regions
// ~~~~~~~~~~~~~~~~~~~~
// Per face, per point the region number.
faceList pointGlobalRegions;
faceList pointLocalRegions;
labelList localToGlobalRegion;
createShellMesh::calcPointRegions
(
mesh.globalData(),
extrudePatch,
nonManifoldEdge,
false, // keep cyclic separated regions apart
pointGlobalRegions,
pointLocalRegions,
localToGlobalRegion
);
// Per local region an originating point
labelList localRegionPoints(localToGlobalRegion.size());
forAll(pointLocalRegions, facei)
{
const face& f = extrudePatch.localFaces()[facei];
const face& pRegions = pointLocalRegions[facei];
forAll(pRegions, fp)
{
localRegionPoints[pRegions[fp]] = f[fp];
}
}
// Calculate region normals by reducing local region normals
pointField localRegionNormals(localToGlobalRegion.size());
{
pointField localSum(localToGlobalRegion.size(), Zero);
forAll(pointLocalRegions, facei)
{
const face& pRegions = pointLocalRegions[facei];
forAll(pRegions, fp)
{
label localRegionI = pRegions[fp];
localSum[localRegionI] += extrudePatch.faceNormals()[facei];
}
}
Map globalSum(2*localToGlobalRegion.size());
forAll(localSum, localRegionI)
{
label globalRegionI = localToGlobalRegion[localRegionI];
globalSum.insert(globalRegionI, localSum[localRegionI]);
}
// Reduce
Pstream::mapCombineReduce(globalSum, plusEqOp());
forAll(localToGlobalRegion, localRegionI)
{
label globalRegionI = localToGlobalRegion[localRegionI];
localRegionNormals[localRegionI] = globalSum[globalRegionI];
}
localRegionNormals /= mag(localRegionNormals);
}
// For debugging: dump hedgehog plot of normals
if (false)
{
OFstream str(runTime.path()/"localRegionNormals.obj");
label vertI = 0;
scalar thickness = model().sumThickness(1);
forAll(pointLocalRegions, facei)
{
const face& f = extrudeFaces[facei];
forAll(f, fp)
{
label region = pointLocalRegions[facei][fp];
const point& pt = extrudePoints[f[fp]];
meshTools::writeOBJ(str, pt);
vertI++;
meshTools::writeOBJ
(
str,
pt+thickness*localRegionNormals[region]
);
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
}
// Use model to create displacements of first layer
vectorField firstDisp(localRegionNormals.size());
forAll(firstDisp, regionI)
{
//const point& regionPt = regionCentres[regionI];
const point& regionPt = extrudePatch.points()
[
extrudePatch.meshPoints()
[
localRegionPoints[regionI]
]
];
const vector& n = localRegionNormals[regionI];
firstDisp[regionI] = model()(regionPt, n, 1) - regionPt;
}
// Create a new mesh
// ~~~~~~~~~~~~~~~~~
createShellMesh extruder
(
extrudePatch,
pointLocalRegions,
localRegionPoints
);
autoPtr shellMap;
fvMesh regionMesh
(
IOobject
(
shellRegionName,
meshInstance,
runTime,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
IOobject::NO_REGISTER
),
Foam::zero{},
false // syncPar
);
// Add the new patches
forAll(regionPatches, patchi)
{
polyPatch* ppPtr = regionPatches[patchi];
regionPatches[patchi] = ppPtr->clone(regionMesh.boundaryMesh()).ptr();
delete ppPtr;
}
#ifdef FULLDEBUG
Pout<< "*** regionPatches : regionPatches:" << nl;
printPatches(Pout, regionPatches);
Pout<< "*** end of regionPatches : regionPatches:" << endl;
#endif
regionMesh.clearOut();
regionMesh.removeFvBoundary();
regionMesh.addFvPatches(regionPatches, true);
{
polyTopoChange meshMod(regionPatches.size());
extruder.setRefinement
(
firstDisp, // first displacement
model().expansionRatio(),
model().nLayers(), // nLayers
extrudeTopPatchID,
extrudeBottomPatchID,
extrudeEdgePatches,
meshMod
);
// Enforce actual point positions according to extrudeModel (model)
// (extruder.setRefinement only does fixed expansionRatio)
// The regionPoints and nLayers are looped in the same way as in
// createShellMesh
DynamicList& newPoints = const_cast&>
(
meshMod.points()
);
label meshPointi = extrudePatch.localPoints().size();
forAll(localRegionPoints, regionI)
{
label pointi = localRegionPoints[regionI];
point pt = extrudePatch.localPoints()[pointi];
const vector& n = localRegionNormals[regionI];
for (label layerI = 1; layerI <= model().nLayers(); layerI++)
{
newPoints[meshPointi++] = model()(pt, n, layerI);
}
}
shellMap = meshMod.changeMesh
(
regionMesh, // mesh to change
false // inflate
);
}
// Necessary?
regionMesh.setInstance(meshInstance);
// Update numbering on extruder.
extruder.updateMesh(shellMap());
// Calculate offsets from shell mesh back to original mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
List topOffsets(zoneNames.size());
List topBbs(zoneNames.size());
List bottomOffsets(zoneNames.size());
List bottomBbs(zoneNames.size());
forAll(regionMesh.boundaryMesh(), patchi)
{
const polyPatch& pp = regionMesh.boundaryMesh()[patchi];
if (isA(pp))
{
if (interRegionTopPatch.found(patchi))
{
label zoneI = interRegionTopPatch.find(patchi);
topOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
topBbs[zoneI] = boundBox(pp.points(), pp.meshPoints(), true);
}
else if (interRegionBottomPatch.found(patchi))
{
label zoneI = interRegionBottomPatch.find(patchi);
bottomOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
bottomBbs[zoneI] = boundBox(pp.points(), pp.meshPoints(), true);
}
}
}
// Change top and bottom boundary conditions on regionMesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
// Correct top patches for offset
setCouplingInfo
(
regionMesh,
interRegionTopPatch,
regionName, // name of main mesh
sampleMode, // sampleMode
topOffsets,
topBbs
);
// Correct bottom patches for offset
setCouplingInfo
(
regionMesh,
interRegionBottomPatch,
regionName,
sampleMode, // sampleMode
bottomOffsets,
bottomBbs
);
// Remove any unused patches
deleteEmptyPatches(regionMesh);
}
// Change top and bottom boundary conditions on main mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (adaptMesh)
{
// Correct top patches for offset
setCouplingInfo
(
mesh,
interMeshTopPatch,
shellRegionName, // name of shell mesh
sampleMode, // sampleMode
-topOffsets,
topBbs
);
// Correct bottom patches for offset
setCouplingInfo
(
mesh,
interMeshBottomPatch,
shellRegionName,
sampleMode,
-bottomOffsets,
bottomBbs
);
}
// Write addressing from region mesh back to originating patch
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelIOList cellToPatchFaceAddressing
(
IOobject
(
"cellToPatchFaceAddressing",
regionMesh.facesInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
extruder.cellToFaceMap()
);
cellToPatchFaceAddressing.note() = "cell to patch face addressing";
labelIOList faceToPatchFaceAddressing
(
IOobject
(
"faceToPatchFaceAddressing",
regionMesh.facesInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
extruder.faceToFaceMap()
);
faceToPatchFaceAddressing.note() =
"front/back face + turning index to patch face addressing";
labelIOList faceToPatchEdgeAddressing
(
IOobject
(
"faceToPatchEdgeAddressing",
regionMesh.facesInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
extruder.faceToEdgeMap()
);
faceToPatchEdgeAddressing.note() =
"side face to patch edge addressing";
labelIOList pointToPatchPointAddressing
(
IOobject
(
"pointToPatchPointAddressing",
regionMesh.facesInstance(),
polyMesh::meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
IOobject::NO_REGISTER
),
extruder.pointToPointMap()
);
pointToPatchPointAddressing.note() =
"point to patch point addressing";
Info<< "Writing mesh " << regionMesh.name()
<< " to " << regionMesh.facesInstance() << nl
<< endl;
bool ok =
regionMesh.write()
&& cellToPatchFaceAddressing.write()
&& faceToPatchFaceAddressing.write()
&& faceToPatchEdgeAddressing.write()
&& pointToPatchPointAddressing.write();
if (!ok)
{
FatalErrorInFunction
<< "Failed writing mesh " << regionMesh.name()
<< " at location " << regionMesh.facesInstance()
<< exit(FatalError);
}
topoSet::removeFiles(regionMesh);
processorMeshes::removeFiles(regionMesh);
// See if we need to extrude coordinates as well
{
autoPtr patchFaceCentresPtr;
IOobject io
(
"patchFaceCentres",
mesh.pointsInstance(),
polyMesh::meshSubDir,
mesh,
IOobject::MUST_READ
);
if (io.typeHeaderOk(true))
{
// Read patchFaceCentres and patchEdgeCentres
Info<< "Reading patch face,edge centres : "
<< io.name() << " and patchEdgeCentres" << endl;
extrudeGeometricProperties
(
mesh,
extrudePatch,
extruder,
regionMesh,
model()
);
}
}
// Insert baffles into original mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr addBafflesMap;
if (adaptMesh)
{
polyTopoChange meshMod(mesh);
// Modify faces to be in bottom (= always coupled) patch
forAll(extrudeMeshFaces, zoneFacei)
{
label meshFacei = extrudeMeshFaces[zoneFacei];
label zoneI = zoneID[zoneFacei];
bool flip = zoneFlipMap[zoneFacei];
const face& f = mesh.faces()[meshFacei];
if (!flip)
{
meshMod.modifyFace
(
f, // modified face
meshFacei, // label of face being modified
mesh.faceOwner()[meshFacei],// owner
-1, // neighbour
false, // face flip
interMeshBottomPatch[zoneI],// patch for face
meshZoneID[zoneI], // zone for face
flip // face flip in zone
);
}
else if (mesh.isInternalFace(meshFacei))
{
meshMod.modifyFace
(
f.reverseFace(), // modified face
meshFacei, // label of modified face
mesh.faceNeighbour()[meshFacei],// owner
-1, // neighbour
true, // face flip
interMeshBottomPatch[zoneI], // patch for face
meshZoneID[zoneI], // zone for face
!flip // face flip in zone
);
}
}
if (zoneShadowNames.size() > 0) //if there is a top faceZone specified
{
forAll(extrudeMeshFaces, zoneFacei)
{
label meshFacei = extrudeMeshShadowFaces[zoneFacei];
label zoneI = zoneShadowID[zoneFacei];
bool flip = zoneShadowFlipMap[zoneFacei];
const face& f = mesh.faces()[meshFacei];
if (!flip)
{
meshMod.modifyFace
(
f, // modified face
meshFacei, // face being modified
mesh.faceOwner()[meshFacei],// owner
-1, // neighbour
false, // face flip
interMeshTopPatch[zoneI], // patch for face
meshZoneID[zoneI], // zone for face
flip // face flip in zone
);
}
else if (mesh.isInternalFace(meshFacei))
{
meshMod.modifyFace
(
f.reverseFace(), // modified face
meshFacei, // label modified face
mesh.faceNeighbour()[meshFacei],// owner
-1, // neighbour
true, // face flip
interMeshTopPatch[zoneI], // patch for face
meshZoneID[zoneI], // zone for face
!flip // face flip in zone
);
}
}
}
else
{
// Add faces (using same points) to be in top patch
forAll(extrudeMeshFaces, zoneFacei)
{
label meshFacei = extrudeMeshFaces[zoneFacei];
label zoneI = zoneID[zoneFacei];
bool flip = zoneFlipMap[zoneFacei];
const face& f = mesh.faces()[meshFacei];
if (!flip)
{
if (mesh.isInternalFace(meshFacei))
{
meshMod.addFace
(
f.reverseFace(), // modified face
mesh.faceNeighbour()[meshFacei],// owner
-1, // neighbour
-1, // master point
-1, // master edge
meshFacei, // master face
true, // flip flux
interMeshTopPatch[zoneI], // patch for face
-1, // zone for face
false //face flip in zone
);
}
}
else
{
meshMod.addFace
(
f, // face
mesh.faceOwner()[meshFacei], // owner
-1, // neighbour
-1, // master point
-1, // master edge
meshFacei, // master face
false, // flip flux
interMeshTopPatch[zoneI], // patch for face
-1, // zone for face
false // zone flip
);
}
}
}
// Change the mesh. Change points directly (no inflation).
addBafflesMap = meshMod.changeMesh(mesh, false);
// Update fields
mesh.updateMesh(addBafflesMap());
//XXXXXX
// Update maps! e.g. faceToPatchFaceAddressing
//XXXXXX
// Move mesh (since morphing might not do this)
if (addBafflesMap().hasMotionPoints())
{
mesh.movePoints(addBafflesMap().preMotionPoints());
}
mesh.setInstance(meshInstance);
// Remove any unused patches
deleteEmptyPatches(mesh);
Info<< "Writing mesh " << mesh.name()
<< " to " << mesh.facesInstance() << nl
<< endl;
if (!mesh.write())
{
FatalErrorInFunction
<< "Failed writing mesh " << mesh.name()
<< " at location " << mesh.facesInstance()
<< exit(FatalError);
}
topoSet::removeFiles(mesh);
processorMeshes::removeFiles(mesh);
}
Info << "End\n" << endl;
return 0;
}
// ************************************************************************* //