songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
6a3571b3e9
openfoam/applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C
mattijs 4d57d3c568 STYLE: extrudeToRegionMesh: comment
2012-03-28 18:10:13 +01:00

3036 lines
85 KiB
C
Raw Blame History

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ 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 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 <http://www.gnu.org/licenses/>.
Description
Extrude faceZones (internal or boundary faces) or faceSets (boundary faces
only) into a separate mesh (as a different region).
- used to e.g. extrude baffles (extrude internal faces) or create
liquid film regions.
- if extruding internal faces:
- create baffles in original mesh with mappedWall patches
- if extruding boundary faces:
- convert boundary faces to mappedWall patches
- extrude edges of faceZone as a \<zone\>_sidePatch
- extrude edges inbetween different faceZones as a
(nonuniformTransform)cyclic \<zoneA\>_\<zoneB\>
- 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"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
//template<class GeoField>
//void addPatchFields(const fvMesh& mesh, const word& patchFieldType)
//{
// HashTable<const GeoField*> flds
// (
// mesh.objectRegistry::lookupClass<GeoField>()
// );
//
// forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
// {
// const GeoField& fld = *iter();
//
// typename GeoField::GeometricBoundaryField& bfld =
// const_cast<typename GeoField::GeometricBoundaryField&>
// (
// fld.boundaryField()
// );
//
// label sz = bfld.size();
//
// for (label i = 0; i < sz; i++)
// {
// bfld.set
// (
// i,
// bfld.clone(GeoField::PatchFieldType::New
// (
// patchFieldType,
// fld.mesh().boundary()[sz],
// fld.dimensionedInternalField()
// )
// );
//
//
//
// Pout<< "fld:" << fld.name() << " had " << sz << " patches." << endl;
// Pout<< "fld before:" << fld << endl;
// Pout<< "adding on patch:" << fld.mesh().boundary()[sz].name() << endl;
//
// bfld.setSize(sz+1);
// bfld.set
// (
// sz,
// GeoField::PatchFieldType::New
// (
// patchFieldType,
// fld.mesh().boundary()[sz],
// fld.dimensionedInternalField()
// )
// );
//
// bfld[sz].operator=(pTraits<typename GeoField::value_type>::zero);
//
// Pout<< "fld:" << fld.name() << " now " << bfld.size() << " patches."
// << endl;
//
// const typename GeoField::PatchFieldType& pfld = bfld[sz];
// Pout<< "pfld:" << pfld << endl;
//
//
// Pout<< "fld value:" << fld << endl;
// }
//}
// Remove last patch field
template<class GeoField>
void trimPatchFields(fvMesh& mesh, const label nPatches)
{
HashTable<const GeoField*> flds
(
mesh.objectRegistry::lookupClass<GeoField>()
);
forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
{
const GeoField& fld = *iter();
const_cast<typename GeoField::GeometricBoundaryField&>
(
fld.boundaryField()
).setSize(nPatches);
}
}
// Reorder patch field
template<class GeoField>
void reorderPatchFields(fvMesh& mesh, const labelList& oldToNew)
{
HashTable<const GeoField*> flds
(
mesh.objectRegistry::lookupClass<GeoField>()
);
forAllConstIter(typename HashTable<const GeoField*>, flds, iter)
{
const GeoField& fld = *iter();
typename GeoField::GeometricBoundaryField& bfld =
const_cast<typename GeoField::GeometricBoundaryField&>
(
fld.boundaryField()
);
bfld.reorder(oldToNew);
}
}
//void addCalculatedPatchFields(const fvMesh& mesh)
//{
// addPatchFields<volScalarField>
// (
// mesh,
// calculatedFvPatchField<scalar>::typeName
// );
// addPatchFields<volVectorField>
// (
// mesh,
// calculatedFvPatchField<vector>::typeName
// );
// addPatchFields<volSphericalTensorField>
// (
// mesh,
// calculatedFvPatchField<sphericalTensor>::typeName
// );
// addPatchFields<volSymmTensorField>
// (
// mesh,
// calculatedFvPatchField<symmTensor>::typeName
// );
// addPatchFields<volTensorField>
// (
// mesh,
// calculatedFvPatchField<tensor>::typeName
// );
//
// // Surface fields
//
// addPatchFields<surfaceScalarField>
// (
// mesh,
// calculatedFvPatchField<scalar>::typeName
// );
// addPatchFields<surfaceVectorField>
// (
// mesh,
// calculatedFvPatchField<vector>::typeName
// );
// addPatchFields<surfaceSphericalTensorField>
// (
// mesh,
// calculatedFvPatchField<sphericalTensor>::typeName
// );
// addPatchFields<surfaceSymmTensorField>
// (
// mesh,
// calculatedFvPatchField<symmTensor>::typeName
// );
// addPatchFields<surfaceTensorField>
// (
// mesh,
// calculatedFvPatchField<tensor>::typeName
// );
//
// // Point fields
//
// addPatchFields<pointScalarField>
// (
// mesh,
// calculatedFvPatchField<scalar>::typeName
// );
// addPatchFields<pointVectorField>
// (
// mesh,
// calculatedFvPatchField<vector>::typeName
// );
//}
//
//
//void addAllPatchFields(fvMesh& mesh, const label insertPatchI)
//{
// polyBoundaryMesh& polyPatches =
// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
//
// label sz = polyPatches.size();
//
// addPatchFields<volScalarField>
// (
// mesh,
// calculatedFvPatchField<scalar>::typeName
// );
// addPatchFields<volVectorField>
// (
// mesh,
// calculatedFvPatchField<vector>::typeName
// );
// addPatchFields<volSphericalTensorField>
// (
// mesh,
// calculatedFvPatchField<sphericalTensor>::typeName
// );
// addPatchFields<volSymmTensorField>
// (
// mesh,
// calculatedFvPatchField<symmTensor>::typeName
// );
// addPatchFields<volTensorField>
// (
// mesh,
// calculatedFvPatchField<tensor>::typeName
// );
//
// // Surface fields
//
// addPatchFields<surfaceScalarField>
// (
// mesh,
// calculatedFvPatchField<scalar>::typeName
// );
// addPatchFields<surfaceVectorField>
// (
// mesh,
// calculatedFvPatchField<vector>::typeName
// );
// addPatchFields<surfaceSphericalTensorField>
// (
// mesh,
// calculatedFvPatchField<sphericalTensor>::typeName
// );
// addPatchFields<surfaceSymmTensorField>
// (
// mesh,
// calculatedFvPatchField<symmTensor>::typeName
// );
// addPatchFields<surfaceTensorField>
// (
// mesh,
// calculatedFvPatchField<tensor>::typeName
// );
//
// // Create reordering list
// // patches before insert position stay as is
// labelList oldToNew(sz);
// for (label i = 0; i < insertPatchI; i++)
// {
// oldToNew[i] = i;
// }
// // patches after insert position move one up
// for (label i = insertPatchI; i < sz-1; i++)
// {
// oldToNew[i] = i+1;
// }
// // appended patch gets moved to insert position
// oldToNew[sz-1] = insertPatchI;
//
// // Shuffle into place
// polyPatches.reorder(oldToNew);
// fvPatches.reorder(oldToNew);
//
// reorderPatchFields<volScalarField>(mesh, oldToNew);
// reorderPatchFields<volVectorField>(mesh, oldToNew);
// reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
// reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
// reorderPatchFields<volTensorField>(mesh, oldToNew);
// reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
// reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
// reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
// reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
// reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
//}
//// Adds patch if not yet there. Returns patchID.
//template<class PatchType>
//label addPatch(fvMesh& mesh, const word& patchName, const dictionary& dict)
//{
// polyBoundaryMesh& polyPatches =
// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
//
// label patchI = polyPatches.findPatchID(patchName);
// if (patchI != -1)
// {
// if (isA<PatchType>(polyPatches[patchI]))
// {
// // Already there
// return patchI;
// }
// else
// {
// FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
// << "Already have patch " << patchName
// << " but of type " << PatchType::typeName
// << exit(FatalError);
// }
// }
//
//
// label insertPatchI = polyPatches.size();
// label startFaceI = mesh.nFaces();
//
// forAll(polyPatches, patchI)
// {
// const polyPatch& pp = polyPatches[patchI];
//
// if (isA<processorPolyPatch>(pp))
// {
// insertPatchI = patchI;
// startFaceI = pp.start();
// break;
// }
// }
//
// dictionary patchDict(dict);
// patchDict.set("type", PatchType::typeName);
// patchDict.set("nFaces", 0);
// patchDict.set("startFace", startFaceI);
//
//
// // Below is all quite a hack. Feel free to change once there is a better
// // mechanism to insert and reorder patches.
//
// // Clear local fields and e.g. polyMesh parallelInfo.
// mesh.clearOut();
//
// label sz = polyPatches.size();
//
// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
//
// // Add polyPatch at the end
// polyPatches.setSize(sz+1);
// polyPatches.set
// (
// sz,
// polyPatch::New
// (
// patchName,
// patchDict,
// insertPatchI,
// polyPatches
// )
// );
// fvPatches.setSize(sz+1);
// fvPatches.set
// (
// sz,
// fvPatch::New
// (
// polyPatches[sz], // point to newly added polyPatch
// mesh.boundary()
// )
// );
//
// addAllPatchFields(mesh, insertPatchI);
//
// return insertPatchI;
//}
//
//
//template<class PatchType>
//label addPatch(fvMesh& mesh, const word& patchName)
//{
//Pout<< "addPatch:" << patchName << endl;
//
// polyBoundaryMesh& polyPatches =
// const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
//
// label patchI = polyPatches.findPatchID(patchName);
// if (patchI != -1)
// {
// if (isA<PatchType>(polyPatches[patchI]))
// {
// // Already there
// return patchI;
// }
// else
// {
// FatalErrorIn("addPatch<PatchType>(fvMesh&, const word&)")
// << "Already have patch " << patchName
// << " but of type " << PatchType::typeName
// << exit(FatalError);
// }
// }
//
//
// label insertPatchI = polyPatches.size();
// label startFaceI = mesh.nFaces();
//
// forAll(polyPatches, patchI)
// {
// const polyPatch& pp = polyPatches[patchI];
//
// if (isA<processorPolyPatch>(pp))
// {
// insertPatchI = patchI;
// startFaceI = pp.start();
// break;
// }
// }
//
// // Below is all quite a hack. Feel free to change once there is a better
// // mechanism to insert and reorder patches.
//
// // Clear local fields and e.g. polyMesh parallelInfo.
// mesh.clearOut();
//
// label sz = polyPatches.size();
//
// fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
//
// // Add polyPatch at the end
// polyPatches.setSize(sz+1);
// polyPatches.set
// (
// sz,
// polyPatch::New
// (
// PatchType::typeName,
// patchName,
// 0, // size
// startFaceI,
// insertPatchI,
// polyPatches
// )
// );
// fvPatches.setSize(sz+1);
// fvPatches.set
// (
// sz,
// fvPatch::New
// (
// polyPatches[sz], // point to newly added polyPatch
// mesh.boundary()
// )
// );
//
// addAllPatchFields(mesh, insertPatchI);
//
// return insertPatchI;
//}
label findPatchID(const List<polyPatch*>& newPatches, const word& name)
{
forAll(newPatches, i)
{
if (newPatches[i]->name() == name)
{
return i;
}
}
return -1;
}
template<class PatchType>
label addPatch
(
const polyBoundaryMesh& patches,
const word& patchName,
DynamicList<polyPatch*>& newPatches
)
{
label patchI = findPatchID(newPatches, patchName);
if (patchI != -1)
{
if (isA<PatchType>(*newPatches[patchI]))
{
// Already there
return patchI;
}
else
{
FatalErrorIn
(
"addPatch<PatchType>(const polyBoundaryMesh&,"
" const word&, DynamicList<polyPatch*>)"
) << "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();
}
newPatches.append
(
polyPatch::New
(
PatchType::typeName,
patchName,
0, // size
startFaceI, // nFaces
patchI,
patches
).ptr()
);
return patchI;
}
template<class PatchType>
label addPatch
(
const polyBoundaryMesh& patches,
const word& patchName,
const dictionary& dict,
DynamicList<polyPatch*>& newPatches
)
{
label patchI = findPatchID(newPatches, patchName);
if (patchI != -1)
{
if (isA<PatchType>(*newPatches[patchI]))
{
// Already there
return patchI;
}
else
{
FatalErrorIn
(
"addPatch<PatchType>(const polyBoundaryMesh&,"
" const word&, DynamicList<polyPatch*>)"
) << "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();
}
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;
}
// Reorder and delete patches.
void reorderPatches
(
fvMesh& mesh,
const labelList& oldToNew,
const label nNewPatches
)
{
polyBoundaryMesh& polyPatches =
const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
// Shuffle into place
polyPatches.reorder(oldToNew);
fvPatches.reorder(oldToNew);
reorderPatchFields<volScalarField>(mesh, oldToNew);
reorderPatchFields<volVectorField>(mesh, oldToNew);
reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
reorderPatchFields<volTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
reorderPatchFields<pointScalarField>(mesh, oldToNew);
reorderPatchFields<pointVectorField>(mesh, oldToNew);
// Remove last.
polyPatches.setSize(nNewPatches);
fvPatches.setSize(nNewPatches);
trimPatchFields<volScalarField>(mesh, nNewPatches);
trimPatchFields<volVectorField>(mesh, nNewPatches);
trimPatchFields<volSphericalTensorField>(mesh, nNewPatches);
trimPatchFields<volSymmTensorField>(mesh, nNewPatches);
trimPatchFields<volTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceScalarField>(mesh, nNewPatches);
trimPatchFields<surfaceVectorField>(mesh, nNewPatches);
trimPatchFields<surfaceSphericalTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceSymmTensorField>(mesh, nNewPatches);
trimPatchFields<surfaceTensorField>(mesh, nNewPatches);
trimPatchFields<pointScalarField>(mesh, nNewPatches);
trimPatchFields<pointVectorField>(mesh, nNewPatches);
}
// Remove zero-sized patches
void deleteEmptyPatches(fvMesh& mesh)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
wordList masterNames;
if (Pstream::master())
{
masterNames = patches.names();
}
Pstream::scatter(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<processorPolyPatch>(patches[patchI]))
{
// Similar named processor patch? Not 'possible'.
if (patches[patchI].size() == 0)
{
Pout<< "Deleting processor patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
}
else
{
oldToNew[patchI] = usedI++;
}
}
else
{
// Common patch.
if (returnReduce(patches[patchI].size(), sumOp<label>()) == 0)
{
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].size() == 0)
{
Pout<< "Deleting processor patch " << patchI
<< " name:" << patches[patchI].name()
<< endl;
oldToNew[patchI] = --notUsedI;
}
else
{
oldToNew[patchI] = usedI++;
}
}
}
reorderPatches(mesh, oldToNew, usedI);
}
void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
{
// Create dummy system/fv*
{
IOobject io
(
"fvSchemes",
mesh.time().system(),
regionName,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
);
Info<< "Testing:" << io.objectPath() << endl;
if (!io.headerOk())
{
Info<< "Writing dummy " << regionName/io.name() << endl;
dictionary dummyDict;
dictionary divDict;
dummyDict.add("divSchemes", divDict);
dictionary gradDict;
dummyDict.add("gradSchemes", gradDict);
dictionary laplDict;
dummyDict.add("laplacianSchemes", laplDict);
IOdictionary(io, dummyDict).regIOobject::write();
}
}
{
IOobject io
(
"fvSolution",
mesh.time().system(),
regionName,
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
);
if (!io.headerOk())
{
Info<< "Writing dummy " << regionName/io.name() << endl;
dictionary dummyDict;
IOdictionary(io, dummyDict).regIOobject::write();
}
}
}
// 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))
{
FatalErrorIn("checkZoneInside(..)")
<< "Zone " << zoneNames[zoneI]
<< " is not consistently all internal or all boundary faces."
<< " Face " << faceI << " at " << mesh.faceCentres()[faceI]
<< " is the first occurrence."
<< exit(FatalError);
}
}
}
// To combineReduce a labelList. Filters out duplicates.
class uniqueEqOp
{
public:
void operator()(labelList& x, const labelList& y) const
{
if (x.empty())
{
if (y.size())
{
x = y;
}
}
else
{
forAll(y, yi)
{
if (findIndex(x, y[yi]) == -1)
{
label sz = x.size();
x.setSize(sz+1);
x[sz] = y[yi];
}
}
}
}
};
// 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)
{
const labelList& eFaces = edgeFaces[edgeI];
// Store pp face and processor as unique tag.
labelList& globalEFaces = globalEdgeFaces[edgeI];
globalEFaces.setSize(eFaces.size());
forAll(eFaces, i)
{
globalEFaces[i] = globalFaces.toGlobal(eFaces[i]);
}
}
// Synchronise across coupled edges.
syncTools::syncEdgeList
(
mesh,
ppMeshEdges,
globalEdgeFaces,
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 UIndirectList<label>& extrudeMeshFaces,// mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
forAll(extrudeMeshFaces, i)
{
extrudeFaceSet.insert(extrudeMeshFaces[i]);
}
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 UIndirectList<label>& extrudeMeshFaces,// mesh faces that are extruded
const label meshEdgeI // mesh edge
)
{
// Make set of extruded faces.
labelHashSet extrudeFaceSet(extrudeMeshFaces.size());
forAll(extrudeMeshFaces, i)
{
extrudeFaceSet.insert(extrudeMeshFaces[i]);
}
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<cyclicPolyPatch>(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] = min(minZoneID[edgeI], zoneI);
maxZoneID[edgeI] = max(maxZoneID[edgeI], zoneI);
}
}
}
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
minZoneID,
minOp<label>(),
labelMax // null value
);
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
maxZoneID,
maxOp<label>(),
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();
WarningIn("countExtrudePatches(..)")
<< "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,
UIndirectList<label>(extrudeMeshFaces, eFaces),
extrudeMeshEdges[edgeI]
);
if (faceI == -1)
{
zoneSidePatch[minZoneID[edgeI]]++;
}
}
}
// Synchronise decistion. Actual numbers are not important, just make
// sure that they're > 0 on all processors.
Pstream::listCombineGather(zoneSidePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneSidePatch);
Pstream::listCombineGather(zoneZonePatch, plusEqOp<label>());
Pstream::listCombineScatter(zoneZonePatch);
}
void addCouplingPatches
(
const fvMesh& mesh,
const word& regionName,
const word& shellRegionName,
const wordList& zoneNames,
const wordList& zoneShadowNames,
const boolList& isInternal,
const labelList& zoneIDs,
DynamicList<polyPatch*>& 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<mappedWallPolyPatch>
(
mesh.boundaryMesh(),
interName + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
(
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<polyPatch>
(
mesh.boundaryMesh(),
zoneNames[zoneI] + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
(
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<mappedWallPolyPatch>
(
mesh.boundaryMesh(),
zoneShadowNames[zoneI] + "_top",
newPatches
);
Pout<< interRegionTopPatch[zoneI]
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->name()
<< '\t' << newPatches[interRegionTopPatch[zoneI]]->type()
<< nl;
interRegionBottomPatch[zoneI] = addPatch<mappedWallPolyPatch>
(
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<polyPatch*>& 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] = min(minProcID[edgeI], procI);
maxProcID[edgeI] = max(maxProcID[edgeI], procI);
}
}
}
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
minProcID,
minOp<label>(),
labelMax // null value
);
syncTools::syncEdgeList
(
mesh,
extrudeMeshEdges,
maxProcID,
maxOp<label>(),
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,
UIndirectList<label>(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
{
FatalErrorIn
(
"void addCoupledPatches"
"("
"const fvMesh&, "
"const primitiveFacePatch&, "
"const labelList&, "
"const labelList&, "
"const mapDistribute&, "
"const labelListList&, "
"labelList&, "
"DynamicList<polyPatch*>&"
")"
) << "Unable to determine coupled patch addressing"
<< abort(FatalError);
}
}
else
{
// Rrocessor patch
word name =
"procBoundary"
+ Foam::name(Pstream::myProcNo())
+ "to"
+ Foam::name(nbrProcI);
sidePatchID[edgeI] = findPatchID(newPatches, name);
if (sidePatchID[edgeI] == -1)
{
dictionary patchDict;
patchDict.add("myProcNo", Pstream::myProcNo());
patchDict.add("neighbProcNo", nbrProcI);
sidePatchID[edgeI] = addPatch<processorPolyPatch>
(
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<polyPatch*>& 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 != word::null)
{
// Reuse single empty patch.
word patchName;
if (oneDPolyPatchType == "empty")
{
patchName = "oneDEmptyPatch";
zoneSidePatch[zoneI] = addPatch<emptyPolyPatch>
(
mesh.boundaryMesh(),
patchName,
newPatches
);
}
else if (oneDPolyPatchType == "wedge")
{
patchName = "oneDWedgePatch";
zoneSidePatch[zoneI] = addPatch<wedgePolyPatch>
(
mesh.boundaryMesh(),
patchName,
newPatches
);
}
else
{
FatalErrorIn("addZoneSidePatches()")
<< "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<polyPatch>
(
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<polyPatch*>& 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<nonuniformTransformCyclicPolyPatch>
(
mesh.boundaryMesh(),
minToMax,
transformDict,
newPatches
);
Pout<< zoneZonePatch_min[index] << '\t' << minToMax
<< nl;
}
{
transformDict.set("neighbourPatch", minToMax);
zoneZonePatch_max[index] =
addPatch<nonuniformTransformCyclicPolyPatch>
(
mesh.boundaryMesh(),
maxToMin,
transformDict,
newPatches
);
Pout<< zoneZonePatch_max[index] << '\t' << maxToMin
<< nl;
}
}
}
}
}
Pout<< "Added " << newPatches.size()-nOldPatches << " inter-zone patches."
<< nl << endl;
}
tmp<pointField> calcOffset
(
const primitiveFacePatch& extrudePatch,
const createShellMesh& extruder,
const polyPatch& pp
)
{
vectorField::subField fc = pp.faceCentres();
tmp<pointField> toffsets(new pointField(fc.size()));
pointField& offsets = toffsets();
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 mappedWallPolyPatch::sampleMode mode,
const List<pointField>& offsets
)
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
List<polyPatch*> newPatches(patches.size(), static_cast<polyPatch*>(NULL));
forAll(zoneToPatch, zoneI)
{
label patchI = zoneToPatch[zoneI];
if (patchI != -1)
{
const polyPatch& pp = patches[patchI];
if (isA<mappedWallPolyPatch>(pp))
{
newPatches[patchI] = new mappedWallPolyPatch
(
pp.name(),
pp.size(),
pp.start(),
patchI,
sampleRegion, // sampleRegion
mode, // sampleMode
pp.name(), // samplePatch
offsets[zoneI], // offset
patches
);
}
}
}
forAll(newPatches, patchI)
{
if (!newPatches[patchI])
{
newPatches[patchI] = patches[patchI].clone(patches).ptr();
}
}
mesh.removeFvBoundary();
mesh.addFvPatches(newPatches, true);
}
// Main program:
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", "name", "specify alternative dictionary");
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedMesh.H"
if (mesh.boundaryMesh().checkParallelSync(true))
{
List<wordList> allNames(Pstream::nProcs());
allNames[Pstream::myProcNo()] = mesh.boundaryMesh().names();
Pstream::gatherList(allNames);
Pstream::scatterList(allNames);
FatalErrorIn(args.executable())
<< "Patches are not synchronised on all processors."
<< " Per processor patches " << allNames
<< exit(FatalError);
}
const word oldInstance = mesh.pointsInstance();
bool overwrite = args.optionFound("overwrite");
const word dictName
(args.optionLookupOrDefault<word>("dict", "extrudeToRegionMeshDict"));
IOdictionary dict
(
IOobject
(
dictName,
runTime.system(),
runTime,
IOobject::MUST_READ_IF_MODIFIED
)
);
// Point generator
autoPtr<extrudeModel> model(extrudeModel::New(dict));
// Region
const word shellRegionName(dict.lookup("region"));
// Faces to extrude - either faceZones or faceSets (boundary faces only)
wordList zoneNames;
wordList zoneShadowNames;
bool hasZones = dict.found("faceZones");
if (hasZones)
{
dict.lookup("faceZones") >> zoneNames;
dict.readIfPresent("faceZonesShadow", zoneShadowNames);
// Check
if (dict.found("faceSets"))
{
FatalIOErrorIn(args.executable().c_str(), dict)
<< "Please supply faces to extrude either through 'faceZones'"
<< " or 'faceSets' entry. Found both."
<< exit(FatalIOError);
}
}
else
{
dict.lookup("faceSets") >> zoneNames;
dict.readIfPresent("faceSetsShadow", zoneShadowNames);
}
mappedPatchBase::sampleMode sampleMode =
mappedPatchBase::sampleModeNames_[dict.lookup("sampleMode")];
const Switch oneD(dict.lookup("oneD"));
const Switch adaptMesh(dict.lookup("adaptMesh"));
if (hasZones)
{
Pout<< "Extruding zones " << zoneNames
<< " on mesh " << regionName
<< " into shell mesh " << shellRegionName
<< endl;
}
else
{
Pout<< "Extruding faceSets " << zoneNames
<< " on mesh " << regionName
<< " into shell mesh " << shellRegionName
<< endl;
}
if (shellRegionName == regionName)
{
FatalIOErrorIn(args.executable().c_str(), dict)
<< "Cannot extrude into same region as mesh." << endl
<< "Mesh region : " << regionName << endl
<< "Shell region : " << shellRegionName
<< exit(FatalIOError);
}
//// Read objects in time directory
//IOobjectList objects(mesh, runTime.timeName());
//
//// Read vol fields.
//
//PtrList<volScalarField> vsFlds;
//ReadFields(mesh, objects, vsFlds);
//
//PtrList<volVectorField> vvFlds;
//ReadFields(mesh, objects, vvFlds);
//
//PtrList<volSphericalTensorField> vstFlds;
//ReadFields(mesh, objects, vstFlds);
//
//PtrList<volSymmTensorField> vsymtFlds;
//ReadFields(mesh, objects, vsymtFlds);
//
//PtrList<volTensorField> vtFlds;
//ReadFields(mesh, objects, vtFlds);
//
//// Read surface fields.
//
//PtrList<surfaceScalarField> ssFlds;
//ReadFields(mesh, objects, ssFlds);
//
//PtrList<surfaceVectorField> svFlds;
//ReadFields(mesh, objects, svFlds);
//
//PtrList<surfaceSphericalTensorField> sstFlds;
//ReadFields(mesh, objects, sstFlds);
//
//PtrList<surfaceSymmTensorField> ssymtFlds;
//ReadFields(mesh, objects, ssymtFlds);
//
//PtrList<surfaceTensorField> stFlds;
//ReadFields(mesh, objects, stFlds);
//
//// Read point fields.
//
//PtrList<pointScalarField> psFlds;
//ReadFields(pointMesh::New(mesh), objects, psFlds);
//
//PtrList<pointVectorField> pvFlds;
//ReadFields(pointMesh::New(mesh), objects, pvFlds);
// Create dummy fv* files
createDummyFvMeshFiles(mesh, shellRegionName);
word meshInstance;
if (!overwrite)
{
runTime++;
meshInstance = runTime.timeName();
}
else
{
meshInstance = oldInstance;
}
Pout<< "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)
{
FatalIOErrorIn(args.executable().c_str(), 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)
{
FatalIOErrorIn(args.executable().c_str(), 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] = zoneShadowIDs[i];
nShadowFaces++;
}
}
}
}
else
{
meshZoneID.setSize(zoneNames.size(), -1);
// Load faceSets
PtrList<faceSet> 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];
forAllConstIter(faceSet, fz, iter)
{
label faceI = iter.key();
if (mesh.isInternalFace(faceI))
{
FatalIOErrorIn(args.executable().c_str(), 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<faceSet> 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;
forAll(shadowZones, i)
{
nShadowFaces += shadowZones[i].size();
}
if (nExtrudeFaces != nShadowFaces)
{
FatalIOErrorIn(args.executable().c_str(), 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];
forAllConstIter(faceSet, fz, iter)
{
label faceI = iter.key();
if (mesh.isInternalFace(faceI))
{
FatalIOErrorIn(args.executable().c_str(), 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(zoneFaces.xfer(), mesh.points());
Pstream::listCombineGather(isInternal, orEqOp<bool>());
Pstream::listCombineScatter(isInternal);
// 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();
Pout<< "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<Map<label> > 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
);
DynamicList<polyPatch*> regionPatches(patches.size());
// Copy all non-local patches since these are used on boundary edges of
// the extrusion
forAll(patches, patchI)
{
if (!isA<processorPolyPatch>(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
// Clone existing patches
DynamicList<polyPatch*> newPatches(patches.size());
forAll(patches, patchI)
{
newPatches.append(patches[patchI].clone(patches).ptr());
}
// Add new patches
addCouplingPatches
(
mesh,
regionName,
shellRegionName,
zoneNames,
zoneShadowNames,
isInternal,
meshZoneID,
newPatches,
interMeshTopPatch,
interMeshBottomPatch
);
// 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(), 0);
// Patch to use for minZone
labelList zoneZonePatch_min(zoneNames.size()*zoneNames.size(), 0);
// Patch to use for maxZone
labelList zoneZonePatch_max(zoneNames.size()*zoneNames.size(), 0);
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 ? dict.lookup("oneDPolyPatchType") : 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
PackedBoolList 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]]];
}
//- 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[edgeI] = 1;
}
}
else if (eFaces.size() == 2)
{
label zone0 = zoneID[eFaces[0]];
label zone1 = zoneID[eFaces[1]];
if (zone0 != zone1) // || (cos(angle) > blabla))
{
label minZone = min(zone0,zone1);
label maxZone = 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[edgeI] = 1;
}
}
else if (sidePatchID[edgeI] != -1)
{
// Coupled extrusion
ePatches.setSize(eFaces.size());
forAll(eFaces, i)
{
ePatches[i] = sidePatchID[edgeI];
}
}
else
{
label faceI = findUncoveredPatchFace
(
mesh,
UIndirectList<label>(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[edgeI] = 1;
}
}
// 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(), vector::zero);
forAll(pointLocalRegions, faceI)
{
const face& pRegions = pointLocalRegions[faceI];
forAll(pRegions, fp)
{
label localRegionI = pRegions[fp];
localSum[localRegionI] += extrudePatch.faceNormals()[faceI];
}
}
Map<point> globalSum(2*localToGlobalRegion.size());
forAll(localSum, localRegionI)
{
label globalRegionI = localToGlobalRegion[localRegionI];
globalSum.insert(globalRegionI, localSum[localRegionI]);
}
// Reduce
Pstream::mapCombineGather(globalSum, plusEqOp<point>());
Pstream::mapCombineScatter(globalSum);
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<mapPolyMesh> shellMap;
fvMesh regionMesh
(
IOobject
(
shellRegionName,
meshInstance,
runTime,
IOobject::NO_READ,
IOobject::AUTO_WRITE,
false
),
xferCopy(pointField()),
xferCopy(faceList()),
xferCopy(labelList()),
xferCopy(labelList()),
false
);
// Add the new patches
forAll(regionPatches, patchI)
{
regionPatches[patchI] = regionPatches[patchI]->clone
(
regionMesh.boundaryMesh()
).ptr();
}
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
);
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<pointField> topOffsets(zoneNames.size());
List<pointField> bottomOffsets(zoneNames.size());
forAll(regionMesh.boundaryMesh(), patchI)
{
const polyPatch& pp = regionMesh.boundaryMesh()[patchI];
if (isA<mappedWallPolyPatch>(pp))
{
if (findIndex(interRegionTopPatch, patchI) != -1)
{
label zoneI = findIndex(interRegionTopPatch, patchI);
topOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
}
else if (findIndex(interRegionBottomPatch, patchI) != -1)
{
label zoneI = findIndex(interRegionBottomPatch, patchI);
bottomOffsets[zoneI] = calcOffset(extrudePatch, extruder, pp);
}
}
}
// Change top and bottom boundary conditions on regionMesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
{
// Correct top patches for offset
setCouplingInfo
(
regionMesh,
interRegionTopPatch,
regionName, // name of main mesh
sampleMode, // sampleMode
topOffsets
);
// Correct bottom patches for offset
setCouplingInfo
(
regionMesh,
interRegionBottomPatch,
regionName,
sampleMode, // sampleMode
bottomOffsets
);
// 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
);
// Correct bottom patches for offset
setCouplingInfo
(
mesh,
interMeshBottomPatch,
shellRegionName,
sampleMode,
-bottomOffsets
);
}
// Write addressing from region mesh back to originating patch
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelIOList cellToPatchFaceAddressing
(
IOobject
(
"cellToPatchFaceAddressing",
regionMesh.facesInstance(),
regionMesh.meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
extruder.cellToFaceMap()
);
cellToPatchFaceAddressing.note() = "cell to patch face addressing";
labelIOList faceToPatchFaceAddressing
(
IOobject
(
"faceToPatchFaceAddressing",
regionMesh.facesInstance(),
regionMesh.meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
extruder.faceToFaceMap()
);
faceToPatchFaceAddressing.note() =
"front/back face + turning index to patch face addressing";
labelIOList faceToPatchEdgeAddressing
(
IOobject
(
"faceToPatchEdgeAddressing",
regionMesh.facesInstance(),
regionMesh.meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
extruder.faceToEdgeMap()
);
faceToPatchEdgeAddressing.note() =
"side face to patch edge addressing";
labelIOList pointToPatchPointAddressing
(
IOobject
(
"pointToPatchPointAddressing",
regionMesh.facesInstance(),
regionMesh.meshSubDir,
regionMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
extruder.pointToPointMap()
);
pointToPatchPointAddressing.note() =
"point to patch point addressing";
Pout<< "Writing mesh " << regionMesh.name()
<< " to " << regionMesh.facesInstance() << nl
<< endl;
bool ok =
regionMesh.write()
&& cellToPatchFaceAddressing.write()
&& faceToPatchFaceAddressing.write()
&& faceToPatchEdgeAddressing.write()
&& pointToPatchPointAddressing.write();
if (!ok)
{
FatalErrorIn(args.executable())
<< "Failed writing mesh " << regionMesh.name()
<< " at location " << regionMesh.facesInstance()
<< exit(FatalError);
}
// Insert baffles into original mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
autoPtr<mapPolyMesh> 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);
Pout<< "Writing mesh " << mesh.name()
<< " to " << mesh.facesInstance() << nl
<< endl;
if (!mesh.write())
{
FatalErrorIn(args.executable())
<< "Failed writing mesh " << mesh.name()
<< " at location " << mesh.facesInstance()
<< exit(FatalError);
}
}
Info << "End\n" << endl;
return 0;
}
// ************************************************************************* //
Reference in New Issue View Git Blame Copy Permalink