From 7faf688709a022bdbbf1fe3cacc1167afae1d92b Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 23 Apr 2008 12:20:12 +0100 Subject: [PATCH] Initial mattijsLib merge --- src/Allwmake | 2 +- src/OpenFOAM/Make/files | 3 +- .../algorithms/MeshWave/FaceCellWave.C | 10 +- .../algorithms/MeshWave/FaceCellWave.H | 6 +- .../containers/Lists/PackedList/PackedList.H | 5 +- .../containers/Lists/PackedList/PackedListI.H | 6 +- src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H | 32 + .../IOstreams/Pstreams/combineGatherScatter.C | 141 +++++ src/OpenFOAM/meshes/meshShapes/face/face.H | 15 +- .../meshes/meshShapes/face/faceIntersection.C | 46 ++ src/OpenFOAM/meshes/meshTools/matchPoints.C | 20 +- .../mapPolyMesh/mapDistribute/mapDistribute.C | 197 +++++++ .../mapPolyMesh/mapDistribute/mapDistribute.H | 203 +++++++ .../mapDistribute/mapDistributeLagrangian.H | 146 +++++ .../mapDistribute/mapDistributePolyMesh.C | 235 ++++++++ .../mapDistribute/mapDistributePolyMesh.H | 291 ++++++++++ .../mapDistribute/mapDistributeTemplates.C | 190 ++++++ .../basic/coupled/coupledPolyPatch.C | 200 ++++--- .../basic/coupled/coupledPolyPatch.H | 41 +- .../constraint/cyclic/cyclicPolyPatch.C | 316 ++++++++-- .../constraint/cyclic/cyclicPolyPatch.H | 6 +- .../constraint/processor/processorPolyPatch.C | 433 ++++++++------ .../constraint/processor/processorPolyPatch.H | 22 +- .../polyMesh/syncTools/syncToolsTemplates.C | 445 +++++++------- .../PrimitivePatch/PrimitivePatch.C | 230 ++++++-- .../PrimitivePatch/PrimitivePatch.H | 61 +- .../PrimitivePatch/PrimitivePatchAddressing.C | 24 +- .../PrimitivePatch/PrimitivePatchBdryPoints.C | 18 +- .../PrimitivePatch/PrimitivePatchCheck.C | 53 +- .../PrimitivePatch/PrimitivePatchClear.C | 50 +- .../PrimitivePatch/PrimitivePatchEdgeLoops.C | 28 +- .../PrimitivePatchLocalPointOrder.C | 19 +- .../PrimitivePatch/PrimitivePatchMeshData.C | 112 ++-- .../PrimitivePatch/PrimitivePatchMeshEdges.C | 24 +- .../PrimitivePatchPointAddressing.C | 54 +- .../PrimitivePatchProjectPoints.C | 67 ++- .../primitiveShapes/triangle/triangle.H | 30 +- .../primitiveShapes/triangle/triangleI.H | 221 ++++--- src/dynamicMesh/fvMeshAdder/fvMeshAdder.C | 4 + .../fvMeshDistribute/fvMeshDistribute.C | 68 ++- .../layerAdditionRemoval/addCellLayer.C | 89 +-- .../layerAdditionRemoval.C | 62 +- .../layerAdditionRemoval.H | 3 + .../layerAdditionRemoval/removeCellLayer.C | 2 +- .../layerAdditionRemoval/setLayerPairing.C | 3 + .../motionSmoother/motionSmoother.C | 37 +- .../motionSmoother/motionSmoother.H | 47 +- .../motionSmoother/motionSmootherCheck.C | 55 ++ .../motionSmoother/motionSmootherTemplates.C | 2 +- .../polyMeshGeometry/polyMeshGeometry.C | 220 ++++++- .../polyMeshGeometry/polyMeshGeometry.H | 10 + src/dynamicMesh/motionSolver/motionSolver.C | 2 +- src/dynamicMesh/motionSolver/motionSolver.H | 3 +- src/dynamicMesh/polyMeshAdder/polyMeshAdder.C | 1 + .../polyTopoChange/addPatchCellLayer.C | 374 ++++++++---- .../polyTopoChange/addPatchCellLayer.H | 104 +++- .../polyTopoChange/combineFaces.C | 125 ++-- .../polyTopoChange/combineFaces.H | 31 +- .../polyTopoChange/polyTopoChange/hexRef8.C | 382 ++++++++++-- .../polyTopoChange/polyTopoChange/hexRef8.H | 41 +- src/meshTools/Make/files | 6 + src/meshTools/PointEdgeWave/PointEdgeWave.C | 270 +++++---- src/meshTools/PointEdgeWave/PointEdgeWave.H | 7 +- .../edgeFaceCirculator/edgeFaceCirculator.C | 42 ++ .../edgeFaceCirculator/edgeFaceCirculator.H | 211 +++++++ .../edgeFaceCirculator/edgeFaceCirculatorI.H | 428 ++++++++++++++ src/meshTools/indexedOctree/indexedOctree.C | 280 ++++----- src/meshTools/indexedOctree/indexedOctree.H | 4 +- src/meshTools/indexedOctree/treeDataCell.C | 250 ++++++++ src/meshTools/indexedOctree/treeDataCell.H | 202 +++++++ src/meshTools/indexedOctree/treeDataEdge.C | 209 +++++++ src/meshTools/indexedOctree/treeDataEdge.H | 191 ++++++ src/meshTools/indexedOctree/treeDataFace.C | 545 ++++++++++++++++++ src/meshTools/indexedOctree/treeDataFace.H | 210 +++++++ src/meshTools/indexedOctree/treeDataPoint.C | 160 +++++ src/meshTools/indexedOctree/treeDataPoint.H | 162 ++++++ .../indexedOctree/treeDataTriSurface.C | 27 +- src/meshTools/octree/treeBoundBox.C | 177 +++--- src/meshTools/octree/treeBoundBox.H | 125 +++- src/meshTools/octree/treeBoundBoxI.H | 162 +++++- src/meshTools/regionSplit/regionSplit.C | 180 ++++-- src/meshTools/regionSplit/regionSplit.H | 21 + .../cellSources/regionToCell/regionToCell.C | 228 ++++++++ .../cellSources/regionToCell/regionToCell.H | 124 ++++ .../searchableSurface/searchableBox.C | 247 ++++++++ .../searchableSurface/searchableBox.H | 161 ++++++ .../searchableSurface/searchableSurface.C | 116 ++++ .../searchableSurface/searchableSurface.H | 257 +++++++++ .../searchableSurface/triSurfaceMesh.C | 271 +++++++++ .../searchableSurface/triSurfaceMesh.H | 194 +++++++ .../triSurfaceMeshes/triSurfaceMeshes.C | 144 ++--- .../triSurfaceMeshes/triSurfaceMeshes.H | 18 - .../triSurfaceSearch/triSurfaceSearch.C | 153 ++--- .../triSurfaceSearch/triSurfaceSearch.H | 13 +- .../triSurface/interfaces/STL/STLtriangle.H | 6 +- .../triSurface/interfaces/STL/STLtriangleI.H | 4 +- src/triSurface/triSurface/triSurface.C | 34 +- src/triSurface/triSurface/triSurface.H | 2 +- 98 files changed, 9511 insertions(+), 1997 deletions(-) create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeLagrangian.H create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.C create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.H create mode 100644 src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C create mode 100644 src/meshTools/edgeFaceCirculator/edgeFaceCirculator.C create mode 100644 src/meshTools/edgeFaceCirculator/edgeFaceCirculator.H create mode 100644 src/meshTools/edgeFaceCirculator/edgeFaceCirculatorI.H create mode 100644 src/meshTools/indexedOctree/treeDataCell.C create mode 100644 src/meshTools/indexedOctree/treeDataCell.H create mode 100644 src/meshTools/indexedOctree/treeDataEdge.C create mode 100644 src/meshTools/indexedOctree/treeDataEdge.H create mode 100644 src/meshTools/indexedOctree/treeDataFace.C create mode 100644 src/meshTools/indexedOctree/treeDataFace.H create mode 100644 src/meshTools/indexedOctree/treeDataPoint.C create mode 100644 src/meshTools/indexedOctree/treeDataPoint.H create mode 100644 src/meshTools/sets/cellSources/regionToCell/regionToCell.C create mode 100644 src/meshTools/sets/cellSources/regionToCell/regionToCell.H create mode 100644 src/meshTools/triSurface/searchableSurface/searchableBox.C create mode 100644 src/meshTools/triSurface/searchableSurface/searchableBox.H create mode 100644 src/meshTools/triSurface/searchableSurface/searchableSurface.C create mode 100644 src/meshTools/triSurface/searchableSurface/searchableSurface.H create mode 100644 src/meshTools/triSurface/searchableSurface/triSurfaceMesh.C create mode 100644 src/meshTools/triSurface/searchableSurface/triSurfaceMesh.H diff --git a/src/Allwmake b/src/Allwmake index 29dfe503c3..1a6c80560b 100755 --- a/src/Allwmake +++ b/src/Allwmake @@ -1,7 +1,7 @@ #!/bin/sh set -x -( cd $FOAM_SRC/other && ./Allwmake ) +#( cd $FOAM_SRC/other && ./Allwmake ) ( cd $FOAM_SRC/OpenFOAM && wmakeLnInclude . ) diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index a42b49e04b..fbb0fecdc6 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -400,7 +400,8 @@ mapPolyMesh = $(polyMesh)/mapPolyMesh $(mapPolyMesh)/mapPolyMesh.C $(mapPolyMesh)/faceMapper/faceMapper.C $(mapPolyMesh)/cellMapper/cellMapper.C -$(mapPolyMesh)/mapDistributePolyMesh/mapDistributePolyMesh.C +$(mapPolyMesh)/mapDistribute/mapDistribute.C +$(mapPolyMesh)/mapDistribute/mapDistributePolyMesh.C $(mapPolyMesh)/mapAddedPolyMesh.C PrimitivePatch = $(primitiveMesh)/PrimitivePatch diff --git a/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C b/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C index be4298a426..41549df18f 100644 --- a/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C +++ b/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.C @@ -527,7 +527,7 @@ void Foam::FaceCellWave::handleProcPatches() { const polyPatch& patch = mesh_.boundaryMesh()[patchI]; - if (Pstream::parRun() && isA(patch)) + if (isA(patch)) { // Allocate buffers label nSendFaces; @@ -580,7 +580,7 @@ void Foam::FaceCellWave::handleProcPatches() { const polyPatch& patch = mesh_.boundaryMesh()[patchI]; - if (Pstream::parRun() && isA(patch)) + if (isA(patch)) { const processorPolyPatch& procPatch = refCast(patch); @@ -812,7 +812,6 @@ Foam::FaceCellWave::FaceCellWave changedCells_(mesh_.nCells()), nChangedCells_(0), hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)), - hasProcPatches_(hasPatchType(processorPolyPatch::typeName)), nEvals_(0), nUnvisitedCells_(mesh_.nCells()), nUnvisitedFaces_(mesh_.nFaces()), @@ -843,7 +842,6 @@ Foam::FaceCellWave::FaceCellWave changedCells_(mesh_.nCells()), nChangedCells_(0), hasCyclicPatches_(hasPatchType(cyclicPolyPatch::typeName)), - hasProcPatches_(hasPatchType(processorPolyPatch::typeName)), nEvals_(0), nUnvisitedCells_(mesh_.nCells()), nUnvisitedFaces_(mesh_.nFaces()), @@ -1031,7 +1029,7 @@ Foam::label Foam::FaceCellWave::cellToFace() // Transfer changed faces across cyclic halves handleCyclicPatches(); } - if (hasProcPatches_) + if (Pstream::parRun()) { // Transfer changed faces from neighbouring processors. handleProcPatches(); @@ -1060,7 +1058,7 @@ Foam::label Foam::FaceCellWave::iterate(const label maxIter) // Transfer changed faces across cyclic halves handleCyclicPatches(); } - if (hasProcPatches_) + if (Pstream::parRun()) { // Transfer changed faces from neighbouring processors. handleProcPatches(); diff --git a/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H b/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H index 2b404ee0d3..a3522b6c94 100644 --- a/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H +++ b/src/OpenFOAM/algorithms/MeshWave/FaceCellWave.H @@ -32,8 +32,7 @@ Description Handles parallel and cyclics and non-parallel cyclics. -Note - Whether to propagate depends on the return value of Type::update + Note: whether to propagate depends on the return value of Type::update which returns true (i.e. propagate) if the value changes by more than a certain tolerance. This tolerance can be very strict for normal face-cell and parallel @@ -106,9 +105,6 @@ class FaceCellWave //- Contains cyclics bool hasCyclicPatches_; - //- Contains processor patches - bool hasProcPatches_; - //- Number of evaluations label nEvals_; diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedList.H b/src/OpenFOAM/containers/Lists/PackedList/PackedList.H index 6a807653ce..4bec222c55 100644 --- a/src/OpenFOAM/containers/Lists/PackedList/PackedList.H +++ b/src/OpenFOAM/containers/Lists/PackedList/PackedList.H @@ -164,9 +164,8 @@ public: //- Get value at index I inline unsigned int get(const label i) const; - - //- Set value at index I - inline void set(const label i, const unsigned int val); + //- Set value at index I. Return true if value changed. + inline bool set(const label i, const unsigned int val); // Member operators diff --git a/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H b/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H index c06567edae..59f00fa5c2 100644 --- a/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H +++ b/src/OpenFOAM/containers/Lists/PackedList/PackedListI.H @@ -151,7 +151,7 @@ inline unsigned int PackedList::operator[](const label i) const // Set value at i template -inline void PackedList::set(const label i, const unsigned int val) +inline bool PackedList::set(const label i, const unsigned int val) { # ifdef DEBUGList checkIndex(i); @@ -175,7 +175,11 @@ inline void PackedList::set(const label i, const unsigned int val) unsigned int& elem = List::operator[](intIndex(i)); + unsigned int oldElem = elem; + elem = (elem & ~shiftedMask) | shiftedVal; + + return elem != oldElem; } diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H b/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H index 9c5b37dc5e..e158ccff24 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H +++ b/src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H @@ -438,6 +438,38 @@ public: template static void listCombineScatter(List& Value); + // Combine variants working on whole map at a time. Container needs to + // have iterators and find() defined. + + template + static void mapCombineGather + ( + const List& comms, + Container& Values, + const CombineOp& cop + ); + + //- Like above but switches between linear/tree communication + template + static void mapCombineGather + ( + Container& Values, + const CombineOp& cop + ); + + //- Scatter data. Reverse of combineGather + template + static void mapCombineScatter + ( + const List& comms, + Container& Values + ); + + //- Like above but switches between linear/tree communication + template + static void mapCombineScatter(Container& Values); + + // Gather/scatter keeping the individual processor data separate. // Values is a List of size Pstream::nProcs() where diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/combineGatherScatter.C b/src/OpenFOAM/db/IOstreams/Pstreams/combineGatherScatter.C index fb3b56dcb7..d8d303513a 100644 --- a/src/OpenFOAM/db/IOstreams/Pstreams/combineGatherScatter.C +++ b/src/OpenFOAM/db/IOstreams/Pstreams/combineGatherScatter.C @@ -409,6 +409,147 @@ void Pstream::listCombineScatter(List& Values) +// Same thing but for sparse list (map) +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + +template +void Pstream::mapCombineGather +( + const List& comms, + Container& Values, + const CombineOp& cop +) +{ + if (Pstream::parRun()) + { + // Get my communication order + const commsStruct& myComm = comms[Pstream::myProcNo()]; + + // Receive from my downstairs neighbours + forAll(myComm.below(), belowI) + { + label belowID = myComm.below()[belowI]; + + IPstream fromBelow(Pstream::scheduled, belowID); + Container receivedValues(fromBelow); + + if (debug & 2) + { + Pout<< " received from " + << belowID << " data:" << receivedValues << endl; + } + + for + ( + typename Container::const_iterator slaveIter = + receivedValues.begin(); + slaveIter != receivedValues.end(); + ++slaveIter + ) + { + typename Container::iterator + masterIter = Values.find(slaveIter.key()); + + if (masterIter != Values.end()) + { + cop(masterIter(), slaveIter()); + } + else + { + Values.insert(slaveIter.key(), slaveIter()); + } + } + } + + // Send up Value + if (myComm.above() != -1) + { + if (debug & 2) + { + Pout<< " sending to " << myComm.above() + << " data:" << Values << endl; + } + + OPstream toAbove(Pstream::scheduled, myComm.above()); + toAbove << Values; + } + } +} + + +template +void Pstream::mapCombineGather(Container& Values, const CombineOp& cop) +{ + if (Pstream::nProcs() < Pstream::nProcsSimpleSum) + { + mapCombineGather(Pstream::linearCommunication(), Values, cop); + } + else + { + mapCombineGather(Pstream::treeCommunication(), Values, cop); + } +} + + +template +void Pstream::mapCombineScatter +( + const List& comms, + Container& Values +) +{ + if (Pstream::parRun()) + { + // Get my communication order + const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()]; + + // Reveive from up + if (myComm.above() != -1) + { + IPstream fromAbove(Pstream::scheduled, myComm.above()); + fromAbove >> Values; + + if (debug & 2) + { + Pout<< " received from " + << myComm.above() << " data:" << Values << endl; + } + } + + // Send to my downstairs neighbours + forAll(myComm.below(), belowI) + { + label belowID = myComm.below()[belowI]; + + if (debug & 2) + { + Pout<< " sending to " << belowID << " data:" << Values << endl; + } + + OPstream toBelow(Pstream::scheduled, belowID); + toBelow << Values; + } + } +} + + +template +void Pstream::mapCombineScatter(Container& Values) +{ + if (Pstream::nProcs() < Pstream::nProcsSimpleSum) + { + mapCombineScatter(Pstream::linearCommunication(), Values); + } + else + { + mapCombineScatter(Pstream::treeCommunication(), Values); + } +} + + + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index 9e7b463fec..f56919647a 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/src/OpenFOAM/meshes/meshShapes/face/face.H @@ -188,7 +188,7 @@ public: //- Return potential intersection with face with a ray starting // at p, direction n (does not need to be normalized) // Does face-center decomposition and returns triangle intersection - // point closest to p. + // point closest to p. Face-center is calculated from point average. // For a hit, the distance is signed. Positive number // represents the point in front of triangle // In case of miss the point is the nearest point on the face @@ -206,6 +206,19 @@ public: const intersection::direction dir = intersection::VECTOR ) const; + //- Fast intersection with a ray. + // For a hit, the pointHit.distance() is the line parameter t : + // intersection=p+t*q. Only defined for FULL_RAY or + // HALF_RAY. + pointHit intersection + ( + const point& p, + const vector& q, + const point& ctr, + const pointField& meshPoints, + const intersection::algorithm alg + ) const; + //- Return nearest point to face pointHit nearestPoint ( diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C index 680183038d..78b9c6e82f 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C @@ -133,6 +133,52 @@ pointHit face::ray } +pointHit face::intersection +( + const point& p, + const vector& q, + const point& ctr, + const pointField& meshPoints, + const intersection::algorithm alg +) const +{ + scalar nearestHitDist = VGREAT; + + // Initialize to miss, distance = GREAT + pointHit nearest(p); + + const labelList& f = *this; + + forAll(f, pI) + { + // Note: for best accuracy, centre point always comes last + pointHit curHit = triPointRef + ( + meshPoints[f[pI]], + meshPoints[f[fcIndex(pI)]], + ctr + ).intersection(p, q, alg); + + if (curHit.hit()) + { + if (Foam::mag(curHit.distance()) < nearestHitDist) + { + nearestHitDist = Foam::mag(curHit.distance()); + nearest.setHit(); + nearest.setPoint(curHit.hitPoint()); + } + } + } + + if (nearest.hit()) + { + nearest.setDistance(nearestHitDist); + } + + return nearest; +} + + pointHit face::nearestPoint ( const point& p, diff --git a/src/OpenFOAM/meshes/meshTools/matchPoints.C b/src/OpenFOAM/meshes/meshTools/matchPoints.C index fd6fcc831d..e60160fb6e 100644 --- a/src/OpenFOAM/meshes/meshTools/matchPoints.C +++ b/src/OpenFOAM/meshes/meshTools/matchPoints.C @@ -64,9 +64,11 @@ bool Foam::matchPoints startI = 0; } - label face1I = -1; - // Go through range of equal mag and find equal vector. + // Go through range of equal mag and find nearest vector. + scalar minDistSqr = VGREAT; + label minFaceI = -1; + for ( label j = startI; @@ -78,17 +80,17 @@ bool Foam::matchPoints ) { label faceI = pts1MagSqr.indices()[j]; - // Compare actual vectors - if (magSqr(pts0[face0I] - pts1[faceI]) <= sqr(matchDist)) - { - face1I = faceI; + scalar distSqr = magSqr(pts0[face0I] - pts1[faceI]); - break; + if (distSqr <= sqr(matchDist) && distSqr < minDistSqr) + { + minDistSqr = distSqr; + minFaceI = faceI; } } - if (face1I == -1) + if (minFaceI == -1) { fullMatch = false; @@ -120,7 +122,7 @@ bool Foam::matchPoints } } - from0To1[face0I] = face1I; + from0To1[face0I] = minFaceI; } return fullMatch; diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C new file mode 100644 index 0000000000..277b18c78c --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C @@ -0,0 +1,197 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +\*---------------------------------------------------------------------------*/ + +#include "mapDistribute.H" +#include "commSchedule.H" +#include "HashSet.H" +#include "ListOps.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::mapDistribute::calcSchedule() const +{ + // Communications: send and receive processor + List allComms; + + { + HashSet > commsSet(Pstream::nProcs()); + + // Find what communication is required + forAll(subMap_, procI) + { + if (procI != Pstream::myProcNo()) + { + if (subMap_[procI].size() > 0) + { + // I need to send to procI + commsSet.insert(labelPair(Pstream::myProcNo(), procI)); + } + if (constructMap_[procI].size() > 0) + { + // I need to receive from procI + commsSet.insert(labelPair(procI, Pstream::myProcNo())); + } + } + } + allComms = commsSet.toc(); + } + + + // Reduce + if (Pstream::master()) + { + // Receive and merge + for + ( + int slave=Pstream::firstSlave(); + slave<=Pstream::lastSlave(); + slave++ + ) + { + IPstream fromSlave(Pstream::blocking, slave); + List nbrData(fromSlave); + + forAll(nbrData, i) + { + if (findIndex(allComms, nbrData[i]) == -1) + { + label sz = allComms.size(); + allComms.setSize(sz+1); + allComms[sz] = nbrData[i]; + } + } + } + // Send back + for + ( + int slave=Pstream::firstSlave(); + slave<=Pstream::lastSlave(); + slave++ + ) + { + OPstream toSlave(Pstream::blocking, slave); + toSlave << allComms; + } + } + else + { + { + OPstream toMaster(Pstream::blocking, Pstream::masterNo()); + toMaster << allComms; + } + { + IPstream fromMaster(Pstream::blocking, Pstream::masterNo()); + fromMaster >> allComms; + } + } + + + // Determine my schedule. + labelList mySchedule + ( + commSchedule + ( + Pstream::nProcs(), + allComms + ).procSchedule()[Pstream::myProcNo()] + ); + + // Processors involved in my schedule + schedulePtr_.reset + ( + new List + ( + IndirectList(allComms, mySchedule) + ) + ); + + + //if (debug) + //{ + // Pout<< "I need to:" << endl; + // const List& comms = schedule(); + // forAll(comms, i) + // { + // const labelPair& twoProcs = comms[i]; + // label sendProc = twoProcs[0]; + // label recvProc = twoProcs[1]; + // + // if (recvProc == Pstream::myProcNo()) + // { + // Pout<< " receive from " << sendProc << endl; + // } + // else + // { + // Pout<< " send to " << recvProc << endl; + // } + // } + //} +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +//- Construct from components +Foam::mapDistribute::mapDistribute +( + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap +) +: + constructSize_(constructSize), + subMap_(subMap), + constructMap_(constructMap), + schedulePtr_() +{} + + +//- (optionally destructively) construct from components +Foam::mapDistribute::mapDistribute +( + const label constructSize, + labelListList& subMap, + labelListList& constructMap, + const bool reUse // clone or reuse +) +: + constructSize_(constructSize), + subMap_(subMap, reUse), + constructMap_(constructMap, reUse), + schedulePtr_() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H new file mode 100644 index 0000000000..6af058a88c --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H @@ -0,0 +1,203 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +Class + Foam::mapDistribute + +Description + Class containing processor-to-processor mapping information. + + We store mapping from the bits-to-send to the complete starting list + (subXXXMap) and from the received bits to their location in the new + list (constructXXXMap). + +Note: + Schedule is a list of processor pairs (one send, one receive. One of + them will be myself) which forms a scheduled (i.e. non-buffered) exchange. + See distribute on how to use it. + + +SourceFiles + mapDistribute.C + +\*---------------------------------------------------------------------------*/ + +#ifndef mapDistribute_H +#define mapDistribute_H + +#include "labelList.H" +#include "labelPair.H" +#include "Pstream.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class mapPolyMesh; + +/*---------------------------------------------------------------------------*\ + Class mapDistribute Declaration +\*---------------------------------------------------------------------------*/ + +class mapDistribute +{ + // Private data + + //- Size of reconstructed data + label constructSize_; + + //- Maps from subsetted data back to original data + labelListList subMap_; + + //- Maps from subsetted data to new reconstructed data + labelListList constructMap_; + + //- Schedule + mutable autoPtr > schedulePtr_; + + + // Private Member Functions + + void calcSchedule() const; + + //- Disallow default bitwise copy construct + mapDistribute(const mapDistribute&); + + //- Disallow default bitwise assignment + void operator=(const mapDistribute&); + + + +public: + + // Constructors + + //- Construct from components + mapDistribute + ( + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap + ); + + //- (optionally destructively) construct from components + mapDistribute + ( + const label constructSize, + labelListList& subMap, + labelListList& constructMap, + const bool reUse // clone or reuse + ); + + + // Member Functions + + // Access + + //- Constructed data size + label constructSize() const + { + return constructSize_; + } + + //- From subsetted data back to original data + const labelListList& subMap() const + { + return subMap_; + } + + //- From subsetted data to new reconstructed data + const labelListList& constructMap() const + { + return constructMap_; + } + + //- Return a schedule. Demand driven. See above. + const List& schedule() const + { + if (!schedulePtr_.valid()) + { + calcSchedule(); + } + return schedulePtr_(); + } + + + // Other + + //- Distribute data. Note:schedule only used for Pstream::scheduled + // for now, all others just use send-to-all, receive-from-all. + template + static void distribute + ( + const Pstream::commsTypes commsType, + const List& schedule, + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap, + List& + ); + + //- Distribute data using scheduling. + template + void distribute(List& fld) const + { + distribute + ( + Pstream::scheduled, + schedule(), + constructSize_, + subMap_, + constructMap_, + fld + ); + } + + //- Correct for topo change. + void updateMesh(const mapPolyMesh&) + { + notImplemented + ( + "mapDistribute::updateMesh(const mapPolyMesh&)" + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "mapDistributeTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeLagrangian.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeLagrangian.H new file mode 100644 index 0000000000..3ed51b365f --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeLagrangian.H @@ -0,0 +1,146 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +Class + Foam::mapDistributeLagrangian + +Description + Class containing mesh-to-mesh mapping information for particles + +SourceFiles + mapDistributeLagrangian.C + +\*---------------------------------------------------------------------------*/ + +#ifndef mapDistributeLagrangian_H +#define mapDistributeLagrangian_H + +#include "mapDistribute.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class mapPolyMesh; + +/*---------------------------------------------------------------------------*\ + Class mapDistributeLagrangian Declaration +\*---------------------------------------------------------------------------*/ + +class mapDistributeLagrangian +{ + // Private data + + //- Map to distribute particles + const mapDistribute particleMap_; + + //- Per element in subsetted mesh the cell label + const labelListList constructCellLabels_; + + +public: + + // Constructors + + //- Construct from components + mapDistributeLagrangian + ( + const label nNewParticles, + const labelListList& subParticleMap, + const labelListList& constructParticleMap, + const labelListList& constructCellLabels + ) + : + particleMap_(nNewParticles, subParticleMap, constructParticleMap), + constructCellLabels_(constructCellLabels) + {} + + //- Construct from components and steal storage + mapDistributeLagrangian + ( + const label nNewParticles, + labelListList& subParticleMap, + labelListList& constructParticleMap, + labelListList& constructCellLabels, + const bool reUse + ) + : + particleMap_ + ( + nNewParticles, + subParticleMap, + constructParticleMap, + reUse + ), + constructCellLabels_(constructCellLabels, reUse) + {} + + + // Member Functions + + // Access + + //- Distribution map + const mapDistribute& particleMap() const + { + return particleMap_; + } + + //- Per received particle the destination cell label + const labelListList& constructCellLabels() const + { + return constructCellLabels_; + } + + + // Edit + + //- distribute list of lagrangian data + template + void distributeLagrangianData(List& lst) const + { + particleMap_.distribute(lst); + } + + //- Correct for topo change. + void updateMesh(const mapPolyMesh&) + { + notImplemented + ( + "mapDistributeLagrangian::updateMesh(const mapPolyMesh&)" + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.C new file mode 100644 index 0000000000..ae733c338c --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.C @@ -0,0 +1,235 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +\*---------------------------------------------------------------------------*/ + +#include "mapDistributePolyMesh.H" +#include "polyMesh.H" + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::mapDistributePolyMesh::calcPatchSizes() +{ + oldPatchSizes_.setSize(oldPatchStarts_.size()); + + // Calculate old patch sizes + for (label patchI = 0; patchI < oldPatchStarts_.size() - 1; patchI++) + { + oldPatchSizes_[patchI] = + oldPatchStarts_[patchI + 1] - oldPatchStarts_[patchI]; + } + + // Set the last one by hand + const label lastPatchID = oldPatchStarts_.size() - 1; + + oldPatchSizes_[lastPatchID] = nOldFaces_ - oldPatchStarts_[lastPatchID]; + + if (min(oldPatchSizes_) < 0) + { + FatalErrorIn("mapDistributePolyMesh::calcPatchSizes()") + << "Calculated negative old patch size:" << oldPatchSizes_ << nl + << "Error in mapping data" << abort(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +//- Construct from components +Foam::mapDistributePolyMesh::mapDistributePolyMesh +( + const polyMesh& mesh, + + // mesh before changes + const label nOldPoints, + const label nOldFaces, + const label nOldCells, + const labelList& oldPatchStarts, + const labelList& oldPatchNMeshPoints, + + // how to subset pieces of mesh to send across + const labelListList& subPointMap, + const labelListList& subFaceMap, + const labelListList& subCellMap, + const labelListList& subPatchMap, + + // how to reconstruct received mesh + const labelListList& constructPointMap, + const labelListList& constructFaceMap, + const labelListList& constructCellMap, + const labelListList& constructPatchMap +) +: + mesh_(mesh), + nOldPoints_(nOldPoints), + nOldFaces_(nOldFaces), + nOldCells_(nOldCells), + oldPatchSizes_(oldPatchStarts.size()), + oldPatchStarts_(oldPatchStarts), + oldPatchNMeshPoints_(oldPatchNMeshPoints), + pointMap_(mesh.nPoints(), subPointMap, constructPointMap), + faceMap_(mesh.nFaces(), subFaceMap, constructFaceMap), + cellMap_(mesh.nCells(), subCellMap, constructCellMap), + patchMap_(mesh.boundaryMesh().size(), subPatchMap, constructPatchMap) +{ + calcPatchSizes(); +} + + +//- (optionally destructively) construct from components +Foam::mapDistributePolyMesh::mapDistributePolyMesh +( + const polyMesh& mesh, + const label nOldPoints, + const label nOldFaces, + const label nOldCells, + labelList& oldPatchStarts, + labelList& oldPatchNMeshPoints, + + labelListList& subPointMap, + labelListList& subFaceMap, + labelListList& subCellMap, + labelListList& subPatchMap, + labelListList& constructPointMap, + labelListList& constructFaceMap, + labelListList& constructCellMap, + labelListList& constructPatchMap, + const bool reUse // clone or reuse +) +: + mesh_(mesh), + nOldPoints_(nOldPoints), + nOldFaces_(nOldFaces), + nOldCells_(nOldCells), + oldPatchSizes_(oldPatchStarts.size()), + oldPatchStarts_(oldPatchStarts, reUse), + oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse), + + pointMap_(mesh.nPoints(), subPointMap, constructPointMap, reUse), + faceMap_(mesh.nFaces(), subFaceMap, constructFaceMap, reUse), + cellMap_(mesh.nCells(), subCellMap, constructCellMap, reUse), + patchMap_(mesh.boundaryMesh().size(), subPatchMap, constructPatchMap, reUse) +{ + calcPatchSizes(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::mapDistributePolyMesh::distributePointIndices(labelList& lst) const +{ + // Construct boolList from selected elements + boolList isSelected + ( + createWithValues + ( + nOldPoints(), + false, + lst, + true + ) + ); + + // Distribute + distributePointData(isSelected); + + // Collect selected elements + lst = findIndices(isSelected, true); +} + + +void Foam::mapDistributePolyMesh::distributeFaceIndices(labelList& lst) const +{ + // Construct boolList from selected elements + boolList isSelected + ( + createWithValues + ( + nOldFaces(), + false, + lst, + true + ) + ); + + // Distribute + distributeFaceData(isSelected); + + // Collect selected elements + lst = findIndices(isSelected, true); +} + + +void Foam::mapDistributePolyMesh::distributeCellIndices(labelList& lst) const +{ + // Construct boolList from selected elements + boolList isSelected + ( + createWithValues + ( + nOldCells(), + false, + lst, + true + ) + ); + + // Distribute + distributeCellData(isSelected); + + // Collect selected elements + lst = findIndices(isSelected, true); +} + + +void Foam::mapDistributePolyMesh::distributePatchIndices(labelList& lst) const +{ + // Construct boolList from selected elements + boolList isSelected + ( + createWithValues + ( + oldPatchStarts().size(), // nOldPatches + false, + lst, + true + ) + ); + + // Distribute + distributePatchData(isSelected); + + // Collect selected elements + lst = findIndices(isSelected, true); +} + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.H new file mode 100644 index 0000000000..ad48d21c0e --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributePolyMesh.H @@ -0,0 +1,291 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +Class + Foam::mapDistributePolyMesh + +Description + Class containing mesh-to-mesh mapping information after a mesh distribution + where we send parts of meshes (using subsetting) to other processors + and receive and reconstruct mesh. + + We store mapping from the bits-to-send to the complete starting mesh + (subXXXMap) and from the received bits to their location in the new + mesh (constructXXXMap). + +SourceFiles + mapDistributePolyMesh.C + +\*---------------------------------------------------------------------------*/ + +#ifndef mapDistributePolyMesh_H +#define mapDistributePolyMesh_H + +#include "mapDistribute.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class mapPolyMesh; +class polyMesh; + +/*---------------------------------------------------------------------------*\ + Class mapDistributePolyMesh Declaration +\*---------------------------------------------------------------------------*/ + +class mapDistributePolyMesh +{ + // Private data + + const polyMesh& mesh_; + + //- Number of old live points + const label nOldPoints_; + + //- Number of old live faces + const label nOldFaces_; + + //- Number of old live cells + const label nOldCells_; + + //- List of the old patch sizes + labelList oldPatchSizes_; + + //- List of the old patch start labels + const labelList oldPatchStarts_; + + //- List of numbers of mesh points per old patch + const labelList oldPatchNMeshPoints_; + + + //- Point distribute map + const mapDistribute pointMap_; + + //- Face distribute map + const mapDistribute faceMap_; + + //- Cell distribute map + const mapDistribute cellMap_; + + //- Patch distribute map + const mapDistribute patchMap_; + + + + // Private Member Functions + + void calcPatchSizes(); + + //- Disallow default bitwise copy construct + mapDistributePolyMesh(const mapDistributePolyMesh&); + + //- Disallow default bitwise assignment + void operator=(const mapDistributePolyMesh&); + + +public: + + // Constructors + + //- Construct from components. Note that mesh has to be changed already + // since uses mesh.nPoints etc as the new size. + mapDistributePolyMesh + ( + const polyMesh& mesh, + + // mesh before changes + const label nOldPoints, + const label nOldFaces, + const label nOldCells, + const labelList& oldPatchStarts, + const labelList& oldPatchNMeshPoints, + + // how to subset pieces of mesh to send across + const labelListList& subPointMap, + const labelListList& subFaceMap, + const labelListList& subCellMap, + const labelListList& subPatchMap, + + // how to reconstruct received mesh + const labelListList& constructPointMap, + const labelListList& constructFaceMap, + const labelListList& constructCellMap, + const labelListList& constructPatchMap + ); + + //- (optionally destructively) construct from components + // Note that mesh has to be changed already! + mapDistributePolyMesh + ( + const polyMesh& mesh, + const label nOldPoints, + const label nOldFaces, + const label nOldCells, + labelList& oldPatchStarts, + labelList& oldPatchNMeshPoints, + + labelListList& subPointMap, + labelListList& subFaceMap, + labelListList& subCellMap, + labelListList& subPatchMap, + labelListList& constructPointMap, + labelListList& constructFaceMap, + labelListList& constructCellMap, + labelListList& constructPatchMap, + const bool reUse // clone or reuse + ); + + + // Member Functions + + // Access + + const polyMesh& mesh() const + { + return mesh_; + } + + //- Number of points in mesh before distribution + label nOldPoints() const + { + return nOldPoints_; + } + + //- Number of faces in mesh before distribution + label nOldFaces() const + { + return nOldFaces_; + } + + //- Number of cells in mesh before distribution + label nOldCells() const + { + return nOldCells_; + } + + //- List of the old patch sizes + const labelList& oldPatchSizes() const + { + return oldPatchSizes_; + } + + //- List of the old patch start labels + const labelList& oldPatchStarts() const + { + return oldPatchStarts_; + } + + //- List of numbers of mesh points per old patch + const labelList& oldPatchNMeshPoints() const + { + return oldPatchNMeshPoints_; + } + + //- Point distribute map + const mapDistribute& pointMap() const + { + return pointMap_; + } + + //- Face distribute map + const mapDistribute& faceMap() const + { + return faceMap_; + } + + //- Cell distribute map + const mapDistribute& cellMap() const + { + return cellMap_; + } + + //- Patch distribute map + const mapDistribute& patchMap() const + { + return patchMap_; + } + + + // Edit + + //- distribute list of point data + template + void distributePointData(List& lst) const + { + pointMap_.distribute(lst); + } + + //- distribute list of face data + template + void distributeFaceData(List& lst) const + { + faceMap_.distribute(lst); + } + + //- distribute list of cell data + template + void distributeCellData(List& lst) const + { + cellMap_.distribute(lst); + } + + //- distribute list of patch data + template + void distributePatchData(List& lst) const + { + patchMap_.distribute(lst); + } + + + //- distribute list of point/face/cell/patch indices. + // (Converts to boolList, distributes boolList and reconstructs) + void distributePointIndices(labelList& pointIDs) const; + + void distributeFaceIndices(labelList& faceIDs) const; + void distributeCellIndices(labelList& cellIDs) const; + void distributePatchIndices(labelList& patchIDs) const; + + + //- Correct for topo change. + void updateMesh(const mapPolyMesh&) + { + notImplemented + ( + "mapDistributePolyMesh::updateMesh(const mapPolyMesh&)" + ); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C new file mode 100644 index 0000000000..9e4cd48054 --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C @@ -0,0 +1,190 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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 + +\*---------------------------------------------------------------------------*/ + +#include "Pstream.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Distribute list. +template +void Foam::mapDistribute::distribute +( + const Pstream::commsTypes commsType, + const List& schedule, + const label constructSize, + const labelListList& subMap, + const labelListList& constructMap, + List& field +) +{ + if (commsType == Pstream::blocking) + { + // Since buffered sending can reuse the field to collect the + // received data. + + // Send sub field to neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + if (domain != Pstream::myProcNo()) + { + OPstream toNbr(Pstream::blocking, domain); + toNbr << IndirectList(field, subMap[domain])(); + } + } + + // Subset myself + List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + + // Receive sub field from myself (subField) + const labelList& map = constructMap[Pstream::myProcNo()]; + + field.setSize(constructSize); + + forAll(map, i) + { + field[map[i]] = subField[i]; + } + + // Receive sub field from neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + if (domain != Pstream::myProcNo()) + { + IPstream fromNbr(Pstream::blocking, domain); + List subField(fromNbr); + + const labelList& map = constructMap[domain]; + + forAll(map, i) + { + field[map[i]] = subField[i]; + } + } + } + } + else if (commsType == Pstream::scheduled) + { + // Need to make sure I don't overwrite field with received data + // since the data might need to be sent to another processor. So + // allocate a new field for the results. + List newField(constructSize); + + // Subset myself + List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + + // Receive sub field from myself (subField) + const labelList& map = constructMap[Pstream::myProcNo()]; + + forAll(map, i) + { + newField[map[i]] = subField[i]; + } + + forAll(schedule, i) + { + const labelPair& twoProcs = schedule[i]; + label sendProc = twoProcs[0]; + label recvProc = twoProcs[1]; + + if (Pstream::myProcNo() == sendProc) + { + // I am sender. Send to recvProc. + OPstream toNbr(Pstream::scheduled, recvProc); + toNbr << IndirectList(field, subMap[recvProc])(); + } + else + { + // I am receiver. Receive from sendProc. + IPstream fromNbr(Pstream::scheduled, sendProc); + List subField(fromNbr); + + const labelList& map = constructMap[sendProc]; + + forAll(map, i) + { + newField[map[i]] = subField[i]; + } + } + } + field.transfer(newField); + } + else if (commsType == Pstream::nonBlocking) + { + List newField(constructSize); + + // Subset myself + List subField(IndirectList(field, subMap[Pstream::myProcNo()])); + + // Receive sub field from myself (subField) + const labelList& map = constructMap[Pstream::myProcNo()]; + + forAll(map, i) + { + newField[map[i]] = subField[i]; + } + + // Send sub field to neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + if (domain != Pstream::myProcNo()) + { + OPstream toNbr(Pstream::nonBlocking, domain); + toNbr << IndirectList(field, subMap[domain])(); + } + } + + + // Receive sub field from neighbour + for (label domain = 0; domain < Pstream::nProcs(); domain++) + { + if (domain != Pstream::myProcNo()) + { + IPstream fromNbr(Pstream::nonBlocking, domain); + List subField(fromNbr); + + const labelList& map = constructMap[domain]; + + forAll(map, i) + { + newField[map[i]] = subField[i]; + } + } + } + OPstream::waitRequests(); + IPstream::waitRequests(); + + field.transfer(newField); + } + else + { + FatalErrorIn("mapDistribute::distribute(..)") + << "Unknown communication schedule " << commsType + << abort(FatalError); + } +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C index a8e1b59711..40fb615e3c 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C @@ -37,7 +37,7 @@ namespace Foam defineTypeNameAndDebug(coupledPolyPatch, 0); } -Foam::scalar Foam::coupledPolyPatch::matchTol_ = 1E-3; +Foam::scalar Foam::coupledPolyPatch::matchTol = 1E-3; // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // @@ -83,7 +83,7 @@ void Foam::coupledPolyPatch::writeOBJ void Foam::coupledPolyPatch::writeOBJ ( const fileName& fName, - const faceList& faces, + const UList& faces, const pointField& points ) { @@ -118,7 +118,7 @@ void Foam::coupledPolyPatch::writeOBJ Foam::pointField Foam::coupledPolyPatch::calcFaceCentres ( - const faceList& faces, + const UList& faces, const pointField& points ) { @@ -135,7 +135,7 @@ Foam::pointField Foam::coupledPolyPatch::calcFaceCentres Foam::pointField Foam::coupledPolyPatch::getAnchorPoints ( - const faceList& faces, + const UList& faces, const pointField& points ) { @@ -189,7 +189,7 @@ Foam::label Foam::coupledPolyPatch::whichPatch Foam::scalarField Foam::coupledPolyPatch::calcFaceTol ( - const faceList& faces, + const UList& faces, const pointField& points, const pointField& faceCentres ) @@ -203,13 +203,13 @@ Foam::scalarField Foam::coupledPolyPatch::calcFaceTol const face& f = faces[faceI]; - scalar maxLen = -GREAT; + scalar maxLenSqr = -GREAT; forAll(f, fp) { - maxLen = max(maxLen, mag(points[f[fp]] - cc)); + maxLenSqr = max(maxLenSqr, magSqr(points[f[fp]] - cc)); } - tols[faceI] = matchTol_ * maxLen; + tols[faceI] = matchTol * Foam::sqrt(maxLenSqr); } return tols; } @@ -249,71 +249,32 @@ Foam::label Foam::coupledPolyPatch::getRotation } -void Foam::coupledPolyPatch::calcTransformTensors -( - const vector& Cf, - const vector& Cr, - const vector& nf, - const vector& nr -) const -{ - if (debug) - { - Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl - << " nf:" << nf << nl - << " nr:" << nr << nl - << " mag(nf&nr):" << mag(nf & nr) << nl - << " Cf:" << Cf << nl - << " Cr:" << Cr << endl; - } - - if (mag(nf & nr) < 1 - SMALL) - { - separation_.setSize(0); - - forwardT_ = tensorField(1, rotationTensor(-nr, nf)); - reverseT_ = tensorField(1, rotationTensor(nf, -nr)); - } - else - { - forwardT_.setSize(0); - reverseT_.setSize(0); - - vector separation = (nf & (Cr - Cf))*nf; - - if (mag(separation) > SMALL) - { - separation_ = vectorField(1, separation); - } - else - { - separation_.setSize(0); - } - } - - if (debug) - { - Pout<< " separation_:" << separation_ << nl - << " forwardT size:" << forwardT_.size() << endl; - } -} - - void Foam::coupledPolyPatch::calcTransformTensors ( const vectorField& Cf, const vectorField& Cr, const vectorField& nf, - const vectorField& nr + const vectorField& nr, + const scalarField& smallDist, + const scalar absTol ) const { if (debug) { Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl - << " size:" << size() << nl + << " (half)size:" << Cf.size() << nl + << " absTol:" << absTol << nl << " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl; } + // Tolerance calculation. + // - normal calculation: assume absTol is the absolute error in a + // single normal/transformation calculation. Consists both of numerical + // precision (on the order of SMALL and of writing precision + // (from e.g. decomposition) + // Then the overall error of summing the normals is sqrt(size())*absTol + // - separation calculation: pass in from the outside an allowable error. + if (size() == 0) { // Dummy geometry. @@ -321,39 +282,104 @@ void Foam::coupledPolyPatch::calcTransformTensors forwardT_ = I; reverseT_ = I; } - else if (sum(mag(nf & nr)) < size() - 1E-12) - { - separation_.setSize(0); - - forwardT_.setSize(size()); - reverseT_.setSize(size()); - - forAll (forwardT_, facei) - { - forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]); - reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]); - } - - if (sum(mag(forwardT_ - forwardT_[0])) < 1E-12) - { - forwardT_.setSize(1); - reverseT_.setSize(1); - } - } else { - forwardT_.setSize(0); - reverseT_.setSize(0); + scalar error = absTol*Foam::sqrt(1.0*Cf.size()); - separation_ = (nf&(Cr - Cf))*nf; - - if (sum(mag(separation_))/size() < 1E-12) + if (debug) { - separation_.setSize(0); + Pout<< " error:" << error << endl; } - else if (sum(mag(separation_ - separation_[0]))/size() < 1E-12) + + if (sum(mag(nf & nr)) < Cf.size()-error) { - separation_.setSize(1); + // Rotation, no separation + + separation_.setSize(0); + + forwardT_.setSize(Cf.size()); + reverseT_.setSize(Cf.size()); + + forAll (forwardT_, facei) + { + forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]); + reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]); + } + + if (debug) + { + Pout<< " sum(mag(forwardT_ - forwardT_[0])):" + << sum(mag(forwardT_ - forwardT_[0])) + << endl; + } + + if (sum(mag(forwardT_ - forwardT_[0])) < error) + { + forwardT_.setSize(1); + reverseT_.setSize(1); + } + } + else + { + forwardT_.setSize(0); + reverseT_.setSize(0); + + separation_ = (nf&(Cr - Cf))*nf; + + // Three situations: + // - separation is zero. No separation. + // - separation is same. Single separation vector. + // - separation differs per face. Separation vectorField. + + // Check for different separation per face + bool sameSeparation = true; + + forAll(separation_, facei) + { + scalar smallSqr = sqr(smallDist[facei]); + + if (magSqr(separation_[facei] - separation_[0]) > smallSqr) + { + if (debug) + { + Pout<< " separation " << separation_[facei] + << " at " << facei + << " differs from separation[0] " << separation_[0] + << " by more than local tolerance " + << smallDist[facei] + << ". Assuming non-uniform separation." << endl; + } + sameSeparation = false; + break; + } + } + + if (sameSeparation) + { + // Check for zero separation (at 0 so everywhere) + if (magSqr(separation_[0]) < sqr(smallDist[0])) + { + if (debug) + { + Pout<< " separation " << mag(separation_[0]) + << " less than local tolerance " << smallDist[0] + << ". Assuming zero separation." << endl; + } + + separation_.setSize(0); + } + else + { + if (debug) + { + Pout<< " separation " << mag(separation_[0]) + << " more than local tolerance " << smallDist[0] + << ". Assuming uniform separation." << endl; + } + + separation_.setSize(1); + } + } } } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H index c5dc127382..5c37ea0528 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H @@ -23,7 +23,7 @@ License Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA Class - Foam::coupledPolyPatch + coupledPolyPatch Description coupledPolyPatch is an abstract base class for patches that couple regions @@ -64,34 +64,29 @@ class coupledPolyPatch //- Neighbour-cell transformation tensor mutable tensorField reverseT_; +public: - // Private static data + // Static data members - //- Relative tolerance (for geometric matching). Is factor of - // maximum edge length per face. - static scalar matchTol_; + //- Relative tolerance (for geometric matching). + static scalar matchTol; protected: // Protected Member Functions - //- Calculate the uniform transformation tensors - void calcTransformTensors - ( - const vector& Cf, - const vector& Cr, - const vector& nf, - const vector& nr - ) const; - //- Calculate the transformation tensors + // smallDist : matching distance per face + // absTol : absolute error in normal void calcTransformTensors ( const vectorField& Cf, const vectorField& Cr, const vectorField& nf, - const vectorField& nr + const vectorField& nr, + const scalarField& smallDist, + const scalar absTol = matchTol ) const; //- Initialise the calculation of the patch geometry @@ -123,7 +118,7 @@ protected: static void writeOBJ ( const fileName&, - const faceList&, + const UList&, const pointField& ); @@ -137,10 +132,18 @@ protected: ); //- Calculate face centres - static pointField calcFaceCentres(const faceList&, const pointField&); + static pointField calcFaceCentres + ( + const UList&, + const pointField& + ); //- Get f[0] for all faces - static pointField getAnchorPoints(const faceList&, const pointField&); + static pointField getAnchorPoints + ( + const UList&, + const pointField& + ); //- Is face (in old face labels) in current patch? bool inPatch @@ -161,7 +164,7 @@ protected: // from face centre to any of the face vertices. static scalarField calcFaceTol ( - const faceList& faces, + const UList& faces, const pointField& points, const pointField& faceCentres ); diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index fffe015533..21b03fdaa6 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -33,6 +33,7 @@ License #include "patchZones.H" #include "matchPoints.H" #include "EdgeMap.H" +#include "Time.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -51,23 +52,108 @@ void Foam::cyclicPolyPatch::calcTransforms() { if (size() > 0) { - const pointField& points = this->points(); + primitivePatch half0 + ( + SubList + ( + *this, + size()/2 + ), + points() + ); + pointField half0Ctrs(calcFaceCentres(half0, half0.points())); + scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs)); - const face& f0 = static_cast(*this)[0]; - const face& fn2 = static_cast(*this)[size()/2]; + primitivePatch half1 + ( + SubList + ( + *this, + size()/2, + size()/2 + ), + points() + ); + pointField half1Ctrs(calcFaceCentres(half1, half1.points())); - vector nf0 = f0.normal(points); - nf0 /= mag(nf0); + // Dump halves + if (debug) + { + fileName casePath(boundaryMesh().mesh().time().path()); - vector nfn2 = fn2.normal(points); - nfn2 /= mag(nfn2); + fileName nm0(casePath/name()+"_half0_faces.obj"); + Pout<< "cyclicPolyPatch::calcTransforms : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0, half0.points()); + fileName nm1(casePath/name()+"_half1_faces.obj"); + Pout<< "cyclicPolyPatch::calcTransforms : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1, half1.points()); + } + + vectorField half0Normals(half0.size()); + vectorField half1Normals(half1.size()); + + for (label facei = 0; facei < size()/2; facei++) + { + half0Normals[facei] = operator[](facei).normal(points()); + label nbrFacei = facei+size()/2; + half1Normals[facei] = operator[](nbrFacei).normal(points()); + + scalar magSf = mag(half0Normals[facei]); + scalar nbrMagSf = mag(half1Normals[facei]); + scalar avSf = (magSf + nbrMagSf)/2.0; + + if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL) + { + // Undetermined normal. Use dummy normal to force separation + // check. (note use of sqrt(VSMALL) since that is how mag + // scales) + half0Normals[facei] = point(1, 0, 0); + half1Normals[facei] = half0Normals[facei]; + } + else if (mag(magSf - nbrMagSf)/avSf > coupledPolyPatch::matchTol) + { + FatalErrorIn + ( + "cyclicPolyPatch::calcTransforms()" + ) << "face " << facei << " area does not match neighbour " + << nbrFacei << " by " + << 100*mag(magSf - nbrMagSf)/avSf + << "% -- possible face ordering problem." << endl + << "patch:" << name() + << " my area:" << magSf + << " neighbour area:" << nbrMagSf + << " matching tolerance:" << coupledPolyPatch::matchTol + << endl + << "Mesh face:" << start()+facei + << " vertices:" + << IndirectList(points(), operator[](facei))() + << endl + << "Neighbour face:" << start()+nbrFacei + << " vertices:" + << IndirectList(points(), operator[](nbrFacei))() + << endl + << "Rerun with cyclic debug flag set" + << " for more information." << exit(FatalError); + } + else + { + half0Normals[facei] /= magSf; + half1Normals[facei] /= nbrMagSf; + } + } + + + // Calculate transformation tensors calcTransformTensors ( - f0.centre(points), - fn2.centre(points), - nf0, - nfn2 + half0Ctrs, + half1Ctrs, + half0Normals, + half1Normals, + half0Tols ); } } @@ -86,14 +172,7 @@ bool Foam::cyclicPolyPatch::getGeometricHalves ) const { // Calculate normals - vectorField normals(pp.size()); - - forAll(pp, faceI) - { - normals[faceI] = pp[faceI].normal(pp.points()); - } - normals /= mag(normals) + VSMALL; - + const vectorField& faceNormals = pp.faceNormals(); // Find edges with sharp angles. boolList regionEdge(pp.nEdges(), false); @@ -106,9 +185,11 @@ bool Foam::cyclicPolyPatch::getGeometricHalves { const labelList& eFaces = edgeFaces[edgeI]; + // Check manifold edges for sharp angle. + // (Non-manifold already handled by patchZones) if (eFaces.size() == 2) { - if ((normals[eFaces[0]] & normals[eFaces[1]])< featureCos_) + if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_) { regionEdge[edgeI] = true; @@ -138,7 +219,11 @@ bool Foam::cyclicPolyPatch::getGeometricHalves for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++) { - OFstream stream(name()+"_zone_"+Foam::name(zoneI)+".obj"); + OFstream stream + ( + boundaryMesh().mesh().time().path() + /name()+"_zone_"+Foam::name(zoneI)+".obj" + ); Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone " << zoneI << " face centres to OBJ file " << stream.name() << endl; @@ -152,9 +237,6 @@ bool Foam::cyclicPolyPatch::getGeometricHalves nZoneFaces[zoneI] = zoneFaces.size(); } - - Pout<< "cyclicPolyPatch::getGeometricHalves : Number of faces per zone:" - << nZoneFaces << endl; } @@ -168,7 +250,7 @@ bool Foam::cyclicPolyPatch::getGeometricHalves if (debug) { Pout<< "cyclicPolyPatch::getGeometricHalves :" - << " falling back to normal comparison" << endl; + << " falling back to face-normal comparison" << endl; } label n0Faces = 0; half0ToPatch.setSize(pp.size()); @@ -177,9 +259,9 @@ bool Foam::cyclicPolyPatch::getGeometricHalves half1ToPatch.setSize(pp.size()); // Compare to face 0 normal. - forAll(normals, faceI) + forAll(faceNormals, faceI) { - if ((normals[faceI] & normals[0]) > 0) + if ((faceNormals[faceI] & faceNormals[0]) > 0) { half0ToPatch[n0Faces++] = faceI; } @@ -191,21 +273,26 @@ bool Foam::cyclicPolyPatch::getGeometricHalves half0ToPatch.setSize(n0Faces); half1ToPatch.setSize(n1Faces); - Pout<< "cyclicPolyPatch::getGeometricHalves :" - << " Number of faces per zone:(" - << n0Faces << ' ' << n1Faces << ')' << endl; + if (debug) + { + Pout<< "cyclicPolyPatch::getGeometricHalves :" + << " Number of faces per zone:(" + << n0Faces << ' ' << n1Faces << ')' << endl; + } } if (half0ToPatch.size() != half1ToPatch.size()) { + fileName casePath(boundaryMesh().mesh().time().path()); + // Dump halves { - fileName nm0(name()+"_half0_faces.obj"); + fileName nm0(casePath/name()+"_half0_faces.obj"); Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0" << " faces to OBJ file " << nm0 << endl; writeOBJ(nm0, IndirectList(pp, half0ToPatch)(), pp.points()); - fileName nm1(name()+"_half1_faces.obj"); + fileName nm1(casePath/name()+"_half1_faces.obj"); Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1" << " faces to OBJ file " << nm1 << endl; writeOBJ(nm1, IndirectList(pp, half1ToPatch)(), pp.points()); @@ -213,7 +300,7 @@ bool Foam::cyclicPolyPatch::getGeometricHalves // Dump face centres { - OFstream str0(name()+"_half0.obj"); + OFstream str0(casePath/name()+"_half0.obj"); Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0" << " face centres to OBJ file " << str0.name() << endl; @@ -222,7 +309,7 @@ bool Foam::cyclicPolyPatch::getGeometricHalves writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points())); } - OFstream str1(name()+"_half1.obj"); + OFstream str1(casePath/name()+"_half1.obj"); Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1" << " face centres to OBJ file " << str1.name() << endl; forAll(half1ToPatch, i) @@ -235,15 +322,13 @@ bool Foam::cyclicPolyPatch::getGeometricHalves ( "cyclicPolyPatch::getGeometricHalves" "(const primitivePatch&, labelList&, labelList&) const" - ) << " patch:" << name() << " : " - << "Patch " << name() << " gets decomposed in two zones of" + ) << "Patch " << name() << " gets decomposed in two zones of" << "inequal size: " << half0ToPatch.size() << " and " << half1ToPatch.size() << endl << "This means that the patch is either not two separate regions" << " or one region where the angle between the different regions" << " is not sufficiently sharp." << endl - << "Please use topological matching or adapt the featureCos" - << " setting" << endl + << "Please adapt the featureCos setting." << endl << "Continuing with incorrect face ordering from now on!" << endl; return false; @@ -278,7 +363,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors anchors0 = getAnchorPoints(half0Faces, pp.points()); half1Ctrs = calcFaceCentres(half1Faces, pp.points()); - if (mag(n0 & n1) < 1-SMALL) + if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol) { if (debug) { @@ -288,6 +373,7 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors // Rotation (around origin) const tensor reverseT(rotationTensor(n0, -n1)); + const tensor forwardT(rotationTensor(-n1, n0)); // Rotation forAll(half0Ctrs, faceI) @@ -393,7 +479,7 @@ bool Foam::cyclicPolyPatch::matchAnchors SeriousErrorIn ( "cyclicPolyPatch::matchAnchors(..)" - ) << " patch:" << name() << " : " + ) << "Patch:" << name() << " : " << "Cannot find point on face " << f << " with vertices:" << IndirectList(pp.points(), f)() @@ -583,16 +669,20 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const if (debug) { - Pout<< "Writing file coupledPoints.obj with coordinates of " - << "coupled points" << endl; - - OFstream str("coupledPoints.obj"); + OFstream str + ( + boundaryMesh().mesh().time().path() + /"coupledPoints.obj" + ); label vertI = 0; + Pout<< "Writing file " << str.name() << " with coordinates of " + << "coupled points" << endl; + forAll(connected, i) { - const point& a = localPoints()[connected[i][0]]; - const point& b = localPoints()[connected[i][1]]; + const point& a = points()[meshPoints()[connected[i][0]]]; + const point& b = points()[meshPoints()[connected[i][1]]]; str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl; str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl; @@ -601,6 +691,10 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledPoints() const str<< "l " << vertI-1 << ' ' << vertI << nl; } } + + // Remove any addressing calculated for the coupled edges calculation + const_cast(static_cast(*this)) + .clearOut(); } return *coupledPointsPtr_; } @@ -637,7 +731,15 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const // Convert edge end points to corresponding points on halfB // side. - edgeMap.insert(edge(aToB[e[0]], aToB[e[1]]), edgeI); + Map