diff options
| author |   Thomas Sachau <tommy@gentoo.org> | 2008-08-31 20:26:05 +0000 |
|---|---|---|
| committer | Thomas Sachau <tommy@gentoo.org> | 2008-08-31 20:26:05 +0000 |
| commit | f56dfadd2c9d2c33867cd214a4b3291d782d1358 (patch) | |
| tree | fd933c84cb7465569c3b06ab154433e933331726 | |
| parent | sci-libs/openfoam-utilities: Drop old version (diff) | |
| download | sunrise-f56dfadd2c9d2c33867cd214a4b3291d782d1358.tar.gz sunrise-f56dfadd2c9d2c33867cd214a4b3291d782d1358.tar.bz2 sunrise-f56dfadd2c9d2c33867cd214a4b3291d782d1358.zip | |
sci-libs/openfoam-utilities: Move big patch to external source
svn path=/sunrise/; revision=6923
4 files changed, 10 insertions, 7367 deletions
diff --git a/sci-libs/openfoam-utilities/ChangeLog b/sci-libs/openfoam-utilities/ChangeLog index 729628176..9e74a08f2 100644 --- a/sci-libs/openfoam-utilities/ChangeLog +++ b/sci-libs/openfoam-utilities/ChangeLog @@ -3,6 +3,10 @@ # $Header: $ 31 Aug 2008; Thomas Sachau (Tommy[D]) <tommy@gentoo.org> + -files/openfoam-utilities-1.4.1_p20080827.patch: + Move big patch to external source + + 31 Aug 2008; Thomas Sachau (Tommy[D]) <tommy@gentoo.org> -openfoam-utilities-1.4.1_p20080328.ebuild, -files/openfoam-utilities-1.4.1_p20080328.patch, -files/openfoam-utilities-compile-1.4.1_p20080328.patch: diff --git a/sci-libs/openfoam-utilities/Manifest b/sci-libs/openfoam-utilities/Manifest index 7f2681a5a..f6cb12b94 100644 --- a/sci-libs/openfoam-utilities/Manifest +++ b/sci-libs/openfoam-utilities/Manifest @@ -1,9 +1,9 @@ AUX OpenFOAM-1.5-compile.patch 15006 RMD160 ba8423526b5244e3c30d9d38830a2fe79e3c2a1a SHA1 7d275039cea1fe8a3c28fafeda1fef3665360f83 SHA256 6cb940b6c559a846ec65184db8f7c7966d1bef105d5bdad6ca4afd3f1b4d5b89 -AUX openfoam-utilities-1.4.1_p20080827.patch 219964 RMD160 3b207a535e5978404d62976f38d26d0b4542443e SHA1 b564ca1dc2531bc15ceac06a31386e625d3f75e2 SHA256 01f6fdac61ca2ddbbf6c77e45b9f23d1d2ced617260b97710da22a5b4f4c04d1 AUX openfoam-utilities-compile-1.4.1_p20080827.patch 996 RMD160 0debfeb112bed547317d8d61a3fbdae1490050ba SHA1 87c9b4ace94e4eb7b0c16a39e041ec72fca7877e SHA256 1e6062ff63e4367067c229cc8348f5bbd5e6e96f2b9f109bbacb73f3f16310d5 DIST OpenFOAM-1.4.1.General.gtgz 148526808 RMD160 e25d8bdfa63f15eeeb7b9f1cef09cc26fb7bef74 SHA1 56bbbf5b33c49d08cda35088a65b24d7dc59014f SHA256 c765b36639b42c737bc9ba1ac13c0f66efe20ee4a9f71a6ef987e86ebd50da28 DIST OpenFOAM-1.5.General.gtgz 117334661 RMD160 bccaa9f8f99d31aa6c791d40b30dd9ad4f534041 SHA1 3577f562dc1f54bb32e1e0ef43f979418212c2f7 SHA256 d4cba2d9475523a53cea80b8d39da70d12bfffb9f46e2d1442946ba4a23efd31 -EBUILD openfoam-utilities-1.4.1_p20080827.ebuild 4310 RMD160 0817be669c8dec4c3fa67f2103bec34c60ecbfd9 SHA1 ae43d8095accd00e97397c5067d60209dd9173f0 SHA256 0779f9ce3996ae8cc02911401520f7e4123abfa41cbbf9d6b7af8177f0f3c862 +DIST openfoam-utilities-1.4.1_p20080827.patch 219964 RMD160 3b207a535e5978404d62976f38d26d0b4542443e SHA1 b564ca1dc2531bc15ceac06a31386e625d3f75e2 SHA256 01f6fdac61ca2ddbbf6c77e45b9f23d1d2ced617260b97710da22a5b4f4c04d1 +EBUILD openfoam-utilities-1.4.1_p20080827.ebuild 4372 RMD160 331f929bdd0067fad5ea04c5c60344f05e14efa7 SHA1 47d257c3b670c74b525551786dcd8702f1e12f63 SHA256 93d43803db750ff010c12d10b6b76fb2d607a417b06e90997b4c59a6d53f3975 EBUILD openfoam-utilities-1.5.ebuild 2194 RMD160 c72a0fd5da83a4c41ae2923bcc3fda49959459bc SHA1 312fc2064f041f315f06ed388eacf87a7d8589cc SHA256 920b09bfc513e9db6dac706bade7b23e1d3216d7a0ba8f5d03bc44ea1a74de02 -MISC ChangeLog 2332 RMD160 d07da05a27d70078ef6e6f3bef9c7532e923f6b3 SHA1 c1cccff4bdaaa5b117ed1f17ddf33c60dc64772e SHA256 00d173c1c7276306638c01509255d2989ad6b66c9bf6c2f1c109ede220e112a8 +MISC ChangeLog 2479 RMD160 271e770579fbd58f44514de55e162358745a828f SHA1 57ae54cf74e92563e54c5665594aa4626cfa7d9b SHA256 403b178da8f510c8c3660004bba3460e104847b6c34cbfc59fc19386c8317f1e MISC metadata.xml 170 RMD160 645927a396fdc21cdeb089fe42c5397332420ea6 SHA1 ac7f48a14fec325926f9ce1be8fbf1f311b4f2e4 SHA256 d797a2ec6f9dc516c9f9c1a758ee87ad3e8c43101b5dc76c2f872d5bd4639b42 diff --git a/sci-libs/openfoam-utilities/files/openfoam-utilities-1.4.1_p20080827.patch b/sci-libs/openfoam-utilities/files/openfoam-utilities-1.4.1_p20080827.patch deleted file mode 100644 index 491f32e62..000000000 --- a/sci-libs/openfoam-utilities/files/openfoam-utilities-1.4.1_p20080827.patch +++ /dev/null @@ -1,7362 +0,0 @@ -Index: postProcessing/dataConversion/foamToVTK/writeFuns.C -=================================================================== ---- applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C (Revision 30) -+++ applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.C (Revision 784) -@@ -229,14 +229,36 @@ - void Foam::writeFuns::insert(const vector& pt, DynamicList<floatScalar>& dest) - { - for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++) - { - dest.append(float(pt[cmpt])); - } - } -+void Foam::writeFuns::insert -+( -+ const symmTensor& pt, -+ DynamicList<floatScalar>& dest -+) -+{ -+ for (direction cmpt = 0; cmpt < symmTensor::nComponents; cmpt++) -+ { -+ dest.append(float(pt[cmpt])); -+ } -+} -+void Foam::writeFuns::insert -+( -+ const sphericalTensor& pt, -+ DynamicList<floatScalar>& dest -+) -+{ -+ for (direction cmpt = 0; cmpt < sphericalTensor::nComponents; cmpt++) -+ { -+ dest.append(float(pt[cmpt])); -+ } -+} - void Foam::writeFuns::insert(const tensor& pt, DynamicList<floatScalar>& dest) - { - for (direction cmpt = 0; cmpt < tensor::nComponents; cmpt++) - { - dest.append(float(pt[cmpt])); - } - } -Index: postProcessing/dataConversion/foamToVTK/writeFuns.H -=================================================================== ---- applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H (Revision 30) -+++ applications/utilities/postProcessing/dataConversion/foamToVTK/writeFuns.H (Revision 784) -@@ -167,14 +167,16 @@ - const label, - const label - ); - - //- convert to VTK and store - static void insert(const scalar&, DynamicList<floatScalar>& dest); - static void insert(const point&, DynamicList<floatScalar>& dest); -+ static void insert(const symmTensor&, DynamicList<floatScalar>& dest); -+ static void insert(const sphericalTensor&, DynamicList<floatScalar>& dest); - static void insert(const tensor&, DynamicList<floatScalar>& dest); - - //- Append elements to DynamicList - static void insert(const labelList&, DynamicList<label>&); - template<class Type> - static void insert(const List<Type>&, DynamicList<floatScalar>&); - - -Index: mesh/manipulation/Optional/Allwmake -=================================================================== ---- applications/utilities/mesh/manipulation/Optional/Allwmake (Revision 30) -+++ applications/utilities/mesh/manipulation/Optional/Allwmake (Revision 784) -@@ -1,11 +1,11 @@ - #!/bin/sh - - READLINE=0 - if [ -f /usr/include/readline/readline.h ]; then - echo "Found readline/readline.h include file. Enabling readline support." - READLINE=1 -- export READLINELINK="-lreadline" -+ export READLINELINK="-lreadline -lcurses" - break - fi - export READLINE - wmake setSet - -Index: mesh/advanced/autoRefineMesh/hexRef8.H -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/hexRef8.H (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/hexRef8.H (Revision 784) -@@ -0,0 +1,534 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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::hexRef8 -+ -+Description -+ Refinement of (split) hexes using polyTopoChange. -+ -+SourceFiles -+ hexRef8.C -+ -+\*---------------------------------------------------------------------------*/ -+ -+#ifndef hexRef8_H -+#define hexRef8_H -+ -+#include "labelIOList.H" -+#include "face.H" -+#include "labelHashSet.H" -+#include "DynamicList.H" -+#include "primitivePatch.H" -+#include "removeFaces.H" -+#include "refinementHistory.H" -+#include "PackedList.H" -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+namespace Foam -+{ -+ -+// Forward declaration of classes -+class polyMesh; -+class polyPatch; -+class polyTopoChange; -+class mapPolyMesh; -+class mapDistributePolyMesh; -+ -+/*---------------------------------------------------------------------------*\ -+ Class hexRef8 Declaration -+\*---------------------------------------------------------------------------*/ -+ -+class hexRef8 -+{ -+ // Private data -+ -+ //- Reference to underlying mesh. -+ const polyMesh& mesh_; -+ -+ //- Per cell the refinement level -+ labelIOList cellLevel_; -+ -+ //- Per point the refinement level -+ labelIOList pointLevel_; -+ -+ //- Typical edge length between unrefined points -+ const scalar level0Edge_; -+ -+ //- Refinement history -+ refinementHistory history_; -+ -+ //- Face remover engine -+ removeFaces faceRemover_; -+ -+ //- Level of saved points -+ Map<label> savedPointLevel_; -+ -+ //- Level of saved cells -+ Map<label> savedCellLevel_; -+ -+ -+ // Private Member Functions -+ -+ //- Reorder according to map. -+ static void reorder -+ ( -+ const labelList& map, -+ const label len, -+ const label null, -+ labelList& elems -+ ); -+ -+ //- Get patch and zone info -+ void getFaceInfo -+ ( -+ const label faceI, -+ label& patchID, -+ label& zoneID, -+ label& zoneFlip -+ ) const; -+ -+ //- Adds a face on top of existing faceI. Reverses if nessecary. -+ label addFace -+ ( -+ polyTopoChange& meshMod, -+ const label faceI, -+ const face& newFace, -+ const label own, -+ const label nei -+ ) const; -+ -+ //- Adds internal face from point. No checks on reversal. -+ label addInternalFace -+ ( -+ polyTopoChange& meshMod, -+ const label meshFaceI, -+ const label meshPointI, -+ const face& newFace, -+ const label own, -+ const label nei -+ ) const; -+ -+ //- Modifies existing faceI for either new owner/neighbour or new face -+ // points. Reverses if nessecary. -+ void modFace -+ ( -+ polyTopoChange& meshMod, -+ const label faceI, -+ const face& newFace, -+ const label own, -+ const label nei -+ ) const; -+ -+ scalar getLevel0EdgeLength() const; -+ -+ //- Get cell added to point of cellI (if any) -+ label getAnchorCell -+ ( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const label cellI, -+ const label faceI, -+ const label pointI -+ ) const; -+ -+ //- Get new owner and neighbour (in unspecified order) of pointI -+ // on faceI. -+ void getFaceNeighbours -+ ( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const label faceI, -+ const label pointI, -+ -+ label& own, -+ label& nei -+ ) const; -+ -+ -+ //- Get index of minimum pointlevel. -+ label findMinLevel(const labelList& f) const; -+ //- Get maximum pointlevel. -+ label findMaxLevel(const labelList& f) const; -+ //- Count number of vertices <= anchorLevel -+ label countAnchors(const labelList&, const label) const; -+ //- Find index of point with wantedLevel, starting from fp. -+ label findLevel -+ ( -+ const face& f, -+ const label startFp, -+ const bool searchForward, -+ const label wantedLevel -+ ) const; -+ //- Gets level such that the face has four points <= level. -+ label getAnchorLevel(const label faceI) const; -+ -+ ////- Print levels of list of points. -+ //void printLevels(Ostream&, const labelList&) const; -+ -+ //- debug:check orientation of added internal face -+ static void checkInternalOrientation -+ ( -+ polyTopoChange& meshMod, -+ const label cellI, -+ const label faceI, -+ const point& ownPt, -+ const point& neiPt, -+ const face& newFace -+ ); -+ -+ //- debug:check orientation of new boundary face -+ static void checkBoundaryOrientation -+ ( -+ polyTopoChange& meshMod, -+ const label cellI, -+ const label faceI, -+ const point& ownPt, -+ const point& boundaryPt, -+ const face& newFace -+ ); -+ -+ //- If p0 and p1 are existing vertices check if edge is split and insert -+ // splitPoint. -+ void insertEdgeSplit -+ ( -+ const labelList& edgeMidPoint, -+ const label p0, -+ const label p1, -+ DynamicList<label>& verts -+ ) const; -+ -+ //- Store in maps correspondence from midpoint to anchors and faces. -+ label storeMidPointInfo -+ ( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const labelList& cellMidPoint, -+ const labelList& edgeMidPoint, -+ const label cellI, -+ const label faceI, -+ const bool faceOrder, -+ const label midPointI, -+ const label anchorPointI, -+ const label faceMidPointI, -+ -+ Map<edge>& midPointToAnchors, -+ Map<edge>& midPointToFaceMids, -+ polyTopoChange& meshMod -+ ) const; -+ -+ //- Create all internal faces from an unsplit face. -+ void createInternalFromSplitFace -+ ( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const labelList& cellMidPoint, -+ const labelList& faceMidPoint, -+ const labelList& edgeMidPoint, -+ const label cellI, -+ const label faceI, -+ -+ Map<edge>& midPointToAnchors, -+ Map<edge>& midPointToFaceMids, -+ polyTopoChange& meshMod, -+ label& nFacesAdded -+ ) const; -+ -+ //- Create all internal faces to split cellI into 8. -+ void createInternalFaces -+ ( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const labelList& cellMidPoint, -+ const labelList& faceMidPoint, -+ const labelList& faceAnchorLevel, -+ const labelList& edgeMidPoint, -+ const label cellI, -+ polyTopoChange& meshMod -+ ) const; -+ -+ //- Store vertices from startFp upto face split point. -+ // Used when splitting face into 4. -+ void walkFaceToMid -+ ( -+ const labelList& edgeMidPoint, -+ const label cLevel, -+ const label faceI, -+ const label startFp, -+ DynamicList<label>& faceVerts -+ ) const; -+ -+ //- Same as walkFaceToMid but now walk back. -+ void walkFaceFromMid -+ ( -+ const labelList& edgeMidPoint, -+ const label cLevel, -+ const label faceI, -+ const label startFp, -+ DynamicList<label>& faceVerts -+ ) const; -+ -+ //- Updates refineCell so consistent 2:1 refinement. Returns local -+ // number of cells changed. -+ label faceConsistentRefinement -+ ( -+ const bool maxSet, -+ PackedList<1>& refineCell -+ ) const; -+ -+ //- Check wanted refinement for 2:1 consistency -+ void checkWantedRefinementLevels(const labelList&) const; -+ -+ -+ -+ //- Disallow default bitwise copy construct -+ hexRef8(const hexRef8&); -+ -+ //- Disallow default bitwise assignment -+ void operator=(const hexRef8&); -+ -+ -+public: -+ -+ //- Runtime type information -+ ClassName("hexRef8"); -+ -+ -+ // Constructors -+ -+ //- Construct from mesh, read_if_present refinement data -+ // (from write below) -+ hexRef8(const polyMesh& mesh); -+ -+ //- Construct from mesh and un/refinement data. -+ hexRef8 -+ ( -+ const polyMesh& mesh, -+ const labelList& cellLevel, -+ const labelList& pointLevel, -+ const refinementHistory& history -+ ); -+ -+ //- Construct from mesh and refinement data. -+ hexRef8 -+ ( -+ const polyMesh& mesh, -+ const labelList& cellLevel, -+ const labelList& pointLevel -+ ); -+ -+ -+ // Member Functions -+ -+ // Access -+ -+ const labelIOList& cellLevel() const -+ { -+ return cellLevel_; -+ } -+ -+ const labelIOList& pointLevel() const -+ { -+ return pointLevel_; -+ } -+ -+ const refinementHistory& history() const -+ { -+ return history_; -+ } -+ -+ //- Typical edge length between unrefined points -+ scalar level0EdgeLength() const -+ { -+ return level0Edge_; -+ } -+ -+ // Refinement -+ -+ //- Helper:get points of a cell without using cellPoints addressing -+ labelList cellPoints(const label cellI) const; -+ -+ //- Given valid mesh and current cell level and proposed -+ // cells to refine calculate any clashes (due to 2:1) and return -+ // ok list of cells to refine. -+ // Either adds cells to refine to set (maxSet = true) or -+ // removes cells to refine (maxSet = false) -+ labelList consistentRefinement -+ ( -+ const labelList& cellsToRefine, -+ const bool maxSet -+ ) const; -+ -+ //- Like consistentRefinement but slower: -+ // - specify number of cells between consecutive refinement levels -+ // (consistentRefinement equivalent to 1) -+ // - specify max level difference between point-connected cells. -+ // (-1 to disable) Note that with normal 2:1 limitation -+ // (maxFaceDiff=1) there can be 8:1 size difference across point -+ // connected cells so maxPointDiff allows you to make that less. -+ // cellsToRefine : cells we're thinking about refining. It will -+ // extend this set. All refinement levels will be -+ // at least maxFaceDiff layers thick. -+ // facesToCheck : additional faces where to implement the -+ // maxFaceDiff thickness (usually only boundary -+ // faces) -+ labelList consistentSlowRefinement -+ ( -+ const label maxFaceDiff, -+ const labelList& cellsToRefine, -+ const labelList& facesToCheck, -+ const label maxPointDiff, -+ const labelList& pointsToCheck -+ ) const; -+ -+ //- Like consistentSlowRefinement but uses different meshWave -+ // (proper distance instead of toplogical count). No point checks -+ // yet. -+ labelList consistentSlowRefinement2 -+ ( -+ const label maxFaceDiff, -+ const labelList& cellsToRefine, -+ const labelList& facesToCheck -+ ) const; -+ -+ //- Insert refinement. All selected cells will be split into 8. -+ // Returns per element in cells the 8 cells they were split into. -+ // Guarantees that the 0th element is the original cell label. -+ // Mapping: -+ // -split cells: 7 new ones get added from original -+ // -split faces: original gets modified; 3 new ones get added -+ // from original -+ // -added internal faces: added from original cell face(if -+ // that was internal) or created out-of-nothing (so will not -+ // get mapped!). Note: could make this inflate from point but -+ // that will allocate interpolation. -+ // -points added to split edge: added from edge start() -+ // -midpoints added: added from cellPoints[0]. -+ labelListList setRefinement -+ ( -+ const labelList& cells, -+ polyTopoChange& -+ ); -+ -+ //- Update local numbering for changed mesh. -+ void updateMesh(const mapPolyMesh&); -+ -+ -+ // Restoring : is where other processes delete and reinsert data. -+ // These callbacks allow this to restore the cellLevel -+ // and pointLevel for reintroduced points. -+ // Is not related to undoing my refinement -+ -+ //- Signal points/face/cells for which to store data -+ void storeData -+ ( -+ const labelList& pointsToStore, -+ const labelList& facesToStore, -+ const labelList& cellsToStore -+ ); -+ -+ //- Update local numbering + undo -+ // Data to restore given as new pointlabel + stored pointlabel -+ // (i.e. what was in pointsToStore) -+ void updateMesh -+ ( -+ const mapPolyMesh&, -+ const Map<label>& pointsToRestore, -+ const Map<label>& facesToRestore, -+ const Map<label>& cellsToRestore -+ ); -+ -+ -+ //- Update local numbering for subsetted mesh. -+ // Gets new-to-old maps. Not compatible with unrefinement. -+ void subset -+ ( -+ const labelList& pointMap, -+ const labelList& faceMap, -+ const labelList& cellMap -+ ); -+ -+ //- Update local numbering for mesh redistribution -+ void distribute(const mapDistributePolyMesh&); -+ -+ //- Debug: Check coupled mesh for correctness -+ void checkMesh() const; -+ -+ //- Debug: Check 2:1 consistency across faces. -+ // maxPointDiff==-1 : only check 2:1 across faces -+ // maxPointDiff!=-1 : check point-connected cells. -+ void checkRefinementLevels -+ ( -+ const label maxPointDiff, -+ const labelList& pointsToCheck -+ ) const; -+ -+ // Unrefinement (undoing refinement, not arbitrary coarsening) -+ -+ //- Return the points at the centre of top-level split cells -+ // that can be unsplit. -+ labelList getSplitPoints() const; -+ -+ //- Given proposed -+ // splitPoints to unrefine according to calculate any clashes -+ // (due to 2:1) and return ok list of points to unrefine. -+ // Either adds points to refine to set (maxSet = true) or -+ // removes points to refine (maxSet = false) -+ labelList consistentUnrefinement -+ ( -+ const labelList& pointsToUnrefine, -+ const bool maxSet -+ ) const; -+ -+ //- Remove some refinement. Needs to be supplied output of -+ // consistentUnrefinement. Only call if undoable set. -+ // All 8 pointCells of a split point will be combined into -+ // the lowest numbered cell of those 8. -+ void setUnrefinement -+ ( -+ const labelList& splitPointLabels, -+ polyTopoChange& -+ ); -+ -+ // Write -+ -+ // Set instance for mesh files -+ void setInstance(const fileName& inst); -+ -+ //- Force writing refinement+history to polyMesh directory. -+ bool write() const; -+ -+}; -+ -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+} // End namespace Foam -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+#endif -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/refinementDistanceDataI.H -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceDataI.H (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceDataI.H (Revision 784) -@@ -0,0 +1,286 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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 "transform.H" -+#include "polyMesh.H" -+ -+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -+ -+// Returns the wanted level -+inline Foam::label Foam::refinementDistanceData::wantedLevel(const point& pt) -+ const -+{ -+ const scalar distSqr = magSqr(pt-origin_); -+ -+ // Get the size at the origin level -+ scalar levelSize = level0Size_/(1<<originLevel_); -+ -+ scalar r = 0; -+ -+ for (label level = originLevel_; level >= 0; --level) -+ { -+ // Current range -+ r += levelSize; -+ -+ // Check if our distance is within influence sphere -+ if (sqr(r) > distSqr) -+ { -+ return level; -+ } -+ -+ // Lower level will have double the size -+ levelSize *= 2; -+ } -+ return 0; -+} -+ -+ -+inline bool Foam::refinementDistanceData::update -+( -+ const point& pos, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+) -+{ -+ if (!valid()) -+ { -+ if (!neighbourInfo.valid()) -+ { -+ FatalErrorIn("refinementDistanceData::update(..)") -+ << "problem" << abort(FatalError); -+ } -+ operator=(neighbourInfo); -+ return true; -+ } -+ -+ // Determine wanted level at current position. -+ label cellLevel = wantedLevel(pos); -+ -+ // Determine wanted level coming through the neighbour -+ label nbrLevel = neighbourInfo.wantedLevel(pos); -+ -+ if (nbrLevel > cellLevel) -+ { -+ operator=(neighbourInfo); -+ return true; -+ } -+ else if (nbrLevel == cellLevel) -+ { -+ scalar myDistSqr = magSqr(pos-origin_); -+ scalar nbrDistSqr = magSqr(pos - neighbourInfo.origin()); -+ scalar diff = myDistSqr - nbrDistSqr; -+ -+ if (diff < 0) -+ { -+ // already nearest -+ return false; -+ } -+ -+ if ((diff < SMALL) || ((myDistSqr > SMALL) && (diff/myDistSqr < tol))) -+ { -+ // don't propagate small changes -+ return false; -+ } -+ else -+ { -+ // update with new values -+ operator=(neighbourInfo); -+ return true; -+ } -+ } -+ else -+ { -+ return false; -+ } -+} -+ -+ -+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -+ -+// Null constructor -+inline Foam::refinementDistanceData::refinementDistanceData() -+: -+ level0Size_(-1) -+{} -+ -+ -+// Construct from components -+inline Foam::refinementDistanceData::refinementDistanceData -+( -+ const scalar level0Size, -+ const point& origin, -+ const label originLevel -+) -+: -+ level0Size_(level0Size), -+ origin_(origin), -+ originLevel_(originLevel) -+{} -+ -+ -+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -+ -+inline bool Foam::refinementDistanceData::valid() const -+{ -+ return level0Size_ != -1; -+} -+ -+ -+// No geometric data so never any problem on cyclics -+inline bool Foam::refinementDistanceData::sameGeometry -+( -+ const polyMesh&, -+ const refinementDistanceData&, -+ const scalar -+) const -+{ -+ return true; -+} -+ -+ -+inline void Foam::refinementDistanceData::leaveDomain -+( -+ const polyMesh&, -+ const polyPatch& patch, -+ const label patchFaceI, -+ const point& faceCentre -+) -+{ -+ origin_ -= faceCentre; -+} -+ -+ -+inline void Foam::refinementDistanceData::transform -+( -+ const polyMesh&, -+ const tensor& rotTensor -+) -+{ -+ origin_ = Foam::transform(rotTensor, origin_); -+} -+ -+ -+// Update absolute geometric quantities. -+inline void Foam::refinementDistanceData::enterDomain -+( -+ const polyMesh&, -+ const polyPatch& patch, -+ const label patchFaceI, -+ const point& faceCentre -+) -+{ -+ // back to absolute form -+ origin_ += faceCentre; -+} -+ -+ -+// Update cell with neighbouring face information -+inline bool Foam::refinementDistanceData::updateCell -+( -+ const polyMesh& mesh, -+ const label thisCellI, -+ const label neighbourFaceI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+) -+{ -+ const point& pos = mesh.cellCentres()[thisCellI]; -+ -+ return update(pos, neighbourInfo, tol); -+} -+ -+ -+// Update face with neighbouring cell information -+inline bool Foam::refinementDistanceData::updateFace -+( -+ const polyMesh& mesh, -+ const label thisFaceI, -+ const label neighbourCellI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+) -+{ -+ const point& pos = mesh.faceCentres()[thisFaceI]; -+ -+ return update(pos, neighbourInfo, tol); -+} -+ -+ -+// Update face with coupled face information -+inline bool Foam::refinementDistanceData::updateFace -+( -+ const polyMesh& mesh, -+ const label thisFaceI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+) -+{ -+ const point& pos = mesh.faceCentres()[thisFaceI]; -+ -+ return update(pos, neighbourInfo, tol); -+} -+ -+ -+// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // -+ -+inline bool Foam::refinementDistanceData::operator== -+( -+ const Foam::refinementDistanceData& rhs -+) -+ const -+{ -+ if (!valid()) -+ { -+ if (!rhs.valid()) -+ { -+ return true; -+ } -+ else -+ { -+ return false; -+ } -+ } -+ else -+ { -+ return -+ level0Size_ == rhs.level0Size_ -+ && origin_ == rhs.origin_ -+ && originLevel_ == rhs.originLevel_; -+ } -+} -+ -+ -+inline bool Foam::refinementDistanceData::operator!= -+( -+ const Foam::refinementDistanceData& rhs -+) -+ const -+{ -+ return !(*this == rhs); -+} -+ -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/surfaceToCell.H -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/surfaceToCell.H (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/surfaceToCell.H (Revision 784) -@@ -0,0 +1,217 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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::surfaceToCell -+ -+Description -+ A topoSetSource to select cells based on relation to surface. -+ -+ Selects: -+ - all cells inside/outside/cut by surface -+ - cells with centre nearer than XXX to surface -+ - cells with centre nearer than XXX to surface @b and with normal -+ at nearest point to centre and cell-corners differing by -+ more than YYY (i.e., point of high curvature) -+ -+SourceFiles -+ surfaceToCell.C -+ -+\*---------------------------------------------------------------------------*/ -+ -+#ifndef surfaceToCell_H -+#define surfaceToCell_H -+ -+#include "topoSetSource.H" -+#include "Map.H" -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+namespace Foam -+{ -+class triSurfaceSearch; -+class triSurface; -+ -+/*---------------------------------------------------------------------------*\ -+ Class surfaceToCell Declaration -+\*---------------------------------------------------------------------------*/ -+ -+class surfaceToCell -+: -+ public topoSetSource -+{ -+ -+ // Private data -+ -+ //- Add usage string -+ static addToUsageTable usage_; -+ -+ //- Name of surface file -+ fileName surfName_; -+ -+ //- Points which are outside -+ pointField outsidePoints_; -+ -+ //- Include cut cells -+ bool includeCut_; -+ -+ //- Include inside cells -+ bool includeInside_; -+ -+ //- Include outside cells -+ bool includeOutside_; -+ -+ //- if > 0 : include cells with distance from cellCentre to surface -+ // less than nearDist. -+ scalar nearDist_; -+ -+ //- if > -1 : include cells with normals at nearest surface points -+ // varying more than curvature_. -+ scalar curvature_; -+ -+ //- triSurface to search on. On pointer since can be external. -+ const triSurface* surfPtr_; -+ -+ //- search engine on surface. -+ const triSurfaceSearch* querySurfPtr_; -+ -+ //- whether I allocated above surface ptrs or whether they are -+ // external. -+ bool IOwnPtrs_; -+ -+ -+ // Private Member Functions -+ -+ //- Find index of nearest triangle to point. Returns triangle or -1 if -+ // not found within search span. -+ // Cache result under pointI. -+ static label getNearest -+ ( -+ const triSurfaceSearch& querySurf, -+ const label pointI, -+ const point& pt, -+ const vector& searchSpan, -+ Map<label>& cache -+ ); -+ -+ //- Return true if surface normal of nearest points to vertices on -+ // cell differ from that on cell centre. Points cached in -+ // pointToNearest. -+ bool differingPointNormals -+ ( -+ const triSurfaceSearch& querySurf, -+ const vector& span, -+ const label cellI, -+ const label cellTriI, -+ Map<label>& pointToNearest -+ ) const; -+ -+ -+ //- Depending on surface add to or delete from cellSet. -+ void combine(topoSet& set, const bool add) const; -+ -+ //- Check values at construction time. -+ void checkSettings() const; -+ -+ const triSurfaceSearch& querySurf() const -+ { -+ return *querySurfPtr_; -+ } -+ -+ -+public: -+ -+ //- Runtime type information -+ TypeName("surfaceToCell"); -+ -+ // Constructors -+ -+ //- Construct from components -+ surfaceToCell -+ ( -+ const polyMesh& mesh, -+ const fileName& surfName, -+ const pointField& outsidePoints, -+ const bool includeCut, -+ const bool includeInside, -+ const bool includeOutside, -+ const scalar nearDist, -+ const scalar curvature -+ ); -+ -+ //- Construct from components (supplied surface, surfaceSearch) -+ surfaceToCell -+ ( -+ const polyMesh& mesh, -+ const fileName& surfName, -+ const triSurface& surf, -+ const triSurfaceSearch& querySurf, -+ const pointField& outsidePoints, -+ const bool includeCut, -+ const bool includeInside, -+ const bool includeOutside, -+ const scalar nearDist, -+ const scalar curvature -+ ); -+ -+ //- Construct from dictionary -+ surfaceToCell -+ ( -+ const polyMesh& mesh, -+ const dictionary& dict -+ ); -+ -+ //- Construct from Istream -+ surfaceToCell -+ ( -+ const polyMesh& mesh, -+ Istream& -+ ); -+ -+ -+ // Destructor -+ -+ virtual ~surfaceToCell(); -+ -+ -+ // Member Functions -+ -+ virtual void applyToSet -+ ( -+ const topoSetSource::setAction action, -+ topoSet& -+ ) const; -+ -+}; -+ -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+} // End namespace Foam -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+#endif -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/refinementDistanceData.C -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceData.C (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceData.C (Revision 784) -@@ -0,0 +1,53 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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 "refinementDistanceData.H" -+ -+// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // -+ -+Foam::Ostream& Foam::operator<< -+( -+ Foam::Ostream& os, -+ const Foam::refinementDistanceData& wDist -+) -+{ -+ return os -+ << wDist.level0Size_ << token::SPACE << wDist.origin_ -+ << token::SPACE << wDist.originLevel_; -+} -+ -+ -+Foam::Istream& Foam::operator>> -+( -+ Foam::Istream& is, -+ Foam::refinementDistanceData& wDist -+) -+{ -+ return is >> wDist.level0Size_ >> wDist.origin_ >> wDist.originLevel_; -+} -+ -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/refinementDistanceData.H -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceData.H (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/refinementDistanceData.H (Revision 784) -@@ -0,0 +1,233 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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::refinementDistanceData -+ -+Description -+ Transfers refinement levels such that slow transition between levels is -+ maintained. Used in FaceCellWave. -+ -+SourceFiles -+ refinementDistanceDataI.H -+ refinementDistanceData.C -+ -+\*---------------------------------------------------------------------------*/ -+ -+#ifndef refinementDistanceData_H -+#define refinementDistanceData_H -+ -+#include "point.H" -+#include "tensor.H" -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+namespace Foam -+{ -+ -+class polyPatch; -+class polyMesh; -+ -+/*---------------------------------------------------------------------------*\ -+ Class refinementDistanceData Declaration -+\*---------------------------------------------------------------------------*/ -+ -+class refinementDistanceData -+{ -+ -+ // Private data -+ -+ //- unrefined (level0) buffer size (nBufferLayers*level0Size) -+ scalar level0Size_; -+ -+ //- nearest point with highest level -+ point origin_; -+ label originLevel_; -+ -+ -+ // Private Member Functions -+ -+ //- Updates with neighbouring data. Returns true if something changed. -+ inline bool update -+ ( -+ const point&, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+ ); -+ -+public: -+ -+ // Constructors -+ -+ //- Construct null -+ inline refinementDistanceData(); -+ -+ //- Construct from count -+ inline refinementDistanceData -+ ( -+ const scalar level0Size, -+ const point& origin, -+ const label level -+ ); -+ -+ -+ // Member Functions -+ -+ // Access -+ -+ inline scalar level0Size() const -+ { -+ return level0Size_; -+ } -+ -+ inline scalar& level0Size() -+ { -+ return level0Size_; -+ } -+ -+ inline const point& origin() const -+ { -+ return origin_; -+ } -+ -+ inline point& origin() -+ { -+ return origin_; -+ } -+ -+ inline label originLevel() const -+ { -+ return originLevel_; -+ } -+ -+ inline label& originLevel() -+ { -+ return originLevel_; -+ } -+ -+ -+ // Other -+ -+ //- Calculates the wanted level at a given point. Walks out from -+ // the origin. -+ inline label wantedLevel(const point& pt) const; -+ -+ -+ // Needed by FaceCellWave -+ -+ //- Check whether origin has been changed at all or -+ // still contains original (invalid) value. -+ inline bool valid() const; -+ -+ //- Check for identical geometrical data. Used for cyclics checking. -+ inline bool sameGeometry -+ ( -+ const polyMesh&, -+ const refinementDistanceData&, -+ const scalar -+ ) const; -+ -+ //- Convert any absolute coordinates into relative to (patch)face -+ // centre -+ inline void leaveDomain -+ ( -+ const polyMesh&, -+ const polyPatch&, -+ const label patchFaceI, -+ const point& faceCentre -+ ); -+ -+ //- Reverse of leaveDomain -+ inline void enterDomain -+ ( -+ const polyMesh&, -+ const polyPatch&, -+ const label patchFaceI, -+ const point& faceCentre -+ ); -+ -+ //- Apply rotation matrix to any coordinates -+ inline void transform -+ ( -+ const polyMesh&, -+ const tensor& -+ ); -+ -+ //- Influence of neighbouring face. -+ inline bool updateCell -+ ( -+ const polyMesh&, -+ const label thisCellI, -+ const label neighbourFaceI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+ ); -+ -+ //- Influence of neighbouring cell. -+ inline bool updateFace -+ ( -+ const polyMesh&, -+ const label thisFaceI, -+ const label neighbourCellI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+ ); -+ -+ //- Influence of different value on same face. -+ inline bool updateFace -+ ( -+ const polyMesh&, -+ const label thisFaceI, -+ const refinementDistanceData& neighbourInfo, -+ const scalar tol -+ ); -+ -+ // Member Operators -+ -+ // Needed for List IO -+ inline bool operator==(const refinementDistanceData&) const; -+ -+ inline bool operator!=(const refinementDistanceData&) const; -+ -+ -+ // IOstream Operators -+ -+ friend Ostream& operator<<(Ostream&, const refinementDistanceData&); -+ friend Istream& operator>>(Istream&, refinementDistanceData&); -+}; -+ -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+} // End namespace Foam -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+#include "refinementDistanceDataI.H" -+ -+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -+ -+#endif -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/hexRef8.C -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/hexRef8.C (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/hexRef8.C (Revision 784) -@@ -0,0 +1,5283 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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 "hexRef8.H" -+ -+#include "polyMesh.H" -+#include "polyTopoChange.H" -+#include "meshTools.H" -+#include "polyAddFace.H" -+#include "polyAddPoint.H" -+#include "polyAddCell.H" -+#include "polyModifyFace.H" -+#include "syncTools.H" -+#include "faceSet.H" -+#include "cellSet.H" -+#include "pointSet.H" -+#include "OFstream.H" -+#include "Time.H" -+#include "FaceCellWave.H" -+#include "mapDistributePolyMesh.H" -+#include "refinementData.H" -+#include "refinementDistanceData.H" -+ -+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -+ -+namespace Foam -+{ -+ defineTypeNameAndDebug(hexRef8, 0); -+ -+ // Reduction class. If x and y are not equal assign value. -+ template< int value > -+ class ifEqEqOp -+ { -+ public: -+ void operator()( label& x, const label& y ) const -+ { -+ x = (x==y) ? x:value; -+ } -+ }; -+} -+ -+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -+ -+void Foam::hexRef8::reorder -+( -+ const labelList& map, -+ const label len, -+ const label null, -+ labelList& elems -+) -+{ -+ labelList newElems(len, null); -+ -+ forAll(elems, i) -+ { -+ label newI = map[i]; -+ -+ if (newI >= len) -+ { -+ FatalErrorIn("hexRef8::reorder") << abort(FatalError); -+ } -+ -+ if (newI >= 0) -+ { -+ newElems[newI] = elems[i]; -+ } -+ } -+ -+ elems.transfer(newElems); -+} -+ -+ -+void Foam::hexRef8::getFaceInfo -+( -+ const label faceI, -+ label& patchID, -+ label& zoneID, -+ label& zoneFlip -+) const -+{ -+ patchID = -1; -+ -+ if (!mesh_.isInternalFace(faceI)) -+ { -+ patchID = mesh_.boundaryMesh().whichPatch(faceI); -+ } -+ -+ zoneID = mesh_.faceZones().whichZone(faceI); -+ -+ zoneFlip = false; -+ -+ if (zoneID >= 0) -+ { -+ const faceZone& fZone = mesh_.faceZones()[zoneID]; -+ -+ zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; -+ } -+} -+ -+ -+// Adds a face on top of existing faceI. -+Foam::label Foam::hexRef8::addFace -+( -+ polyTopoChange& meshMod, -+ const label faceI, -+ const face& newFace, -+ const label own, -+ const label nei -+) const -+{ -+ label patchID, zoneID, zoneFlip; -+ -+ getFaceInfo(faceI, patchID, zoneID, zoneFlip); -+ -+ label newFaceI = -1; -+ -+ if ((nei == -1) || (own < nei)) -+ { -+ // Ordering ok. -+ newFaceI = meshMod.setAction -+ ( -+ polyAddFace -+ ( -+ newFace, // face -+ own, // owner -+ nei, // neighbour -+ -1, // master point -+ -1, // master edge -+ faceI, // master face for addition -+ false, // flux flip -+ patchID, // patch for face -+ zoneID, // zone for face -+ zoneFlip // face zone flip -+ ) -+ ); -+ } -+ else -+ { -+ // Reverse owner/neighbour -+ newFaceI = meshMod.setAction -+ ( -+ polyAddFace -+ ( -+ newFace.reverseFace(), // face -+ nei, // owner -+ own, // neighbour -+ -1, // master point -+ -1, // master edge -+ faceI, // master face for addition -+ false, // flux flip -+ patchID, // patch for face -+ zoneID, // zone for face -+ zoneFlip // face zone flip -+ ) -+ ); -+ } -+ return newFaceI; -+} -+ -+ -+// Adds an internal face from an edge. Assumes orientation correct. -+// Problem is that the face is between four new vertices. So what do we provide -+// as master? The only existing mesh item we have is the edge we have split. -+// Have to be careful in only using it if it has internal faces since otherwise -+// polyMeshMorph will complain (because it cannot generate a sensible mapping -+// for the face) -+Foam::label Foam::hexRef8::addInternalFace -+( -+ polyTopoChange& meshMod, -+ const label meshFaceI, -+ const label meshPointI, -+ const face& newFace, -+ const label own, -+ const label nei -+) const -+{ -+ if (mesh_.isInternalFace(meshFaceI)) -+ { -+ return meshMod.setAction -+ ( -+ polyAddFace -+ ( -+ newFace, // face -+ own, // owner -+ nei, // neighbour -+ -1, // master point -+ -1, // master edge -+ meshFaceI, // master face for addition -+ false, // flux flip -+ -1, // patch for face -+ -1, // zone for face -+ false // face zone flip -+ ) -+ ); -+ } -+ else -+ { -+ // Two choices: -+ // - append (i.e. create out of nothing - will not be mapped) -+ // problem: field does not get mapped. -+ // - inflate from point. -+ // problem: does interpolative mapping which constructs full -+ // volPointInterpolation! -+ -+ // For now create out of nothing -+ -+ return meshMod.setAction -+ ( -+ polyAddFace -+ ( -+ newFace, // face -+ own, // owner -+ nei, // neighbour -+ -1, // master point -+ -1, // master edge -+ -1, // master face for addition -+ false, // flux flip -+ -1, // patch for face -+ -1, // zone for face -+ false // face zone flip -+ ) -+ ); -+ -+ -+ ////- Inflate-from-point: -+ //// Check if point has any internal faces we can use. -+ //label masterPointI = -1; -+ // -+ //const labelList& pFaces = mesh_.pointFaces()[meshPointI]; -+ // -+ //forAll(pFaces, i) -+ //{ -+ // if (mesh_.isInternalFace(pFaces[i])) -+ // { -+ // // meshPoint uses internal faces so ok to inflate from it -+ // masterPointI = meshPointI; -+ // -+ // break; -+ // } -+ //} -+ // -+ //return meshMod.setAction -+ //( -+ // polyAddFace -+ // ( -+ // newFace, // face -+ // own, // owner -+ // nei, // neighbour -+ // masterPointI, // master point -+ // -1, // master edge -+ // -1, // master face for addition -+ // false, // flux flip -+ // -1, // patch for face -+ // -1, // zone for face -+ // false // face zone flip -+ // ) -+ //); -+ } -+} -+ -+ -+// Modifies existing faceI for either new owner/neighbour or new face points. -+void Foam::hexRef8::modFace -+( -+ polyTopoChange& meshMod, -+ const label faceI, -+ const face& newFace, -+ const label own, -+ const label nei -+) const -+{ -+ label patchID, zoneID, zoneFlip; -+ -+ getFaceInfo(faceI, patchID, zoneID, zoneFlip); -+ -+ if -+ ( -+ (own != mesh_.faceOwner()[faceI]) -+ || ( -+ mesh_.isInternalFace(faceI) -+ && (nei != mesh_.faceNeighbour()[faceI]) -+ ) -+ || (newFace != mesh_.faces()[faceI]) -+ ) -+ { -+ if ((nei == -1) || (own < nei)) -+ { -+ meshMod.setAction -+ ( -+ polyModifyFace -+ ( -+ newFace, // modified face -+ faceI, // label of face being modified -+ own, // owner -+ nei, // neighbour -+ false, // face flip -+ patchID, // patch for face -+ false, // remove from zone -+ zoneID, // zone for face -+ zoneFlip // face flip in zone -+ ) -+ ); -+ } -+ else -+ { -+ meshMod.setAction -+ ( -+ polyModifyFace -+ ( -+ newFace.reverseFace(), // modified face -+ faceI, // label of face being modified -+ nei, // owner -+ own, // neighbour -+ false, // face flip -+ patchID, // patch for face -+ false, // remove from zone -+ zoneID, // zone for face -+ zoneFlip // face flip in zone -+ ) -+ ); -+ } -+ } -+} -+ -+ -+// Bit complex way to determine the unrefined edge length. -+Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const -+{ -+ // Determine minimum edge length per refinement level -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ const scalar GREAT2 = sqr(GREAT); -+ -+ label nLevels = gMax(cellLevel_)+1; -+ -+ scalarField typEdgeLenSqr(nLevels, GREAT2); -+ -+ -+ // 1. Look only at edges surrounded by cellLevel cells only. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ { -+ // Per edge the cellLevel of connected cells. -1 if not set, -+ // labelMax if different levels, otherwise levels of connected cells. -+ labelList edgeLevel(mesh_.nEdges(), -1); -+ -+ forAll(cellLevel_, cellI) -+ { -+ const label cLevel = cellLevel_[cellI]; -+ -+ const labelList& cEdges = mesh_.cellEdges()[cellI]; -+ -+ forAll(cEdges, i) -+ { -+ label edgeI = cEdges[i]; -+ -+ if (edgeLevel[edgeI] == -1) -+ { -+ edgeLevel[edgeI] = cLevel; -+ } -+ else if (edgeLevel[edgeI] == labelMax) -+ { -+ // Already marked as on different cellLevels -+ } -+ else if (edgeLevel[edgeI] != cLevel) -+ { -+ edgeLevel[edgeI] = labelMax; -+ } -+ } -+ } -+ -+ // Make sure that edges with different levels on different processors -+ // are also marked. Do the same test (edgeLevel != cLevel) on coupled -+ // edges. -+ syncTools::syncEdgeList -+ ( -+ mesh_, -+ edgeLevel, -+ ifEqEqOp<labelMax>(), -+ labelMin, -+ false // no separation -+ ); -+ -+ // Now use the edgeLevel with a valid value to determine the -+ // length per level. -+ forAll(edgeLevel, edgeI) -+ { -+ const label eLevel = edgeLevel[edgeI]; -+ -+ if (eLevel >= 0 && eLevel < labelMax) -+ { -+ const edge& e = mesh_.edges()[edgeI]; -+ -+ scalar edgeLenSqr = magSqr(e.vec(mesh_.points())); -+ -+ typEdgeLenSqr[eLevel] = min(typEdgeLenSqr[eLevel], edgeLenSqr); -+ } -+ } -+ } -+ -+ // Get the minimum per level over all processors. Note minimum so if -+ // cells are not cubic we use the smallest edge side. -+ Pstream::listCombineGather(typEdgeLenSqr, minEqOp<scalar>()); -+ Pstream::listCombineScatter(typEdgeLenSqr); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::getLevel0EdgeLength() :" -+ << " After phase1: Edgelengths (squared) per refinementlevel:" -+ << typEdgeLenSqr << endl; -+ } -+ -+ -+ // 2. For any levels where we haven't determined a valid length yet -+ // use any surrounding cell level. Here we use the max so we don't -+ // pick up levels between celllevel and higher celllevel (will have -+ // edges sized according to highest celllevel) -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ scalarField maxEdgeLenSqr(nLevels, -GREAT2); -+ -+ forAll(cellLevel_, cellI) -+ { -+ const label cLevel = cellLevel_[cellI]; -+ -+ const labelList& cEdges = mesh_.cellEdges()[cellI]; -+ -+ forAll(cEdges, i) -+ { -+ const edge& e = mesh_.edges()[cEdges[i]]; -+ -+ scalar edgeLenSqr = magSqr(e.vec(mesh_.points())); -+ -+ maxEdgeLenSqr[cLevel] = max(maxEdgeLenSqr[cLevel], edgeLenSqr); -+ } -+ } -+ -+ Pstream::listCombineGather(maxEdgeLenSqr, maxEqOp<scalar>()); -+ Pstream::listCombineScatter(maxEdgeLenSqr); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::getLevel0EdgeLength() :" -+ << " Crappy Edgelengths (squared) per refinementlevel:" -+ << maxEdgeLenSqr << endl; -+ } -+ -+ -+ // 3. Combine the two sets of lengths -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ forAll(typEdgeLenSqr, levelI) -+ { -+ if (typEdgeLenSqr[levelI] == GREAT2) -+ { -+ typEdgeLenSqr[levelI] = maxEdgeLenSqr[levelI]; -+ } -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::getLevel0EdgeLength() :" -+ << " Final Edgelengths (squared) per refinementlevel:" -+ << typEdgeLenSqr << endl; -+ } -+ -+ // Find lowest level present -+ scalar level0Size = -1; -+ -+ forAll(typEdgeLenSqr, levelI) -+ { -+ scalar lenSqr = typEdgeLenSqr[levelI]; -+ -+ if (lenSqr < GREAT2) -+ { -+ level0Size = Foam::sqrt(lenSqr)/(1<<levelI); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::getLevel0EdgeLength() :" -+ << " For level:" << levelI -+ << " found edgeLen:" << level0Size -+ << endl; -+ } -+ break; -+ } -+ } -+ -+ if (level0Size == -1) -+ { -+ FatalErrorIn("hexRef8::getLevel0EdgeLength()") -+ << "Problem : typEdgeLenSqr:" << typEdgeLenSqr << abort(FatalError); -+ } -+ -+ return level0Size; -+} -+ -+ -+// Check whether pointI is an anchor on cellI. -+// If it is not check whether any other point on the face is an anchor cell. -+Foam::label Foam::hexRef8::getAnchorCell -+( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const label cellI, -+ const label faceI, -+ const label pointI -+) const -+{ -+ if (cellAnchorPoints[cellI].size() > 0) -+ { -+ label index = findIndex(cellAnchorPoints[cellI], pointI); -+ -+ if (index != -1) -+ { -+ return cellAddedCells[cellI][index]; -+ } -+ -+ -+ // pointI is not an anchor cell. -+ // Maybe we are already a refined face so check all the face -+ // vertices. -+ const face& f = mesh_.faces()[faceI]; -+ -+ forAll(f, fp) -+ { -+ label index = findIndex(cellAnchorPoints[cellI], f[fp]); -+ -+ if (index != -1) -+ { -+ return cellAddedCells[cellI][index]; -+ } -+ } -+ -+ // Problem. -+ -+ // Pick up points of the cell -+ const labelList cPoints(cellPoints(cellI)); -+ -+ Perr<< "cell:" << cellI << " points:" << endl; -+ forAll(cPoints, i) -+ { -+ label pointI = cPoints[i]; -+ -+ Perr<< " " << pointI << " coord:" << mesh_.points()[pointI] -+ << nl; -+ } -+ Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI] -+ << endl; -+ -+ FatalErrorIn("hexRef8::getAnchorCell") -+ << "Could not find point " << pointI -+ << " in the anchorPoints for cell " << cellI << endl -+ << "Does your original mesh obey the 2:1 constraint and" -+ << " did you use consistentRefinement to make your cells to refine" -+ << " obey this constraint as well?" -+ << abort(FatalError); -+ -+ return -1; -+ } -+ else -+ { -+ return cellI; -+ } -+} -+ -+ -+// Get new owner and neighbour -+void Foam::hexRef8::getFaceNeighbours -+( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const label faceI, -+ const label pointI, -+ -+ label& own, -+ label& nei -+) const -+{ -+ // Is owner split? -+ own = getAnchorCell -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ mesh_.faceOwner()[faceI], -+ faceI, -+ pointI -+ ); -+ -+ if (mesh_.isInternalFace(faceI)) -+ { -+ nei = getAnchorCell -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ mesh_.faceNeighbour()[faceI], -+ faceI, -+ pointI -+ ); -+ } -+ else -+ { -+ nei = -1; -+ } -+} -+ -+ -+// Get point with the lowest pointLevel -+Foam::label Foam::hexRef8::findMinLevel(const labelList& f) const -+{ -+ label minLevel = labelMax; -+ label minFp = -1; -+ -+ forAll(f, fp) -+ { -+ label level = pointLevel_[f[fp]]; -+ -+ if (level < minLevel) -+ { -+ minLevel = level; -+ minFp = fp; -+ } -+ } -+ -+ return minFp; -+} -+ -+ -+// Get point with the highest pointLevel -+Foam::label Foam::hexRef8::findMaxLevel(const labelList& f) const -+{ -+ label maxLevel = labelMin; -+ label maxFp = -1; -+ -+ forAll(f, fp) -+ { -+ label level = pointLevel_[f[fp]]; -+ -+ if (level > maxLevel) -+ { -+ maxLevel = level; -+ maxFp = fp; -+ } -+ } -+ -+ return maxFp; -+} -+ -+ -+Foam::label Foam::hexRef8::countAnchors -+( -+ const labelList& f, -+ const label anchorLevel -+) const -+{ -+ label nAnchors = 0; -+ -+ forAll(f, fp) -+ { -+ if (pointLevel_[f[fp]] <= anchorLevel) -+ { -+ nAnchors++; -+ } -+ } -+ return nAnchors; -+} -+ -+ -+// Find point with certain pointLevel. Skip any higher levels. -+Foam::label Foam::hexRef8::findLevel -+( -+ const face& f, -+ const label startFp, -+ const bool searchForward, -+ const label wantedLevel -+) const -+{ -+ label fp = startFp; -+ -+ forAll(f, i) -+ { -+ label pointI = f[fp]; -+ -+ if (pointLevel_[pointI] < wantedLevel) -+ { -+ FatalErrorIn("hexRef8::findLevel") -+ << "face:" << f -+ << " level:" << IndirectList<label>(pointLevel_, f)() -+ << " startFp:" << startFp -+ << " wantedLevel:" << wantedLevel -+ << abort(FatalError); -+ } -+ else if (pointLevel_[pointI] == wantedLevel) -+ { -+ return fp; -+ } -+ -+ if (searchForward) -+ { -+ fp = f.fcIndex(fp); -+ } -+ else -+ { -+ fp = f.rcIndex(fp); -+ } -+ } -+ -+ FatalErrorIn("hexRef8::findLevel") -+ << "face:" << f -+ << " level:" << IndirectList<label>(pointLevel_, f)() -+ << " startFp:" << startFp -+ << " wantedLevel:" << wantedLevel -+ << abort(FatalError); -+ -+ return -1; -+} -+ -+ -+// Gets cell level such that the face has four points <= level. -+Foam::label Foam::hexRef8::getAnchorLevel(const label faceI) const -+{ -+ const face& f = mesh_.faces()[faceI]; -+ -+ if (f.size() <= 4) -+ { -+ return pointLevel_[f[findMaxLevel(f)]]; -+ } -+ else -+ { -+ label ownLevel = cellLevel_[mesh_.faceOwner()[faceI]]; -+ -+ if (countAnchors(f, ownLevel) == 4) -+ { -+ return ownLevel; -+ } -+ else if (countAnchors(f, ownLevel+1) == 4) -+ { -+ return ownLevel+1; -+ } -+ else -+ { -+ //FatalErrorIn("hexRef8::getAnchorLevel(const label) const") -+ // << "face:" << faceI -+ // << " centre:" << mesh_.faceCentres()[faceI] -+ // << " verts:" << f -+ // << " point levels:" << IndirectList<label>(pointLevel_, f)() -+ // << " own:" << mesh_.faceOwner()[faceI] -+ // << " ownLevel:" << cellLevel_[mesh_.faceOwner()[faceI]] -+ // << abort(FatalError); -+ -+ return -1; -+ } -+ } -+} -+ -+ -+void Foam::hexRef8::checkInternalOrientation -+( -+ polyTopoChange& meshMod, -+ const label cellI, -+ const label faceI, -+ const point& ownPt, -+ const point& neiPt, -+ const face& newFace -+) -+{ -+ face compactFace(identity(newFace.size())); -+ pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)()); -+ -+ vector n(compactFace.normal(compactPoints)); -+ -+ vector dir(neiPt - ownPt); -+ -+ if ((dir & n) < 0) -+ { -+ FatalErrorIn("checkInternalOrientation") -+ << "cell:" << cellI << " old face:" << faceI -+ << " newFace:" << newFace << endl -+ << " coords:" << compactPoints -+ << " ownPt:" << ownPt -+ << " neiPt:" << neiPt -+ << abort(FatalError); -+ } -+ -+ vector fcToOwn(compactFace.centre(compactPoints) - ownPt); -+ -+ scalar s = (fcToOwn&n) / (dir&n); -+ -+ if (s < 0.1 || s > 0.9) -+ { -+ FatalErrorIn("checkInternalOrientation") -+ << "cell:" << cellI << " old face:" << faceI -+ << " newFace:" << newFace << endl -+ << " coords:" << compactPoints -+ << " ownPt:" << ownPt -+ << " neiPt:" << neiPt -+ << " s:" << s -+ << abort(FatalError); -+ } -+} -+ -+ -+void Foam::hexRef8::checkBoundaryOrientation -+( -+ polyTopoChange& meshMod, -+ const label cellI, -+ const label faceI, -+ const point& ownPt, -+ const point& boundaryPt, -+ const face& newFace -+) -+{ -+ face compactFace(identity(newFace.size())); -+ pointField compactPoints(IndirectList<point>(meshMod.points(), newFace)()); -+ -+ vector n(compactFace.normal(compactPoints)); -+ -+ vector dir(boundaryPt - ownPt); -+ -+ if ((dir & n) < 0) -+ { -+ FatalErrorIn("checkBoundaryOrientation") -+ << "cell:" << cellI << " old face:" << faceI -+ << " newFace:" << newFace -+ << " coords:" << compactPoints -+ << " ownPt:" << ownPt -+ << " boundaryPt:" << boundaryPt -+ << abort(FatalError); -+ } -+ -+ vector fcToOwn(compactFace.centre(compactPoints) - ownPt); -+ -+ scalar s = (fcToOwn&dir) / magSqr(dir); -+ -+ if (s < 0.7 || s > 1.3) -+ { -+ WarningIn("checkBoundaryOrientation") -+ << "cell:" << cellI << " old face:" << faceI -+ << " newFace:" << newFace -+ << " coords:" << compactPoints -+ << " ownPt:" << ownPt -+ << " boundaryPt:" << boundaryPt -+ << " s:" << s -+ << endl; -+ //<< abort(FatalError); -+ } -+} -+ -+ -+// If p0 and p1 are existing vertices check if edge is split and insert -+// splitPoint. -+void Foam::hexRef8::insertEdgeSplit -+( -+ const labelList& edgeMidPoint, -+ const label p0, -+ const label p1, -+ DynamicList<label>& verts -+) const -+{ -+ if (p0 < mesh_.nPoints() && p1 < mesh_.nPoints()) -+ { -+ label edgeI = meshTools::findEdge(mesh_, p0, p1); -+ -+ if (edgeI != -1 && edgeMidPoint[edgeI] != -1) -+ { -+ verts.append(edgeMidPoint[edgeI]); -+ } -+ } -+} -+ -+ -+// Internal faces are one per edge between anchor points. So one per midPoint -+// between the anchor points. Here we store the information on the midPoint -+// and if we have enough information: -+// - two anchors -+// - two face mid points -+// we add the face. Note that this routine can get called anywhere from -+// two times (two unrefined faces) to four times (two refined faces) so -+// the first call that adds the information creates the face. -+Foam::label Foam::hexRef8::storeMidPointInfo -+( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const labelList& cellMidPoint, -+ const labelList& edgeMidPoint, -+ const label cellI, -+ const label faceI, -+ const bool faceOrder, -+ const label edgeMidPointI, -+ const label anchorPointI, -+ const label faceMidPointI, -+ -+ Map<edge>& midPointToAnchors, -+ Map<edge>& midPointToFaceMids, -+ polyTopoChange& meshMod -+) const -+{ -+ // See if need to store anchors. -+ -+ bool changed = false; -+ bool haveTwoAnchors = false; -+ -+ Map<edge>::iterator edgeMidFnd = midPointToAnchors.find(edgeMidPointI); -+ -+ if (edgeMidFnd == midPointToAnchors.end()) -+ { -+ midPointToAnchors.insert(edgeMidPointI, edge(anchorPointI, -1)); -+ } -+ else -+ { -+ edge& e = edgeMidFnd(); -+ -+ if (anchorPointI != e[0]) -+ { -+ if (e[1] == -1) -+ { -+ e[1] = anchorPointI; -+ changed = true; -+ } -+ } -+ -+ if (e[0] != -1 && e[1] != -1) -+ { -+ haveTwoAnchors = true; -+ } -+ } -+ -+ bool haveTwoFaceMids = false; -+ -+ Map<edge>::iterator faceMidFnd = midPointToFaceMids.find(edgeMidPointI); -+ -+ if (faceMidFnd == midPointToFaceMids.end()) -+ { -+ midPointToFaceMids.insert(edgeMidPointI, edge(faceMidPointI, -1)); -+ } -+ else -+ { -+ edge& e = faceMidFnd(); -+ -+ if (faceMidPointI != e[0]) -+ { -+ if (e[1] == -1) -+ { -+ e[1] = faceMidPointI; -+ changed = true; -+ } -+ } -+ -+ if (e[0] != -1 && e[1] != -1) -+ { -+ haveTwoFaceMids = true; -+ } -+ } -+ -+ // Check if this call of storeMidPointInfo is the one that completed all -+ // the nessecary information. -+ -+ if (changed && haveTwoAnchors && haveTwoFaceMids) -+ { -+ const edge& anchors = midPointToAnchors[edgeMidPointI]; -+ const edge& faceMids = midPointToFaceMids[edgeMidPointI]; -+ -+ label otherFaceMidPointI = faceMids.otherVertex(faceMidPointI); -+ -+ // Create face consistent with anchorI being the owner. -+ // Note that the edges between the edge mid point and the face mids -+ // might be marked for splitting. Note that these edge splits cannot -+ // be between cellMid and face mids. -+ -+ DynamicList<label> newFaceVerts(4); -+ if (faceOrder == (mesh_.faceOwner()[faceI] == cellI)) -+ { -+ newFaceVerts.append(faceMidPointI); -+ -+ // Check & insert edge split if any -+ insertEdgeSplit -+ ( -+ edgeMidPoint, -+ faceMidPointI, // edge between faceMid and -+ edgeMidPointI, // edgeMid -+ newFaceVerts -+ ); -+ -+ newFaceVerts.append(edgeMidPointI); -+ -+ insertEdgeSplit -+ ( -+ edgeMidPoint, -+ edgeMidPointI, -+ otherFaceMidPointI, -+ newFaceVerts -+ ); -+ -+ newFaceVerts.append(otherFaceMidPointI); -+ newFaceVerts.append(cellMidPoint[cellI]); -+ } -+ else -+ { -+ newFaceVerts.append(otherFaceMidPointI); -+ -+ insertEdgeSplit -+ ( -+ edgeMidPoint, -+ otherFaceMidPointI, -+ edgeMidPointI, -+ newFaceVerts -+ ); -+ -+ newFaceVerts.append(edgeMidPointI); -+ -+ insertEdgeSplit -+ ( -+ edgeMidPoint, -+ edgeMidPointI, -+ faceMidPointI, -+ newFaceVerts -+ ); -+ -+ newFaceVerts.append(faceMidPointI); -+ newFaceVerts.append(cellMidPoint[cellI]); -+ } -+ -+ face newFace; -+ newFace.transfer(newFaceVerts.shrink()); -+ newFaceVerts.clear(); -+ -+ label anchorCell0 = getAnchorCell -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ cellI, -+ faceI, -+ anchorPointI -+ ); -+ label anchorCell1 = getAnchorCell -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ cellI, -+ faceI, -+ anchors.otherVertex(anchorPointI) -+ ); -+ -+ -+ label own, nei; -+ point ownPt, neiPt; -+ -+ if (anchorCell0 < anchorCell1) -+ { -+ own = anchorCell0; -+ nei = anchorCell1; -+ -+ ownPt = mesh_.points()[anchorPointI]; -+ neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)]; -+ -+ } -+ else -+ { -+ own = anchorCell1; -+ nei = anchorCell0; -+ newFace = newFace.reverseFace(); -+ -+ ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)]; -+ neiPt = mesh_.points()[anchorPointI]; -+ } -+ -+ if (debug) -+ { -+ point ownPt, neiPt; -+ -+ if (anchorCell0 < anchorCell1) -+ { -+ ownPt = mesh_.points()[anchorPointI]; -+ neiPt = mesh_.points()[anchors.otherVertex(anchorPointI)]; -+ } -+ else -+ { -+ ownPt = mesh_.points()[anchors.otherVertex(anchorPointI)]; -+ neiPt = mesh_.points()[anchorPointI]; -+ } -+ -+ checkInternalOrientation -+ ( -+ meshMod, -+ cellI, -+ faceI, -+ ownPt, -+ neiPt, -+ newFace -+ ); -+ } -+ -+ return addInternalFace -+ ( -+ meshMod, -+ faceI, -+ anchorPointI, -+ newFace, -+ own, -+ nei -+ ); -+ } -+ else -+ { -+ return -1; -+ } -+} -+ -+ -+// Creates all the 12 internal faces for cellI. -+void Foam::hexRef8::createInternalFaces -+( -+ const labelListList& cellAnchorPoints, -+ const labelListList& cellAddedCells, -+ const labelList& cellMidPoint, -+ const labelList& faceMidPoint, -+ const labelList& faceAnchorLevel, -+ const labelList& edgeMidPoint, -+ const label cellI, -+ -+ polyTopoChange& meshMod -+) const -+{ -+ // Find in every face the cellLevel+1 points (from edge subdivision) -+ // and the anchor points. -+ -+ const cell& cFaces = mesh_.cells()[cellI]; -+ const label cLevel = cellLevel_[cellI]; -+ -+ // From edge mid to anchor points -+ Map<edge> midPointToAnchors(24); -+ // From edge mid to face mids -+ Map<edge> midPointToFaceMids(24); -+ -+ // Running count of number of internal faces added so far. -+ label nFacesAdded = 0; -+ -+ forAll(cFaces, i) -+ { -+ label faceI = cFaces[i]; -+ -+ const face& f = mesh_.faces()[faceI]; -+ const labelList& fEdges = mesh_.faceEdges()[faceI]; -+ -+ // We are on the cellI side of face f. The face will have 1 or 4 -+ // cLevel points and lots of higher numbered ones. -+ -+ label faceMidPointI = -1; -+ -+ label nAnchors = countAnchors(f, cLevel); -+ -+ if (nAnchors == 1) -+ { -+ // Only one anchor point. So the other side of the face has already -+ // been split using cLevel+1 and cLevel+2 points. -+ -+ // Find the one anchor. -+ label anchorFp = -1; -+ -+ forAll(f, fp) -+ { -+ if (pointLevel_[f[fp]] <= cLevel) -+ { -+ anchorFp = fp; -+ break; -+ } -+ } -+ -+ // Now the face mid point is the second cLevel+1 point -+ label edgeMid = findLevel(f, f.fcIndex(anchorFp), true, cLevel+1); -+ label faceMid = findLevel(f, f.fcIndex(edgeMid), true, cLevel+1); -+ -+ faceMidPointI = f[faceMid]; -+ } -+ else if (nAnchors == 4) -+ { -+ // There is no face middle yet but the face will be marked for -+ // splitting. -+ -+ faceMidPointI = faceMidPoint[faceI]; -+ } -+ else -+ { -+ FatalErrorIn("createInternalFaces") -+ << "nAnchors:" << nAnchors -+ << " faceI:" << faceI -+ << abort(FatalError); -+ } -+ -+ -+ -+ // Now loop over all the anchors (might be just one) and store -+ // the edge mids connected to it. storeMidPointInfo will collect -+ // all the info and combine it all. -+ -+ forAll(f, fp0) -+ { -+ label point0 = f[fp0]; -+ -+ if (pointLevel_[point0] <= cLevel) -+ { -+ // Anchor. -+ -+ // Walk forward -+ // ~~~~~~~~~~~~ -+ // to cLevel+1 or edgeMidPoint of this level. -+ -+ -+ label edgeMidPointI = -1; -+ -+ label fp1 = f.fcIndex(fp0); -+ -+ if (pointLevel_[f[fp1]] <= cLevel) -+ { -+ // Anchor. Edge will be split. -+ label edgeI = fEdges[fp0]; -+ -+ edgeMidPointI = edgeMidPoint[edgeI]; -+ -+ if (edgeMidPointI == -1) -+ { -+ const labelList cPoints(cellPoints(cellI)); -+ -+ FatalErrorIn("createInternalFaces") -+ << "cell:" << cellI << " cLevel:" << cLevel -+ << " cell points:" << cPoints -+ << " pointLevel:" -+ << IndirectList<label>(pointLevel_, cPoints)() -+ << " face:" << faceI -+ << " f:" << f -+ << " pointLevel:" -+ << IndirectList<label>(pointLevel_, f)() -+ << " faceAnchorLevel:" << faceAnchorLevel[faceI] -+ << " faceMidPoint:" << faceMidPoint[faceI] -+ << " faceMidPointI:" << faceMidPointI -+ << " fp:" << fp0 -+ << abort(FatalError); -+ } -+ } -+ else -+ { -+ // Search forward in face to clevel+1 -+ label edgeMid = findLevel(f, fp1, true, cLevel+1); -+ -+ edgeMidPointI = f[edgeMid]; -+ } -+ -+ label newFaceI = storeMidPointInfo -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ cellMidPoint, -+ edgeMidPoint, -+ -+ cellI, -+ faceI, -+ true, // mid point after anchor -+ edgeMidPointI, // edgemid -+ point0, // anchor -+ faceMidPointI, -+ -+ midPointToAnchors, -+ midPointToFaceMids, -+ meshMod -+ ); -+ -+ if (newFaceI != -1) -+ { -+ nFacesAdded++; -+ -+ if (nFacesAdded == 12) -+ { -+ break; -+ } -+ } -+ -+ -+ -+ // Walk backward -+ // ~~~~~~~~~~~~~ -+ -+ label fpMin1 = f.rcIndex(fp0); -+ -+ if (pointLevel_[f[fpMin1]] <= cLevel) -+ { -+ // Anchor. Edge will be split. -+ label edgeI = fEdges[fpMin1]; -+ -+ edgeMidPointI = edgeMidPoint[edgeI]; -+ -+ if (edgeMidPointI == -1) -+ { -+ const labelList cPoints(cellPoints(cellI)); -+ -+ FatalErrorIn("createInternalFaces") -+ << "cell:" << cellI << " cLevel:" << cLevel -+ << " cell points:" << cPoints -+ << " pointLevel:" -+ << IndirectList<label>(pointLevel_, cPoints)() -+ << " face:" << faceI -+ << " f:" << f -+ << " pointLevel:" -+ << IndirectList<label>(pointLevel_, f)() -+ << " faceAnchorLevel:" << faceAnchorLevel[faceI] -+ << " faceMidPoint:" << faceMidPoint[faceI] -+ << " faceMidPointI:" << faceMidPointI -+ << " fp:" << fp0 -+ << abort(FatalError); -+ } -+ } -+ else -+ { -+ // Search back to clevel+1 -+ label edgeMid = findLevel(f, fpMin1, false, cLevel+1); -+ -+ edgeMidPointI = f[edgeMid]; -+ } -+ -+ newFaceI = storeMidPointInfo -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ cellMidPoint, -+ edgeMidPoint, -+ -+ cellI, -+ faceI, -+ false, // mid point before anchor -+ edgeMidPointI, // edgemid -+ point0, // anchor -+ faceMidPointI, -+ -+ midPointToAnchors, -+ midPointToFaceMids, -+ meshMod -+ ); -+ -+ if (newFaceI != -1) -+ { -+ nFacesAdded++; -+ -+ if (nFacesAdded == 12) -+ { -+ break; -+ } -+ } -+ } // done anchor -+ } // done face -+ -+ if (nFacesAdded == 12) -+ { -+ break; -+ } -+ } -+} -+ -+ -+void Foam::hexRef8::walkFaceToMid -+( -+ const labelList& edgeMidPoint, -+ const label cLevel, -+ const label faceI, -+ const label startFp, -+ DynamicList<label>& faceVerts -+) const -+{ -+ const face& f = mesh_.faces()[faceI]; -+ const labelList& fEdges = mesh_.faceEdges()[faceI]; -+ -+ label fp = startFp; -+ -+ // Starting from fp store all (1 or 2) vertices until where the face -+ // gets split -+ -+ while (true) -+ { -+ if (edgeMidPoint[fEdges[fp]] >= 0) -+ { -+ faceVerts.append(edgeMidPoint[fEdges[fp]]); -+ } -+ -+ fp = f.fcIndex(fp); -+ -+ if (pointLevel_[f[fp]] <= cLevel) -+ { -+ // Next anchor. Have already append split point on edge in code -+ // above. -+ return; -+ } -+ else if (pointLevel_[f[fp]] == cLevel+1) -+ { -+ // Mid level -+ faceVerts.append(f[fp]); -+ -+ return; -+ } -+ else if (pointLevel_[f[fp]] == cLevel+2) -+ { -+ // Store and continue to cLevel+1. -+ faceVerts.append(f[fp]); -+ } -+ } -+} -+ -+ -+// Same as walkFaceToMid but now walk back. -+void Foam::hexRef8::walkFaceFromMid -+( -+ const labelList& edgeMidPoint, -+ const label cLevel, -+ const label faceI, -+ const label startFp, -+ DynamicList<label>& faceVerts -+) const -+{ -+ const face& f = mesh_.faces()[faceI]; -+ const labelList& fEdges = mesh_.faceEdges()[faceI]; -+ -+ label fp = f.rcIndex(startFp); -+ -+ while (true) -+ { -+ if (pointLevel_[f[fp]] <= cLevel) -+ { -+ // anchor. -+ break; -+ } -+ else if (pointLevel_[f[fp]] == cLevel+1) -+ { -+ // Mid level -+ faceVerts.append(f[fp]); -+ break; -+ } -+ else if (pointLevel_[f[fp]] == cLevel+2) -+ { -+ // Continue to cLevel+1. -+ } -+ fp = f.rcIndex(fp); -+ } -+ -+ // Store -+ while (true) -+ { -+ if (edgeMidPoint[fEdges[fp]] >= 0) -+ { -+ faceVerts.append(edgeMidPoint[fEdges[fp]]); -+ } -+ -+ fp = f.fcIndex(fp); -+ -+ if (fp == startFp) -+ { -+ break; -+ } -+ faceVerts.append(f[fp]); -+ } -+} -+ -+ -+// Updates refineCell (cells marked for refinement) so across all faces -+// there will be 2:1 consistency after refinement. -+Foam::label Foam::hexRef8::faceConsistentRefinement -+( -+ const bool maxSet, -+ PackedList<1>& refineCell -+) const -+{ -+ label nChanged = 0; -+ -+ // Internal faces. -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own] + refineCell.get(own); -+ -+ label nei = mesh_.faceNeighbour()[faceI]; -+ label neiLevel = cellLevel_[nei] + refineCell.get(nei); -+ -+ if (ownLevel > (neiLevel+1)) -+ { -+ if (maxSet) -+ { -+ refineCell.set(nei, 1); -+ } -+ else -+ { -+ refineCell.set(own, 0); -+ } -+ nChanged++; -+ } -+ else if (neiLevel > (ownLevel+1)) -+ { -+ if (maxSet) -+ { -+ refineCell.set(own, 1); -+ } -+ else -+ { -+ refineCell.set(nei, 0); -+ } -+ nChanged++; -+ } -+ } -+ -+ -+ // Coupled faces. Swap owner level to get neighbouring cell level. -+ // (only boundary faces of neiLevel used) -+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ -+ neiLevel[i] = cellLevel_[own] + refineCell.get(own); -+ } -+ -+ // Swap to neighbour -+ syncTools::swapBoundaryFaceList(mesh_, neiLevel, false); -+ -+ // Now we have neighbour value see which cells need refinement -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ label ownLevel = cellLevel_[own] + refineCell.get(own); -+ -+ if (ownLevel > (neiLevel[i]+1)) -+ { -+ if (!maxSet) -+ { -+ refineCell.set(own, 0); -+ nChanged++; -+ } -+ } -+ else if (neiLevel[i] > (ownLevel+1)) -+ { -+ if (maxSet) -+ { -+ refineCell.set(own, 1); -+ nChanged++; -+ } -+ } -+ } -+ -+ return nChanged; -+} -+ -+ -+// Debug: check if wanted refinement is compatible with 2:1 -+void Foam::hexRef8::checkWantedRefinementLevels -+( -+ const labelList& cellsToRefine -+) const -+{ -+ PackedList<1> refineCell(mesh_.nCells(), 0); -+ forAll(cellsToRefine, i) -+ { -+ refineCell.set(cellsToRefine[i], 1); -+ } -+ -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own] + refineCell.get(own); -+ -+ label nei = mesh_.faceNeighbour()[faceI]; -+ label neiLevel = cellLevel_[nei] + refineCell.get(nei); -+ -+ if (mag(ownLevel-neiLevel) > 1) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const labelList&)" -+ ) << "cell:" << own -+ << " current level:" << cellLevel_[own] -+ << " level after refinement:" << ownLevel -+ << nl -+ << "neighbour cell:" << nei -+ << " current level:" << cellLevel_[nei] -+ << " level after refinement:" << neiLevel -+ << nl -+ << "which does not satisfy 2:1 constraints anymore." -+ << abort(FatalError); -+ } -+ } -+ -+ // Coupled faces. Swap owner level to get neighbouring cell level. -+ // (only boundary faces of neiLevel used) -+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ -+ neiLevel[i] = cellLevel_[own] + refineCell.get(own); -+ } -+ -+ // Swap to neighbour -+ syncTools::swapBoundaryFaceList(mesh_, neiLevel, false); -+ -+ // Now we have neighbour value see which cells need refinement -+ forAll(neiLevel, i) -+ { -+ label faceI = i + mesh_.nInternalFaces(); -+ -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own] + refineCell.get(own); -+ -+ if (mag(ownLevel - neiLevel[i]) > 1) -+ { -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const labelList&)" -+ ) << "Celllevel does not satisfy 2:1 constraint." -+ << " On coupled face " -+ << faceI -+ << " on patch " << patchI << " " -+ << mesh_.boundaryMesh()[patchI].name() -+ << " owner cell " << own -+ << " current level:" << cellLevel_[own] -+ << " level after refinement:" << ownLevel -+ << nl -+ << " (coupled) neighbour cell will get refinement " -+ << neiLevel[i] -+ << abort(FatalError); -+ } -+ } -+} -+ -+ -+// Set instance for mesh files -+void Foam::hexRef8::setInstance(const fileName& inst) -+{ -+ if (debug) -+ { -+ Pout<< "hexRef8::setInstance(const fileName& inst) : " -+ << "Resetting file instance to " << inst << endl; -+ } -+ -+ cellLevel_.instance() = inst; -+ pointLevel_.instance() = inst; -+ history_.instance() = inst; -+} -+ -+ -+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -+ -+// Construct from mesh, read refinement data -+Foam::hexRef8::hexRef8(const polyMesh& mesh) -+: -+ mesh_(mesh), -+ cellLevel_ -+ ( -+ IOobject -+ ( -+ "cellLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::READ_IF_PRESENT, -+ IOobject::AUTO_WRITE -+ ), -+ labelList(mesh_.nCells(), 0) -+ ), -+ pointLevel_ -+ ( -+ IOobject -+ ( -+ "pointLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::READ_IF_PRESENT, -+ IOobject::AUTO_WRITE -+ ), -+ labelList(mesh_.nPoints(), 0) -+ ), -+ level0Edge_(getLevel0EdgeLength()), -+ history_ -+ ( -+ IOobject -+ ( -+ "refinementHistory", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::READ_IF_PRESENT, -+ IOobject::AUTO_WRITE -+ ), -+ mesh_.nCells() // All cells visible if could not be read -+ ), -+ faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible -+ savedPointLevel_(0), -+ savedCellLevel_(0) -+{ -+ if (history_.active() && history_.visibleCells().size() != mesh_.nCells()) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::hexRef8(const polyMesh&)" -+ ) << "History enabled but number of visible cells in it " -+ << history_.visibleCells().size() -+ << " is not equal to the number of cells in the mesh " -+ << mesh_.nCells() -+ << abort(FatalError); -+ } -+ -+ // Check refinement levels for consistency -+ checkRefinementLevels(-1, labelList(0)); -+ -+ -+ // Check initial mesh for consistency -+ -+ //if (debug) -+ { -+ checkMesh(); -+ } -+} -+ -+ -+// Construct from components -+Foam::hexRef8::hexRef8 -+( -+ const polyMesh& mesh, -+ const labelList& cellLevel, -+ const labelList& pointLevel, -+ const refinementHistory& history -+) -+: -+ mesh_(mesh), -+ cellLevel_ -+ ( -+ IOobject -+ ( -+ "cellLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ cellLevel -+ ), -+ pointLevel_ -+ ( -+ IOobject -+ ( -+ "pointLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ pointLevel -+ ), -+ level0Edge_(getLevel0EdgeLength()), -+ history_ -+ ( -+ IOobject -+ ( -+ "refinementHistory", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ history -+ ), -+ faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible -+ savedPointLevel_(0), -+ savedCellLevel_(0) -+{ -+ if (history_.active() && history_.visibleCells().size() != mesh_.nCells()) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::hexRef8(const polyMesh&, const labelList&" -+ ", const labelList&, const refinementHistory&)" -+ ) << "History enabled but number of visible cells in it " -+ << history_.visibleCells().size() -+ << " is not equal to the number of cells in the mesh " -+ << mesh_.nCells() -+ << abort(FatalError); -+ } -+ -+ // Check refinement levels for consistency -+ checkRefinementLevels(-1, labelList(0)); -+ -+ -+ // Check initial mesh for consistency -+ -+ //if (debug) -+ { -+ checkMesh(); -+ } -+} -+ -+ -+// Construct from components -+Foam::hexRef8::hexRef8 -+( -+ const polyMesh& mesh, -+ const labelList& cellLevel, -+ const labelList& pointLevel -+) -+: -+ mesh_(mesh), -+ cellLevel_ -+ ( -+ IOobject -+ ( -+ "cellLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ cellLevel -+ ), -+ pointLevel_ -+ ( -+ IOobject -+ ( -+ "pointLevel", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ pointLevel -+ ), -+ level0Edge_(getLevel0EdgeLength()), -+ history_ -+ ( -+ IOobject -+ ( -+ "refinementHistory", -+ mesh_.facesInstance(), -+ polyMesh::meshSubDir, -+ mesh_, -+ IOobject::NO_READ, -+ IOobject::AUTO_WRITE -+ ), -+ List<refinementHistory::splitCell8>(0), -+ labelList(0) -+ ), -+ faceRemover_(mesh_, GREAT), // merge boundary faces wherever possible -+ savedPointLevel_(0), -+ savedCellLevel_(0) -+{ -+ // Check refinement levels for consistency -+ checkRefinementLevels(-1, labelList(0)); -+ -+ // Check initial mesh for consistency -+ -+ //if (debug) -+ { -+ checkMesh(); -+ } -+} -+ -+ -+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -+ -+//- Get points of a cell (without using cellPoints addressing) -+Foam::labelList Foam::hexRef8::cellPoints(const label cellI) const -+{ -+ // Pick up points of the cell -+ const cell& cFaces = mesh_.cells()[cellI]; -+ -+ labelHashSet cPoints(4*cFaces.size()); -+ -+ forAll(cFaces, i) -+ { -+ const face& f = mesh_.faces()[cFaces[i]]; -+ -+ forAll(f, fp) -+ { -+ cPoints.insert(f[fp]); -+ } -+ } -+ return cPoints.toc(); -+} -+ -+ -+Foam::labelList Foam::hexRef8::consistentRefinement -+( -+ const labelList& cellsToRefine, -+ const bool maxSet -+) const -+{ -+ // Loop, modifying cellsToRefine, until no more changes to due to 2:1 -+ // conflicts. -+ // maxSet = false : unselect cells to refine -+ // maxSet = true : select cells to refine -+ -+ // Go to straight boolList. -+ PackedList<1> refineCell(mesh_.nCells(), 0); -+ forAll(cellsToRefine, i) -+ { -+ refineCell.set(cellsToRefine[i], 1); -+ } -+ -+ while (true) -+ { -+ label nChanged = faceConsistentRefinement(maxSet, refineCell); -+ -+ reduce(nChanged, sumOp<label>()); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentRefinement : Changed " << nChanged -+ << " refinement levels due to 2:1 conflicts." -+ << endl; -+ } -+ -+ if (nChanged == 0) -+ { -+ break; -+ } -+ } -+ -+ -+ // Convert back to labelList. -+ label nRefined = 0; -+ -+ forAll(refineCell, cellI) -+ { -+ if (refineCell.get(cellI) == 1) -+ { -+ nRefined++; -+ } -+ } -+ -+ labelList newCellsToRefine(nRefined); -+ nRefined = 0; -+ -+ forAll(refineCell, cellI) -+ { -+ if (refineCell.get(cellI) == 1) -+ { -+ newCellsToRefine[nRefined++] = cellI; -+ } -+ } -+ -+ if (debug) -+ { -+ checkWantedRefinementLevels(newCellsToRefine); -+ } -+ -+ return newCellsToRefine; -+} -+ -+ -+// Given a list of cells to refine determine additional cells to refine -+// such that the overall refinement: -+// - satisfies maxFaceDiff (e.g. 2:1) across neighbouring faces -+// - satisfies maxPointDiff (e.g. 4:1) across selected point connected -+// cells. This is used to ensure that e.g. cells on the surface are not -+// point connected to cells which are 8 times smaller. -+Foam::labelList Foam::hexRef8::consistentSlowRefinement -+( -+ const label maxFaceDiff, -+ const labelList& cellsToRefine, -+ const labelList& facesToCheck, -+ const label maxPointDiff, -+ const labelList& pointsToCheck -+) const -+{ -+ const labelList& faceOwner = mesh_.faceOwner(); -+ const labelList& faceNeighbour = mesh_.faceNeighbour(); -+ -+ -+ if (maxFaceDiff <= 0) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement" -+ "(const label, const labelList&, const labelList&" -+ ", const label, const labelList&)" -+ ) << "Illegal maxFaceDiff " << maxFaceDiff << nl -+ << "Value should be >= 1" << exit(FatalError); -+ } -+ -+ -+ // Bit tricky. Say we want a distance of three cells between two -+ // consecutive refinement levels. This is done by using FaceCellWave to -+ // transport out the new refinement level. It gets decremented by one -+ // every cell it crosses so if we initialize it to maxFaceDiff -+ // we will get a field everywhere that tells us whether an unselected cell -+ // needs refining as well. -+ -+ -+ // Initial information about (distance to) cellLevel on all cells -+ List<refinementData> allCellInfo(mesh_.nCells()); -+ -+ // Initial information about (distance to) cellLevel on all faces -+ List<refinementData> allFaceInfo(mesh_.nFaces()); -+ -+ forAll(allCellInfo, cellI) -+ { -+ // maxFaceDiff since refinementData counts both -+ // faces and cells. -+ allCellInfo[cellI] = refinementData -+ ( -+ maxFaceDiff*(cellLevel_[cellI]+1),// when cell is to be refined -+ maxFaceDiff*cellLevel_[cellI] // current level -+ ); -+ } -+ -+ // Cells to be refined will have cellLevel+1 -+ forAll(cellsToRefine, i) -+ { -+ label cellI = cellsToRefine[i]; -+ -+ allCellInfo[cellI].count() = allCellInfo[cellI].refinementCount(); -+ } -+ -+ -+ // Labels of seed faces -+ DynamicList<label> seedFaces(mesh_.nFaces()/100); -+ // refinementLevel data on seed faces -+ DynamicList<refinementData> seedFacesInfo(mesh_.nFaces()/100); -+ -+ -+ // Additional buffer layer thickness by changing initial count. Usually -+ // this happens on boundary faces. Bit tricky. Use allFaceInfo to mark -+ // off thus marked faces so they're skipped in the next loop. -+ forAll(facesToCheck, i) -+ { -+ label faceI = facesToCheck[i]; -+ -+ if (allFaceInfo[faceI].valid()) -+ { -+ // Can only occur if face has already gone through loop below. -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement" -+ "(const label, const labelList&, const labelList&" -+ ", const label, const labelList&)" -+ ) << "Argument facesToCheck seems to have duplicate entries!" -+ << endl -+ << "face:" << faceI << " occurs at positions " -+ << findIndices(facesToCheck, faceI) -+ << abort(FatalError); -+ } -+ -+ -+ const refinementData& ownData = allCellInfo[faceOwner[faceI]]; -+ -+ if (mesh_.isInternalFace(faceI)) -+ { -+ // Seed face if neighbouring cell (after possible refinement) -+ // will be refined one more than the current owner or neighbour. -+ -+ const refinementData& neiData = allCellInfo[faceNeighbour[faceI]]; -+ -+ label faceCount; -+ label faceRefineCount; -+ if (neiData.count() > ownData.count()) -+ { -+ faceCount = neiData.count() + maxFaceDiff; -+ faceRefineCount = faceCount + maxFaceDiff; -+ } -+ else -+ { -+ faceCount = ownData.count() + maxFaceDiff; -+ faceRefineCount = faceCount + maxFaceDiff; -+ } -+ -+ seedFaces.append(faceI); -+ seedFacesInfo.append -+ ( -+ refinementData -+ ( -+ faceRefineCount, -+ faceCount -+ ) -+ ); -+ allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1]; -+ } -+ else -+ { -+ label faceCount = ownData.count() + maxFaceDiff; -+ label faceRefineCount = faceCount + maxFaceDiff; -+ -+ seedFaces.append(faceI); -+ seedFacesInfo.append -+ ( -+ refinementData -+ ( -+ faceRefineCount, -+ faceCount -+ ) -+ ); -+ allFaceInfo[faceI] = seedFacesInfo[seedFacesInfo.size()-1]; -+ } -+ } -+ -+ -+ // Just seed with all faces inbetween different refinement levels for now -+ // (alternatively only seed faces on cellsToRefine but that gives problems -+ // if no cells to refine) -+ forAll(faceNeighbour, faceI) -+ { -+ // Check if face already handled in loop above -+ if (!allFaceInfo[faceI].valid()) -+ { -+ label own = faceOwner[faceI]; -+ label nei = faceNeighbour[faceI]; -+ -+ // Seed face with transported data from highest cell. -+ -+ if (allCellInfo[own].count() > allCellInfo[nei].count()) -+ { -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ own, -+ allCellInfo[own], -+ FaceCellWave<refinementData>::propagationTol() -+ ); -+ seedFaces.append(faceI); -+ seedFacesInfo.append(allFaceInfo[faceI]); -+ } -+ else if (allCellInfo[own].count() < allCellInfo[nei].count()) -+ { -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ nei, -+ allCellInfo[nei], -+ FaceCellWave<refinementData>::propagationTol() -+ ); -+ seedFaces.append(faceI); -+ seedFacesInfo.append(allFaceInfo[faceI]); -+ } -+ } -+ } -+ -+ // Seed all boundary faces with owner value. This is to make sure that -+ // they are visited (probably only important for coupled faces since -+ // these need to be visited from both sides) -+ for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++) -+ { -+ // Check if face already handled in loop above -+ if (!allFaceInfo[faceI].valid()) -+ { -+ label own = faceOwner[faceI]; -+ -+ // Seed face with transported data from owner. -+ refinementData faceData; -+ faceData.updateFace -+ ( -+ mesh_, -+ faceI, -+ own, -+ allCellInfo[own], -+ FaceCellWave<refinementData>::propagationTol() -+ ); -+ seedFaces.append(faceI); -+ seedFacesInfo.append(faceData); -+ } -+ } -+ -+ -+ // face-cell-face transport engine -+ FaceCellWave<refinementData> levelCalc -+ ( -+ mesh_, -+ allFaceInfo, -+ allCellInfo -+ ); -+ -+ while(true) -+ { -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentSlowRefinement : Seeded " -+ << seedFaces.size() << " faces between cells with different" -+ << " refinement level." << endl; -+ } -+ -+ // Set seed faces -+ levelCalc.setFaceInfo(seedFaces.shrink(), seedFacesInfo.shrink()); -+ seedFaces.clear(); -+ seedFacesInfo.clear(); -+ -+ // Iterate until no change. Now 2:1 face difference should be satisfied -+ levelCalc.iterate(mesh_.globalData().nTotalFaces()); -+ -+ -+ // Now check point-connected cells (face-connected cells already ok): -+ // - get per point max of connected cells -+ // - sync across coupled points -+ // - check cells against above point max -+ -+ if (maxPointDiff == -1) -+ { -+ // No need to do any point checking. -+ break; -+ } -+ -+ // Determine per point the max cell level. (done as count, not -+ // as cell level purely for ease) -+ labelList maxPointCount(mesh_.nPoints(), 0); -+ -+ const labelListList& pointCells = mesh_.pointCells(); -+ -+ forAll(pointCells, pointI) -+ { -+ label& pLevel = maxPointCount[pointI]; -+ -+ const labelList& pCells = pointCells[pointI]; -+ -+ forAll(pCells, i) -+ { -+ pLevel = max(pLevel, allCellInfo[pCells[i]].count()); -+ } -+ } -+ -+ // Sync maxPointCount to neighbour -+ syncTools::syncPointList -+ ( -+ mesh_, -+ maxPointCount, -+ maxEqOp<label>(), -+ labelMin, // null value -+ false // no separation -+ ); -+ -+ // Update allFaceInfo from maxPointCount for all points to check -+ // (usually on boundary faces) -+ -+ // Per face the new refinement data -+ Map<refinementData> changedFacesInfo(pointsToCheck.size()); -+ -+ forAll(pointsToCheck, i) -+ { -+ label pointI = pointsToCheck[i]; -+ -+ // Loop over all cells using the point and check whether their -+ // refinement level is much less than the maximum. -+ -+ const labelList& pCells = pointCells[pointI]; -+ -+ forAll(pCells, pCellI) -+ { -+ label cellI = pCells[pCellI]; -+ -+ refinementData& cellInfo = allCellInfo[cellI]; -+ -+ if -+ ( -+ !cellInfo.isRefined() -+ && ( -+ maxPointCount[pointI] -+ > cellInfo.count() + maxFaceDiff*maxPointDiff -+ ) -+ ) -+ { -+ // Mark cell for refinement -+ cellInfo.count() = cellInfo.refinementCount(); -+ -+ // Insert faces of cell as seed faces. -+ const cell& cFaces = mesh_.cells()[cellI]; -+ -+ forAll(cFaces, cFaceI) -+ { -+ label faceI = cFaces[cFaceI]; -+ -+ refinementData faceData; -+ faceData.updateFace -+ ( -+ mesh_, -+ faceI, -+ cellI, -+ cellInfo, -+ FaceCellWave<refinementData>::propagationTol() -+ ); -+ -+ if (faceData.count() > allFaceInfo[faceI].count()) -+ { -+ changedFacesInfo.insert(faceI, faceData); -+ } -+ } -+ } -+ } -+ } -+ -+ label nChanged = changedFacesInfo.size(); -+ reduce(nChanged, sumOp<label>()); -+ -+ if (nChanged == 0) -+ { -+ break; -+ } -+ -+ -+ // Transfer into seedFaces, seedFacesInfo -+ seedFaces.setSize(changedFacesInfo.size()); -+ seedFacesInfo.setSize(changedFacesInfo.size()); -+ -+ forAllConstIter(Map<refinementData>, changedFacesInfo, iter) -+ { -+ seedFaces.append(iter.key()); -+ seedFacesInfo.append(iter()); -+ } -+ } -+ -+ -+ if (debug) -+ { -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = -+ cellLevel_[own] -+ + (allCellInfo[own].isRefined() ? 1 : 0); -+ -+ label nei = mesh_.faceNeighbour()[faceI]; -+ label neiLevel = -+ cellLevel_[nei] -+ + (allCellInfo[nei].isRefined() ? 1 : 0); -+ -+ if (mag(ownLevel-neiLevel) > 1) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement" -+ ) << "cell:" << own -+ << " current level:" << cellLevel_[own] -+ << " current refData:" << allCellInfo[own] -+ << " level after refinement:" << ownLevel -+ << nl -+ << "neighbour cell:" << nei -+ << " current level:" << cellLevel_[nei] -+ << " current refData:" << allCellInfo[nei] -+ << " level after refinement:" << neiLevel -+ << nl -+ << "which does not satisfy 2:1 constraints anymore." << nl -+ << "face:" << faceI << " faceRefData:" << allFaceInfo[faceI] -+ << abort(FatalError); -+ } -+ } -+ -+ -+ // Coupled faces. Swap owner level to get neighbouring cell level. -+ // (only boundary faces of neiLevel used) -+ -+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces()); -+ labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ neiLevel[i] = cellLevel_[own]; -+ neiCount[i] = allCellInfo[own].count(); -+ neiRefCount[i] = allCellInfo[own].refinementCount(); -+ } -+ -+ // Swap to neighbour -+ syncTools::swapBoundaryFaceList(mesh_, neiLevel, false); -+ syncTools::swapBoundaryFaceList(mesh_, neiCount, false); -+ syncTools::swapBoundaryFaceList(mesh_, neiRefCount, false); -+ -+ // Now we have neighbour value see which cells need refinement -+ forAll(neiLevel, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = -+ cellLevel_[own] -+ + (allCellInfo[own].isRefined() ? 1 : 0); -+ -+ label nbrLevel = -+ neiLevel[i] -+ + ((neiCount[i] >= neiRefCount[i]) ? 1 : 0); -+ -+ if (mag(ownLevel - nbrLevel) > 1) -+ { -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const labelList&)" -+ ) << "Celllevel does not satisfy 2:1 constraint." -+ << " On coupled face " -+ << faceI -+ << " refData:" << allFaceInfo[faceI] -+ << " on patch " << patchI << " " -+ << mesh_.boundaryMesh()[patchI].name() << nl -+ << "owner cell " << own -+ << " current level:" << cellLevel_[own] -+ << " current count:" << allCellInfo[own].count() -+ << " current refCount:" -+ << allCellInfo[own].refinementCount() -+ << " level after refinement:" << ownLevel -+ << nl -+ << "(coupled) neighbour cell" -+ << " has current level:" << neiLevel[i] -+ << " current count:" << neiCount[i] -+ << " current refCount:" << neiRefCount[i] -+ << " level after refinement:" << nbrLevel -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ // Convert back to labelList of cells to refine. -+ -+ label nRefined = 0; -+ -+ forAll(allCellInfo, cellI) -+ { -+ if (allCellInfo[cellI].isRefined()) -+ { -+ nRefined++; -+ } -+ } -+ -+ // Updated list of cells to refine -+ labelList newCellsToRefine(nRefined); -+ nRefined = 0; -+ -+ forAll(allCellInfo, cellI) -+ { -+ if (allCellInfo[cellI].isRefined()) -+ { -+ newCellsToRefine[nRefined++] = cellI; -+ } -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentSlowRefinement : From " -+ << cellsToRefine.size() << " to " << newCellsToRefine.size() -+ << " cells to refine." << endl; -+ } -+ -+ return newCellsToRefine; -+} -+ -+ -+Foam::labelList Foam::hexRef8::consistentSlowRefinement2 -+( -+ const label maxFaceDiff, -+ const labelList& cellsToRefine, -+ const labelList& facesToCheck -+) const -+{ -+ const labelList& faceOwner = mesh_.faceOwner(); -+ const labelList& faceNeighbour = mesh_.faceNeighbour(); -+ -+ if (maxFaceDiff <= 0) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement2" -+ "(const label, const labelList&, const labelList&)" -+ ) << "Illegal maxFaceDiff " << maxFaceDiff << nl -+ << "Value should be >= 1" << exit(FatalError); -+ } -+ -+ const scalar level0Size = 2*maxFaceDiff*level0Edge_; -+ -+ -+ // Bit tricky. Say we want a distance of three cells between two -+ // consecutive refinement levels. This is done by using FaceCellWave to -+ // transport out the 'refinement shell'. Anything inside the refinement -+ // shell (given by a distance) gets marked for refinement. -+ -+ // Initial information about (distance to) cellLevel on all cells -+ List<refinementDistanceData> allCellInfo(mesh_.nCells()); -+ -+ // Initial information about (distance to) cellLevel on all faces -+ List<refinementDistanceData> allFaceInfo(mesh_.nFaces()); -+ -+ -+ // Mark cells with wanted refinement level -+ forAll(cellsToRefine, i) -+ { -+ label cellI = cellsToRefine[i]; -+ -+ allCellInfo[cellI] = refinementDistanceData -+ ( -+ level0Size, -+ mesh_.cellCentres()[cellI], -+ cellLevel_[cellI]+1 // wanted refinement -+ ); -+ } -+ // Mark all others with existing refinement level -+ forAll(allCellInfo, cellI) -+ { -+ if (!allCellInfo[cellI].valid()) -+ { -+ allCellInfo[cellI] = refinementDistanceData -+ ( -+ level0Size, -+ mesh_.cellCentres()[cellI], -+ cellLevel_[cellI] // wanted refinement -+ ); -+ } -+ } -+ -+ -+ // Labels of seed faces -+ DynamicList<label> seedFaces(mesh_.nFaces()/100); -+ // refinementLevel data on seed faces -+ DynamicList<refinementDistanceData> seedFacesInfo(mesh_.nFaces()/100); -+ -+ -+ const pointField& cc = mesh_.cellCentres(); -+ -+ forAll(facesToCheck, i) -+ { -+ label faceI = facesToCheck[i]; -+ -+ if (allFaceInfo[faceI].valid()) -+ { -+ // Can only occur if face has already gone through loop below. -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement2" -+ "(const label, const labelList&, const labelList&)" -+ ) << "Argument facesToCheck seems to have duplicate entries!" -+ << endl -+ << "face:" << faceI << " occurs at positions " -+ << findIndices(facesToCheck, faceI) -+ << abort(FatalError); -+ } -+ -+ label own = faceOwner[faceI]; -+ -+ label ownLevel = -+ ( -+ allCellInfo[own].valid() -+ ? allCellInfo[own].originLevel() -+ : cellLevel_[own] -+ ); -+ -+ if (!mesh_.isInternalFace(faceI)) -+ { -+ // Do as if boundary face would have neighbour with one higher -+ // refinement level. -+ const point& fc = mesh_.faceCentres()[faceI]; -+ -+ refinementDistanceData neiData -+ ( -+ level0Size, -+ 2*fc - cc[own], // est'd cell centre -+ ownLevel+1 -+ ); -+ -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ own, // not used (should be nei) -+ neiData, -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ } -+ else -+ { -+ label nei = faceNeighbour[faceI]; -+ -+ label neiLevel = -+ ( -+ allCellInfo[nei].valid() -+ ? allCellInfo[nei].originLevel() -+ : cellLevel_[nei] -+ ); -+ -+ if (ownLevel == neiLevel) -+ { -+ // Fake as if nei>own or own>nei (whichever one 'wins') -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ nei, -+ refinementDistanceData(level0Size, cc[nei], neiLevel+1), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ own, -+ refinementDistanceData(level0Size, cc[own], ownLevel+1), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ } -+ else -+ { -+ // Difference in level anyway. -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ nei, -+ refinementDistanceData(level0Size, cc[nei], neiLevel), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ own, -+ refinementDistanceData(level0Size, cc[own], ownLevel), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ } -+ } -+ seedFaces.append(faceI); -+ seedFacesInfo.append(allFaceInfo[faceI]); -+ } -+ -+ -+ // Create some initial seeds to start walking from. This is only if there -+ // are no facesToCheck. -+ // Just seed with all faces inbetween different refinement levels for now -+ forAll(faceNeighbour, faceI) -+ { -+ // Check if face already handled in loop above -+ if (!allFaceInfo[faceI].valid()) -+ { -+ label own = faceOwner[faceI]; -+ -+ label ownLevel = -+ ( -+ allCellInfo[own].valid() -+ ? allCellInfo[own].originLevel() -+ : cellLevel_[own] -+ ); -+ -+ label nei = faceNeighbour[faceI]; -+ -+ label neiLevel = -+ ( -+ allCellInfo[nei].valid() -+ ? allCellInfo[nei].originLevel() -+ : cellLevel_[nei] -+ ); -+ -+ if (ownLevel > neiLevel) -+ { -+ // Set face to owner data. (since face not yet would be copy) -+ seedFaces.append(faceI); -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ own, -+ refinementDistanceData(level0Size, cc[own], ownLevel), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ seedFacesInfo.append(allFaceInfo[faceI]); -+ } -+ else if (neiLevel > ownLevel) -+ { -+ seedFaces.append(faceI); -+ allFaceInfo[faceI].updateFace -+ ( -+ mesh_, -+ faceI, -+ nei, -+ refinementDistanceData(level0Size, cc[nei], neiLevel), -+ FaceCellWave<refinementDistanceData>::propagationTol() -+ ); -+ seedFacesInfo.append(allFaceInfo[faceI]); -+ } -+ } -+ } -+ -+ seedFaces.shrink(); -+ seedFacesInfo.shrink(); -+ -+ // face-cell-face transport engine -+ FaceCellWave<refinementDistanceData> levelCalc -+ ( -+ mesh_, -+ seedFaces, -+ seedFacesInfo, -+ allFaceInfo, -+ allCellInfo, -+ mesh_.globalData().nTotalCells() -+ ); -+ -+ -+ //if (debug) -+ //{ -+ // // Dump wanted level -+ // volScalarField wantedLevel -+ // ( -+ // IOobject -+ // ( -+ // "wantedLevel", -+ // fMesh.time().timeName(), -+ // fMesh, -+ // IOobject::NO_READ, -+ // IOobject::AUTO_WRITE, -+ // false -+ // ), -+ // fMesh, -+ // dimensionedScalar("zero", dimless, 0) -+ // ); -+ // -+ // forAll(wantedLevel, cellI) -+ // { -+ // wantedLevel[cellI] = allCellInfo[cellI].wantedLevel(cc[cellI]); -+ // } -+ // -+ // Pout<< "Writing " << wantedLevel.objectPath() << endl; -+ // wantedLevel.write(); -+ //} -+ -+ -+ // Convert back to labelList of cells to refine. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ // 1. Force original refinement cells to be picked up by setting the -+ // originLevel of input cells to be a very large level (but within range -+ // of 1<< shift inside refinementDistanceData::wantedLevel) -+ forAll(cellsToRefine, i) -+ { -+ label cellI = cellsToRefine[i]; -+ -+ allCellInfo[cellI].originLevel() = sizeof(label)*8-2; -+ allCellInfo[cellI].origin() = cc[cellI]; -+ } -+ -+ // 2. Extend to 2:1. I don't understand yet why this is not done -+ PackedList<1> refineCell(mesh_.nCells(), 0); -+ forAll(allCellInfo, cellI) -+ { -+ label wanted = allCellInfo[cellI].wantedLevel(cc[cellI]); -+ -+ if (wanted > cellLevel_[cellI]+1) -+ { -+ refineCell.set(cellI, 1); -+ } -+ } -+ faceConsistentRefinement(true, refineCell); -+ -+ // 3. Convert back to labelList. -+ label nRefined = 0; -+ -+ forAll(refineCell, cellI) -+ { -+ if (refineCell.get(cellI) == 1) -+ { -+ nRefined++; -+ } -+ } -+ -+ labelList newCellsToRefine(nRefined); -+ nRefined = 0; -+ -+ forAll(refineCell, cellI) -+ { -+ if (refineCell.get(cellI) == 1) -+ { -+ newCellsToRefine[nRefined++] = cellI; -+ } -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentSlowRefinement2 : From " -+ << cellsToRefine.size() << " to " << newCellsToRefine.size() -+ << " cells to refine." << endl; -+ -+//XXXX -+ // Check that newCellsToRefine obeys at least 2:1. -+ -+ { -+ cellSet cellsIn(mesh_, "cellsToRefineIn", cellsToRefine); -+ Pout<< "hexRef8::consistentSlowRefinement2 : writing " -+ << cellsIn.size() << " to cellSet " -+ << cellsIn.objectPath() << endl; -+ cellsIn.write(); -+ } -+ { -+ cellSet cellsOut(mesh_, "cellsToRefineOut", newCellsToRefine); -+ Pout<< "hexRef8::consistentSlowRefinement2 : writing " -+ << cellsOut.size() << " to cellSet " -+ << cellsOut.objectPath() << endl; -+ cellsOut.write(); -+ } -+ -+ // Extend to 2:1 -+ PackedList<1> refineCell(mesh_.nCells(), 0); -+ forAll(newCellsToRefine, i) -+ { -+ refineCell.set(newCellsToRefine[i], 1); -+ } -+ const PackedList<1> savedRefineCell(refineCell); -+ -+ label nChanged = faceConsistentRefinement(true, refineCell); -+ -+ { -+ cellSet cellsOut2 -+ ( -+ mesh_, "cellsToRefineOut2", newCellsToRefine.size() -+ ); -+ forAll(refineCell, cellI) -+ { -+ if (refineCell.get(cellI) == 1) -+ { -+ cellsOut2.insert(cellI); -+ } -+ } -+ Pout<< "hexRef8::consistentSlowRefinement2 : writing " -+ << cellsOut2.size() << " to cellSet " -+ << cellsOut2.objectPath() << endl; -+ cellsOut2.write(); -+ } -+ -+ if (nChanged > 0) -+ { -+ forAll(refineCell, cellI) -+ { -+ if -+ ( -+ refineCell.get(cellI) == 1 -+ && savedRefineCell.get(cellI) == 0 -+ ) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::consistentSlowRefinement2" -+ "(const label, const labelList&, const labelList&)" -+ ) << "Cell:" << cellI << " cc:" -+ << mesh_.cellCentres()[cellI] -+ << " was not marked for refinement but does not obey" -+ << " 2:1 constraints." -+ << abort(FatalError); -+ } -+ } -+ } -+//XXXX -+ } -+ -+ return newCellsToRefine; -+} -+ -+ -+// Top level driver to insert topo changes to do all refinement. -+Foam::labelListList Foam::hexRef8::setRefinement -+( -+ const labelList& cellLabels, -+ polyTopoChange& meshMod -+) -+{ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Checking initial mesh just to make sure" << endl; -+ -+ checkMesh(); -+ checkRefinementLevels(-1, labelList(0)); -+ } -+ -+ // Clear any saved point/cell data. -+ savedPointLevel_.clear(); -+ savedCellLevel_.clear(); -+ -+ -+ // New point/cell level. Copy of pointLevel for existing points. -+ DynamicList<label> newCellLevel(cellLevel_.size()); -+ forAll(cellLevel_, cellI) -+ { -+ newCellLevel.append(cellLevel_[cellI]); -+ } -+ DynamicList<label> newPointLevel(pointLevel_.size()); -+ forAll(pointLevel_, pointI) -+ { -+ newPointLevel.append(pointLevel_[pointI]); -+ } -+ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Allocating " << cellLabels.size() << " cell midpoints." -+ << endl; -+ } -+ -+ -+ // Mid point per refined cell. -+ // -1 : not refined -+ // >=0: label of mid point. -+ labelList cellMidPoint(mesh_.nCells(), -1); -+ -+ forAll(cellLabels, i) -+ { -+ label cellI = cellLabels[i]; -+ -+ label anchorPointI = mesh_.faces()[mesh_.cells()[cellI][0]][0]; -+ -+ cellMidPoint[cellI] = meshMod.setAction -+ ( -+ polyAddPoint -+ ( -+ mesh_.cellCentres()[cellI], // point -+ anchorPointI, // master point -+ -1, // zone for point -+ true // supports a cell -+ ) -+ ); -+ -+ newPointLevel(cellMidPoint[cellI]) = cellLevel_[cellI]+1; -+ } -+ -+ -+ if (debug) -+ { -+ cellSet splitCells(mesh_, "splitCells", cellLabels.size()); -+ -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ splitCells.insert(cellI); -+ } -+ } -+ -+ Pout<< "hexRef8::setRefinement : Dumping " << splitCells.size() -+ << " cells to split to cellSet " << splitCells.objectPath() -+ << endl; -+ -+ splitCells.write(); -+ } -+ -+ -+ -+ // Split edges -+ // ~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Allocating edge midpoints." -+ << endl; -+ } -+ -+ // Unrefined edges are ones between cellLevel or lower points. -+ // If any cell using this edge gets split then the edge needs to be split. -+ -+ // -1 : no need to split edge -+ // >=0 : label of introduced mid point -+ labelList edgeMidPoint(mesh_.nEdges(), -1); -+ -+ // Note: Loop over cells to be refined or edges? -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ const labelList& cEdges = mesh_.cellEdges()[cellI]; -+ -+ forAll(cEdges, i) -+ { -+ label edgeI = cEdges[i]; -+ -+ const edge& e = mesh_.edges()[edgeI]; -+ -+ if -+ ( -+ pointLevel_[e[0]] <= cellLevel_[cellI] -+ && pointLevel_[e[1]] <= cellLevel_[cellI] -+ ) -+ { -+ edgeMidPoint[edgeI] = 12345; // mark need for splitting -+ } -+ } -+ } -+ } -+ -+ // Synchronize edgeMidPoint across coupled patches. Take max so that -+ // any split takes precedence. -+ syncTools::syncEdgeList -+ ( -+ mesh_, -+ edgeMidPoint, -+ maxEqOp<label>(), -+ labelMin, -+ false // no separation -+ ); -+ -+ -+ // Introduce edge points -+ // ~~~~~~~~~~~~~~~~~~~~~ -+ -+ { -+ // Phase 1: calculate midpoints and sync. -+ // This needs doing for if people do not write binary and we slowly -+ // get differences. -+ -+ pointField edgeMids(mesh_.nEdges(), point(-GREAT, -GREAT, -GREAT)); -+ -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ // Edge marked to be split. -+ edgeMids[edgeI] = mesh_.edges()[edgeI].centre(mesh_.points()); -+ } -+ } -+ syncTools::syncEdgeList -+ ( -+ mesh_, -+ edgeMids, -+ maxEqOp<vector>(), -+ point(-GREAT, -GREAT, -GREAT), -+ true // apply separation -+ ); -+ -+ -+ // Phase 2: introduce points at the synced locations. -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ // Edge marked to be split. Replace edgeMidPoint with actual -+ // point label. -+ -+ const edge& e = mesh_.edges()[edgeI]; -+ -+ edgeMidPoint[edgeI] = meshMod.setAction -+ ( -+ polyAddPoint -+ ( -+ edgeMids[edgeI], // point -+ e[0], // master point -+ -1, // zone for point -+ true // supports a cell -+ ) -+ ); -+ -+ newPointLevel(edgeMidPoint[edgeI]) = -+ max -+ ( -+ pointLevel_[e[0]], -+ pointLevel_[e[1]] -+ ) -+ + 1; -+ } -+ } -+ } -+ -+ if (debug) -+ { -+ OFstream str(mesh_.time().path()/"edgeMidPoint.obj"); -+ -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ const edge& e = mesh_.edges()[edgeI]; -+ -+ meshTools::writeOBJ(str, e.centre(mesh_.points())); -+ } -+ } -+ -+ Pout<< "hexRef8::setRefinement :" -+ << " Dumping edge centres to split to file " << str.name() << endl; -+ } -+ -+ -+ // Calculate face level -+ // ~~~~~~~~~~~~~~~~~~~~ -+ // (after splitting) -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Allocating face midpoints." -+ << endl; -+ } -+ -+ // Face anchor level. There are guaranteed 4 points with level -+ // <= anchorLevel. These are the corner points. -+ labelList faceAnchorLevel(mesh_.nFaces()); -+ -+ for (label faceI = 0; faceI < mesh_.nFaces(); faceI++) -+ { -+ faceAnchorLevel[faceI] = getAnchorLevel(faceI); -+ } -+ -+ // -1 : no need to split face -+ // >=0 : label of introduced mid point -+ labelList faceMidPoint(mesh_.nFaces(), -1); -+ -+ -+ // Internal faces: look at cells on both sides. Uniquely determined since -+ // face itself guaranteed to be same level as most refined neighbour. -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ if (faceAnchorLevel[faceI] >= 0) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own]; -+ label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0); -+ -+ label nei = mesh_.faceNeighbour()[faceI]; -+ label neiLevel = cellLevel_[nei]; -+ label newNeiLevel = neiLevel + (cellMidPoint[nei] >= 0 ? 1 : 0); -+ -+ if -+ ( -+ newOwnLevel > faceAnchorLevel[faceI] -+ || newNeiLevel > faceAnchorLevel[faceI] -+ ) -+ { -+ faceMidPoint[faceI] = 12345; // mark to be split -+ } -+ } -+ } -+ -+ // Coupled patches handled like internal faces except now all information -+ // from neighbour comes from across processor. -+ // Boundary faces are more complicated since the boundary face can -+ // be more refined than its owner (or neighbour for coupled patches) -+ // (does not happen if refining/unrefining only, but does e.g. when -+ // refinining and subsetting) -+ -+ { -+ labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(newNeiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ label ownLevel = cellLevel_[own]; -+ label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0); -+ -+ newNeiLevel[i] = newOwnLevel; -+ } -+ -+ // Swap. -+ syncTools::swapBoundaryFaceList(mesh_, newNeiLevel, false); -+ -+ // So now we have information on the neighbour. -+ -+ forAll(newNeiLevel, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ if (faceAnchorLevel[faceI] >= 0) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own]; -+ label newOwnLevel = ownLevel + (cellMidPoint[own] >= 0 ? 1 : 0); -+ -+ if -+ ( -+ newOwnLevel > faceAnchorLevel[faceI] -+ || newNeiLevel[i] > faceAnchorLevel[faceI] -+ ) -+ { -+ faceMidPoint[faceI] = 12345; // mark to be split -+ } -+ } -+ } -+ } -+ -+ -+ // Synchronize faceMidPoint across coupled patches. (logical or) -+ syncTools::syncFaceList -+ ( -+ mesh_, -+ faceMidPoint, -+ maxEqOp<label>(), -+ false -+ ); -+ -+ -+ -+ // Introduce face points -+ // ~~~~~~~~~~~~~~~~~~~~~ -+ -+ { -+ // Phase 1: determine mid points and sync. See comment for edgeMids -+ // above -+ pointField bFaceMids -+ ( -+ mesh_.nFaces()-mesh_.nInternalFaces(), -+ point(-GREAT, -GREAT, -GREAT) -+ ); -+ -+ forAll(bFaceMids, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ if (faceMidPoint[faceI] >= 0) -+ { -+ bFaceMids[i] = mesh_.faceCentres()[faceI]; -+ } -+ } -+ syncTools::syncBoundaryFaceList -+ ( -+ mesh_, -+ bFaceMids, -+ maxEqOp<vector>(), -+ true // apply separation -+ ); -+ -+ forAll(faceMidPoint, faceI) -+ { -+ if (faceMidPoint[faceI] >= 0) -+ { -+ // Face marked to be split. Replace faceMidPoint with actual -+ // point label. -+ -+ const face& f = mesh_.faces()[faceI]; -+ -+ faceMidPoint[faceI] = meshMod.setAction -+ ( -+ polyAddPoint -+ ( -+ ( -+ faceI < mesh_.nInternalFaces() -+ ? mesh_.faceCentres()[faceI] -+ : bFaceMids[faceI-mesh_.nInternalFaces()] -+ ), // point -+ f[0], // master point -+ -1, // zone for point -+ true // supports a cell -+ ) -+ ); -+ -+ // Determine the level of the corner points and midpoint will -+ // be one higher. -+ newPointLevel(faceMidPoint[faceI]) = faceAnchorLevel[faceI]+1; -+ } -+ } -+ } -+ -+ if (debug) -+ { -+ faceSet splitFaces(mesh_, "splitFaces", cellLabels.size()); -+ -+ forAll(faceMidPoint, faceI) -+ { -+ if (faceMidPoint[faceI] >= 0) -+ { -+ splitFaces.insert(faceI); -+ } -+ } -+ -+ Pout<< "hexRef8::setRefinement : Dumping " << splitFaces.size() -+ << " faces to split to faceSet " << splitFaces.objectPath() << endl; -+ -+ splitFaces.write(); -+ } -+ -+ -+ // Information complete -+ // ~~~~~~~~~~~~~~~~~~~~ -+ // At this point we have all the information we need. We should no -+ // longer reference the cellLabels to refine. All the information is: -+ // - cellMidPoint >= 0 : cell needs to be split -+ // - faceMidPoint >= 0 : face needs to be split -+ // - edgeMidPoint >= 0 : edge needs to be split -+ -+ -+ -+ // Get the corner/anchor points -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Finding cell anchorPoints (8 per cell)" -+ << endl; -+ } -+ -+ // There will always be 8 points on the hex that have were introduced -+ // with the hex and will have the same or lower refinement level. -+ -+ // Per cell the 8 corner points. -+ labelListList cellAnchorPoints(mesh_.nCells()); -+ -+ { -+ labelList nAnchorPoints(mesh_.nCells(), 0); -+ -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ cellAnchorPoints[cellI].setSize(8); -+ } -+ } -+ -+ forAll(pointLevel_, pointI) -+ { -+ const labelList& pCells = mesh_.pointCells()[pointI]; -+ -+ forAll(pCells, pCellI) -+ { -+ label cellI = pCells[pCellI]; -+ -+ if -+ ( -+ cellMidPoint[cellI] >= 0 -+ && pointLevel_[pointI] <= cellLevel_[cellI] -+ ) -+ { -+ if (nAnchorPoints[cellI] == 8) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::setRefinement(const labelList&" -+ ", polyTopoChange&)" -+ ) << "cell " << cellI -+ << " of level " << cellLevel_[cellI] -+ << " uses more than 8 points of equal or" -+ << " lower level" << nl -+ << "Points so far:" << cellAnchorPoints[cellI] -+ << abort(FatalError); -+ } -+ -+ cellAnchorPoints[cellI][nAnchorPoints[cellI]++] = pointI; -+ } -+ } -+ } -+ -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ if (nAnchorPoints[cellI] != 8) -+ { -+ const labelList cPoints(cellPoints(cellI)); -+ -+ FatalErrorIn -+ ( -+ "hexRef8::setRefinement(const labelList&" -+ ", polyTopoChange&)" -+ ) << "cell " << cellI -+ << " of level " << cellLevel_[cellI] -+ << " does not seem to have 8 points of equal or" -+ << " lower level" << endl -+ << "cellPoints:" << cPoints << endl -+ << "pointLevels:" -+ << IndirectList<label>(pointLevel_, cPoints)() << endl -+ << abort(FatalError); -+ } -+ } -+ } -+ } -+ -+ -+ // Add the cells -+ // ~~~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Adding cells (1 per anchorPoint)" -+ << endl; -+ } -+ -+ // Per cell the 7 added cells (+ original cell) -+ labelListList cellAddedCells(mesh_.nCells()); -+ -+ forAll(cellAnchorPoints, cellI) -+ { -+ const labelList& cAnchors = cellAnchorPoints[cellI]; -+ -+ if (cAnchors.size() == 8) -+ { -+ labelList& cAdded = cellAddedCells[cellI]; -+ cAdded.setSize(8); -+ -+ // Original cell at 0 -+ cAdded[0] = cellI; -+ -+ // Update cell level -+ newCellLevel[cellI] = cellLevel_[cellI]+1; -+ -+ -+ for (label i = 1; i < 8; i++) -+ { -+ cAdded[i] = meshMod.setAction -+ ( -+ polyAddCell -+ ( -+ -1, // master point -+ -1, // master edge -+ -1, // master face -+ cellI, // master cell -+ mesh_.cellZones().whichZone(cellI) // zone for cell -+ ) -+ ); -+ -+ newCellLevel(cAdded[i]) = cellLevel_[cellI]+1; -+ } -+ } -+ } -+ -+ -+ // Faces -+ // ~~~~~ -+ // 1. existing faces that get split (into four always) -+ // 2. existing faces that do not get split but only edges get split -+ // 3. existing faces that do not get split but get new owner/neighbour -+ // 4. new internal faces inside split cells. -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Marking faces to be handled" -+ << endl; -+ } -+ -+ // Get all affected faces. -+ PackedList<1> affectedFace(mesh_.nFaces(), 0); -+ -+ { -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ const cell& cFaces = mesh_.cells()[cellI]; -+ -+ forAll(cFaces, i) -+ { -+ affectedFace.set(cFaces[i], 1); -+ } -+ } -+ } -+ -+ forAll(faceMidPoint, faceI) -+ { -+ if (faceMidPoint[faceI] >= 0) -+ { -+ affectedFace.set(faceI, 1); -+ } -+ } -+ -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ const labelList& eFaces = mesh_.edgeFaces()[edgeI]; -+ -+ forAll(eFaces, i) -+ { -+ affectedFace.set(eFaces[i], 1); -+ } -+ } -+ } -+ } -+ -+ -+ // 1. Faces that get split -+ // ~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement : Splitting faces" << endl; -+ } -+ -+ forAll(faceMidPoint, faceI) -+ { -+ if (faceMidPoint[faceI] >= 0 && affectedFace.get(faceI) == 1) -+ { -+ // Face needs to be split and hasn't yet been done in some way -+ // (affectedFace - is impossible since this is first change but -+ // just for completeness) -+ -+ const face& f = mesh_.faces()[faceI]; -+ -+ // Has original faceI been used (three faces added, original gets -+ // modified) -+ bool modifiedFace = false; -+ -+ label anchorLevel = faceAnchorLevel[faceI]; -+ -+ face newFace(4); -+ -+ forAll(f, fp) -+ { -+ label pointI = f[fp]; -+ -+ if (pointLevel_[pointI] <= anchorLevel) -+ { -+ // point is anchor. Start collecting face. -+ -+ DynamicList<label> faceVerts(4); -+ -+ faceVerts.append(pointI); -+ -+ // Walk forward to mid point. -+ // - if next is +2 midpoint is +1 -+ // - if next is +1 it is midpoint -+ // - if next is +0 there has to be edgeMidPoint -+ -+ walkFaceToMid -+ ( -+ edgeMidPoint, -+ anchorLevel, -+ faceI, -+ fp, -+ faceVerts -+ ); -+ -+ faceVerts.append(faceMidPoint[faceI]); -+ -+ walkFaceFromMid -+ ( -+ edgeMidPoint, -+ anchorLevel, -+ faceI, -+ fp, -+ faceVerts -+ ); -+ -+ // Convert dynamiclist to face. -+ newFace.transfer(faceVerts.shrink()); -+ faceVerts.clear(); -+ -+ //Pout<< "Split face:" << faceI << " verts:" << f -+ // << " into quad:" << newFace << endl; -+ -+ // Get new owner/neighbour -+ label own, nei; -+ getFaceNeighbours -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ faceI, -+ pointI, // Anchor point -+ -+ own, -+ nei -+ ); -+ -+ -+ if (debug) -+ { -+ if (mesh_.isInternalFace(faceI)) -+ { -+ label oldOwn = mesh_.faceOwner()[faceI]; -+ label oldNei = mesh_.faceNeighbour()[faceI]; -+ -+ checkInternalOrientation -+ ( -+ meshMod, -+ oldOwn, -+ faceI, -+ mesh_.cellCentres()[oldOwn], -+ mesh_.cellCentres()[oldNei], -+ newFace -+ ); -+ } -+ else -+ { -+ label oldOwn = mesh_.faceOwner()[faceI]; -+ -+ checkBoundaryOrientation -+ ( -+ meshMod, -+ oldOwn, -+ faceI, -+ mesh_.cellCentres()[oldOwn], -+ mesh_.faceCentres()[faceI], -+ newFace -+ ); -+ } -+ } -+ -+ -+ if (!modifiedFace) -+ { -+ modifiedFace = true; -+ -+ modFace(meshMod, faceI, newFace, own, nei); -+ } -+ else -+ { -+ addFace(meshMod, faceI, newFace, own, nei); -+ } -+ } -+ } -+ -+ // Mark face as having been handled -+ affectedFace.set(faceI, 0); -+ } -+ } -+ -+ -+ // 2. faces that do not get split but use edges that get split -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Adding edge splits to unsplit faces" -+ << endl; -+ } -+ -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ // Split edge. Check that face not already handled above. -+ -+ const labelList& eFaces = mesh_.edgeFaces()[edgeI]; -+ -+ forAll(eFaces, i) -+ { -+ label faceI = eFaces[i]; -+ -+ if (faceMidPoint[faceI] < 0 && affectedFace.get(faceI) == 1) -+ { -+ // Unsplit face. Add edge splits to face. -+ -+ const face& f = mesh_.faces()[faceI]; -+ const labelList& fEdges = mesh_.faceEdges()[faceI]; -+ -+ DynamicList<label> newFaceVerts(f.size()); -+ -+ forAll(f, fp) -+ { -+ newFaceVerts.append(f[fp]); -+ -+ label edgeI = fEdges[fp]; -+ -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ newFaceVerts.append(edgeMidPoint[edgeI]); -+ } -+ } -+ -+ face newFace; -+ newFace.transfer(newFaceVerts.shrink()); -+ -+ -+ // The point with the lowest level should be an anchor -+ // point of the neighbouring cells. -+ label anchorFp = findMinLevel(f); -+ -+ label own, nei; -+ getFaceNeighbours -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ faceI, -+ f[anchorFp], // Anchor point -+ -+ own, -+ nei -+ ); -+ -+ -+ if (debug) -+ { -+ if (mesh_.isInternalFace(faceI)) -+ { -+ label oldOwn = mesh_.faceOwner()[faceI]; -+ label oldNei = mesh_.faceNeighbour()[faceI]; -+ -+ checkInternalOrientation -+ ( -+ meshMod, -+ oldOwn, -+ faceI, -+ mesh_.cellCentres()[oldOwn], -+ mesh_.cellCentres()[oldNei], -+ newFace -+ ); -+ } -+ else -+ { -+ label oldOwn = mesh_.faceOwner()[faceI]; -+ -+ checkBoundaryOrientation -+ ( -+ meshMod, -+ oldOwn, -+ faceI, -+ mesh_.cellCentres()[oldOwn], -+ mesh_.faceCentres()[faceI], -+ newFace -+ ); -+ } -+ } -+ -+ modFace(meshMod, faceI, newFace, own, nei); -+ -+ // Mark face as having been handled -+ affectedFace.set(faceI, 0); -+ } -+ } -+ } -+ } -+ -+ -+ // 3. faces that do not get split but whose owner/neighbour change -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Changing owner/neighbour for otherwise unaffected faces" -+ << endl; -+ } -+ -+ forAll(affectedFace, faceI) -+ { -+ if (affectedFace.get(faceI) == 1) -+ { -+ const face& f = mesh_.faces()[faceI]; -+ -+ // The point with the lowest level should be an anchor -+ // point of the neighbouring cells. -+ label anchorFp = findMinLevel(f); -+ -+ label own, nei; -+ getFaceNeighbours -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ faceI, -+ f[anchorFp], // Anchor point -+ -+ own, -+ nei -+ ); -+ -+ modFace(meshMod, faceI, f, own, nei); -+ -+ // Mark face as having been handled -+ affectedFace.set(faceI, 0); -+ } -+ } -+ -+ -+ // 4. new internal faces inside split cells. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ -+ // This is the hard one. We have to find the splitting points between -+ // the anchor points. But the edges between the anchor points might have -+ // been split (into two,three or four edges). -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Create new internal faces for split cells" -+ << endl; -+ } -+ -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ createInternalFaces -+ ( -+ cellAnchorPoints, -+ cellAddedCells, -+ cellMidPoint, -+ faceMidPoint, -+ faceAnchorLevel, -+ edgeMidPoint, -+ cellI, -+ meshMod -+ ); -+ } -+ } -+ -+ // Extend pointLevels and cellLevels for the new cells. Could also be done -+ // in updateMesh but saves passing cellAddedCells out of this routine. -+ -+ // Check -+ if (debug) -+ { -+ label minPointI = labelMax; -+ label maxPointI = labelMin; -+ -+ forAll(cellMidPoint, cellI) -+ { -+ if (cellMidPoint[cellI] >= 0) -+ { -+ minPointI = min(minPointI, cellMidPoint[cellI]); -+ maxPointI = max(maxPointI, cellMidPoint[cellI]); -+ } -+ } -+ forAll(faceMidPoint, faceI) -+ { -+ if (faceMidPoint[faceI] >= 0) -+ { -+ minPointI = min(minPointI, faceMidPoint[faceI]); -+ maxPointI = max(maxPointI, faceMidPoint[faceI]); -+ } -+ } -+ forAll(edgeMidPoint, edgeI) -+ { -+ if (edgeMidPoint[edgeI] >= 0) -+ { -+ minPointI = min(minPointI, edgeMidPoint[edgeI]); -+ maxPointI = max(maxPointI, edgeMidPoint[edgeI]); -+ } -+ } -+ -+ if (minPointI != labelMax && minPointI != mesh_.nPoints()) -+ { -+ FatalErrorIn("hexRef8::setRefinement") -+ << "Added point labels not consecutive to existing mesh points." -+ << nl -+ << "mesh_.nPoints():" << mesh_.nPoints() -+ << " minPointI:" << minPointI -+ << " maxPointI:" << maxPointI -+ << abort(FatalError); -+ } -+ } -+ -+ pointLevel_.transfer(newPointLevel.shrink()); -+ newPointLevel.clear(); -+ cellLevel_.transfer(newCellLevel.shrink()); -+ newCellLevel.clear(); -+ -+ // Mark files as changed -+ setInstance(mesh_.facesInstance()); -+ -+ -+ // Update the live split cells tree. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ // New unrefinement structure -+ if (history_.active()) -+ { -+ if (debug) -+ { -+ Pout<< "hexRef8::setRefinement :" -+ << " Updating refinement history to " << cellLevel_.size() -+ << " cells" << endl; -+ } -+ -+ // Extend refinement history for new cells -+ history_.resize(cellLevel_.size()); -+ -+ forAll(cellAddedCells, cellI) -+ { -+ const labelList& addedCells = cellAddedCells[cellI]; -+ -+ if (addedCells.size() > 0) -+ { -+ // Cell was split. -+ history_.storeSplit(cellI, addedCells); -+ } -+ } -+ } -+ -+ // Compact cellAddedCells. -+ -+ labelListList refinedCells(cellLabels.size()); -+ -+ forAll(cellLabels, i) -+ { -+ label cellI = cellLabels[i]; -+ -+ refinedCells[i].transfer(cellAddedCells[cellI]); -+ } -+ -+ return refinedCells; -+} -+ -+ -+void Foam::hexRef8::storeData -+( -+ const labelList& pointsToStore, -+ const labelList& facesToStore, -+ const labelList& cellsToStore -+) -+{ -+ savedPointLevel_.resize(2*pointsToStore.size()); -+ forAll(pointsToStore, i) -+ { -+ label pointI = pointsToStore[i]; -+ savedPointLevel_.insert(pointI, pointLevel_[pointI]); -+ } -+ -+ savedCellLevel_.resize(2*cellsToStore.size()); -+ forAll(cellsToStore, i) -+ { -+ label cellI = cellsToStore[i]; -+ savedCellLevel_.insert(cellI, cellLevel_[cellI]); -+ } -+} -+ -+ -+// Gets called after the mesh change. setRefinement will already have made -+// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we -+// only need to account for reordering. -+void Foam::hexRef8::updateMesh(const mapPolyMesh& map) -+{ -+ Map<label> dummyMap(0); -+ -+ updateMesh(map, dummyMap, dummyMap, dummyMap); -+} -+ -+ -+// Gets called after the mesh change. setRefinement will already have made -+// sure the pointLevel_ and cellLevel_ are the size of the new mesh so we -+// only need to account for reordering. -+void Foam::hexRef8::updateMesh -+( -+ const mapPolyMesh& map, -+ const Map<label>& pointsToRestore, -+ const Map<label>& facesToRestore, -+ const Map<label>& cellsToRestore -+) -+{ -+ // Update celllevel -+ if (debug) -+ { -+ Pout<< "hexRef8::updateMesh :" -+ << " Updating various lists" -+ << endl; -+ } -+ -+ { -+ const labelList& reverseCellMap = map.reverseCellMap(); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::updateMesh :" -+ << " reverseCellMap:" << map.reverseCellMap().size() -+ << " cellMap:" << map.cellMap().size() -+ << " nCells:" << mesh_.nCells() -+ << " nOldCells:" << map.nOldCells() -+ << " cellLevel_:" << cellLevel_.size() -+ << " reversePointMap:" << map.reversePointMap().size() -+ << " pointMap:" << map.pointMap().size() -+ << " nPoints:" << mesh_.nPoints() -+ << " nOldPoints:" << map.nOldPoints() -+ << " pointLevel_:" << pointLevel_.size() -+ << endl; -+ } -+ -+ if (reverseCellMap.size() == cellLevel_.size()) -+ { -+ // Assume it is after hexRef8 that this routine is called. -+ // Just account for reordering. We cannot use cellMap since -+ // then cells created from cells would get cellLevel_ of -+ // cell they were created from. -+ reorder(reverseCellMap, mesh_.nCells(), -1, cellLevel_); -+ } -+ else -+ { -+ // Map data -+ const labelList& cellMap = map.cellMap(); -+ -+ labelList newCellLevel(cellMap.size()); -+ forAll(cellMap, newCellI) -+ { -+ label oldCellI = cellMap[newCellI]; -+ -+ if (oldCellI == -1) -+ { -+ //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ // << "Problem : cell " << newCellI -+ // << " at " << mesh_.cellCentres()[newCellI] -+ // << " does not originate from another cell" -+ // << " (i.e. is inflated)." << nl -+ // << "Hence we cannot determine the new cellLevel" -+ // << " for it." << abort(FatalError); -+ newCellLevel[newCellI] = -1; -+ } -+ else -+ { -+ newCellLevel[newCellI] = cellLevel_[oldCellI]; -+ } -+ } -+ cellLevel_.transfer(newCellLevel); -+ } -+ -+ // See if any cells to restore. This will be for some new cells -+ // the corresponding old cell. -+ forAllConstIter(Map<label>, cellsToRestore, iter) -+ { -+ label newCellI = iter.key(); -+ label storedCellI = iter(); -+ -+ Map<label>::iterator fnd = savedCellLevel_.find(storedCellI); -+ -+ if (fnd == savedCellLevel_.end()) -+ { -+ FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ << "Problem : trying to restore old value for new cell " -+ << newCellI << " but cannot find old cell " << storedCellI -+ << " in map of stored values " << savedCellLevel_ -+ << abort(FatalError); -+ } -+ cellLevel_[newCellI] = fnd(); -+ } -+ -+ //if (findIndex(cellLevel_, -1) != -1) -+ //{ -+ // WarningIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ // << "Problem : " -+ // << "cellLevel_ contains illegal value -1 after mapping at cell " -+ // << findIndex(cellLevel_, -1) << endl -+ // << "This means that another program has inflated cells" -+ // << " (created cells out-of-nothing) and hence we don't know" -+ // << " their cell level. Continuing with illegal value." -+ // << abort(FatalError); -+ //} -+ } -+ -+ -+ // Update pointlevel -+ { -+ const labelList& reversePointMap = map.reversePointMap(); -+ -+ if (reversePointMap.size() == pointLevel_.size()) -+ { -+ // Assume it is after hexRef8 that this routine is called. -+ reorder(reversePointMap, mesh_.nPoints(), -1, pointLevel_); -+ } -+ else -+ { -+ // Map data -+ const labelList& pointMap = map.pointMap(); -+ -+ labelList newPointLevel(pointMap.size()); -+ -+ forAll(pointMap, newPointI) -+ { -+ label oldPointI = pointMap[newPointI]; -+ -+ if (oldPointI == -1) -+ { -+ //FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ // << "Problem : point " << newPointI -+ // << " at " << mesh_.points()[newPointI] -+ // << " does not originate from another point" -+ // << " (i.e. is inflated)." << nl -+ // << "Hence we cannot determine the new pointLevel" -+ // << " for it." << abort(FatalError); -+ newPointLevel[newPointI] = -1; -+ } -+ else -+ { -+ newPointLevel[newPointI] = pointLevel_[oldPointI]; -+ } -+ } -+ pointLevel_.transfer(newPointLevel); -+ } -+ -+ // See if any points to restore. This will be for some new points -+ // the corresponding old point (the one from the call to storeData) -+ forAllConstIter(Map<label>, pointsToRestore, iter) -+ { -+ label newPointI = iter.key(); -+ label storedPointI = iter(); -+ -+ Map<label>::iterator fnd = savedPointLevel_.find(storedPointI); -+ -+ if (fnd == savedPointLevel_.end()) -+ { -+ FatalErrorIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ << "Problem : trying to restore old value for new point " -+ << newPointI << " but cannot find old point " -+ << storedPointI -+ << " in map of stored values " << savedPointLevel_ -+ << abort(FatalError); -+ } -+ pointLevel_[newPointI] = fnd(); -+ } -+ -+ //if (findIndex(pointLevel_, -1) != -1) -+ //{ -+ // WarningIn("hexRef8::updateMesh(const mapPolyMesh&)") -+ // << "Problem : " -+ // << "pointLevel_ contains illegal value -1 after mapping" -+ // << " at point" << findIndex(pointLevel_, -1) << endl -+ // << "This means that another program has inflated points" -+ // << " (created points out-of-nothing) and hence we don't know" -+ // << " their point level. Continuing with illegal value." -+ // //<< abort(FatalError); -+ //} -+ } -+ -+ // Update refinement tree -+ if (history_.active()) -+ { -+ history_.updateMesh(map); -+ } -+ -+ // Mark files as changed -+ setInstance(mesh_.facesInstance()); -+ -+ // Update face removal engine -+ faceRemover_.updateMesh(map); -+} -+ -+ -+// Gets called after mesh subsetting. Maps are from new back to old. -+void Foam::hexRef8::subset -+( -+ const labelList& pointMap, -+ const labelList& faceMap, -+ const labelList& cellMap -+) -+{ -+ // Update celllevel -+ if (debug) -+ { -+ Pout<< "hexRef8::subset :" -+ << " Updating various lists" -+ << endl; -+ } -+ -+ if (history_.active()) -+ { -+ WarningIn -+ ( -+ "hexRef8::subset(const labelList&, const labelList&" -+ ", const labelList&)" -+ ) << "Subsetting will not work in combination with unrefinement." -+ << nl -+ << "Proceed at your own risk." << endl; -+ } -+ -+ -+ // Update celllevel -+ { -+ labelList newCellLevel(cellMap.size()); -+ -+ forAll(cellMap, newCellI) -+ { -+ newCellLevel[newCellI] = cellLevel_[cellMap[newCellI]]; -+ } -+ -+ cellLevel_.transfer(newCellLevel); -+ -+ if (findIndex(cellLevel_, -1) != -1) -+ { -+ FatalErrorIn("hexRef8::subset") -+ << "Problem : " -+ << "cellLevel_ contains illegal value -1 after mapping:" -+ << cellLevel_ -+ << abort(FatalError); -+ } -+ } -+ -+ // Update pointlevel -+ { -+ labelList newPointLevel(pointMap.size()); -+ -+ forAll(pointMap, newPointI) -+ { -+ newPointLevel[newPointI] = pointLevel_[pointMap[newPointI]]; -+ } -+ -+ pointLevel_.transfer(newPointLevel); -+ -+ if (findIndex(pointLevel_, -1) != -1) -+ { -+ FatalErrorIn("hexRef8::subset") -+ << "Problem : " -+ << "pointLevel_ contains illegal value -1 after mapping:" -+ << pointLevel_ -+ << abort(FatalError); -+ } -+ } -+ -+ // Update refinement tree -+ if (history_.active()) -+ { -+ history_.subset(pointMap, faceMap, cellMap); -+ } -+ -+ // Mark files as changed -+ setInstance(mesh_.facesInstance()); -+ -+ // Nothing needs doing to faceRemover. -+ //faceRemover_.subset(pointMap, faceMap, cellMap); -+} -+ -+ -+// Gets called after the mesh distribution -+void Foam::hexRef8::distribute(const mapDistributePolyMesh& map) -+{ -+ if (debug) -+ { -+ Pout<< "hexRef8::distribute :" -+ << " Updating various lists" -+ << endl; -+ } -+ -+ // Update celllevel -+ map.distributeCellData(cellLevel_); -+ // Update pointlevel -+ map.distributePointData(pointLevel_); -+ -+ // Update refinement tree -+ if (history_.active()) -+ { -+ history_.distribute(map); -+ } -+ -+ // Update face removal engine -+ faceRemover_.distribute(map); -+} -+ -+ -+void Foam::hexRef8::checkMesh() const -+{ -+ const boundBox& meshBb = mesh_.globalData().bb(); -+ const scalar smallDim = 1E-6*mag(meshBb.max() - meshBb.min()); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::checkMesh : Using matching tolerance smallDim:" -+ << smallDim << endl; -+ } -+ -+ // Check owner on coupled faces. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ // There should be only one coupled face between two cells. Why? Since -+ // otherwise mesh redistribution might cause multiple faces between two -+ // cells -+ { -+ labelList nei(mesh_.nFaces()-mesh_.nInternalFaces()); -+ forAll(nei, i) -+ { -+ nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ } -+ -+ // Replace data on coupled patches with their neighbour ones. -+ syncTools::swapBoundaryFaceList(mesh_, nei, false); -+ -+ const polyBoundaryMesh& patches = mesh_.boundaryMesh(); -+ -+ forAll(patches, patchI) -+ { -+ const polyPatch& pp = patches[patchI]; -+ -+ if (pp.coupled()) -+ { -+ // Check how many faces between owner and neighbour. Should -+ // be only one. -+ HashTable<label, labelPair, labelPair::Hash<> > -+ cellToFace(2*pp.size()); -+ -+ label faceI = pp.start(); -+ -+ forAll(pp, i) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label bFaceI = faceI-mesh_.nInternalFaces(); -+ -+ if (!cellToFace.insert(labelPair(own, nei[bFaceI]), faceI)) -+ { -+ FatalErrorIn("hexRef8::checkMesh()") -+ << "Faces do not seem to be correct across coupled" -+ << " boundaries" << endl -+ << "Coupled face " << faceI -+ << " between owner " << own -+ << " on patch " << pp.name() -+ << " and coupled neighbour " << nei[bFaceI] -+ << " has two faces connected to it:" -+ << faceI << " and " -+ << cellToFace[labelPair(own, nei[bFaceI])] -+ << abort(FatalError); -+ } -+ -+ faceI++; -+ } -+ } -+ } -+ } -+ -+ // Check face areas. -+ // ~~~~~~~~~~~~~~~~~ -+ -+ { -+ scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces()); -+ forAll(neiFaceAreas, i) -+ { -+ neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]); -+ } -+ -+ // Replace data on coupled patches with their neighbour ones. -+ syncTools::swapBoundaryFaceList(mesh_, neiFaceAreas, false); -+ -+ forAll(neiFaceAreas, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ const scalar magArea = mag(mesh_.faceAreas()[faceI]); -+ -+ if (mag(magArea - neiFaceAreas[i]) > smallDim) -+ { -+ const face& f = mesh_.faces()[faceI]; -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn("hexRef8::checkMesh()") -+ << "Faces do not seem to be correct across coupled" -+ << " boundaries" << endl -+ << "Coupled face " << faceI -+ << " on patch " << patchI -+ << " " << mesh_.boundaryMesh()[patchI].name() -+ << " coords:" << IndirectList<point>(mesh_.points(), f)() -+ << " has face area:" << magArea -+ << " (coupled) neighbour face area differs:" -+ << neiFaceAreas[i] -+ << " to within tolerance " << smallDim -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ -+ // Check number of points on faces. -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ { -+ labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(nVerts, i) -+ { -+ nVerts[i] = mesh_.faces()[i+mesh_.nInternalFaces()].size(); -+ } -+ -+ // Replace data on coupled patches with their neighbour ones. -+ syncTools::swapBoundaryFaceList(mesh_, nVerts, false); -+ -+ forAll(nVerts, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ const face& f = mesh_.faces()[faceI]; -+ -+ if (f.size() != nVerts[i]) -+ { -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn("hexRef8::checkMesh()") -+ << "Faces do not seem to be correct across coupled" -+ << " boundaries" << endl -+ << "Coupled face " << faceI -+ << " on patch " << patchI -+ << " " << mesh_.boundaryMesh()[patchI].name() -+ << " coords:" << IndirectList<point>(mesh_.points(), f)() -+ << " has size:" << f.size() -+ << " (coupled) neighbour face has size:" -+ << nVerts[i] -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ -+ // Check points of face -+ // ~~~~~~~~~~~~~~~~~~~~ -+ { -+ // Anchor points. -+ pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(anchorPoints, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ const point& fc = mesh_.faceCentres()[faceI]; -+ const face& f = mesh_.faces()[faceI]; -+ const vector anchorVec(mesh_.points()[f[0]] - fc); -+ -+ anchorPoints[i] = anchorVec; -+ } -+ -+ // Replace data on coupled patches with their neighbour ones. Apply -+ // rotation transformation (but not separation since is relative vector -+ // to point on same face. -+ syncTools::swapBoundaryFaceList(mesh_, anchorPoints, false); -+ -+ forAll(anchorPoints, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ const point& fc = mesh_.faceCentres()[faceI]; -+ const face& f = mesh_.faces()[faceI]; -+ const vector anchorVec(mesh_.points()[f[0]] - fc); -+ -+ if (mag(anchorVec - anchorPoints[i]) > smallDim) -+ { -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn("hexRef8::checkMesh()") -+ << "Faces do not seem to be correct across coupled" -+ << " boundaries" << endl -+ << "Coupled face " << faceI -+ << " on patch " << patchI -+ << " " << mesh_.boundaryMesh()[patchI].name() -+ << " coords:" << IndirectList<point>(mesh_.points(), f)() -+ << " has anchor vector:" << anchorVec -+ << " (coupled) neighbour face anchor vector differs:" -+ << anchorPoints[i] -+ << " to within tolerance " << smallDim -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::checkMesh : Returning" << endl; -+ } -+} -+ -+ -+void Foam::hexRef8::checkRefinementLevels -+( -+ const label maxPointDiff, -+ const labelList& pointsToCheck -+) const -+{ -+ if (debug) -+ { -+ Pout<< "hexRef8::checkRefinementLevels :" -+ << " Checking 2:1 refinement level" << endl; -+ } -+ -+ if -+ ( -+ cellLevel_.size() != mesh_.nCells() -+ || pointLevel_.size() != mesh_.nPoints() -+ ) -+ { -+ FatalErrorIn("hexRef8::checkRefinementLevels(const label)") -+ << "cellLevel size should be number of cells" -+ << " and pointLevel size should be number of points."<< nl -+ << "cellLevel:" << cellLevel_.size() -+ << " mesh.nCells():" << mesh_.nCells() << nl -+ << "pointLevel:" << pointLevel_.size() -+ << " mesh.nPoints():" << mesh_.nPoints() -+ << abort(FatalError); -+ } -+ -+ -+ // Check 2:1 consistency. -+ // ~~~~~~~~~~~~~~~~~~~~~~ -+ -+ { -+ // Internal faces. -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label nei = mesh_.faceNeighbour()[faceI]; -+ -+ if (mag(cellLevel_[own] - cellLevel_[nei]) > 1) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const label)" -+ ) << "Celllevel does not satisfy 2:1 constraint." << nl -+ << "On face " << faceI << " owner cell " << own -+ << " has refinement " << cellLevel_[own] -+ << " neighbour cell " << nei << " has refinement " -+ << cellLevel_[nei] -+ << abort(FatalError); -+ } -+ } -+ -+ // Coupled faces. Get neighbouring value -+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ -+ neiLevel[i] = cellLevel_[own]; -+ } -+ -+ // No separation -+ syncTools::swapBoundaryFaceList(mesh_, neiLevel, false); -+ -+ forAll(neiLevel, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ -+ label own = mesh_.faceOwner()[faceI]; -+ -+ if (mag(cellLevel_[own] - neiLevel[i]) > 1) -+ { -+ label patchI = mesh_.boundaryMesh().whichPatch(faceI); -+ -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const label)" -+ ) << "Celllevel does not satisfy 2:1 constraint." -+ << " On coupled face " << faceI -+ << " on patch " << patchI << " " -+ << mesh_.boundaryMesh()[patchI].name() -+ << " owner cell " << own << " has refinement " -+ << cellLevel_[own] -+ << " (coupled) neighbour cell has refinement " -+ << neiLevel[i] -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ -+ // Check pointLevel is synchronized -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ { -+ labelList syncPointLevel(pointLevel_); -+ -+ // Get min level -+ syncTools::syncPointList -+ ( -+ mesh_, -+ syncPointLevel, -+ minEqOp<label>(), -+ labelMax, -+ false // no separation -+ ); -+ -+ -+ forAll(syncPointLevel, pointI) -+ { -+ if (pointLevel_[pointI] != syncPointLevel[pointI]) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const label)" -+ ) << "PointLevel is not consistent across coupled patches." -+ << endl -+ << "point:" << pointI << " coord:" << mesh_.points()[pointI] -+ << " has level " << pointLevel_[pointI] -+ << " whereas the coupled point has level " -+ << syncPointLevel[pointI] -+ << abort(FatalError); -+ } -+ } -+ } -+ -+ -+ // Check 2:1 across points (instead of faces) -+ if (maxPointDiff != -1) -+ { -+ const labelListList& pointCells = mesh_.pointCells(); -+ -+ // Determine per point the max cell level. -+ labelList maxPointLevel(mesh_.nPoints(), 0); -+ -+ forAll(pointCells, pointI) -+ { -+ const labelList& pCells = pointCells[pointI]; -+ -+ label& pLevel = maxPointLevel[pointI]; -+ -+ forAll(pCells, i) -+ { -+ pLevel = max(pLevel, cellLevel_[pCells[i]]); -+ } -+ } -+ -+ // Sync maxPointLevel to neighbour -+ syncTools::syncPointList -+ ( -+ mesh_, -+ maxPointLevel, -+ maxEqOp<label>(), -+ labelMin, // null value -+ false // no separation -+ ); -+ -+ // Check 2:1 across boundary points -+ forAll(pointsToCheck, i) -+ { -+ label pointI = pointsToCheck[i]; -+ -+ const labelList& pCells = pointCells[pointI]; -+ -+ forAll(pCells, i) -+ { -+ label cellI = pCells[i]; -+ -+ if -+ ( -+ mag(cellLevel_[cellI]-maxPointLevel[pointI]) -+ > maxPointDiff -+ ) -+ { -+ FatalErrorIn -+ ( -+ "hexRef8::checkRefinementLevels(const label)" -+ ) << "Too big a difference between" -+ << " point-connected cells." << nl -+ << "cell:" << cellI -+ << " cellLevel:" << cellLevel_[cellI] -+ << " uses point:" << pointI -+ << " coord:" << mesh_.points()[pointI] -+ << " which is also used by a cell with level:" -+ << maxPointLevel[pointI] -+ << abort(FatalError); -+ } -+ } -+ } -+ } -+} -+ -+ -+ -+// -+// Unrefinement -+// ~~~~~~~~~~~~ -+// -+ -+ -+Foam::labelList Foam::hexRef8::getSplitPoints() const -+{ -+ if (debug) -+ { -+ checkRefinementLevels(-1, labelList(0)); -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::getSplitPoints :" -+ << " Calculating unrefineable points" << endl; -+ } -+ -+ -+ if (!history_.active()) -+ { -+ FatalErrorIn("hexRef8::getSplitPoints()") -+ << "Only call if constructed with history capability" -+ << abort(FatalError); -+ } -+ -+ // Master cell -+ // -1 undetermined -+ // -2 certainly not split point -+ // >= label of master cell -+ labelList splitMaster(mesh_.nPoints(), -1); -+ labelList splitMasterLevel(mesh_.nPoints(), 0); -+ -+ // Unmark all with not 8 cells -+ const labelListList& pointCells = mesh_.pointCells(); -+ -+ forAll(pointCells, pointI) -+ { -+ const labelList& pCells = pointCells[pointI]; -+ -+ if (pCells.size() != 8) -+ { -+ splitMaster[pointI] = -2; -+ } -+ } -+ -+ // Unmark all with different master cells -+ const labelList& visibleCells = history_.visibleCells(); -+ -+ forAll(visibleCells, cellI) -+ { -+ //const labelList& cPoints = mesh_.cellPoints()[cellI]; -+ const labelList cPoints(cellPoints(cellI)); -+ -+ if (visibleCells[cellI] != -1 && history_.parentIndex(cellI) >= 0) -+ { -+ label parentIndex = history_.parentIndex(cellI); -+ -+ // Check same master. -+ forAll(cPoints, i) -+ { -+ label pointI = cPoints[i]; -+ -+ label masterCellI = splitMaster[pointI]; -+ -+ if (masterCellI == -1) -+ { -+ // First time visit of point. Store parent cell and -+ // level of the parent cell (with respect to cellI). This -+ // is additional guarantee that we're referring to the -+ // same master at the same refinement level. -+ -+ splitMaster[pointI] = parentIndex; -+ splitMasterLevel[pointI] = cellLevel_[cellI] - 1; -+ } -+ else if (masterCellI == -2) -+ { -+ // Already decided that point is not splitPoint -+ } -+ else if -+ ( -+ (masterCellI != parentIndex) -+ || (splitMasterLevel[pointI] != cellLevel_[cellI] - 1) -+ ) -+ { -+ // Different masters so point is on two refinement -+ // patterns -+ splitMaster[pointI] = -2; -+ } -+ } -+ } -+ else -+ { -+ // Either not visible or is unrefined cell -+ forAll(cPoints, i) -+ { -+ label pointI = cPoints[i]; -+ -+ splitMaster[pointI] = -2; -+ } -+ } -+ } -+ -+ // Unmark boundary faces -+ for -+ ( -+ label faceI = mesh_.nInternalFaces(); -+ faceI < mesh_.nFaces(); -+ faceI++ -+ ) -+ { -+ const face& f = mesh_.faces()[faceI]; -+ -+ forAll(f, fp) -+ { -+ splitMaster[f[fp]] = -2; -+ } -+ } -+ -+ -+ // Collect into labelList -+ -+ label nSplitPoints = 0; -+ -+ forAll(splitMaster, pointI) -+ { -+ if (splitMaster[pointI] >= 0) -+ { -+ nSplitPoints++; -+ } -+ } -+ -+ labelList splitPoints(nSplitPoints); -+ nSplitPoints = 0; -+ -+ forAll(splitMaster, pointI) -+ { -+ if (splitMaster[pointI] >= 0) -+ { -+ splitPoints[nSplitPoints++] = pointI; -+ } -+ } -+ -+ return splitPoints; -+} -+ -+ -+//void Foam::hexRef8::markIndex -+//( -+// const label maxLevel, -+// const label level, -+// const label index, -+// const label markValue, -+// labelList& indexValues -+//) const -+//{ -+// if (level < maxLevel && indexValues[index] == -1) -+// { -+// // Mark -+// indexValues[index] = markValue; -+// -+// // Mark parent -+// const splitCell8& split = history_.splitCells()[index]; -+// -+// if (split.parent_ >= 0) -+// { -+// markIndex -+// ( -+// maxLevel, level+1, split.parent_, markValue, indexValues); -+// ) -+// } -+// } -+//} -+// -+// -+//// Get all cells which (down to level) originate from the same cell. -+//// level=0 returns cell only, level=1 returns the 8 cells this cell -+//// originates from, level=2 returns 64 cells etc. -+//// If the cell does not originate from refinement returns just itself. -+//void Foam::hexRef8::markCellClusters -+//( -+// const label maxLevel, -+// labelList& cluster -+//) const -+//{ -+// cluster.setSize(mesh_.nCells()); -+// cluster = -1; -+// -+// const DynamicList<splitCell8>& splitCells = history_.splitCells(); -+// -+// // Mark all splitCells down to level maxLevel with a cell originating from -+// // it. -+// -+// labelList indexLevel(splitCells.size(), -1); -+// -+// forAll(visibleCells, cellI) -+// { -+// label index = visibleCells[cellI]; -+// -+// if (index >= 0) -+// { -+// markIndex(maxLevel, 0, index, cellI, indexLevel); -+// } -+// } -+// -+// // Mark cells with splitCell -+//} -+ -+ -+Foam::labelList Foam::hexRef8::consistentUnrefinement -+( -+ const labelList& pointsToUnrefine, -+ const bool maxSet -+) const -+{ -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentUnrefinement :" -+ << " Determining 2:1 consistent unrefinement" << endl; -+ } -+ -+ if (maxSet) -+ { -+ FatalErrorIn("hexRef8::consistentUnrefinement") -+ << "maxSet not implemented yet." -+ << abort(FatalError); -+ } -+ -+ // Loop, modifying pointsToUnrefine, until no more changes to due to 2:1 -+ // conflicts. -+ // maxSet = false : unselect points to refine -+ // maxSet = true: select points to refine -+ -+ // Maintain boolList for pointsToUnrefine and cellsToUnrefine -+ PackedList<1> unrefinePoint(mesh_.nPoints(), 0); -+ -+ forAll(pointsToUnrefine, i) -+ { -+ label pointI = pointsToUnrefine[i]; -+ -+ unrefinePoint.set(pointI, 1); -+ } -+ -+ -+ while (true) -+ { -+ // Construct cells to unrefine -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ PackedList<1> unrefineCell(mesh_.nCells(), 0); -+ -+ forAll(unrefinePoint, pointI) -+ { -+ if (unrefinePoint.get(pointI) == 1) -+ { -+ const labelList& pCells = mesh_.pointCells()[pointI]; -+ -+ forAll(pCells, j) -+ { -+ unrefineCell.set(pCells[j], 1); -+ } -+ } -+ } -+ -+ -+ label nChanged = 0; -+ -+ -+ // Check 2:1 consistency taking refinement into account -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ // Internal faces. -+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) -+ { -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own] - unrefineCell.get(own); -+ -+ label nei = mesh_.faceNeighbour()[faceI]; -+ label neiLevel = cellLevel_[nei] - unrefineCell.get(nei); -+ -+ if (ownLevel < (neiLevel-1)) -+ { -+ // Since was 2:1 this can only occur if own is marked for -+ // unrefinement. -+ -+ if (maxSet) -+ { -+ unrefineCell.set(nei, 1); -+ } -+ else -+ { -+ if (unrefineCell.get(own) == 0) -+ { -+ FatalErrorIn("hexRef8::consistentUnrefinement") -+ << "problem" << abort(FatalError); -+ } -+ -+ unrefineCell.set(own, 0); -+ } -+ nChanged++; -+ } -+ else if (neiLevel < (ownLevel-1)) -+ { -+ if (maxSet) -+ { -+ unrefineCell.set(own, 1); -+ } -+ else -+ { -+ if (unrefineCell.get(nei) == 0) -+ { -+ FatalErrorIn("hexRef8::consistentUnrefinement") -+ << "problem" << abort(FatalError); -+ } -+ -+ unrefineCell.set(nei, 0); -+ } -+ nChanged++; -+ } -+ } -+ -+ -+ // Coupled faces. Swap owner level to get neighbouring cell level. -+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces()); -+ -+ forAll(neiLevel, i) -+ { -+ label own = mesh_.faceOwner()[i+mesh_.nInternalFaces()]; -+ -+ neiLevel[i] = cellLevel_[own] - unrefineCell.get(own); -+ } -+ -+ // Swap to neighbour -+ syncTools::swapBoundaryFaceList(mesh_, neiLevel, false); -+ -+ forAll(neiLevel, i) -+ { -+ label faceI = i+mesh_.nInternalFaces(); -+ label own = mesh_.faceOwner()[faceI]; -+ label ownLevel = cellLevel_[own] - unrefineCell.get(own); -+ -+ if (ownLevel < (neiLevel[i]-1)) -+ { -+ if (!maxSet) -+ { -+ if (unrefineCell.get(own) == 0) -+ { -+ FatalErrorIn("hexRef8::consistentUnrefinement") -+ << "problem" << abort(FatalError); -+ } -+ -+ unrefineCell.set(own, 0); -+ nChanged++; -+ } -+ } -+ else if (neiLevel[i] < (ownLevel-1)) -+ { -+ if (maxSet) -+ { -+ if (unrefineCell.get(own) == 1) -+ { -+ FatalErrorIn("hexRef8::consistentUnrefinement") -+ << "problem" << abort(FatalError); -+ } -+ -+ unrefineCell.set(own, 1); -+ nChanged++; -+ } -+ } -+ } -+ -+ reduce(nChanged, sumOp<label>()); -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::consistentUnrefinement :" -+ << " Changed " << nChanged -+ << " refinement levels due to 2:1 conflicts." -+ << endl; -+ } -+ -+ if (nChanged == 0) -+ { -+ break; -+ } -+ -+ -+ // Convert cellsToUnrefine back into points to unrefine -+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -+ -+ // Knock out any point whose cell neighbour cannot be unrefined. -+ forAll(unrefinePoint, pointI) -+ { -+ if (unrefinePoint.get(pointI) == 1) -+ { -+ const labelList& pCells = mesh_.pointCells()[pointI]; -+ -+ forAll(pCells, j) -+ { -+ if (unrefineCell.get(pCells[j]) == 0) -+ { -+ unrefinePoint.set(pointI, 0); -+ break; -+ } -+ } -+ } -+ } -+ } -+ -+ -+ // Convert back to labelList. -+ label nSet = 0; -+ -+ forAll(unrefinePoint, pointI) -+ { -+ if (unrefinePoint.get(pointI) == 1) -+ { -+ nSet++; -+ } -+ } -+ -+ labelList newPointsToUnrefine(nSet); -+ nSet = 0; -+ -+ forAll(unrefinePoint, pointI) -+ { -+ if (unrefinePoint.get(pointI) == 1) -+ { -+ newPointsToUnrefine[nSet++] = pointI; -+ } -+ } -+ -+ return newPointsToUnrefine; -+} -+ -+ -+void Foam::hexRef8::setUnrefinement -+( -+ const labelList& splitPointLabels, -+ polyTopoChange& meshMod -+) -+{ -+ if (!history_.active()) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement()") -+ << "Only call if constructed with history capability" -+ << abort(FatalError); -+ } -+ -+ if (debug) -+ { -+ Pout<< "hexRef8::setUnrefinement :" -+ << " Checking initial mesh just to make sure" << endl; -+ -+ checkMesh(); -+ -+ forAll(cellLevel_, cellI) -+ { -+ if (cellLevel_[cellI] < 0) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement()") -+ << "Illegal cell level " << cellLevel_[cellI] -+ << " for cell " << cellI -+ << abort(FatalError); -+ } -+ } -+ -+ -+ // Write to sets. -+ pointSet pSet(mesh_, "splitPoints", splitPointLabels); -+ pSet.write(); -+ -+ cellSet cSet(mesh_, "splitPointCells", splitPointLabels.size()); -+ -+ forAll(splitPointLabels, i) -+ { -+ const labelList& pCells = mesh_.pointCells()[splitPointLabels[i]]; -+ -+ forAll(pCells, j) -+ { -+ cSet.insert(pCells[j]); -+ } -+ } -+ cSet.write(); -+ -+ Pout<< "hexRef8::setRefinement : Dumping " << pSet.size() -+ << " points and " -+ << cSet.size() << " cells for unrefinement to" << nl -+ << " pointSet " << pSet.objectPath() << nl -+ << " cellSet " << cSet.objectPath() -+ << endl; -+ } -+ -+ -+ labelList cellRegion; -+ labelList cellRegionMaster; -+ labelList facesToRemove; -+ -+ { -+ labelHashSet splitFaces(12*splitPointLabels.size()); -+ -+ forAll(splitPointLabels, i) -+ { -+ const labelList& pFaces = mesh_.pointFaces()[splitPointLabels[i]]; -+ -+ forAll(pFaces, j) -+ { -+ splitFaces.insert(pFaces[j]); -+ } -+ } -+ -+ // Check with faceRemover what faces will get removed. Note that this -+ // can be more (but never less) than splitFaces provided. -+ faceRemover_.compatibleRemoves -+ ( -+ splitFaces.toc(), // pierced faces -+ cellRegion, // per cell -1 or region it is merged into -+ cellRegionMaster, // per region the master cell -+ facesToRemove // new faces to be removed. -+ ); -+ -+ if (facesToRemove.size() != splitFaces.size()) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement") -+ << "Ininitial set of split points to unrefine does not" -+ << " seem to be consistent or not mid points of refined cells" -+ << abort(FatalError); -+ } -+ } -+ -+ // Redo the region master so it is consistent with our master. -+ // This will guarantee that the new cell (for which faceRemover uses -+ // the region master) is already compatible with our refinement structure. -+ -+ forAll(splitPointLabels, i) -+ { -+ label pointI = splitPointLabels[i]; -+ -+ // Get original cell label -+ -+ const labelList& pCells = mesh_.pointCells()[pointI]; -+ -+ // Check -+ if (pCells.size() != 8) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement") -+ << "splitPoint " << pointI -+ << " should have 8 cells using it. It has " << pCells -+ << abort(FatalError); -+ } -+ -+ -+ // Check that the lowest numbered pCells is the master of the region -+ // (should be guaranteed by directRemoveFaces) -+ //if (debug) -+ { -+ label masterCellI = min(pCells); -+ -+ forAll(pCells, j) -+ { -+ label cellI = pCells[j]; -+ -+ label region = cellRegion[cellI]; -+ -+ if (region == -1) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement") -+ << "Ininitial set of split points to unrefine does not" -+ << " seem to be consistent or not mid points" -+ << " of refined cells" << nl -+ << "cell:" << cellI << " on splitPoint " << pointI -+ << " has no region to be merged into" -+ << abort(FatalError); -+ } -+ -+ if (masterCellI != cellRegionMaster[region]) -+ { -+ FatalErrorIn("hexRef8::setUnrefinement") -+ << "cell:" << cellI << " on splitPoint:" << pointI -+ << " in region " << region -+ << " has master:" << cellRegionMaster[region] -+ << " which is not the lowest numbered cell" -+ << " among the pointCells:" << pCells -+ << abort(FatalError); -+ } -+ } -+ } -+ } -+ -+ // Insert all commands to combine cells. Never fails so don't have to -+ // test for success. -+ faceRemover_.setRefinement -+ ( -+ facesToRemove, -+ cellRegion, -+ cellRegionMaster, -+ meshMod -+ ); -+ -+ // Remove the 8 cells that originated from merging around the split point -+ // and adapt cell levels (not that pointLevels stay the same since points -+ // either get removed or stay at the same position. -+ forAll(splitPointLabels, i) -+ { -+ label pointI = splitPointLabels[i]; -+ -+ const labelList& pCells = mesh_.pointCells()[pointI]; -+ -+ label masterCellI = min(pCells); -+ -+ forAll(pCells, j) -+ { -+ cellLevel_[pCells[j]]--; -+ } -+ -+ history_.combineCells(masterCellI, pCells); -+ } -+ -+ // Mark files as changed -+ setInstance(mesh_.facesInstance()); -+ -+ // history_.updateMesh will take care of truncating. -+} -+ -+ -+// Write refinement to polyMesh directory. -+bool Foam::hexRef8::write() const -+{ -+ bool writeOk = cellLevel_.write() && pointLevel_.write(); -+ -+ if (history_.active()) -+ { -+ writeOk = writeOk && history_.write(); -+ } -+ -+ return writeOk; -+} -+ -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/surfaceToCell.C -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/surfaceToCell.C (Revision 0) -+++ applications/utilities/mesh/advanced/autoRefineMesh/surfaceToCell.C (Revision 784) -@@ -0,0 +1,486 @@ -+/*---------------------------------------------------------------------------*\ -+ ========= | -+ \\ / 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 -+ -+Description -+ -+\*---------------------------------------------------------------------------*/ -+ -+#include "surfaceToCell.H" -+#include "polyMesh.H" -+#include "meshSearch.H" -+#include "triSurface.H" -+#include "triSurfaceSearch.H" -+#include "octree.H" -+#include "octreeDataTriSurface.H" -+#include "cellClassification.H" -+ -+#include "addToRunTimeSelectionTable.H" -+ -+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -+ -+namespace Foam -+{ -+ -+defineTypeNameAndDebug(surfaceToCell, 0); -+ -+addToRunTimeSelectionTable(topoSetSource, surfaceToCell, word); -+ -+addToRunTimeSelectionTable(topoSetSource, surfaceToCell, istream); -+ -+} -+ -+ -+Foam::topoSetSource::addToUsageTable Foam::surfaceToCell::usage_ -+( -+ surfaceToCell::typeName, -+ "\n Usage: surfaceToCell" -+ "<surface> <outsidePoints> <cut> <inside> <outside> <near> <curvature>\n\n" -+ " <surface> name of triSurface\n" -+ " <outsidePoints> list of points that define outside\n" -+ " <cut> boolean whether to include cells cut by surface\n" -+ " <inside> ,, ,, inside surface\n" -+ " <outside> ,, ,, outside surface\n" -+ " <near> scalar; include cells with centre <= near to surface\n" -+ " <curvature> scalar; include cells close to strong curvature" -+ " on surface\n" -+ " (curvature defined as difference in surface normal at nearest" -+ " point on surface for each vertex of cell)\n\n" -+); -+ -+ -+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -+ -+Foam::label Foam::surfaceToCell::getNearest -+( -+ const triSurfaceSearch& querySurf, -+ const label pointI, -+ const point& pt, -+ const vector& span, -+ Map<label>& cache -+) -+{ -+ Map<label>::const_iterator iter = cache.find(pointI); -+ -+ if (iter != cache.end()) -+ { -+ // Found cached answer -+ return iter(); -+ } -+ else -+ { -+ pointIndexHit inter = querySurf.nearest(pt, span); -+ -+ // Triangle label (can be -1) -+ label triI = inter.index(); -+ -+ // Store triangle on point -+ cache.insert(pointI, triI); -+ -+ return triI; -+ } -+} -+ -+ -+// Return true if nearest surface to points on cell makes largish angle -+// with nearest surface to cell centre. Returns false otherwise. Points visited -+// are cached in pointToNearest -+bool Foam::surfaceToCell::differingPointNormals -+( -+ const triSurfaceSearch& querySurf, -+ -+ const vector& span, // current search span -+ const label cellI, -+ const label cellTriI, // nearest (to cell centre) surface triangle -+ -+ Map<label>& pointToNearest // cache for nearest triangle to point -+) const -+{ -+ const triSurface& surf = querySurf.surface(); -+ const vectorField& normals = surf.faceNormals(); -+ -+ const faceList& faces = mesh().faces(); -+ const pointField& points = mesh().points(); -+ -+ const labelList& cFaces = mesh().cells()[cellI]; -+ -+ forAll(cFaces, cFaceI) -+ { -+ const face& f = faces[cFaces[cFaceI]]; -+ -+ forAll(f, fp) -+ { -+ label pointI = f[fp]; -+ -+ label pointTriI = -+ getNearest -+ ( -+ querySurf, -+ pointI, -+ points[pointI], -+ span, -+ pointToNearest -+ ); -+ -+ if (pointTriI != -1 && pointTriI != cellTriI) -+ { -+ scalar cosAngle = normals[pointTriI] & normals[cellTriI]; -+ -+ if (cosAngle < 0.9) -+ { -+ return true; -+ } -+ } -+ } -+ } -+ return false; -+} -+ -+ -+void Foam::surfaceToCell::combine(topoSet& set, const bool add) const -+{ -+ cpuTime timer; -+ -+ if (includeCut_ || includeInside_ || includeOutside_) -+ { -+ // -+ // Cut cells with surface and classify cells -+ // -+ -+ -+ // Construct search engine on mesh -+ -+ meshSearch queryMesh(mesh_, true); -+ -+ -+ // Check all 'outside' points -+ forAll(outsidePoints_, outsideI) -+ { -+ const point& outsidePoint = outsidePoints_[outsideI]; -+ -+ // Find cell point is in. Linear search. -+ if (queryMesh.findCell(outsidePoint, -1, false) == -1) -+ { -+ FatalErrorIn("surfaceToCell::combine(topoSet&, const bool)") -+ << "outsidePoint " << outsidePoint -+ << " is not inside any cell" -+ << exit(FatalError); -+ } -+ } -+ -+ // Cut faces with surface and classify cells -+ -+ cellClassification cellType -+ ( -+ mesh_, -+ queryMesh, -+ querySurf(), -+ outsidePoints_ -+ ); -+ -+ -+ Info<< " Marked inside/outside in = " -+ << timer.cpuTimeIncrement() << " s" << endl << endl; -+ -+ -+ forAll(cellType, cellI) -+ { -+ label cType = cellType[cellI]; -+ -+ if -+ ( -+ ( -+ includeCut_ -+ && (cType == cellClassification::CUT) -+ ) -+ || ( -+ includeInside_ -+ && (cType == cellClassification::INSIDE) -+ ) -+ || ( -+ includeOutside_ -+ && (cType == cellClassification::OUTSIDE) -+ ) -+ ) -+ { -+ addOrDelete(set, cellI, add); -+ } -+ } -+ } -+ -+ -+ if (nearDist_ > 0) -+ { -+ // -+ // Determine distance to surface -+ // -+ -+ const pointField& ctrs = mesh_.cellCentres(); -+ -+ // Box dimensions to search in octree. -+ const vector span(nearDist_, nearDist_, nearDist_); -+ -+ -+ if (curvature_ < -1) -+ { -+ Info<< " Selecting cells with cellCentre closer than " -+ << nearDist_ << " to surface" << endl; -+ -+ // No need to test curvature. Insert near cells into set. -+ -+ forAll(ctrs, cellI) -+ { -+ const point& c = ctrs[cellI]; -+ -+ pointIndexHit inter = querySurf().nearest(c, span); -+ -+ if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_)) -+ { -+ addOrDelete(set, cellI, add); -+ } -+ } -+ -+ Info<< " Determined nearest surface point in = " -+ << timer.cpuTimeIncrement() << " s" << endl << endl; -+ -+ } -+ else -+ { -+ // Test near cells for curvature -+ -+ Info<< " Selecting cells with cellCentre closer than " -+ << nearDist_ << " to surface and curvature factor" -+ << " less than " << curvature_ << endl; -+ -+ // Cache for nearest surface triangle for a point -+ Map<label> pointToNearest(mesh_.nCells()/10); -+ -+ forAll(ctrs, cellI) -+ { -+ const point& c = ctrs[cellI]; -+ -+ pointIndexHit inter = querySurf().nearest(c, span); -+ -+ if (inter.hit() && (mag(inter.hitPoint() - c) < nearDist_)) -+ { -+ if -+ ( -+ differingPointNormals -+ ( -+ querySurf(), -+ span, -+ cellI, -+ inter.index(), // nearest surface triangle -+ pointToNearest -+ ) -+ ) -+ { -+ addOrDelete(set, cellI, add); -+ } -+ } -+ } -+ -+ Info<< " Determined nearest surface point in = " -+ << timer.cpuTimeIncrement() << " s" << endl << endl; -+ } -+ } -+} -+ -+ -+void Foam::surfaceToCell::checkSettings() const -+{ -+ if -+ ( -+ (nearDist_ < 0) -+ && (curvature_ < -1) -+ && ( -+ (includeCut_ && includeInside_ && includeOutside_) -+ || (!includeCut_ && !includeInside_ && !includeOutside_) -+ ) -+ ) -+ { -+ FatalErrorIn -+ ( -+ "surfaceToCell:checkSettings()" -+ ) << "Illegal include cell specification." -+ << " Result would be either all or no cells." << endl -+ << "Please set one of includeCut, includeInside, includeOutside" -+ << " to true, set nearDistance to a value > 0" -+ << " or set curvature to a value -1 .. 1." -+ << exit(FatalError); -+ } -+} -+ -+ -+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -+ -+// Construct from components -+Foam::surfaceToCell::surfaceToCell -+( -+ const polyMesh& mesh, -+ const fileName& surfName, -+ const pointField& outsidePoints, -+ const bool includeCut, -+ const bool includeInside, -+ const bool includeOutside, -+ const scalar nearDist, -+ const scalar curvature -+) -+: -+ topoSetSource(mesh), -+ surfName_(surfName), -+ outsidePoints_(outsidePoints), -+ includeCut_(includeCut), -+ includeInside_(includeInside), -+ includeOutside_(includeOutside), -+ nearDist_(nearDist), -+ curvature_(curvature), -+ surfPtr_(new triSurface(surfName_)), -+ querySurfPtr_(new triSurfaceSearch(*surfPtr_)), -+ IOwnPtrs_(true) -+{ -+ checkSettings(); -+} -+ -+ -+// Construct from components. Externally supplied surface. -+Foam::surfaceToCell::surfaceToCell -+( -+ const polyMesh& mesh, -+ const fileName& surfName, -+ const triSurface& surf, -+ const triSurfaceSearch& querySurf, -+ const pointField& outsidePoints, -+ const bool includeCut, -+ const bool includeInside, -+ const bool includeOutside, -+ const scalar nearDist, -+ const scalar curvature -+) -+: -+ topoSetSource(mesh), -+ surfName_(surfName), -+ outsidePoints_(outsidePoints), -+ includeCut_(includeCut), -+ includeInside_(includeInside), -+ includeOutside_(includeOutside), -+ nearDist_(nearDist), -+ curvature_(curvature), -+ surfPtr_(&surf), -+ querySurfPtr_(&querySurf), -+ IOwnPtrs_(false) -+{ -+ checkSettings(); -+} -+ -+ -+// Construct from dictionary -+Foam::surfaceToCell::surfaceToCell -+( -+ const polyMesh& mesh, -+ const dictionary& dict -+) -+: -+ topoSetSource(mesh), -+ surfName_(dict.lookup("file")), -+ outsidePoints_(dict.lookup("outsidePoints")), -+ includeCut_(readBool(dict.lookup("includeCut"))), -+ includeInside_(readBool(dict.lookup("includeInside"))), -+ includeOutside_(readBool(dict.lookup("includeOutside"))), -+ nearDist_(readScalar(dict.lookup("nearDistance"))), -+ curvature_(readScalar(dict.lookup("curvature"))), -+ surfPtr_(new triSurface(surfName_)), -+ querySurfPtr_(new triSurfaceSearch(*surfPtr_)), -+ IOwnPtrs_(true) -+{ -+ checkSettings(); -+} -+ -+ -+// Construct from Istream -+Foam::surfaceToCell::surfaceToCell -+( -+ const polyMesh& mesh, -+ Istream& is -+) -+: -+ topoSetSource(mesh), -+ surfName_(checkIs(is)), -+ outsidePoints_(checkIs(is)), -+ includeCut_(readBool(checkIs(is))), -+ includeInside_(readBool(checkIs(is))), -+ includeOutside_(readBool(checkIs(is))), -+ nearDist_(readScalar(checkIs(is))), -+ curvature_(readScalar(checkIs(is))), -+ surfPtr_(new triSurface(surfName_)), -+ querySurfPtr_(new triSurfaceSearch(*surfPtr_)), -+ IOwnPtrs_(true) -+{ -+ checkSettings(); -+} -+ -+ -+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // -+ -+Foam::surfaceToCell::~surfaceToCell() -+{ -+ if (IOwnPtrs_) -+ { -+ if (surfPtr_) -+ { -+ delete surfPtr_; -+ } -+ if (querySurfPtr_) -+ { -+ delete querySurfPtr_; -+ } -+ } -+} -+ -+ -+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -+ -+void Foam::surfaceToCell::applyToSet -+( -+ const topoSetSource::setAction action, -+ topoSet& set -+) const -+{ -+ if ( (action == topoSetSource::NEW) || (action == topoSetSource::ADD)) -+ { -+ Info<< " Adding cells in relation to surface " << surfName_ -+ << " ..." << endl; -+ -+ combine(set, true); -+ } -+ else if (action == topoSetSource::DELETE) -+ { -+ Info<< " Removing cells in relation to surface " << surfName_ -+ << " ..." << endl; -+ -+ combine(set, false); -+ } -+} -+ -+ -+// ************************************************************************* // -Index: mesh/advanced/autoRefineMesh/Make/files -=================================================================== ---- applications/utilities/mesh/advanced/autoRefineMesh/Make/files (Revision 30) -+++ applications/utilities/mesh/advanced/autoRefineMesh/Make/files (Revision 784) -@@ -1,4 +1,8 @@ -+refinementDistanceData.C -+/* FaceCellWaveName.C */ -+hexRef8.C -+surfaceToCell.C - autoRefineMesh.C - - EXE = $(FOAM_APPBIN)/autoRefineMesh - - -Index: mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C -=================================================================== ---- applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C (Revision 30) -+++ applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C (Revision 784) -@@ -276,26 +276,26 @@ - is.getLine(line); - } - else if (feID == 171) - { - // Rod. Skip. - is.getLine(line); - } -- else if (feID == 41) -+ else if (feID == 41 || feID == 91) - { - // Triangle. Save - used for patching later on. - is.getLine(line); - - face cVerts(3); - IStringStream lineStr(line); - lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2]; - boundaryFaces.append(cVerts); - boundaryFaceIndices.append(cellI); - } -- else if (feID == 94) -+ else if (feID == 44 || feID == 94) - { - // Quad. Save - used for patching later on. - is.getLine(line); - - face cVerts(4); - IStringStream lineStr(line); - lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]; -@@ -389,29 +389,38 @@ - faceIndices.setSize(nFaces); - label faceI = 0; - - while (faceI < faceIndices.size()) - { - is.getLine(line); - IStringStream lineStr(line); -- label typeCode1, tag1, nodeLeaf1, component1, -- typeCode2, tag2, nodeLeaf2, component2; -- lineStr -- >> typeCode1 >> tag1 >> nodeLeaf1 >> component1 -- >> typeCode2 >> tag2 >> nodeLeaf2 >> component2; - -- if (typeCode1 != 8 || typeCode2 != 8) -+ // Read one (for last face) or two entries from line. -+ label nRead = 2; -+ if (faceI == faceIndices.size()-1) - { -- FatalErrorIn("readPatches") -- << "When reading patches expect Entity Type Code 8" << nl -- << "At line " << is.lineNumber() << exit(FatalError); -+ nRead = 1; - } - -- faceIndices[faceI++] = tag1; -- faceIndices[faceI++] = tag2; -+ for (label i = 0; i < nRead; i++) -+ { -+ label typeCode, tag, nodeLeaf, component; -+ -+ lineStr >> typeCode >> tag >> nodeLeaf >> component; -+ -+ if (typeCode != 8) -+ { -+ FatalErrorIn("readPatches") -+ << "When reading patches expect Entity Type Code 8" -+ << nl << "At line " << is.lineNumber() -+ << exit(FatalError); -+ } -+ -+ faceIndices[faceI++] = tag; -+ } - } - } - - patchNames.shrink(); - patchFaceIndices.shrink(); - } - -@@ -688,15 +697,15 @@ - // sets (dofVertIndices). - - List<faceList> patchFaceVerts; - - - if (dofVertIndices.size() > 0) - { -- // Use the vertex constraints to patch. Is of course bit dodge since -+ // Use the vertex constraints to patch. Is of course bit dodgy since - // face goes on patch if all its vertices are on a constraint. - // Note: very inefficient since goes through all faces (including - // internal ones) twice. Should do a construct without patches - // and then repatchify. - - Info<< "Using " << dofVertIndices.size() - << " DOF sets to detect boundary faces."<< endl; - -Index: mesh/generation/blockMesh/curvedEdges/polySplineEdge.C -=================================================================== ---- applications/utilities/mesh/generation/blockMesh/curvedEdges/polySplineEdge.C (Revision 30) -+++ applications/utilities/mesh/generation/blockMesh/curvedEdges/polySplineEdge.C (Revision 784) -@@ -68,15 +68,15 @@ - - label nSize(nsize(otherknots.size(), nbetweenKnots)); - - pointField ans(nSize); - - label N = spl.nKnots(); - scalar init = 1.0/(N - 1); -- scalar interval = (N - 3)/N; -+ scalar interval = (N - scalar(3))/N; - interval /= otherknots.size() + 1; - interval /= nbetweenKnots + 1; - - ans[0] = points_[start_]; - - register scalar index(init); - for (register label i=1; i<nSize-1; i++) -@@ -121,15 +121,15 @@ - // Construct from Istream - polySplineEdge::polySplineEdge - ( - const pointField& points, - Istream& is - ) - : -- curvedEdge(points, readLabel(is), readLabel(is)), -+ curvedEdge(points, is), - polyLine(pointField(0)), - otherKnots_(is) - { - label nInterKnots(20); - vector fstend(is); - vector sndend(is); - - diff --git a/sci-libs/openfoam-utilities/openfoam-utilities-1.4.1_p20080827.ebuild b/sci-libs/openfoam-utilities/openfoam-utilities-1.4.1_p20080827.ebuild index 9e57b2f96..7ae2c0834 100644 --- a/sci-libs/openfoam-utilities/openfoam-utilities-1.4.1_p20080827.ebuild +++ b/sci-libs/openfoam-utilities/openfoam-utilities-1.4.1_p20080827.ebuild @@ -10,7 +10,8 @@ MY_P="${MY_PN}-${MY_PV}" DESCRIPTION="OpenFOAM - Utilities" HOMEPAGE="http://www.opencfd.co.uk/openfoam/" -SRC_URI="mirror://sourceforge/foam/${MY_P}.General.gtgz" +SRC_URI="mirror://sourceforge/foam/${MY_P}.General.gtgz + http://dev.gentooexperimental.org/~tommy/distfiles/${P}.patch" LICENSE="GPL-2" SLOT="0" @@ -40,7 +41,7 @@ src_unpack() { unpack ./${MY_P}.General.tgz cd "${S}" - epatch "${FILESDIR}"/${P}.patch + epatch "${DISTDIR}"/${P}.patch epatch "${FILESDIR}"/${PN}-compile-${PV}.patch } |