From aa217c28678d84cedad3544054f42bc9c257b4e5 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 5 Jan 2010 22:13:45 +0000 Subject: [PATCH 01/24] Added combine function to combine two pointConstraints --- .../pointConstraint/pointConstraint.H | 3 ++ .../pointConstraint/pointConstraintI.H | 41 +++++++++++++++++++ 2 files changed, 44 insertions(+) diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H index 5e28ac1e14..10aa63f2c9 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H @@ -77,6 +77,9 @@ public: //- Apply and accumulate the effect of the given constraint direction inline void applyConstraint(const vector& cd); + //- Combine constraints + inline void combine(const pointConstraint&); + //- Return the accumulated constraint transformation tensor inline tensor constraintTransformation() const; }; diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H index a0af05b5af..c488bf603e 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H @@ -69,6 +69,47 @@ void Foam::pointConstraint::applyConstraint(const vector& cd) } +void Foam::pointConstraint::combine(const pointConstraint& pc) +{ + if (first() == 0) + { + operator=(pc); + } + else if (first() == 1) + { + // Save single normal + vector n = second(); + // Apply to supplied point constaint + operator=(pc); + applyConstraint(n); + } + else if (first() == 2) + { + if (pc.first() == 0) + {} + else if (pc.first() == 1) + { + applyConstraint(pc.second()); + } + else if (pc.first() == 2) + { + // Both constrained to line. Same (+-)direction? + if (mag(second() & pc.second()) <= (1.0-1e-3)) + { + // Different directions + first() = 3; + second() = vector::zero; + } + } + else + { + first() = 3; + second() = vector::zero; + } + } +} + + Foam::tensor Foam::pointConstraint::constraintTransformation() const { if (first() == 0) From 210f2305b6085906d6a4ce80b152ac2412bee92c Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 13 Jan 2010 19:06:11 +0000 Subject: [PATCH 02/24] Moved dumping of graph up to ease debugging --- .../scotchDecomp/scotchDecomp.C | 90 +++++++++---------- 1 file changed, 45 insertions(+), 45 deletions(-) diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C index 5375fb529b..c14e7e39f1 100644 --- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C +++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C @@ -126,6 +126,51 @@ Foam::label Foam::scotchDecomp::decompose List& finalDecomp ) { + // Dump graph + if (decompositionDict_.found("scotchCoeffs")) + { + const dictionary& scotchCoeffs = + decompositionDict_.subDict("scotchCoeffs"); + + if (scotchCoeffs.found("writeGraph")) + { + Switch writeGraph(scotchCoeffs.lookup("writeGraph")); + + if (writeGraph) + { + OFstream str(mesh_.time().path() / mesh_.name() + ".grf"); + + Info<< "Dumping Scotch graph file to " << str.name() << endl + << "Use this in combination with gpart." << endl; + + label version = 0; + str << version << nl; + // Numer of vertices + str << xadj.size()-1 << ' ' << adjncy.size() << nl; + // Numbering starts from 0 + label baseval = 0; + // Has weights? + label hasEdgeWeights = 0; + label hasVertexWeights = 0; + label numericflag = 10*hasEdgeWeights+hasVertexWeights; + str << baseval << ' ' << numericflag << nl; + for (label cellI = 0; cellI < xadj.size()-1; cellI++) + { + label start = xadj[cellI]; + label end = xadj[cellI+1]; + str << end-start; + + for (label i = start; i < end; i++) + { + str << ' ' << adjncy[i]; + } + str << nl; + } + } + } + } + + // Strategy // ~~~~~~~~ @@ -206,51 +251,6 @@ Foam::label Foam::scotchDecomp::decompose check(SCOTCH_graphCheck(&grafdat), "SCOTCH_graphCheck"); - // Dump graph - if (decompositionDict_.found("scotchCoeffs")) - { - const dictionary& scotchCoeffs = - decompositionDict_.subDict("scotchCoeffs"); - - if (scotchCoeffs.found("writeGraph")) - { - Switch writeGraph(scotchCoeffs.lookup("writeGraph")); - - if (writeGraph) - { - OFstream str(mesh_.time().path() / mesh_.name() + ".grf"); - - Info<< "Dumping Scotch graph file to " << str.name() << endl - << "Use this in combination with gpart." << endl; - - label version = 0; - str << version << nl; - // Numer of vertices - str << xadj.size()-1 << ' ' << adjncy.size() << nl; - // Numbering starts from 0 - label baseval = 0; - // Has weights? - label hasEdgeWeights = 0; - label hasVertexWeights = 0; - label numericflag = 10*hasEdgeWeights+hasVertexWeights; - str << baseval << ' ' << numericflag << nl; - for (label cellI = 0; cellI < xadj.size()-1; cellI++) - { - label start = xadj[cellI]; - label end = xadj[cellI+1]; - str << end-start; - - for (label i = start; i < end; i++) - { - str << ' ' << adjncy[i]; - } - str << nl; - } - } - } - } - - // Architecture // ~~~~~~~~~~~~ // (fully connected network topology since using switch) From 735de0c758c6a4be4831a418a2d223dae769a144 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 13 Jan 2010 19:07:04 +0000 Subject: [PATCH 03/24] Assign points even when no merging; remove points used by illegal triangles --- src/triSurface/triSurface/stitchTriangles.C | 93 ++++++++++++++++----- 1 file changed, 72 insertions(+), 21 deletions(-) diff --git a/src/triSurface/triSurface/stitchTriangles.C b/src/triSurface/triSurface/stitchTriangles.C index 7a441ac679..70fe304808 100644 --- a/src/triSurface/triSurface/stitchTriangles.C +++ b/src/triSurface/triSurface/stitchTriangles.C @@ -26,6 +26,7 @@ License #include "triSurface.H" #include "mergePoints.H" +#include "PackedBoolList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -46,46 +47,44 @@ bool triSurface::stitchTriangles pointField newPoints; bool hasMerged = mergePoints(rawPoints, tol, verbose, pointMap, newPoints); + pointField& ps = storedPoints(); + // Set the coordinates to the merged ones + ps.transfer(newPoints); if (hasMerged) { if (verbose) { Pout<< "stitchTriangles : Merged from " << rawPoints.size() - << " points down to " << newPoints.size() << endl; + << " points down to " << ps.size() << endl; } - pointField& ps = storedPoints(); - - // Set the coordinates to the merged ones - ps = newPoints; - // Reset the triangle point labels to the unique points array label newTriangleI = 0; forAll(*this, i) { - label newA = pointMap[operator[](i)[0]]; - label newB = pointMap[operator[](i)[1]]; - label newC = pointMap[operator[](i)[2]]; + const labelledTri& tri = operator[](i); + labelledTri newTri + ( + pointMap[tri[0]], + pointMap[tri[1]], + pointMap[tri[2]], + tri.region() + ); - if ((newA != newB) && (newA != newC) && (newB != newC)) + if ((newTri[0] != newTri[1]) && (newTri[0] != newTri[2]) && (newTri[1] != newTri[2])) { - operator[](newTriangleI)[0] = newA; - operator[](newTriangleI)[1] = newB; - operator[](newTriangleI)[2] = newC; - operator[](newTriangleI).region() = operator[](i).region(); - newTriangleI++; + operator[](newTriangleI++) = newTri; } else if (verbose) { Pout<< "stitchTriangles : " - << "Removing triangle " << i << " with non-unique vertices." - << endl - << " vertices :" << newA << ' ' << newB << ' ' << newC - << endl - << " coordinates:" << ps[newA] << ' ' << ps[newB] - << ' ' << ps[newC] << endl; + << "Removing triangle " << i + << " with non-unique vertices." << endl + << " vertices :" << newTri << endl + << " coordinates:" << newTri.points(ps) + << endl; } } @@ -98,6 +97,58 @@ bool triSurface::stitchTriangles << " triangles" << endl; } setSize(newTriangleI); + + // And possibly compact out any unused points (since used only + // by triangles that have just been deleted) + // Done in two passes to save memory (pointField) + + // 1. Detect only + PackedBoolList pointIsUsed(ps.size()); + + label nPoints = 0; + + forAll(*this, i) + { + const labelledTri& tri = operator[](i); + + forAll(tri, fp) + { + label pointI = tri[fp]; + if (pointIsUsed.set(pointI, 1)) + { + nPoints++; + } + } + } + + if (nPoints != ps.size()) + { + // 2. Compact. + pointMap.setSize(ps.size()); + label newPointI = 0; + forAll(pointIsUsed, pointI) + { + if (pointIsUsed[pointI]) + { + ps[newPointI] = ps[pointI]; + pointMap[pointI] = newPointI++; + } + } + ps.setSize(newPointI); + + newTriangleI = 0; + forAll(*this, i) + { + const labelledTri& tri = operator[](i); + operator[](newTriangleI++) = labelledTri + ( + pointMap[tri[0]], + pointMap[tri[1]], + pointMap[tri[2]], + tri.region() + ); + } + } } } From 1e78a1486552d9ab5e6c83de2fa4cc2878985038 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 13 Jan 2010 19:08:01 +0000 Subject: [PATCH 04/24] Handle tracking on face or edge --- src/meshTools/indexedOctree/indexedOctree.C | 54 +++++++++++---------- 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/src/meshTools/indexedOctree/indexedOctree.C b/src/meshTools/indexedOctree/indexedOctree.C index 128e6bf600..1d444eff3d 100644 --- a/src/meshTools/indexedOctree/indexedOctree.C +++ b/src/meshTools/indexedOctree/indexedOctree.C @@ -957,6 +957,7 @@ Foam::point Foam::indexedOctree::pushPointIntoFace } else if (nFaces == 1) { + // Point is on a single face keepFaceID = faceIndices[0]; } else @@ -1782,16 +1783,6 @@ Foam::pointIndexHit Foam::indexedOctree::findLine label i = 0; for (; i < 100000; i++) { - if (verbose) - { - Pout<< "iter:" << i - << " at startPoint:" << hitInfo.rawPoint() << endl - << " node:" << nodeI - << " octant:" << octant - << " bb:" << subBbox(nodeI, octant) << endl; - } - - // Ray-trace to end of current node. Updates point (either on triangle // in case of hit or on node bounding box in case of miss) @@ -1808,6 +1799,19 @@ Foam::pointIndexHit Foam::indexedOctree::findLine ) ); + if (verbose) + { + Pout<< "iter:" << i + << " at current:" << hitInfo.rawPoint() + << " (perturbed:" << startPoint << ")" << endl + << " node:" << nodeI + << " octant:" << octant + << " bb:" << subBbox(nodeI, octant) << endl; + } + + + + // Faces of current bounding box current point is on direction hitFaceID = 0; @@ -1833,12 +1837,23 @@ Foam::pointIndexHit Foam::indexedOctree::findLine break; } - if (hitFaceID == 0) + if (hitFaceID == 0 || hitInfo.rawPoint() == treeEnd) { // endpoint inside the tree. Return miss. break; } + // Create a point on other side of face. + point perturbedPoint + ( + pushPoint + ( + octantBb, + hitFaceID, + hitInfo.rawPoint(), + false // push outside of octantBb + ) + ); if (verbose) { @@ -1848,14 +1863,7 @@ Foam::pointIndexHit Foam::indexedOctree::findLine << " node:" << nodeI << " octant:" << octant << " bb:" << subBbox(nodeI, octant) << nl - << " walking to neighbour containing:" - << pushPoint - ( - octantBb, - hitFaceID, - hitInfo.rawPoint(), - false // push outside of octantBb - ) + << " walking to neighbour containing:" << perturbedPoint << endl; } @@ -1866,13 +1874,7 @@ Foam::pointIndexHit Foam::indexedOctree::findLine bool ok = walkToNeighbour ( - pushPoint - ( - octantBb, - hitFaceID, - hitInfo.rawPoint(), - false // push outside of octantBb - ), + perturbedPoint, hitFaceID, // face(s) that hitInfo is on nodeI, From 0496e181236503bd21f8121bf0dbced60177a514 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 13 Jan 2010 19:08:47 +0000 Subject: [PATCH 05/24] Modified tolerances to take truncation error into account --- .../polyPatches/basic/coupled/coupledPolyPatch.C | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C index ee2d966836..8e18f68c77 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C @@ -199,13 +199,20 @@ Foam::scalarField Foam::coupledPolyPatch::calcFaceTol const face& f = faces[faceI]; + // 1. calculate a typical size of the face. Use maximum distance + // to face centre scalar maxLenSqr = -GREAT; + // 2. as measure of truncation error when comparing two coordinates + // use SMALL * maximum component + scalar maxCmpt = -GREAT; forAll(f, fp) { - maxLenSqr = max(maxLenSqr, magSqr(points[f[fp]] - cc)); + const point& pt = points[f[fp]]; + maxLenSqr = max(maxLenSqr, magSqr(pt - cc)); + maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt))); } - tols[faceI] = matchTol * Foam::sqrt(maxLenSqr); + tols[faceI] = max(SMALL*maxCmpt, matchTol*Foam::sqrt(maxLenSqr)); } return tols; } From 21048d96005bb8371806eaa2ccc93a222a184d78 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 13 Jan 2010 19:10:54 +0000 Subject: [PATCH 06/24] Disabled writing zero-sized faceZones since upset tecio library --- .../foamToTecplot360/foamToTecplot360.C | 145 ++++++++++-------- 1 file changed, 77 insertions(+), 68 deletions(-) diff --git a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C index 6c1065a489..44384b7efe 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C +++ b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C @@ -1077,83 +1077,92 @@ int main(int argc, char *argv[]) { const faceZone& pp = zones[zoneI]; - const indirectPrimitivePatch ipp - ( - IndirectList(mesh.faces(), pp), - mesh.points() - ); + if (pp.size() > 0) + { + const indirectPrimitivePatch ipp + ( + IndirectList(mesh.faces(), pp), + mesh.points() + ); - writer.writePolygonalZone - ( - pp.name(), - strandID++, //1+patchIDs.size()+zoneI, //strandID, - ipp, - allVarLocation - ); + writer.writePolygonalZone + ( + pp.name(), + strandID++, //1+patchIDs.size()+zoneI, //strandID, + ipp, + allVarLocation + ); - // Write coordinates - writer.writeField(ipp.localPoints().component(0)()); - writer.writeField(ipp.localPoints().component(1)()); - writer.writeField(ipp.localPoints().component(2)()); + // Write coordinates + writer.writeField(ipp.localPoints().component(0)()); + writer.writeField(ipp.localPoints().component(1)()); + writer.writeField(ipp.localPoints().component(2)()); - // Write all volfields - forAll(vsf, i) - { - writer.writeField - ( - writer.getFaceField + // Write all volfields + forAll(vsf, i) + { + writer.writeField ( - linearInterpolate(vsf[i])(), - pp - )() - ); - } - forAll(vvf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vsf[i])(), + pp + )() + ); + } + forAll(vvf, i) + { + writer.writeField ( - linearInterpolate(vvf[i])(), - pp - )() - ); - } - forAll(vSpheretf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vvf[i])(), + pp + )() + ); + } + forAll(vSpheretf, i) + { + writer.writeField ( - linearInterpolate(vSpheretf[i])(), - pp - )() - ); - } - forAll(vSymmtf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vSpheretf[i])(), + pp + )() + ); + } + forAll(vSymmtf, i) + { + writer.writeField ( - linearInterpolate(vSymmtf[i])(), - pp - )() - ); - } - forAll(vtf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vSymmtf[i])(), + pp + )() + ); + } + forAll(vtf, i) + { + writer.writeField ( - linearInterpolate(vtf[i])(), - pp - )() - ); - } + writer.getFaceField + ( + linearInterpolate(vtf[i])(), + pp + )() + ); + } - writer.writeConnectivity(ipp); + writer.writeConnectivity(ipp); + } + else + { + Info<< " Skipping zero sized faceZone " << zoneI + << "\t" << pp.name() + << nl << endl; + } } } From c9033fc36fc2a8474b0e7dccb85fee933e775700 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:14:06 +0000 Subject: [PATCH 07/24] Added optional settings at comment --- .../utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index c7f2d5b6f3..783469ee7d 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -40,6 +40,10 @@ geometry { type triSurfaceMesh; + //tolerance 1E-6; // optional:non-default tolerance on intersections + //maxTreeDepth 10; // optional:depth of octree. Decrease only in case + // of memory limitations. + // Per region the patchname. If not provided will be _. regions { From 7f64e0271d9aeafaab4ad482be508fedb9532cd2 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:14:42 +0000 Subject: [PATCH 08/24] Changed comment in header --- .../scotchDecomp/scotchDecomp.C | 111 +++++++++++++----- 1 file changed, 79 insertions(+), 32 deletions(-) diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C index c14e7e39f1..53931689cb 100644 --- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C +++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C @@ -22,40 +22,87 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA - Default strategy: - Strat=b - { - job=t, - map=t, - poli=S, - sep= - ( - m + From scotch forum: + + By: Francois PELLEGRINI RE: Graph mapping 'strategy' string [ reply ] + 2008-08-22 10:09 Strategy handling in Scotch is a bit tricky. In order + not to be confused, you must have a clear view of how they are built. + Here are some rules: + + 1- Strategies are made up of "methods" which are combined by means of + "operators". + + 2- A method is of the form "m{param=value,param=value,...}", where "m" + is a single character (this is your first error: "f" is a method name, + not a parameter name). + + 3- There exist different sort of strategies : bipartitioning strategies, + mapping strategies, ordering strategies, which cannot be mixed. For + instance, you cannot build a bipartitioning strategy and feed it to a + mapping method (this is your second error). + + To use the "mapCompute" routine, you must create a mapping strategy, not + a bipartitioning one, and so use stratGraphMap() and not + stratGraphBipart(). Your mapping strategy should however be based on the + "recursive bipartitioning" method ("b"). For instance, a simple (and + hence not very efficient) mapping strategy can be : + + "b{sep=f}" + + which computes mappings with the recursive bipartitioning method "b", + this latter using the Fiduccia-Mattheyses method "f" to compute its + separators. + + If you want an exact partition (see your previous post), try + "b{sep=fx}". + + However, these strategies are not the most efficient, as they do not + make use of the multi-level framework. + + To use the multi-level framework, try for instance: + + "b{sep=m{vert=100,low=h,asc=f}x}" + + The current default mapping strategy in Scotch can be seen by using the + "-vs" option of program gmap. It is, to date: + + b + { + job=t, + map=t, + poli=S, + sep= + ( + m + { + asc=b { - asc=b - { - bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005}, - org=f{move=80,pass=-1,bal=0.005},width=3 - }, - low=h{pass=10}f{move=80,pass=-1,bal=0.0005}, - type=h, - vert=80, - rat=0.8 - } - | m + bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005}, + org=f{move=80,pass=-1,bal=0.005}, + width=3 + }, + low=h{pass=10}f{move=80,pass=-1,bal=0.0005}, + type=h, + vert=80, + rat=0.8 + } + | m + { + asc=b { - asc=b - { - bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005}, - org=f{move=80,pass=-1,bal=0.005}, - width=3}, - low=h{pass=10}f{move=80,pass=-1,bal=0.0005}, - type=h, - vert=80, - rat=0.8 - } - ) - } + bnd=d{pass=40,dif=1,rem=1}f{move=80,pass=-1,bal=0.005}, + org=f{move=80,pass=-1,bal=0.005}, + width=3 + }, + low=h{pass=10}f{move=80,pass=-1,bal=0.0005}, + type=h, + vert=80, + rat=0.8 + } + ) + } + + \*---------------------------------------------------------------------------*/ #include "scotchDecomp.H" From b7ba2d273a0464895dc4b9f1b7c721199169a200 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:15:42 +0000 Subject: [PATCH 09/24] Corrected printing to be on master only --- .../mesh/manipulation/setSet/setSet.C | 44 +++++++++---------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C index 7833b4ec24..cf24a1c447 100644 --- a/applications/utilities/mesh/manipulation/setSet/setSet.C +++ b/applications/utilities/mesh/manipulation/setSet/setSet.C @@ -94,7 +94,7 @@ void backup { if (fromSet.size()) { - Pout<< " Backing up " << fromName << " into " << toName << endl; + Info<< " Backing up " << fromName << " into " << toName << endl; topoSet::New(setType, mesh, toName, fromSet)().write(); } @@ -525,7 +525,7 @@ bool doCommand { topoSet& currentSet = currentSetPtr(); - Pout<< " Set:" << currentSet.name() + Info<< " Set:" << currentSet.name() << " Size:" << currentSet.size() << " Action:" << actionName << endl; @@ -622,7 +622,7 @@ bool doCommand + ".vtk" ); - Pout<< " Writing " << currentSet.name() + Info<< " Writing " << currentSet.name() << " (size " << currentSet.size() << ") to " << currentSet.instance()/currentSet.local() /currentSet.name() @@ -634,7 +634,7 @@ bool doCommand } else { - Pout<< " Writing " << currentSet.name() + Info<< " Writing " << currentSet.name() << " (size " << currentSet.size() << ") to " << currentSet.instance()/currentSet.local() /currentSet.name() << endl << endl; @@ -683,7 +683,7 @@ enum commandStatus void printMesh(const Time& runTime, const polyMesh& mesh) { - Pout<< "Time:" << runTime.timeName() + Info<< "Time:" << runTime.timeName() << " cells:" << mesh.nCells() << " faces:" << mesh.nFaces() << " points:" << mesh.nPoints() @@ -703,19 +703,19 @@ commandStatus parseType { if (setType.empty()) { - Pout<< "Type 'help' for usage information" << endl; + Info<< "Type 'help' for usage information" << endl; return INVALID; } else if (setType == "help") { - printHelp(Pout); + printHelp(Info); return INVALID; } else if (setType == "list") { - printAllSets(mesh, Pout); + printAllSets(mesh, Info); return INVALID; } @@ -726,7 +726,7 @@ commandStatus parseType label nearestIndex = Time::findClosestTimeIndex(Times, requestedTime); - Pout<< "Changing time from " << runTime.timeName() + Info<< "Changing time from " << runTime.timeName() << " to " << Times[nearestIndex].name() << endl; @@ -737,24 +737,24 @@ commandStatus parseType { case polyMesh::UNCHANGED: { - Pout<< " mesh not changed." << endl; + Info<< " mesh not changed." << endl; break; } case polyMesh::POINTS_MOVED: { - Pout<< " points moved; topology unchanged." << endl; + Info<< " points moved; topology unchanged." << endl; break; } case polyMesh::TOPO_CHANGE: { - Pout<< " topology changed; patches unchanged." << nl + Info<< " topology changed; patches unchanged." << nl << " "; printMesh(runTime, mesh); break; } case polyMesh::TOPO_PATCH_CHANGE: { - Pout<< " topology changed and patches changed." << nl + Info<< " topology changed and patches changed." << nl << " "; printMesh(runTime, mesh); @@ -773,7 +773,7 @@ commandStatus parseType } else if (setType == "quit") { - Pout<< "Quitting ..." << endl; + Info<< "Quitting ..." << endl; return QUIT; } @@ -864,7 +864,7 @@ int main(int argc, char *argv[]) printMesh(runTime, mesh); // Print current sets - printAllSets(mesh, Pout); + printAllSets(mesh, Info); @@ -874,7 +874,7 @@ int main(int argc, char *argv[]) { fileName batchFile(args.option("batch")); - Pout<< "Reading commands from file " << batchFile << endl; + Info<< "Reading commands from file " << batchFile << endl; // we cannot handle .gz files if (!isFile(batchFile, false)) @@ -888,11 +888,11 @@ int main(int argc, char *argv[]) #if READLINE != 0 else if (!read_history(historyFile)) { - Pout<< "Successfully read history from " << historyFile << endl; + Info<< "Successfully read history from " << historyFile << endl; } #endif - Pout<< "Please type 'help', 'quit' or a set command after prompt." << endl; + Info<< "Please type 'help', 'quit' or a set command after prompt." << endl; bool ok = true; @@ -916,7 +916,7 @@ int main(int argc, char *argv[]) { if (!fileStreamPtr->good()) { - Pout<< "End of batch file" << endl; + Info<< "End of batch file" << endl; break; } @@ -924,7 +924,7 @@ int main(int argc, char *argv[]) if (rawLine.size()) { - Pout<< "Doing:" << rawLine << endl; + Info<< "Doing:" << rawLine << endl; } } else @@ -945,7 +945,7 @@ int main(int argc, char *argv[]) } # else { - Pout<< "Command>" << flush; + Info<< "Command>" << flush; std::getline(std::cin, rawLine); } # endif @@ -992,7 +992,7 @@ int main(int argc, char *argv[]) delete fileStreamPtr; } - Pout<< "\nEnd\n" << endl; + Info<< "\nEnd\n" << endl; return 0; } From d0b95d6949990ba8040b8e18565e0580526d2e0a Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:16:07 +0000 Subject: [PATCH 10/24] Added mesh region option --- .../utilities/mesh/manipulation/stitchMesh/stitchMesh.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C index 70927ca8f3..553ea07138 100644 --- a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C +++ b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C @@ -96,6 +96,7 @@ void checkPatch(const polyBoundaryMesh& bMesh, const word& name) int main(int argc, char *argv[]) { argList::noParallel(); +# include "addRegionOption.H" argList::validArgs.append("masterPatch"); argList::validArgs.append("slavePatch"); @@ -109,7 +110,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" runTime.functionObjects().off(); -# include "createMesh.H" +# include "createNamedMesh.H" const word oldInstance = mesh.pointsInstance(); From 06b1695962e6f0c3dc0619226f3ef006e7673741 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:17:05 +0000 Subject: [PATCH 11/24] Added conversion of zero-sized zones --- .../foamFormatConvert/foamFormatConvert.C | 100 +++++++++++++++++- 1 file changed, 97 insertions(+), 3 deletions(-) diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C index 7e4363ffb0..cd51599ca1 100644 --- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C +++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C @@ -31,6 +31,14 @@ Description Mainly used to convert binary mesh/field files to ASCII. + Problem: any zero-size List written binary gets written as '0'. When + reading the file as a dictionary this is interpreted as a label. This + is (usually) not a problem when doing patch fields since these get the + 'uniform', 'nonuniform' prefix. However zone contents are labelLists + not labelFields and these go wrong. For now hacked a solution where + we detect the keywords in zones and redo the dictionary entries + to be labelLists. + \*---------------------------------------------------------------------------*/ #include "argList.H" @@ -56,6 +64,82 @@ namespace Foam } +// Hack to do zones which have Lists in them. See above. +bool writeZones(const word& name, Time& runTime) +{ + IOobject io + ( + name, + runTime.timeName(), + polyMesh::meshSubDir, + runTime, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ); + + bool writeOk = false; + + if (io.headerOk()) + { + Info<< " Reading " << io.headerClassName() + << " : " << name << endl; + + // Switch off type checking (for reading e.g. faceZones as + // generic list of dictionaries). + const word oldTypeName = IOPtrList::typeName; + const_cast(IOPtrList::typeName) = word::null; + + IOPtrList meshObject(io); + + forAll(meshObject, i) + { + if (meshObject[i].isDict()) + { + dictionary& d = meshObject[i].dict(); + + if (d.found("faceLabels")) + { + d.set("faceLabels", labelList(d.lookup("faceLabels"))); + } + + if (d.found("flipMap")) + { + d.set("flipMap", boolList(d.lookup("flipMap"))); + } + + if (d.found("cellLabels")) + { + d.set("cellLabels", labelList(d.lookup("cellLabels"))); + } + + if (d.found("pointLabels")) + { + d.set("pointLabels", labelList(d.lookup("pointLabels"))); + } + } + } + + const_cast(IOPtrList::typeName) = oldTypeName; + // Fake type back to what was in field + const_cast(meshObject.type()) = io.headerClassName(); + + Info<< " Writing " << name << endl; + + // Force writing as ascii + writeOk = meshObject.regIOobject::writeObject + ( + IOstream::ASCII, + IOstream::currentVersion, + runTime.writeCompression() + ); + } + + return writeOk; +} + + + // Main program: int main(int argc, char *argv[]) @@ -76,9 +160,19 @@ int main(int argc, char *argv[]) writeMeshObject("neighbour", runTime); writeMeshObject("faces", runTime); writeMeshObject("points", runTime); - writeMeshObject >("cellZones", runTime); - writeMeshObject >("faceZones", runTime); - writeMeshObject >("pointZones", runTime); + writeMeshObject("pointProcAddressing", runTime); + writeMeshObject("faceProcAddressing", runTime); + writeMeshObject("cellProcAddressing", runTime); + writeMeshObject("boundaryProcAddressing", runTime); + + if (runTime.writeFormat() == IOstream::ASCII) + { + // Only do zones when converting from binary to ascii + // The other way gives problems since working on dictionary level. + writeZones("cellZones", runTime); + writeZones("faceZones", runTime); + writeZones("pointZones", runTime); + } // Get list of objects from the database IOobjectList objects(runTime, runTime.timeName()); From 86743ace66941c73a8595f3fbc266fa7268a78f1 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:17:33 +0000 Subject: [PATCH 12/24] Corrected spelling in comment --- .../pointPatches/derived/coupled/coupledFacePointPatch.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H index a77af7af30..ec5ae840e4 100644 --- a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H +++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H @@ -132,7 +132,7 @@ public: // associated with any faces virtual const labelList& loneMeshPoints() const; - //- Return point unit normals. Not impelemented. + //- Return point unit normals. Not implemented. virtual const vectorField& pointNormals() const; }; From bca38c7fde143c2c6d30067e28a4a82de1401744 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:17:58 +0000 Subject: [PATCH 13/24] Removed unused comment --- src/dynamicMesh/meshCut/directions/directions.C | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C index b50069d606..b77d2cefa4 100644 --- a/src/dynamicMesh/meshCut/directions/directions.C +++ b/src/dynamicMesh/meshCut/directions/directions.C @@ -446,19 +446,4 @@ Foam::directions::directions } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // - - // ************************************************************************* // From 5a87a542d06eebdd1b52620956ffb4f5368595d6 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:18:29 +0000 Subject: [PATCH 14/24] Loosened tolerances for single precision running --- .../motionSmoother/polyMeshGeometry/polyMeshGeometry.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index f42f2a4ef4..f777081d03 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -824,7 +824,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness // (i.e. treat as if mirror cell on other side) vector faceNormal = faceAreas[faceI]; - faceNormal /= mag(faceNormal) + VSMALL; + faceNormal /= mag(faceNormal) + ROOTVSMALL; vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]]; @@ -834,7 +834,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness scalar skewness = mag(faceCentres[faceI] - faceIntersection) - /(2*mag(dWall) + VSMALL); + /(2*mag(dWall) + ROOTVSMALL); // Check if the skewness vector is greater than the PN vector. // This does not cause trouble but is a good indication of a poor From 00b4267a07b985f2979c7958912232ffa28df0ff Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:18:54 +0000 Subject: [PATCH 15/24] Corrected comment --- .../fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H index 062ffa1d86..5e179da0f9 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H @@ -142,7 +142,7 @@ public: // Access - //- Retirn local reference cast into the cyclic patch + //- Return local reference cast into the cyclic patch const cyclicFvPatch& cyclicPatch() const { return cyclicPatch_; From 09aba41a1989ee5b31e686a143b5dc3e0f4dea80 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:19:14 +0000 Subject: [PATCH 16/24] Corrected indentation --- .../extendedStencil/cellToCell/fullStencils/cellToCellStencil.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H index 7fa48ccecb..f8b94c8d5f 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H @@ -95,7 +95,7 @@ protected: class unionEqOp { public: - void operator()( labelList& x, const labelList& y ) const; + void operator()(labelList& x, const labelList& y) const; }; //- Collect cell neighbours of faces in global numbering From 15d0becbb15fc471d8e98a5fe17b047c3087c49b Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:20:42 +0000 Subject: [PATCH 17/24] Changed field writing to have zeroGradient bc for ease of postprocessing --- src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C index 909cf3e683..c4ea71f434 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C @@ -54,6 +54,7 @@ License #include "Random.H" #include "searchableSurfaces.H" #include "treeBoundBox.H" +#include "zeroGradientFvPatchFields.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -2128,7 +2129,8 @@ void Foam::meshRefinement::dumpRefinementLevel() const false ), mesh_, - dimensionedScalar("zero", dimless, 0) + dimensionedScalar("zero", dimless, 0), + zeroGradientFvPatchScalarField::typeName ); const labelList& cellLevel = meshCutter_.cellLevel(); From c182daf90ac8007ffdaf0a857e084fd75dec7a7e Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:21:20 +0000 Subject: [PATCH 18/24] Corrected usage of treeBoundBox::extend --- src/meshTools/meshSearch/meshSearch.C | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C index eb5ec14d8f..d634603649 100644 --- a/src/meshTools/meshSearch/meshSearch.C +++ b/src/meshTools/meshSearch/meshSearch.C @@ -463,7 +463,7 @@ const Foam::indexedOctree& Foam::meshSearch::boundaryTree() treeBoundBox overallBb(mesh_.points()); Random rndGen(123456); - overallBb.extend(rndGen, 1E-4); + overallBb = overallBb.extend(rndGen, 1E-4); overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); @@ -497,7 +497,7 @@ const Foam::indexedOctree& Foam::meshSearch::cellTree() treeBoundBox overallBb(mesh_.points()); Random rndGen(123456); - overallBb.extend(rndGen, 1E-4); + overallBb = overallBb.extend(rndGen, 1E-4); overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); @@ -531,7 +531,7 @@ const Foam::indexedOctree& treeBoundBox overallBb(mesh_.cellCentres()); Random rndGen(123456); - overallBb.extend(rndGen, 1E-4); + overallBb = overallBb.extend(rndGen, 1E-4); overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); From c5a5079582e57399767459e2f52621566baad477 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:21:42 +0000 Subject: [PATCH 19/24] Changed printing --- .../autoHexMeshDriver/autoLayerDriver.C | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index cc49e48939..acd1673392 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -2490,13 +2490,18 @@ void Foam::autoLayerDriver::mergePatchFacesUndo << "---------------------------" << nl << " - which are on the same patch" << nl << " - which make an angle < " << layerParams.featureAngle() + << "- which are on the same patch" << nl + << "- which make an angle < " << layerParams.featureAngle() << " degrees" << nl << " (cos:" << minCos << ')' << nl << " - as long as the resulting face doesn't become concave" + << " (cos:" << minCos << ')' << nl + << "- as long as the resulting face doesn't become concave" << " by more than " << layerParams.concaveAngle() << " degrees" << nl << " (0=straight, 180=fully concave)" << nl + << " (0=straight, 180=fully concave)" << nl << endl; label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict); @@ -2546,11 +2551,6 @@ void Foam::autoLayerDriver::addLayers ); - // Undistorted edge length - const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); - const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); - - // Point-wise extrusion data // ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -2612,6 +2612,10 @@ void Foam::autoLayerDriver::addLayers // Disable extrusion on warped faces // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Undistorted edge length + const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); + const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); + handleWarpedFaces ( pp, @@ -2651,6 +2655,10 @@ void Foam::autoLayerDriver::addLayers } + // Undistorted edge length + const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); + const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); + // Determine (wanted) point-wise layer thickness and expansion ratio scalarField thickness(pp().nPoints()); scalarField minThickness(pp().nPoints()); From 4d4276fd6d9430fe7a14e67addc80ed77dc4db0e Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:33:12 +0000 Subject: [PATCH 20/24] Adapted for new api of searchableSurface --- .../refinementSurfaces/refinementSurfaces.C | 251 ++++++++++++------ .../refinementSurfaces/refinementSurfaces.H | 4 +- 2 files changed, 168 insertions(+), 87 deletions(-) diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C index c910b4c1b5..0f81dff2f5 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C @@ -496,23 +496,63 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const } -// Count number of triangles per surface region -Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s) -{ - const geometricSurfacePatchList& regions = s.patches(); +// // Count number of triangles per surface region +// Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s) +// { +// const geometricSurfacePatchList& regions = s.patches(); +// +// labelList nTris(regions.size(), 0); +// +// forAll(s, triI) +// { +// nTris[s[triI].region()]++; +// } +// return nTris; +// } - labelList nTris(regions.size(), 0); - forAll(s, triI) - { - nTris[s[triI].region()]++; - } - return nTris; -} +// // Pre-calculate the refinement level for every element +// void Foam::refinementSurfaces::wantedRefinementLevel +// ( +// const shellSurfaces& shells, +// const label surfI, +// const List& info, // Indices +// const pointField& ctrs, // Representative coordinate +// labelList& minLevelField +// ) const +// { +// const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; +// +// // Get per element the region +// labelList region; +// geom.getRegion(info, region); +// +// // Initialise fields to region wise minLevel +// minLevelField.setSize(ctrs.size()); +// minLevelField = -1; +// +// forAll(minLevelField, i) +// { +// if (info[i].hit()) +// { +// minLevelField[i] = minLevel(surfI, region[i]); +// } +// } +// +// // Find out if triangle inside shell with higher level +// // What level does shell want to refine fc to? +// labelList shellLevel; +// shells.findHigherLevel(ctrs, minLevelField, shellLevel); +// +// forAll(minLevelField, i) +// { +// minLevelField[i] = max(minLevelField[i], shellLevel[i]); +// } +// } // Precalculate the refinement level for every element of the searchable -// surface. This is currently hardcoded for triSurfaceMesh only. +// surface. void Foam::refinementSurfaces::setMinLevelFields ( const shellSurfaces& shells @@ -522,58 +562,55 @@ void Foam::refinementSurfaces::setMinLevelFields { const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; - if (isA(geom)) + // Precalculation only makes sense if there are different regions + // (so different refinement levels possible) and there are some + // elements. Possibly should have 'enough' elements to have fine + // enough resolution but for now just make sure we don't catch e.g. + // searchableBox (size=6) + if (geom.regions().size() > 1 && geom.globalSize() > 10) { - const triSurfaceMesh& triMesh = refCast(geom); + // Representative local coordinates + const pointField ctrs = geom.coordinates(); - autoPtr minLevelFieldPtr - ( - new triSurfaceLabelField - ( - IOobject - ( - "minLevel", - triMesh.objectRegistry::time().timeName(), // instance - "triSurface", // local - triMesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - triMesh, - dimless - ) - ); - triSurfaceLabelField& minLevelField = minLevelFieldPtr(); - - const triSurface& s = static_cast(triMesh); - - // Initialise fields to region wise minLevel - forAll(s, triI) + labelList minLevelField(ctrs.size(), -1); { - minLevelField[triI] = minLevel(surfI, s[triI].region()); + // Get the element index in a roundabout way. Problem is e.g. + // distributed surface where local indices differ from global + // ones (needed for getRegion call) + List info; + geom.findNearest + ( + ctrs, + scalarField(ctrs.size(), sqr(GREAT)), + info + ); + + // Get per element the region + labelList region; + geom.getRegion(info, region); + + // From the region get the surface-wise refinement level + forAll(minLevelField, i) + { + if (info[i].hit()) + { + minLevelField[i] = minLevel(surfI, region[i]); + } + } } // Find out if triangle inside shell with higher level - pointField fc(s.size()); - forAll(s, triI) - { - fc[triI] = s[triI].centre(s.points()); - } // What level does shell want to refine fc to? labelList shellLevel; - shells.findHigherLevel(fc, minLevelField, shellLevel); + shells.findHigherLevel(ctrs, minLevelField, shellLevel); - forAll(minLevelField, triI) + forAll(minLevelField, i) { - minLevelField[triI] = max - ( - minLevelField[triI], - shellLevel[triI] - ); + minLevelField[i] = max(minLevelField[i], shellLevel[i]); } - // Store field on triMesh - minLevelFieldPtr.ptr()->store(); + // Store minLevelField on surface + const_cast(geom).setField(minLevelField); } } } @@ -601,44 +638,89 @@ void Foam::refinementSurfaces::findHigherIntersection return; } + if (surfaces_.size() == 1) + { + // Optimisation: single segmented surface. No need to duplicate + // point storage. + + label surfI = 0; + + const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; + + // Do intersection test + List intersectionInfo(start.size()); + geom.findLineAny(start, end, intersectionInfo); + + // See if a cached level field available + labelList minLevelField; + geom.getField(intersectionInfo, minLevelField); + bool haveLevelField = + ( + returnReduce(minLevelField.size(), sumOp