/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
C | 394 lines | 276 code | 70 blank | 48 comment | 34 complexity | 75d7f2dadc3a732f4d1f302f771dc035 MD5 | raw file
- /*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2013-2015 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/>.
- Application
- orientFaceZone
- Description
- Corrects orientation of faceZone.
- - correct in parallel - excludes coupled faceZones from walk
- - correct for non-manifold faceZones - restarts walk
- \*---------------------------------------------------------------------------*/
- #include "argList.H"
- #include "Time.H"
- #include "syncTools.H"
- #include "patchFaceOrientation.H"
- #include "PatchEdgeFaceWave.H"
- #include "orientedSurface.H"
- #include "globalIndex.H"
- using namespace Foam;
- // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
- int main(int argc, char *argv[])
- {
- #include "addRegionOption.H"
- argList::validArgs.append("faceZone");
- argList::validArgs.append("outsidePoint");
- #include "setRootCase.H"
- #include "createTime.H"
- #include "createNamedPolyMesh.H"
- const word zoneName = args[1];
- const point outsidePoint = args.argRead<point>(2);
- Info<< "Orienting faceZone " << zoneName
- << " such that " << outsidePoint << " is outside"
- << nl << endl;
- const faceZone& fZone = mesh.faceZones()[zoneName];
- if (fZone.checkParallelSync())
- {
- FatalErrorIn(args.executable())
- << "Face zone " << fZone.name()
- << " is not parallel synchronised."
- << " Any coupled face also needs its coupled version to be included"
- << " and with opposite flipMap."
- << exit(FatalError);
- }
- const labelList& faceLabels = fZone;
- const indirectPrimitivePatch patch
- (
- IndirectList<face>(mesh.faces(), faceLabels),
- mesh.points()
- );
- const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
- // Data on all edges and faces
- List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
- List<patchFaceOrientation> allFaceInfo(patch.size());
- // Make sure we don't walk through
- // - slaves of coupled faces
- // - non-manifold edges
- {
- const polyBoundaryMesh& bm = mesh.boundaryMesh();
- label nProtected = 0;
- forAll(faceLabels, faceI)
- {
- const label meshFaceI = faceLabels[faceI];
- const label patchI = bm.whichPatch(meshFaceI);
- if
- (
- patchI != -1
- && bm[patchI].coupled()
- && !isMasterFace[meshFaceI]
- )
- {
- // Slave side. Mark so doesn't get visited.
- allFaceInfo[faceI] = orientedSurface::NOFLIP;
- nProtected++;
- }
- }
- Info<< "Protected from visiting "
- << returnReduce(nProtected, sumOp<label>())
- << " slaves of coupled faces" << nl << endl;
- }
- {
- // Number of (master)faces per edge
- labelList nMasterFaces(patch.nEdges(), 0);
- forAll(faceLabels, faceI)
- {
- const label meshFaceI = faceLabels[faceI];
- if (isMasterFace[meshFaceI])
- {
- const labelList& fEdges = patch.faceEdges()[faceI];
- forAll(fEdges, fEdgeI)
- {
- nMasterFaces[fEdges[fEdgeI]]++;
- }
- }
- }
- syncTools::syncEdgeList
- (
- mesh,
- patch.meshEdges(mesh.edges(), mesh.pointEdges()),
- nMasterFaces,
- plusEqOp<label>(),
- label(0)
- );
- label nProtected = 0;
- forAll(nMasterFaces, edgeI)
- {
- if (nMasterFaces[edgeI] > 2)
- {
- allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
- nProtected++;
- }
- }
- Info<< "Protected from visiting "
- << returnReduce(nProtected, sumOp<label>())
- << " non-manifold edges" << nl << endl;
- }
- DynamicList<label> changedEdges;
- DynamicList<patchFaceOrientation> changedInfo;
- const scalar tol = PatchEdgeFaceWave
- <
- indirectPrimitivePatch,
- patchFaceOrientation
- >::propagationTol();
- int dummyTrackData;
- globalIndex globalFaces(patch.size());
- while (true)
- {
- // Pick an unset face
- label unsetFaceI = labelMax;
- forAll(allFaceInfo, faceI)
- {
- if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
- {
- unsetFaceI = globalFaces.toGlobal(faceI);
- break;
- }
- }
- reduce(unsetFaceI, minOp<label>());
- if (unsetFaceI == labelMax)
- {
- break;
- }
- label procI = globalFaces.whichProcID(unsetFaceI);
- label seedFaceI = globalFaces.toLocal(procI, unsetFaceI);
- Info<< "Seeding from processor " << procI << " face " << seedFaceI
- << endl;
- if (procI == Pstream::myProcNo())
- {
- // Determine orientation of seedFace
- vector d = outsidePoint-patch.faceCentres()[seedFaceI];
- const vector& fn = patch.faceNormals()[seedFaceI];
- // Set information to correct orientation
- patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];
- faceInfo = orientedSurface::NOFLIP;
- if ((fn&d) < 0)
- {
- faceInfo.flip();
- Pout<< "Face " << seedFaceI << " at "
- << patch.faceCentres()[seedFaceI]
- << " with normal " << fn
- << " needs to be flipped." << endl;
- }
- else
- {
- Pout<< "Face " << seedFaceI << " at "
- << patch.faceCentres()[seedFaceI]
- << " with normal " << fn
- << " points in positive direction (cos = " << (fn&d)/mag(d)
- << ")" << endl;
- }
- const labelList& fEdges = patch.faceEdges()[seedFaceI];
- forAll(fEdges, fEdgeI)
- {
- label edgeI = fEdges[fEdgeI];
- patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
- if
- (
- edgeInfo.updateEdge<int>
- (
- mesh,
- patch,
- edgeI,
- seedFaceI,
- faceInfo,
- tol,
- dummyTrackData
- )
- )
- {
- changedEdges.append(edgeI);
- changedInfo.append(edgeInfo);
- }
- }
- }
- if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
- {
- break;
- }
- // Walk
- PatchEdgeFaceWave
- <
- indirectPrimitivePatch,
- patchFaceOrientation
- > calc
- (
- mesh,
- patch,
- changedEdges,
- changedInfo,
- allEdgeInfo,
- allFaceInfo,
- returnReduce(patch.nEdges(), sumOp<label>())
- );
- }
- // Push master zone info over to slave (since slave faces never visited)
- {
- const polyBoundaryMesh& bm = mesh.boundaryMesh();
- labelList neiStatus
- (
- mesh.nFaces()-mesh.nInternalFaces(),
- orientedSurface::UNVISITED
- );
- forAll(faceLabels, i)
- {
- const label meshFaceI = faceLabels[i];
- if (!mesh.isInternalFace(meshFaceI))
- {
- neiStatus[meshFaceI-mesh.nInternalFaces()] =
- allFaceInfo[i].flipStatus();
- }
- }
- syncTools::swapBoundaryFaceList(mesh, neiStatus);
- forAll(faceLabels, i)
- {
- const label meshFaceI = faceLabels[i];
- const label patchI = bm.whichPatch(meshFaceI);
- if
- (
- patchI != -1
- && bm[patchI].coupled()
- && !isMasterFace[meshFaceI]
- )
- {
- // Slave side. Take flipped from neighbour
- label bFaceI = meshFaceI-mesh.nInternalFaces();
- if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
- {
- allFaceInfo[i] = orientedSurface::FLIP;
- }
- else if (neiStatus[bFaceI] == orientedSurface::FLIP)
- {
- allFaceInfo[i] = orientedSurface::NOFLIP;
- }
- else
- {
- FatalErrorIn(args.executable())
- << "Incorrect status for face " << meshFaceI
- << abort(FatalError);
- }
- }
- }
- }
- // Convert to flipmap and adapt faceZones
- boolList newFlipMap(allFaceInfo.size(), false);
- label nChanged = 0;
- forAll(allFaceInfo, faceI)
- {
- if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
- {
- newFlipMap[faceI] = false;
- }
- else if (allFaceInfo[faceI] == orientedSurface::FLIP)
- {
- newFlipMap[faceI] = true;
- }
- else
- {
- FatalErrorIn(args.executable())
- << "Problem : unvisited face " << faceI
- << " centre:" << mesh.faceCentres()[faceLabels[faceI]]
- << abort(FatalError);
- }
- if (fZone.flipMap()[faceI] != newFlipMap[faceI])
- {
- nChanged++;
- }
- }
- reduce(nChanged, sumOp<label>());
- if (nChanged > 0)
- {
- Info<< "Flipping " << nChanged << " out of "
- << globalFaces.size() << " faces." << nl << endl;
- mesh.faceZones()[zoneName].resetAddressing(faceLabels, newFlipMap);
- if (!mesh.faceZones().write())
- {
- FatalErrorIn(args.executable())
- << "Failed writing faceZones" << exit(FatalError);
- }
- }
- Info<< "\nEnd\n" << endl;
- return 0;
- }
- // ************************************************************************* //