From 76bb575f45d388f8fc43e269b8ee0e6944de04c2 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 13 Apr 2009 12:28:43 +0100 Subject: [PATCH 1/4] updated comment --- .../meshes/primitiveShapes/triangle/triangle.H | 11 ++++------- .../polyTopoChange/polyTopoChange/addPatchCellLayer.H | 5 +++-- 2 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 1165046636..3bc3ea38a8 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -166,14 +166,11 @@ public: const intersection::direction dir = intersection::VECTOR ) const; - //- Fast intersection with a ray. Distance is normalised distance - // to the intersection (normalised with the vector magnitude). - // For a hit, the distance is signed. Positive number - // represents the point in front of triangle. In case of miss - // pointHit position is undefined. - // Only defined for VISIBLE, FULL_RAY or + //- Fast intersection with a ray. + // For a hit, the pointHit.distance() is the line parameter t : + // intersection=p+t*q. Only defined for VISIBLE, FULL_RAY or // HALF_RAY. tol increases the virtual size of the triangle - // by a relative factor. + // by a relative factor. inline pointHit intersection ( const point& p, diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H index ea6a7674e4..f03d23e214 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H @@ -294,9 +294,10 @@ public: // - nPointLayers : number of layers per (patch)point // - nFaceLayers : number of layers per (patch) face // - firstDisplacement : displacement per point for first - // layer of points. If zero do not add point. + // layer of points (i.e. nearest to original mesh). If zero + // do not add point. // Layer thicknesses are calculated to constant geometric - // expansion. Use expansionRatio for constant size. + // expansion. Use expansionRatio 1 for constant size. // Sets addedPoints_ which is per pp point a list of points // added. // Note: firstDisplacement has to be parallel synchronised before From f0d3088592804310574f11180907e28aa58e8c94 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 13 Apr 2009 12:31:04 +0100 Subject: [PATCH 2/4] added tolerance --- .../searchableSurface/triSurfaceMesh.C | 53 ++++++++++++++++++- .../searchableSurface/triSurfaceMesh.H | 8 +++ 2 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.C b/src/meshTools/searchableSurface/triSurfaceMesh.C index 877f264341..950cd13285 100644 --- a/src/meshTools/searchableSurface/triSurfaceMesh.C +++ b/src/meshTools/searchableSurface/triSurfaceMesh.C @@ -237,6 +237,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s) ) ), triSurface(s), + tolerance_(indexedOctree::perturbTol()), surfaceClosed_(-1) {} @@ -279,6 +280,7 @@ Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io) searchableSurface::objectPath() ) ), + tolerance_(indexedOctree::perturbTol()), surfaceClosed_(-1) {} @@ -324,6 +326,7 @@ Foam::triSurfaceMesh::triSurfaceMesh searchableSurface::objectPath() ) ), + tolerance_(indexedOctree::perturbTol()), surfaceClosed_(-1) { scalar scaleFactor = 0; @@ -332,8 +335,18 @@ Foam::triSurfaceMesh::triSurfaceMesh // eg, CAD geometries are often done in millimeters if (dict.readIfPresent("scale", scaleFactor) && scaleFactor > 0) { + Info<< searchableSurface::name() << " : using scale " << scaleFactor + << endl; triSurface::scalePoints(scaleFactor); } + + + // Have optional non-standard search tolerance for gappy surfaces. + if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0) + { + Info<< searchableSurface::name() << " : using intersection tolerance " + << tolerance_ << endl; + } } @@ -380,6 +393,9 @@ const Foam::indexedOctree& bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + tree_.reset ( new indexedOctree @@ -391,6 +407,8 @@ const Foam::indexedOctree& 3.0 // duplicity ) ); + + indexedOctree::perturbTol() = oldTol; } return tree_(); @@ -491,6 +509,9 @@ void Foam::triSurfaceMesh::findNearest info.setSize(samples.size()); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + forAll(samples, i) { static_cast(info[i]) = octree.findNearest @@ -499,6 +520,8 @@ void Foam::triSurfaceMesh::findNearest nearestDistSqr[i] ); } + + indexedOctree::perturbTol() = oldTol; } @@ -513,6 +536,9 @@ void Foam::triSurfaceMesh::findLine info.setSize(start.size()); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + forAll(start, i) { static_cast(info[i]) = octree.findLine @@ -521,6 +547,8 @@ void Foam::triSurfaceMesh::findLine end[i] ); } + + indexedOctree::perturbTol() = oldTol; } @@ -535,6 +563,9 @@ void Foam::triSurfaceMesh::findLineAny info.setSize(start.size()); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + forAll(start, i) { static_cast(info[i]) = octree.findLineAny @@ -543,6 +574,8 @@ void Foam::triSurfaceMesh::findLineAny end[i] ); } + + indexedOctree::perturbTol() = oldTol; } @@ -557,6 +590,9 @@ void Foam::triSurfaceMesh::findLineAll info.setSize(start.size()); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + // Work array DynamicList hits; @@ -570,7 +606,7 @@ void Foam::triSurfaceMesh::findLineAll const scalarField magSqrDirVec(magSqr(dirVec)); const vectorField smallVec ( - Foam::sqrt(SMALL)*dirVec + indexedOctree::perturbTol()*dirVec + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL) ); @@ -613,6 +649,8 @@ void Foam::triSurfaceMesh::findLineAll info[pointI].clear(); } } + + indexedOctree::perturbTol() = oldTol; } @@ -649,7 +687,13 @@ void Foam::triSurfaceMesh::getNormal { if (info[i].hit()) { - normal[i] = faceNormals()[info[i].index()]; + label triI = info[i].index(); + //- Cached: + //normal[i] = faceNormals()[triI]; + + //- Uncached + normal[i] = triSurface::operator[](triI).normal(points()); + normal[i] /= mag(normal[i]) + VSMALL; } else { @@ -691,6 +735,9 @@ void Foam::triSurfaceMesh::getVolumeType { volType.setSize(points.size()); + scalar oldTol = indexedOctree::perturbTol(); + indexedOctree::perturbTol() = tolerance_; + forAll(points, pointI) { const point& pt = points[pointI]; @@ -699,6 +746,8 @@ void Foam::triSurfaceMesh::getVolumeType // - cheat conversion since same values volType[pointI] = static_cast(tree().getVolumeType(pt)); } + + indexedOctree::perturbTol() = oldTol; } diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H index b2ec61c14a..4ae415fafa 100644 --- a/src/meshTools/searchableSurface/triSurfaceMesh.H +++ b/src/meshTools/searchableSurface/triSurfaceMesh.H @@ -28,6 +28,11 @@ Class Description IOoject and searching on triSurface + Note: when constructing from dictionary has optional parameters: + - scale : scaling factor. + - tolerance : relative tolerance for doing intersections + (see triangle::intersection) + SourceFiles triSurfaceMesh.C @@ -63,6 +68,9 @@ private: // Private member data + //- Optional tolerance to use in searches + scalar tolerance_; + //- Search tree (triangles) mutable autoPtr > tree_; From 79d4641b374bb253fab5bc04a044887fc257b03b Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 14 Apr 2009 12:30:49 +0100 Subject: [PATCH 3/4] cleanup; introduced UIndirectList --- .../mapDistribute/mapDistributeTemplates.C | 66 +++---------------- 1 file changed, 8 insertions(+), 58 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C index 08e12b876e..330f48be00 100644 --- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C +++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C @@ -52,24 +52,13 @@ void Foam::mapDistribute::distribute if (domain != Pstream::myProcNo() && map.size()) { - List subField(map.size()); - forAll(map, i) - { - subField[i] = field[map[i]]; - } OPstream toNbr(Pstream::blocking, domain); - toNbr << subField; + toNbr << UIndirectList(field, map); } } // Subset myself - const labelList& mySubMap = subMap[Pstream::myProcNo()]; - - List subField(mySubMap.size()); - forAll(mySubMap, i) - { - subField[i] = field[mySubMap[i]]; - } + UIndirectList subField(field, subMap[Pstream::myProcNo()]); // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -126,13 +115,7 @@ void Foam::mapDistribute::distribute List newField(constructSize); // Subset myself - const labelList& mySubMap = subMap[Pstream::myProcNo()]; - - List subField(mySubMap.size()); - forAll(mySubMap, i) - { - subField[i] = field[mySubMap[i]]; - } + UIndirectList subField(field, subMap[Pstream::myProcNo()]); // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -152,16 +135,8 @@ void Foam::mapDistribute::distribute if (Pstream::myProcNo() == sendProc) { // I am sender. Send to recvProc. - const labelList& map = subMap[recvProc]; - - List subField(map.size()); - forAll(map, i) - { - subField[i] = field[map[i]]; - } - OPstream toNbr(Pstream::scheduled, recvProc); - toNbr << subField; + toNbr << UIndirectList(field, subMap[recvProc]); } else { @@ -374,24 +349,13 @@ void Foam::mapDistribute::distribute if (domain != Pstream::myProcNo() && map.size()) { - List subField(map.size()); - forAll(map, i) - { - subField[i] = field[map[i]]; - } OPstream toNbr(Pstream::blocking, domain); - toNbr << subField; + toNbr << UIndirectList(field, map); } } // Subset myself - const labelList& mySubMap = subMap[Pstream::myProcNo()]; - - List subField(mySubMap.size()); - forAll(mySubMap, i) - { - subField[i] = field[mySubMap[i]]; - } + UIndirectList subField(field, subMap[Pstream::myProcNo()]); // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -449,13 +413,7 @@ void Foam::mapDistribute::distribute List newField(constructSize, nullValue); // Subset myself - const labelList& mySubMap = subMap[Pstream::myProcNo()]; - - List subField(mySubMap.size()); - forAll(mySubMap, i) - { - subField[i] = field[mySubMap[i]]; - } + UIndirectList subField(field, subMap[Pstream::myProcNo()]); // Receive sub field from myself (subField) const labelList& map = constructMap[Pstream::myProcNo()]; @@ -475,16 +433,8 @@ void Foam::mapDistribute::distribute if (Pstream::myProcNo() == sendProc) { // I am sender. Send to recvProc. - const labelList& map = subMap[recvProc]; - - List subField(map.size()); - forAll(map, i) - { - subField[i] = field[map[i]]; - } - OPstream toNbr(Pstream::scheduled, recvProc); - toNbr << subField; + toNbr << UIndirectList(field, subMap[recvProc]); } else { From 0825825008c2f4add88b584b211191f9f3bfc4a4 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 14 Apr 2009 12:31:38 +0100 Subject: [PATCH 4/4] relaxed mesh quality constraints; absolute layer size --- .../autoHexMeshDriver/autoLayerDriver.C | 292 ++++++++++++------ .../autoHexMeshDriver/autoLayerDriver.H | 7 +- .../layerParameters/layerParameters.C | 31 +- .../layerParameters/layerParameters.H | 34 +- 4 files changed, 236 insertions(+), 128 deletions(-) diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index 0859e37d33..ee6cea3e12 100644 --- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -45,6 +45,7 @@ Description #include "OFstream.H" #include "layerParameters.H" #include "combineFaces.H" +#include "IOmanip.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -1414,8 +1415,11 @@ void Foam::autoLayerDriver::calculateLayerThickness const indirectPrimitivePatch& pp, const labelList& patchIDs, const scalarField& patchExpansionRatio, - const scalarField& patchFinalLayerRatio, - const scalarField& patchRelMinThickness, + + const bool relativeSizes, + const scalarField& patchFinalLayerThickness, + const scalarField& patchMinThickness, + const labelList& cellLevel, const labelList& patchNLayers, const scalar edge0Len, @@ -1428,22 +1432,100 @@ void Foam::autoLayerDriver::calculateLayerThickness const fvMesh& mesh = meshRefiner_.mesh(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); - if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2) + + // Rework patch-wise layer parameters into minimum per point + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // Reuse input fields + expansionRatio.setSize(pp.nPoints()); + expansionRatio = GREAT; + thickness.setSize(pp.nPoints()); + thickness = GREAT; + minThickness.setSize(pp.nPoints()); + minThickness = GREAT; + + forAll(patchIDs, i) { - FatalErrorIn("calculateLayerThickness(..)") - << "Thickness should be factor of local undistorted cell size." - << " Valid values are [0..2]." << nl - << " minThickness:" << patchRelMinThickness - << exit(FatalError); + label patchI = patchIDs[i]; + + const labelList& meshPoints = patches[patchI].meshPoints(); + + forAll(meshPoints, patchPointI) + { + label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]]; + + expansionRatio[ppPointI] = min + ( + expansionRatio[ppPointI], + patchExpansionRatio[patchI] + ); + thickness[ppPointI] = min + ( + thickness[ppPointI], + patchFinalLayerThickness[patchI] + ); + minThickness[ppPointI] = min + ( + minThickness[ppPointI], + patchMinThickness[patchI] + ); + } } + syncTools::syncPointList + ( + mesh, + pp.meshPoints(), + expansionRatio, + minEqOp(), + GREAT, // null value + false // no separation + ); + syncTools::syncPointList + ( + mesh, + pp.meshPoints(), + thickness, + minEqOp(), + GREAT, // null value + false // no separation + ); + syncTools::syncPointList + ( + mesh, + pp.meshPoints(), + minThickness, + minEqOp(), + GREAT, // null value + false // no separation + ); - // Per point the max cell level of connected cells - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - labelList maxPointLevel(pp.nPoints(), labelMin); + // Now the thicknesses are set according to the minimum of connected + // patches. + + // Rework relative thickness into absolute + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // by multiplying with the internal cell size. + + if (relativeSizes) { + if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2) + { + FatalErrorIn("calculateLayerThickness(..)") + << "Thickness should be factor of local undistorted cell size." + << " Valid values are [0..2]." << nl + << " minThickness:" << patchMinThickness + << exit(FatalError); + } + + + // Determine per point the max cell level of connected cells + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + labelList maxPointLevel(pp.nPoints(), labelMin); + forAll(pp, i) { label ownLevel = cellLevel[mesh.faceOwner()[pp.addressing()[i]]]; @@ -1465,113 +1547,44 @@ void Foam::autoLayerDriver::calculateLayerThickness labelMin, // null value false // no separation ); - } - // Rework patch-wise layer parameters into minimum per point - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - expansionRatio.setSize(pp.nPoints()); - expansionRatio = GREAT; - scalarField finalLayerRatio(pp.nPoints(), GREAT); - scalarField relMinThickness(pp.nPoints(), GREAT); - - { - forAll(patchIDs, i) + forAll(maxPointLevel, pointI) { - label patchI = patchIDs[i]; - - const labelList& meshPoints = patches[patchI].meshPoints(); - - forAll(meshPoints, patchPointI) - { - label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]]; - - expansionRatio[ppPointI] = min - ( - expansionRatio[ppPointI], - patchExpansionRatio[patchI] - ); - finalLayerRatio[ppPointI] = min - ( - finalLayerRatio[ppPointI], - patchFinalLayerRatio[patchI] - ); - relMinThickness[ppPointI] = min - ( - relMinThickness[ppPointI], - patchRelMinThickness[patchI] - ); - } + // Find undistorted edge size for this level. + scalar edgeLen = edge0Len/(1<(), - GREAT, // null value - false // no separation - ); - syncTools::syncPointList - ( - mesh, - pp.meshPoints(), - finalLayerRatio, - minEqOp(), - GREAT, // null value - false // no separation - ); - syncTools::syncPointList - ( - mesh, - pp.meshPoints(), - relMinThickness, - minEqOp(), - GREAT, // null value - false // no separation - ); } - // Per mesh point the expansion parameters - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Rework thickness (of final layer) into overall thickness of all layers + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - thickness.setSize(pp.nPoints()); - minThickness.setSize(pp.nPoints()); - - forAll(maxPointLevel, pointI) + forAll(thickness, pointI) { - // Find undistorted edge size for this level. - scalar edgeLen = edge0Len/(1< 0) + { + if (expansionRatio[ppPointI] == 1) + { + sumNearWallThickness += thickness[ppPointI]/nLay; + } + else + { + scalar s = + (1.0-pow(expansionRatio[ppPointI], nLay)) + / (1.0-expansionRatio[ppPointI]); + sumNearWallThickness += thickness[ppPointI]/s; + } + } + } + + label totNPoints = returnReduce(meshPoints.size(), sumOp