diff --git a/applications/utilities/mesh/advanced/PDRMesh/Make/files b/applications/utilities/mesh/advanced/PDRMesh/Make/files new file mode 100644 index 0000000000..998d5028f2 --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/Make/files @@ -0,0 +1,3 @@ +PDRMesh.C + +EXE = $(FOAM_APPBIN)/PDRMesh diff --git a/applications/utilities/mesh/advanced/PDRMesh/Make/options b/applications/utilities/mesh/advanced/PDRMesh/Make/options new file mode 100644 index 0000000000..6ab51e876a --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/Make/options @@ -0,0 +1,11 @@ +EXE_INC = \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lmeshTools \ + -ldynamicMesh \ + -lfiniteVolume \ + -lcompressibleRASModels + diff --git a/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C new file mode 100644 index 0000000000..c69e102bc2 --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C @@ -0,0 +1,1181 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + Mesh and field preparation utility for PDR type simulations. + + Reads + - cellSet giving blockedCells + - faceSets giving blockedFaces and the patch they should go into + + NOTE: To avoid exposing wrong fields values faceSets should include + faces contained in the blockedCells cellset. + + - coupledFaces reads coupledFacesSet to introduces mixe-coupled baffles + + Subsets out the blocked cells and splits the blockedFaces and updates + fields. + + The face splitting is done by duplicating the faces. No duplication of + points for now (so checkMesh will show a lot of 'duplicate face' messages) + +\*---------------------------------------------------------------------------*/ + +#include "fvMeshSubset.H" +#include "argList.H" +#include "cellSet.H" +#include "IOobjectList.H" +#include "volFields.H" +#include "mapPolyMesh.H" +#include "faceSet.H" +#include "cellSet.H" +#include "syncTools.H" +#include "polyTopoChange.H" +#include "polyModifyFace.H" +#include "polyAddFace.H" +#include "regionSplit.H" +#include "Tuple2.H" +#include "cyclicFvPatch.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void modifyOrAddFace +( + polyTopoChange& meshMod, + const face& f, + const label faceI, + const label own, + const bool flipFaceFlux, + const label newPatchI, + const label zoneID, + const bool zoneFlip, + + PackedBoolList& modifiedFace +) +{ + if (!modifiedFace[faceI]) + { + // First usage of face. Modify. + meshMod.setAction + ( + polyModifyFace + ( + f, // modified face + faceI, // label of face + own, // owner + -1, // neighbour + flipFaceFlux, // face flip + newPatchI, // patch for face + false, // remove from zone + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + modifiedFace[faceI] = 1; + } + else + { + // Second or more usage of face. Add. + meshMod.setAction + ( + polyAddFace + ( + f, // modified face + own, // owner + -1, // neighbour + -1, // master point + -1, // master edge + faceI, // master face + flipFaceFlux, // face flip + newPatchI, // patch for face + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + } +} + + +template +void subsetVolFields +( + const fvMeshSubset& subsetter, + const IOobjectList& objectsList, + const label patchI, + const Type& exposedValue, + const word GeomVolType, + PtrList >& subFields +) +{ + const fvMesh& baseMesh = subsetter.baseMesh(); + + label i = 0; + + forAllConstIter(IOobjectList , objectsList, iter) + { + if (iter()->headerClassName() == GeomVolType) + { + const word fieldName = iter()->name(); + + Info<< "Subsetting field " << fieldName << endl; + + GeometricField volField + ( + *iter(), + baseMesh + ); + + subFields.set(i, subsetter.interpolate(volField)); + + // Explicitly set exposed faces (in patchI) to exposedValue. + if (patchI >= 0) + { + fvPatchField& fld = + subFields[i++].boundaryField()[patchI]; + + label newStart = fld.patch().patch().start(); + + label oldPatchI = subsetter.patchMap()[patchI]; + + if (oldPatchI == -1) + { + // New patch. Reset whole value. + fld = exposedValue; + } + else + { + // Reset those faces that originate from different patch + // or internal faces. + label oldSize = volField.boundaryField()[oldPatchI].size(); + label oldStart = volField.boundaryField() + [ + oldPatchI + ].patch().patch().start(); + + forAll(fld, j) + { + label oldFaceI = subsetter.faceMap()[newStart+j]; + + if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize) + { + fld[j] = exposedValue; + } + } + } + } + } + } +} + + +template +void subsetSurfaceFields +( + const fvMeshSubset& subsetter, + const IOobjectList& objectsList, + const label patchI, + const Type& exposedValue, + const word GeomSurfType, + PtrList >& subFields +) +{ + const fvMesh& baseMesh = subsetter.baseMesh(); + + label i(0); + + forAllConstIter(IOobjectList , objectsList, iter) + { + if (iter()->headerClassName() == GeomSurfType) + { + const word& fieldName = iter.key(); + + Info<< "Subsetting field " << fieldName << endl; + + GeometricField volField + ( + *iter(), + baseMesh + ); + + subFields.set(i, subsetter.interpolate(volField)); + + + // Explicitly set exposed faces (in patchI) to exposedValue. + if (patchI >= 0) + { + fvsPatchField& fld = + subFields[i++].boundaryField()[patchI]; + + label newStart = fld.patch().patch().start(); + + label oldPatchI = subsetter.patchMap()[patchI]; + + if (oldPatchI == -1) + { + // New patch. Reset whole value. + fld = exposedValue; + } + else + { + // Reset those faces that originate from different patch + // or internal faces. + label oldSize = volField.boundaryField()[oldPatchI].size(); + label oldStart = volField.boundaryField() + [ + oldPatchI + ].patch().patch().start(); + + forAll(fld, j) + { + label oldFaceI = subsetter.faceMap()[newStart+j]; + + if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize) + { + fld[j] = exposedValue; + } + } + } + } + } + } +} + + +// Faces introduced into zero-sized patches don't get a value at all. +// This is hack to set them to an initial value. +template +void initCreatedPatches +( + const fvMesh& mesh, + const mapPolyMesh& map, + const typename GeoField::value_type initValue +) +{ + HashTable fields + ( + mesh.objectRegistry::lookupClass() + ); + + for + ( + typename HashTable:: + iterator fieldIter = fields.begin(); + fieldIter != fields.end(); + ++fieldIter + ) + { + GeoField& field = const_cast(*fieldIter()); + + forAll(field.boundaryField(), patchi) + { + if (map.oldPatchSizes()[patchi] == 0) + { + // Not mapped. + field.boundaryField()[patchi] = initValue; + + if (field.boundaryField()[patchi].fixesValue()) + { + field.boundaryField()[patchi] == initValue; + } + } + } + } +} + + +void createCoupledBaffles +( + fvMesh& mesh, + const labelList& coupledWantedPatch, + polyTopoChange& meshMod, + PackedBoolList& modifiedFace +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + + forAll(coupledWantedPatch, faceI) + { + if (coupledWantedPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + // Use owner side of face + modifyOrAddFace + ( + meshMod, + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + false, // face flip + coupledWantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + // Use neighbour side of face + modifyOrAddFace + ( + meshMod, + f.reverseFace(), // modified face + faceI, // label of face + mesh.faceNeighbour()[faceI],// owner + false, // face flip + coupledWantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + } + } + } +} + + +void createCyclicCoupledBaffles +( + fvMesh& mesh, + const labelList& cyclicMasterPatch, + const labelList& cyclicSlavePatch, + polyTopoChange& meshMod, + PackedBoolList& modifiedFace +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + + forAll(cyclicMasterPatch, faceI) + { + if (cyclicMasterPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + modifyOrAddFace + ( + meshMod, + f.reverseFace(), // modified face + faceI, // label of face + mesh.faceNeighbour()[faceI], // owner + false, // face flip + cyclicMasterPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add + ); + } + } + + forAll(cyclicSlavePatch, faceI) + { + if (cyclicSlavePatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + // Use owner side of face + modifyOrAddFace + ( + meshMod, + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + false, // face flip + cyclicSlavePatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + } + } + } +} + + +void createBaffles +( + fvMesh& mesh, + const labelList& wantedPatch, + polyTopoChange& meshMod +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + Info << "faceZone:createBaffle " << faceZones << endl; + forAll(wantedPatch, faceI) + { + if (wantedPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + meshMod.setAction + ( + polyModifyFace + ( + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + -1, // neighbour + false, // face flip + wantedPatch[faceI], // patch for face + false, // remove from zone + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + meshMod.setAction + ( + polyAddFace + ( + f.reverseFace(), // modified face + mesh.faceNeighbour()[faceI],// owner + -1, // neighbour + -1, // masterPointID + -1, // masterEdgeID + faceI, // masterFaceID, + false, // face flip + wantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + } + } + } +} + + +// Wrapper around find patch. Also makes sure same patch in parallel. +label findPatch(const polyBoundaryMesh& patches, const word& patchName) +{ + label patchI = patches.findPatchID(patchName); + + if (patchI == -1) + { + FatalErrorIn("findPatch(const polyBoundaryMesh&, const word&)") + << "Illegal patch " << patchName + << nl << "Valid patches are " << patches.names() + << exit(FatalError); + } + + // Check same patch for all procs + { + label newPatch = patchI; + reduce(newPatch, minOp