zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

mesh driver rewrite; initial distributed surfaces

Browse Source
This commit is contained in:
mattijs
2008-07-07 18:41:09 +01:00
parent b7bb04e81b
commit 76f2d62bde
42 changed files with 7945 additions and 4825 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/autoMesh/Make/files
View File

@ -1,15 +1,21 @@
autoHexMesh = autoHexMesh
autoHexMeshDriver = $(autoHexMesh)/autoHexMeshDriver
$(autoHexMeshDriver)/autoHexMeshDriver.C
$(autoHexMeshDriver)/autoHexMeshDriverLayers.C
$(autoHexMeshDriver)/autoHexMeshDriverShrink.C
$(autoHexMeshDriver)/autoHexMeshDriverSnap.C
$(autoHexMeshDriver)/layerParameters/layerParameters.C
$(autoHexMeshDriver)/refinementParameters/refinementParameters.C
$(autoHexMeshDriver)/snapParameters/snapParameters.C
$(autoHexMeshDriver)/pointData/pointData.C
$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriver.C
$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverLayers.C
$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverShrink.C
$(autoHexMesh)/autoHexMeshDriver/autoHexMeshDriverSnap.C
$(autoHexMesh)/autoHexMeshDriver/pointData/pointData.C
$(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
$(autoHexMesh)/meshRefinement/meshRefinement.C
$(autoHexMesh)/meshRefinement/meshRefinementMerge.C
$(autoHexMesh)/meshRefinement/meshRefinementRefine.C
$(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
$(autoHexMesh)/shellSurfaces/shellSurfaces.C
$(autoHexMesh)/trackedParticle/trackedParticle.C
$(autoHexMesh)/trackedParticle/trackedParticleCloud.C

1230
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C
View File

File diff suppressed because it is too large Load Diff

145
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.H
View File

@ -49,8 +49,9 @@ SourceFiles
#include "wallPoint.H"
#include "indirectPrimitivePatch.H"
#include "featureEdgeMesh.H"
#include "searchableSurface.H"
#include "searchableSurfaces.H"
#include "refinementSurfaces.H"
#include "shellSurfaces.H"
#include "meshRefinement.H"
#include "decompositionMethod.H"
#include "fvMeshDistribute.H"
@ -70,6 +71,9 @@ class pointData;
class faceSet;
class addPatchCellLayer;
class mapDistributePolyMesh;
class snapParameters;
class layerParameters;
class refinementParameters;
/*---------------------------------------------------------------------------*\
Class autoHexMeshDriver Declaration
@ -79,15 +83,13 @@ class autoHexMeshDriver
{
// Static data members
//- Default angle for faces to be convcave
static const scalar defaultConcaveAngle;
//- Extrusion controls
enum extrudeMode
{
NOEXTRUDE, /*!< Do not extrude. No layers added. */
EXTRUDE, /*!< Extrude */
EXTRUDEREMOVE /*!< Extrude but afterwards remove added faces locally */
EXTRUDEREMOVE /*!< Extrude but afterwards remove added */
/*!< faces locally */
};
@ -141,42 +143,24 @@ class autoHexMeshDriver
//- Debug level
const label debug_;
//- Total number of cells
const label maxGlobalCells_;
//- Per processor max number of cells
const label maxLocalCells_;
//- When to stop refining
const label minRefineCells_;
//- Curvature
const scalar curvature_;
//- Number of layers between different refinement levels
const label nBufferLayers_;
//- Areas to keep
const pointField keepPoints_;
//- Merge distance
const scalar mergeDist_;
// Refinement only
//- Explicit features
PtrList<featureEdgeMesh> featureMeshes_;
//- Per feature the refinement level
labelList featureLevels_;
//- Explicit features
PtrList<featureEdgeMesh> featureMeshes_;
//- Per feature the refinement level
labelList featureLevels_;
//- Shells
PtrList<searchableSurface> shells_;
//- Per shell the refinement level
labelList shellLevels_;
//- Per shell whether to refine inside or outside
boolList shellRefineInside_;
//- All surface based geometry
autoPtr<searchableSurfaces> allGeometryPtr_;
//- Surfaces with refinement levels built-in
//- Shells (geometry for inside/outside refinement)
autoPtr<shellSurfaces> shellsPtr_;
//- Surfaces (geometry for intersection based refinement)
autoPtr<refinementSurfaces> surfacesPtr_;
//- Per refinement surface region the patch
@ -200,10 +184,6 @@ class autoHexMeshDriver
//- Calculate merge distance. Check against writing tolerance.
scalar getMergeDistance(const scalar mergeTol) const;
// Return per keeppoint -1 or the local cell label the point is in.
// Guaranteed to be only on one processor.
labelList findCells(const pointField&) const;
static void orientOutside(PtrList<searchableSurface>&);
//- Read feature edges
@ -371,6 +351,7 @@ class autoHexMeshDriver
// layers per surface.
void setNumLayers
(
const labelList& patchToNLayers,
const labelList& patchIDs,
const indirectPrimitivePatch& pp,
pointField& patchDisp,
@ -390,17 +371,22 @@ class autoHexMeshDriver
//- Calculate pointwise wanted and minimum thickness.
// thickness: wanted thickness
// minthickness: when to give up and not extrude
// Gets per patch parameters and determine pp pointwise
// parameters.
void calculateLayerThickness
(
const scalar expansionRatio,
const scalar finalLayerRatio,
const scalar relMinThickness,
const indirectPrimitivePatch& pp,
const labelList& patchIDs,
const scalarField& patchExpansionRatio,
const scalarField& patchFinalLayerRatio,
const scalarField& patchRelMinThickness,
const labelList& cellLevel,
const labelList& patchNLayers,
const scalar edge0Len,
scalarField& thickness,
scalarField& minThickness
scalarField& minThickness,
scalarField& expansionRatio
) const;
@ -638,13 +624,13 @@ public:
// Constructors
//- Construct from dictionary and mesh to modify
autoHexMeshDriver
(
fvMesh& mesh,
const dictionary& meshDict,
const dictionary& decomposeDict
);
// //- Construct from dictionary and mesh to modify
// autoHexMeshDriver
// (
// fvMesh& mesh,
// const dictionary& meshDict,
// const dictionary& decomposeDict
// );
//- Alternative constructor from top-level controldictionary and
// refinement specific dictionary
@ -652,6 +638,7 @@ public:
(
fvMesh& mesh,
const dictionary& controlDict,
const dictionary& geometryDict,
const dictionary& refineDict,
const dictionary& decomposeDict
);
@ -677,6 +664,12 @@ public:
return surfacesPtr_();
}
//- Surfaces to volume refinement on
const shellSurfaces& shells() const
{
return shellsPtr_();
}
//- Per refinementsurface, per region the patch
const labelList& globalToPatch() const
{
@ -689,32 +682,53 @@ public:
//- Refine around explicit feature edges
label featureEdgeRefine
(
const refinementParameters& refineParams,
const PtrList<dictionary>& featDicts,
const label maxIter,
const label minRefine
);
//- Refine at surface intersections
label surfaceOnlyRefine(const label maxIter);
label surfaceOnlyRefine
(
const refinementParameters& refineParams,
const label maxIter
);
//- Remove cells not reachable from keep points
void removeInsideCells(const label nBufferLayers);
void removeInsideCells
(
const refinementParameters& refineParams,
const label nBufferLayers
);
//- Refine volume regions (interior of shells)
label shellRefine(const label maxIter);
label shellRefine
(
const refinementParameters& refineParams,
const label maxIter
);
//- Introduce baffles; remove unreachable mesh parts
// handleSnapProblems : whether to remove free floating cells
void baffleAndSplitMesh(const bool handleSnapProblems);
void baffleAndSplitMesh
(
const refinementParameters& refineParams,
const bool handleSnapProblems
);
//- Move cells to zones
void zonify();
void zonify(const refinementParameters&);
//- Split and recombine baffles to get rid of single face baffles.
void splitAndMergeBaffles(const bool handleSnapProblems);
void splitAndMergeBaffles
(
const refinementParameters& refineParams,
const bool handleSnapProblems
);
//- Merge multiple boundary faces on single cell
void mergePatchFaces();
void mergePatchFaces(const refinementParameters&);
//- Redecompose according to cell count
// keepZoneFaces : find all faceZones from zoned surfaces and keep
@ -742,7 +756,7 @@ public:
//- Calculate edge length per patch point.
scalarField calcSnapDistance
(
const dictionary& snapDict,
const snapParameters& snapParams,
const indirectPrimitivePatch&
) const;
@ -753,12 +767,19 @@ public:
// of surface points (on castellated mesh) w.r.t. surface.
void preSmoothPatch
(
const dictionary& snapDict,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother&
) const;
//- Get points both on patch and facezone.
labelList getZoneSurfacePoints
(
const indirectPrimitivePatch&,
const word& zoneName
) const;
//- Per patch point calculate point on nearest surface. Set as
// boundary conditions of motionSmoother displacement field. Return
// displacement of patch points.
@ -771,7 +792,7 @@ public:
//- Smooth the displacement field to the internal.
void smoothDisplacement
(
const dictionary& snapDict,
const snapParameters& snapParams,
motionSmoother&
) const;
@ -779,7 +800,7 @@ public:
// locally relax the displacement.
void scaleMesh
(
const dictionary& snapDict,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother&
@ -791,7 +812,7 @@ public:
//- Merge patch faces on same cell.
void mergePatchFacesUndo
(
const dictionary& shrinkDict,
const layerParameters& layerParams,
const dictionary& motionDict
);
@ -801,7 +822,7 @@ public:
//- Add cell layers
void addLayers
(
const dictionary& shrinkDict,
const layerParameters& layerParams,
const dictionary& motionDict,
const label nAllowableErrors,
motionSmoother&

337
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverLayers.C
View File

@ -27,7 +27,6 @@ Description
\*----------------------------------------------------------------------------*/
#include "ListOps.H"
#include "autoHexMeshDriver.H"
#include "fvMesh.H"
#include "Time.H"
@ -43,11 +42,8 @@ Description
#include "mapPolyMesh.H"
#include "addPatchCellLayer.H"
#include "mapDistributePolyMesh.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::autoHexMeshDriver::defaultConcaveAngle = 90;
#include "OFstream.H"
#include "layerParameters.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -1191,6 +1187,7 @@ void Foam::autoHexMeshDriver::handleWarpedFaces
// No extrusion on faces with differing number of layers for points
void Foam::autoHexMeshDriver::setNumLayers
(
const labelList& patchToNLayers,
const labelList& patchIDs,
const indirectPrimitivePatch& pp,
pointField& patchDisp,
@ -1201,18 +1198,6 @@ void Foam::autoHexMeshDriver::setNumLayers
Info<< nl << "Handling points with inconsistent layer specification ..."
<< endl;
const labelList& nSurfLayers = surfaces().numLayers();
// Build map from patch to layers
Map<label> patchToNLayers(nSurfLayers.size());
forAll(nSurfLayers, region)
{
if (globalToPatch_[region] != -1)
{
patchToNLayers.insert(globalToPatch_[region], nSurfLayers[region]);
}
}
// Get for every point (really only nessecary on patch external points)
// the max and min of any patch faces using it.
labelList maxLayers(patchNLayers.size(), labelMin);
@ -1257,7 +1242,7 @@ void Foam::autoHexMeshDriver::setNumLayers
// Unmark any point with different min and max
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nConflicts = 0;
//label nConflicts = 0;
forAll(maxLayers, i)
{
@ -1294,11 +1279,11 @@ void Foam::autoHexMeshDriver::setNumLayers
}
}
reduce(nConflicts, sumOp<label>());
Info<< "Set displacement to zero for " << nConflicts
<< " points due to points being on multiple regions"
<< " with inconsistent nLayers specification." << endl;
//reduce(nConflicts, sumOp<label>());
//
//Info<< "Set displacement to zero for " << nConflicts
// << " points due to points being on multiple regions"
// << " with inconsistent nLayers specification." << endl;
}
@ -1369,59 +1354,137 @@ void Foam::autoHexMeshDriver::growNoExtrusion
}
// Calculate pointwise wanted and minimum thickness.
// thickness: wanted thickness per point
void Foam::autoHexMeshDriver::calculateLayerThickness
(
const scalar expansionRatio,
const scalar finalLayerRatio,
const scalar relMinThickness,
const indirectPrimitivePatch& pp,
const labelList& patchIDs,
const scalarField& patchExpansionRatio,
const scalarField& patchFinalLayerRatio,
const scalarField& patchRelMinThickness,
const labelList& cellLevel,
const labelList& patchNLayers,
const scalar edge0Len,
scalarField& thickness,
scalarField& minThickness
scalarField& minThickness,
scalarField& expansionRatio
) const
{
if (relMinThickness < 0 || relMinThickness > 2)
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
if (min(patchRelMinThickness) < 0 || max(patchRelMinThickness) > 2)
{
FatalErrorIn("calculateLayerThickness(..)")
<< "Thickness should be factor of local undistorted cell size."
<< " Valid values are [0..2]." << nl
<< " minThickness:" << relMinThickness
<< " minThickness:" << patchRelMinThickness
<< exit(FatalError);
}
thickness.setSize(pp.nPoints());
minThickness.setSize(pp.nPoints());
// 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]]];
const face& f = pp.localFaces()[i];
forAll(f, fp)
forAll(pp, i)
{
maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
label ownLevel = cellLevel[mesh_.faceOwner()[pp.addressing()[i]]];
const face& f = pp.localFaces()[i];
forAll(f, fp)
{
maxPointLevel[f[fp]] = max(maxPointLevel[f[fp]], ownLevel);
}
}
syncTools::syncPointList
(
mesh_,
pp.meshPoints(),
maxPointLevel,
maxEqOp<label>(),
labelMin, // null value
false // no separation
);
}
syncTools::syncPointList
(
mesh_,
pp.meshPoints(),
maxPointLevel,
maxEqOp<label>(),
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)
{
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]
);
}
}
syncTools::syncPointList
(
mesh_,
pp.meshPoints(),
expansionRatio,
minEqOp<scalar>(),
GREAT, // null value
false // no separation
);
syncTools::syncPointList
(
mesh_,
pp.meshPoints(),
finalLayerRatio,
minEqOp<scalar>(),
GREAT, // null value
false // no separation
);
syncTools::syncPointList
(
mesh_,
pp.meshPoints(),
relMinThickness,
minEqOp<scalar>(),
GREAT, // null value
false // no separation
);
}
// Per mesh point the expansion parameters
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
thickness.setSize(pp.nPoints());
minThickness.setSize(pp.nPoints());
forAll(maxPointLevel, pointI)
{
@ -1430,19 +1493,24 @@ void Foam::autoHexMeshDriver::calculateLayerThickness
// Calculate layer thickness based on expansion ratio
// and final layer height
if (expansionRatio == 1)
if (expansionRatio[pointI] == 1)
{
thickness[pointI] = finalLayerRatio*patchNLayers[pointI]*edgeLen;
minThickness[pointI] = relMinThickness*edgeLen;
thickness[pointI] =
finalLayerRatio[pointI]
* patchNLayers[pointI]
* edgeLen;
minThickness[pointI] = relMinThickness[pointI]*edgeLen;
}
else
{
scalar invExpansion = 1.0 / expansionRatio;
scalar invExpansion = 1.0 / expansionRatio[pointI];
label nLay = patchNLayers[pointI];
thickness[pointI] = finalLayerRatio*edgeLen
* (1.0 - pow(invExpansion, nLay))
/ (1.0 - invExpansion);
minThickness[pointI] = relMinThickness*edgeLen;
thickness[pointI] =
finalLayerRatio[pointI]
* edgeLen
* (1.0 - pow(invExpansion, nLay))
/ (1.0 - invExpansion);
minThickness[pointI] = relMinThickness[pointI]*edgeLen;
}
}
@ -2309,31 +2377,34 @@ void Foam::autoHexMeshDriver::getLayerCellsFaces
void Foam::autoHexMeshDriver::mergePatchFacesUndo
(
const dictionary& shrinkDict,
const layerParameters& layerParams,
const dictionary& motionDict
)
{
const scalar featureAngle(readScalar(shrinkDict.lookup("featureAngle")));
scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.0);
scalar minCos = Foam::cos
(
layerParams.featureAngle()
* mathematicalConstant::pi/180.0
);
scalar concaveAngle = defaultConcaveAngle;
if (shrinkDict.found("concaveAngle"))
{
concaveAngle = readScalar(shrinkDict.lookup("concaveAngle"));
}
scalar concaveCos = Foam::cos(concaveAngle*mathematicalConstant::pi/180.0);
scalar concaveCos = Foam::cos
(
layerParams.concaveAngle()
* mathematicalConstant::pi/180.0
);
Info<< nl
<< "Merging all faces of a cell" << nl
<< "---------------------------" << nl
<< " - which are on the same patch" << nl
<< " - which make an angle < " << featureAngle << " degrees"
<< " - which make an angle < " << layerParams.featureAngle()
<< " degrees"
<< nl
<< " (cos:" << minCos << ')' << nl
<< " - as long as the resulting face doesn't become concave"
<< " by more than "
<< concaveAngle << " degrees (0=straight, 180=fully concave)" << nl
<< layerParams.concaveAngle()
<< " degrees (0=straight, 180=fully concave)" << nl
<< endl;
label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict);
@ -2344,80 +2415,18 @@ void Foam::autoHexMeshDriver::mergePatchFacesUndo
void Foam::autoHexMeshDriver::addLayers
(
const dictionary& shrinkDict,
const layerParameters& layerParams,
const dictionary& motionDict,
const label nAllowableErrors,
motionSmoother& meshMover
)
{
// Read some more dictionary settings
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Min thickness per cell
const scalar relMinThickness(readScalar(shrinkDict.lookup("minThickness")));
// Warped faces recoginition
const scalar maxFaceThicknessRatio
(
readScalar(shrinkDict.lookup("maxFaceThicknessRatio"))
);
const scalar featureAngle(readScalar(shrinkDict.lookup("featureAngle")));
// Feature angle
const scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.);
// Number of layers for to grow non-extrusion region by
const label nGrow(readLabel(shrinkDict.lookup("nGrow")));
//const label nSmoothDisp(readLabel(shrinkDict.lookup("nSmoothDispl")));
// Snapping iterations
const label nSnap(readLabel(shrinkDict.lookup("nSnap")));
// Mesh termination and medial axis smoothing quantities
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const scalar minCosLayerTermination =
Foam::cos(0.5*featureAngle*mathematicalConstant::pi/180.);
const scalar expansionRatio
(
readScalar(shrinkDict.lookup("expansionRatio"))
);
const scalar finalLayerRatio
(
readScalar(shrinkDict.lookup("finalLayerRatio"))
);
const label nBufferCellsNoExtrude
(
readLabel(shrinkDict.lookup("nBufferCellsNoExtrude"))
);
const label nSmoothSurfaceNormals
(
readLabel(shrinkDict.lookup("nSmoothSurfaceNormals"))
);
const label nSmoothNormals(readLabel(shrinkDict.lookup("nSmoothNormals")));
const label nSmoothThickness
(
readLabel(shrinkDict.lookup("nSmoothThickness"))
);
const scalar maxThicknessToMedialRatio
(
readScalar(shrinkDict.lookup("maxThicknessToMedialRatio"))
);
// Medial axis setup
const scalar minMedianAxisAngle
(
readScalar(shrinkDict.lookup("minMedianAxisAngle"))
);
const scalar minMedianAxisAngleCos =
Foam::cos(minMedianAxisAngle*mathematicalConstant::pi/180.);
// Precalculate mesh edge labels for patch edges
const indirectPrimitivePatch& pp = meshMover.patch();
const labelList& meshPoints = pp.meshPoints();
// Precalculate mesh edge labels for patch edges
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList meshEdges(pp.nEdges());
forAll(pp.edges(), edgeI)
{
@ -2449,8 +2458,9 @@ void Foam::autoHexMeshDriver::addLayers
setNumLayers
(
meshMover.adaptPatchIDs(),
pp,
layerParams.numLayers(), // per patch the num layers
meshMover.adaptPatchIDs(), // patches that are being moved
pp, // indirectpatch for all faces moving
patchDisp,
patchNLayers,
@ -2477,7 +2487,7 @@ void Foam::autoHexMeshDriver::addLayers
(
pp,
meshEdges,
minCos,
layerParams.featureAngle()*mathematicalConstant::pi/180.0,
patchDisp,
patchNLayers,
@ -2494,7 +2504,7 @@ void Foam::autoHexMeshDriver::addLayers
handleWarpedFaces
(
pp,
maxFaceThicknessRatio,
layerParams.maxFaceThicknessRatio(),
edge0Len,
cellLevel,
@ -2517,7 +2527,7 @@ void Foam::autoHexMeshDriver::addLayers
// Grow out region of non-extrusion
for (label i = 0; i < nGrow; i++)
for (label i = 0; i < layerParams.nGrow(); i++)
{
growNoExtrusion
(
@ -2528,21 +2538,24 @@ void Foam::autoHexMeshDriver::addLayers
);
}
// Determine point-wise layer thickness
// Determine (wanted) point-wise layer thickness and expansion ratio
scalarField thickness(pp.nPoints());
scalarField minThickness(pp.nPoints());
scalarField expansionRatio(pp.nPoints());
calculateLayerThickness
(
expansionRatio,
finalLayerRatio,
relMinThickness,
pp,
meshMover.adaptPatchIDs(),
layerParams.expansionRatio(),
layerParams.finalLayerRatio(),
layerParams.minThickness(),
cellLevel,
patchNLayers,
edge0Len,
thickness,
minThickness
minThickness,
expansionRatio
);
// Calculate wall to medial axis distance for smoothing displacement
@ -2602,9 +2615,9 @@ void Foam::autoHexMeshDriver::addLayers
medialAxisSmoothingInfo
(
meshMover,
nSmoothNormals,
nSmoothSurfaceNormals,
minMedianAxisAngleCos,
layerParams.nSmoothNormals(),
layerParams.nSmoothSurfaceNormals(),
layerParams.minMedianAxisAngleCos(),
dispVec,
medialRatio,
@ -2657,8 +2670,8 @@ void Foam::autoHexMeshDriver::addLayers
// debug,
// meshMover,
// -patchDisp, // Shrink in opposite direction of addedPoints
// nSmoothDisp,
// nSnap
// layerParams.nSmoothDisp(),
// layerParams.nSnap()
//);
// Medial axis based shrinking
@ -2666,11 +2679,11 @@ void Foam::autoHexMeshDriver::addLayers
(
meshMover,
nSmoothThickness,
maxThicknessToMedialRatio,
layerParams.nSmoothThickness(),
layerParams.maxThicknessToMedialRatio(),
nAllowableErrors,
nSnap,
minCosLayerTermination,
layerParams.nSnap(),
layerParams.layerTerminationCos(),
thickness,
minThickness,
@ -2735,7 +2748,7 @@ void Foam::autoHexMeshDriver::addLayers
meshMover,
patchNLayers,
extrudeStatus,
nBufferCellsNoExtrude,
layerParams.nBufferCellsNoExtrude(),
nPatchPointLayers,
nPatchFaceLayers
);
@ -2747,7 +2760,7 @@ void Foam::autoHexMeshDriver::addLayers
{
if (patchNLayers[i] > 0)
{
if (expansionRatio == 1.0)
if (expansionRatio[i] == 1.0)
{
firstDisp[i] = patchDisp[i]/nPatchPointLayers[i];
}
@ -2755,15 +2768,15 @@ void Foam::autoHexMeshDriver::addLayers
{
label nLay = nPatchPointLayers[i];
scalar h =
pow(expansionRatio,nLay - 1)
* (mag(patchDisp[i])*(1.0 - expansionRatio))
/ (1.0 - pow(expansionRatio, nLay));
pow(expansionRatio[i], nLay - 1)
* (mag(patchDisp[i])*(1.0 - expansionRatio[i]))
/ (1.0 - pow(expansionRatio[i], nLay));
firstDisp[i] = h/mag(patchDisp[i])*patchDisp[i];
}
}
}
scalar invExpansionRatio = 1.0 / expansionRatio;
scalarField invExpansionRatio = 1.0 / expansionRatio;
// Add topo regardless of whether extrudeStatus is extruderemove.
// Not add layer if patchDisp is zero.

224
src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriverSnap.C
View File

@ -37,6 +37,7 @@ Description
#include "pointEdgePoint.H"
#include "PointEdgeWave.H"
#include "mergePoints.H"
#include "snapParameters.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -95,8 +96,7 @@ Foam::Map<Foam::label> Foam::autoHexMeshDriver::getZoneBafflePatches
label patchI = globalToPatch_[surfaces().globalRegion(surfI, 0)];
Info<< "For surface "
<< surfaces()[surfI].IOobject::name()
//<< surfaces().names()[surfI]
<< surfaces().names()[surfI]
<< " found faceZone " << fZone.name()
<< " and patch " << mesh_.boundaryMesh()[patchI].name()
<< endl;
@ -712,8 +712,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::createZoneBaffles
if (debug_)
{
const_cast<Time&>(mesh_.time())++;
Pout<< "Writing baffled mesh to time " << mesh_.time().timeName()
<< endl;
Pout<< "Writing baffled mesh to time "
<< mesh_.time().timeName() << endl;
mesh_.write();
}
}
@ -755,14 +755,10 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoHexMeshDriver::mergeZoneBaffles
Foam::scalarField Foam::autoHexMeshDriver::calcSnapDistance
(
const dictionary& snapDict,
const snapParameters& snapParams,
const indirectPrimitivePatch& pp
) const
{
// When to snap
scalar snapTol(readScalar(snapDict.lookup("snapTol")));
const edgeList& edges = pp.edges();
const labelListList& pointEdges = pp.pointEdges();
const pointField& localPoints = pp.localPoints();
@ -793,7 +789,7 @@ Foam::scalarField Foam::autoHexMeshDriver::calcSnapDistance
false // no separation
);
return snapTol*maxEdgeLen;
return snapParams.snapTol()*maxEdgeLen;
}
@ -801,7 +797,7 @@ Foam::scalarField Foam::autoHexMeshDriver::calcSnapDistance
Foam::labelList Foam::autoHexMeshDriver::getSurfacePatches() const
{
// Set of patches originating from surface
labelHashSet surfacePatchSet(surfaces().size());
labelHashSet surfacePatchSet(globalToPatch_.size());
forAll(globalToPatch_, i)
{
@ -826,21 +822,21 @@ Foam::labelList Foam::autoHexMeshDriver::getSurfacePatches() const
void Foam::autoHexMeshDriver::preSmoothPatch
(
const dictionary& snapDict,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother& meshMover
) const
{
// Smoothing iterations
label nSmoothPatch(readLabel(snapDict.lookup("nSmoothPatch")));
// Snapping iterations
label nSnap(readLabel(snapDict.lookup("nSnap")));
labelList checkFaces;
Info<< "Smoothing patch points ..." << endl;
for (label smoothIter = 0; smoothIter < nSmoothPatch; smoothIter++)
for
(
label smoothIter = 0;
smoothIter < snapParams.nSmoothPatch();
smoothIter++
)
{
Info<< "Smoothing iteration " << smoothIter << endl;
checkFaces.setSize(mesh_.nFaces());
@ -857,11 +853,11 @@ void Foam::autoHexMeshDriver::preSmoothPatch
scalar oldErrorReduction = -1;
for (label snapIter = 0; snapIter < 2*nSnap; snapIter++)
for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
{
Info<< nl << "Scaling iteration " << snapIter << endl;
if (snapIter == nSnap)
if (snapIter == snapParams.nSnap())
{
Info<< "Displacement scaling for error reduction set to 0."
<< endl;
@ -901,6 +897,56 @@ void Foam::autoHexMeshDriver::preSmoothPatch
}
// Get (pp-local) indices of points that are both on zone and on patched surface
Foam::labelList Foam::autoHexMeshDriver::getZoneSurfacePoints
(
const indirectPrimitivePatch& pp,
const word& zoneName
) const
{
label zoneI = mesh_.faceZones().findZoneID(zoneName);
if (zoneI == -1)
{
FatalErrorIn
(
"autoHexMeshDriver::getZoneSurfacePoints"
"(const indirectPrimitivePatch&, const word&)"
) << "Cannot find zone " << zoneName
<< exit(FatalError);
}
const faceZone& fZone = mesh_.faceZones()[zoneI];
// Could use PrimitivePatch & localFaces to extract points but might just
// as well do it ourselves.
boolList pointOnZone(pp.nPoints(), false);
forAll(fZone, i)
{
const face& f = mesh_.faces()[fZone[i]];
forAll(f, fp)
{
label meshPointI = f[fp];
Map<label>::const_iterator iter =
pp.meshPointMap().find(meshPointI);
if (iter != pp.meshPointMap().end())
{
label pointI = iter();
pointOnZone[pointI] = true;
}
}
}
return findIndices(pointOnZone, true);
}
Foam::vectorField Foam::autoHexMeshDriver::calcNearestSurface
(
const scalarField& snapDist,
@ -925,33 +971,40 @@ Foam::vectorField Foam::autoHexMeshDriver::calcNearestSurface
// 1. All points to non-interface surfaces
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(localPoints, pointI)
{
pointIndexHit pHit;
label surfI = surfaces().findNearest
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces().findNearest
(
unzonedSurfaces,
localPoints[pointI],
sqr(4*snapDist[pointI]), // sqr of attract distance
pHit
localPoints,
sqr(4*snapDist), // sqr of attract distance
hitSurface,
hitInfo
);
if (surfI != -1)
forAll(hitInfo, pointI)
{
patchDisp[pointI] = pHit.hitPoint() - localPoints[pointI];
if (hitInfo[pointI].hit())
{
patchDisp[pointI] =
hitInfo[pointI].hitPoint()
- localPoints[pointI];
}
//else
//{
// WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
// << "For point:" << pointI
// << " coordinate:" << localPoints[pointI]
// << " did not find any surface within:"
// << 4*snapDist[pointI]
// << " meter." << endl;
//}
}
//else
//{
// WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
// << "For point:" << pointI
// << " coordinate:" << localPoints[pointI]
// << " did not find any surface within:" << 4*snapDist[pointI]
// << " meter." << endl;
//}
}
// 2. All points on zones to their respective surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -964,50 +1017,50 @@ Foam::vectorField Foam::autoHexMeshDriver::calcNearestSurface
const labelList surfacesToTest(1, zoneSurfI);
label zoneI = mesh_.faceZones().findZoneID(faceZoneNames[zoneSurfI]);
// Get indices of points both on faceZone and on pp.
labelList zonePointIndices
(
getZoneSurfacePoints
(
pp,
faceZoneNames[zoneSurfI]
)
);
const faceZone& fZone = mesh_.faceZones()[zoneI];
forAll(fZone, i)
pointField zonePoints(zonePointIndices.size());
forAll(zonePointIndices, i)
{
const face& f = mesh_.faces()[fZone[i]];
zonePoints[i] = localPoints[zonePointIndices[i]];
}
forAll(f, fp)
// Find nearest for points both on faceZone and pp.
List<pointIndexHit> hitInfo;
labelList hitSurface;
surfaces().findNearest
(
labelList(1, zoneSurfI),
zonePoints,
sqr(4*snapDist),
hitSurface,
hitInfo
);
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
label meshPointI = f[fp];
Map<label>::const_iterator iter =
pp.meshPointMap().find(meshPointI);
if (iter != pp.meshPointMap().end())
{
label pointI = iter();
pointIndexHit pHit;
label surfI = surfaces().findNearest
(
surfacesToTest,
localPoints[pointI],
sqr(4*snapDist[pointI]), // sqr of attract distance
pHit
);
if (surfI != -1)
{
patchDisp[pointI] =
pHit.hitPoint() - localPoints[pointI];
}
else
{
WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
<< "For point:" << pointI
<< " coordinate:" << localPoints[pointI]
<< " did not find any surface within:"
<< 4*snapDist[pointI]
<< " meter." << endl;
}
}
patchDisp[pointI] =
hitInfo[pointI].hitPoint()
- localPoints[pointI];
}
else
{
WarningIn("autoHexMeshDriver::calcNearestSurface(..)")
<< "For point:" << pointI
<< " coordinate:" << localPoints[pointI]
<< " did not find any surface within:"
<< 4*snapDist[pointI]
<< " meter." << endl;
}
}
}
@ -1078,7 +1131,7 @@ Foam::vectorField Foam::autoHexMeshDriver::calcNearestSurface
void Foam::autoHexMeshDriver::smoothDisplacement
(
const dictionary& snapDict,
const snapParameters& snapParams,
motionSmoother& meshMover
) const
{
@ -1087,9 +1140,6 @@ void Foam::autoHexMeshDriver::smoothDisplacement
Info<< "Smoothing displacement ..." << endl;
// Smoothing iterations
label nSmoothDisp(readLabel(snapDict.lookup("nSmoothDispl")));
// Set edge diffusivity as inverse of distance to patch
scalarField edgeGamma(1.0/(edgePatchDist(pMesh, pp) + SMALL));
//scalarField edgeGamma(mesh_.nEdges(), 1.0);
@ -1098,7 +1148,7 @@ void Foam::autoHexMeshDriver::smoothDisplacement
// Get displacement field
pointVectorField& disp = meshMover.displacement();
for (label iter = 0; iter < nSmoothDisp; iter++)
for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
{
if ((iter % 10) == 0)
{
@ -1140,16 +1190,12 @@ void Foam::autoHexMeshDriver::smoothDisplacement
void Foam::autoHexMeshDriver::scaleMesh
(
const dictionary& snapDict,
const snapParameters& snapParams,
const label nInitErrors,
const List<labelPair>& baffles,
motionSmoother& meshMover
)
{
// Snapping iterations
label nSnap(readLabel(snapDict.lookup("nSnap")));
// Relax displacement until correct mesh
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList checkFaces(identity(mesh_.nFaces()));
@ -1157,11 +1203,11 @@ void Foam::autoHexMeshDriver::scaleMesh
scalar oldErrorReduction = -1;
Info<< "Moving mesh ..." << endl;
for (label iter = 0; iter < 2*nSnap; iter++)
for (label iter = 0; iter < 2*snapParams.nSnap(); iter++)
{
Info<< nl << "Iteration " << iter << endl;
if (iter == nSnap)
if (iter == snapParams.nSnap())
{
Info<< "Displacement scaling for error reduction set to 0." << endl;
oldErrorReduction = meshMover.setErrorReduction(0.0);

213
src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C Normal file
View File

@ -0,0 +1,213 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "layerParameters.H"
#include "polyBoundaryMesh.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::layerParameters::defaultConcaveAngle = 90;
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::layerParameters::layerParameters
(
const dictionary& dict,
const labelList& numLayers
)
:
numLayers_(numLayers),
expansionRatio_
(
numLayers_.size(),
readScalar(dict.lookup("expansionRatio"))
),
finalLayerRatio_
(
numLayers_.size(),
readScalar(dict.lookup("finalLayerRatio"))
),
minThickness_
(
numLayers_.size(),
readScalar(dict.lookup("minThickness"))
),
featureAngle_(readScalar(dict.lookup("featureAngle"))),
concaveAngle_
(
dict.found("concaveAngle")
? readScalar(dict.lookup("concaveAngle"))
: defaultConcaveAngle
),
nGrow_(readLabel(dict.lookup("nGrow"))),
nSmoothSurfaceNormals_
(
readLabel(dict.lookup("nSmoothSurfaceNormals"))
),
nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
maxFaceThicknessRatio_
(
readScalar(dict.lookup("maxFaceThicknessRatio"))
),
layerTerminationCos_
(
Foam::cos
(
0.5
* featureAngle_
* mathematicalConstant::pi/180.
)
),
maxThicknessToMedialRatio_
(
readScalar(dict.lookup("maxThicknessToMedialRatio"))
),
minMedianAxisAngleCos_
(
Foam::cos(readScalar(dict.lookup("minMedianAxisAngle")))
* mathematicalConstant::pi/180.
),
nBufferCellsNoExtrude_
(
readLabel(dict.lookup("nBufferCellsNoExtrude"))
),
nSnap_(readLabel(dict.lookup("nSnap")))
{}
// Construct from dictionary
Foam::layerParameters::layerParameters
(
const dictionary& dict,
const polyBoundaryMesh& boundaryMesh
)
:
numLayers_(boundaryMesh.size(), 0),
expansionRatio_
(
boundaryMesh.size(),
readScalar(dict.lookup("expansionRatio"))
),
finalLayerRatio_
(
boundaryMesh.size(),
readScalar(dict.lookup("finalLayerRatio"))
),
minThickness_
(
boundaryMesh.size(),
readScalar(dict.lookup("minThickness"))
),
featureAngle_(readScalar(dict.lookup("featureAngle"))),
concaveAngle_
(
dict.found("concaveAngle")
? readScalar(dict.lookup("concaveAngle"))
: defaultConcaveAngle
),
nGrow_(readLabel(dict.lookup("nGrow"))),
nSmoothSurfaceNormals_
(
readLabel(dict.lookup("nSmoothSurfaceNormals"))
),
nSmoothNormals_(readLabel(dict.lookup("nSmoothNormals"))),
nSmoothThickness_(readLabel(dict.lookup("nSmoothThickness"))),
maxFaceThicknessRatio_
(
readScalar(dict.lookup("maxFaceThicknessRatio"))
),
layerTerminationCos_
(
Foam::cos
(
0.5
* featureAngle_
* mathematicalConstant::pi/180.
)
),
maxThicknessToMedialRatio_
(
readScalar(dict.lookup("maxThicknessToMedialRatio"))
),
minMedianAxisAngleCos_
(
Foam::cos(readScalar(dict.lookup("minMedianAxisAngle")))
* mathematicalConstant::pi/180.
),
nBufferCellsNoExtrude_
(
readLabel(dict.lookup("nBufferCellsNoExtrude"))
),
nSnap_(readLabel(dict.lookup("nRelaxIter")))
{
const dictionary& layersDict = dict.subDict("layers");
forAllConstIter(dictionary, layersDict, iter)
{
const word& key = iter().keyword();
if (layersDict.isDict(key))
{
label patchI = boundaryMesh.findPatchID(key);
if (patchI == -1)
{
FatalErrorIn
(
"layerParameters::layerParameters"
"(const dictionary&, const polyBoundaryMesh&)"
) << "Specified illegal patch " << key
<< " in layer dictionary." << endl
<< "Valid patch names are " << boundaryMesh.names()
<< exit(FatalError);
}
const dictionary& layerDict = layersDict.subDict(key);
numLayers_[patchI] =
readLabel(layerDict.lookup("nSurfaceLayers"));
//- Patch-wise layer parameters disabled for now. Just remove
// settings in initialiser list and uncomment below.
//expansionRatio_[patchI] =
// readScalar(layerDict.lookup("expansionRatio"));
//finalLayerRatio_[patchI] =
// readScalar(layerDict.lookup("finalLayerRatio"));
//minThickness_[patchI] =
// readScalar(layerDict.lookup("minThickness"));
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

246
src/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H Normal file
View File

@ -0,0 +1,246 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::layerParameters
Description
Simple container to keep together layer specific information.
SourceFiles
layerParameters.C
\*---------------------------------------------------------------------------*/
#ifndef layerParameters_H
#define layerParameters_H
#include "dictionary.H"
#include "scalarField.H"
#include "labelList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class polyBoundaryMesh;
/*---------------------------------------------------------------------------*\
Class layerParameters Declaration
\*---------------------------------------------------------------------------*/
class layerParameters
{
// Static data members
//- Default angle for faces to be convcave
static const scalar defaultConcaveAngle;
// Private data
// Per patch (not region!) information
//- How many layers to add.
labelList numLayers_;
scalarField expansionRatio_;
//- Wanted thickness of final added cell layer. If multiple layers
// is the thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
scalarField finalLayerRatio_;
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
scalarField minThickness_;
scalar featureAngle_;
scalar concaveAngle_;
label nGrow_;
label nSmoothSurfaceNormals_;
label nSmoothNormals_;
label nSmoothThickness_;
scalar maxFaceThicknessRatio_;
scalar layerTerminationCos_;
scalar maxThicknessToMedialRatio_;
scalar minMedianAxisAngleCos_;
label nBufferCellsNoExtrude_;
label nSnap_;
// Private Member Functions
//- Disallow default bitwise copy construct
layerParameters(const layerParameters&);
//- Disallow default bitwise assignment
void operator=(const layerParameters&);
public:
// Constructors
//- Construct from dictionary - old syntax
layerParameters(const dictionary& dict, const labelList& numLayers);
//- Construct from dictionary - new syntax
layerParameters(const dictionary& dict, const polyBoundaryMesh&);
// Member Functions
// Access
// Per patch information
//- How many layers to add.
const labelList& numLayers() const
{
return numLayers_;
}
// Expansion factor for layer mesh
const scalarField& expansionRatio() const
{
return expansionRatio_;
}
//- Wanted thickness of final added cell layer. If multiple
// layers
// is the thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
const scalarField& finalLayerRatio() const
{
return finalLayerRatio_;
}
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
const scalarField& minThickness() const
{
return minThickness_;
}
scalar featureAngle() const
{
return featureAngle_;
}
scalar concaveAngle() const
{
return concaveAngle_;
}
//- If points get not extruded do nGrow layers of connected faces
// that are not grown. Is used to not do layers at all close to
// features.
label nGrow() const
{
return nGrow_;
}
// Number of smoothing iterations of surface normals
label nSmoothSurfaceNormals() const
{
return nSmoothSurfaceNormals_;
}
// Number of smoothing iterations of interior mesh movement
// direction
label nSmoothNormals() const
{
return nSmoothNormals_;
}
// Smooth layer thickness over surface patches
label nSmoothThickness() const
{
return nSmoothThickness_;
}
// Stop layer growth on highly warped cells
scalar maxFaceThicknessRatio() const
{
return maxFaceThicknessRatio_;
}
scalar layerTerminationCos() const
{
return layerTerminationCos_;
}
// Reduce layer growth where ratio thickness to medial
// distance is large
scalar maxThicknessToMedialRatio() const
{
return maxThicknessToMedialRatio_;
}
// Angle used to pick up medial axis points
scalar minMedianAxisAngleCos() const
{
return minMedianAxisAngleCos_;
}
// Create buffer region for new layer terminations
label nBufferCellsNoExtrude() const
{
return nBufferCellsNoExtrude_;
}
label nSnap() const
{
return nSnap_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

123
src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C Normal file
View File

@ -0,0 +1,123 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "refinementParameters.H"
#include "mathematicalConstants.H"
#include "polyMesh.H"
#include "globalIndex.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::refinementParameters::refinementParameters
(
const dictionary& dict,
const label dummy
)
:
maxGlobalCells_(readLabel(dict.lookup("cellLimit"))),
maxLocalCells_(readLabel(dict.lookup("procCellLimit"))),
minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
curvature_(readScalar(dict.lookup("curvature"))),
nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
keepPoints_(dict.lookup("keepPoints"))
{}
Foam::refinementParameters::refinementParameters(const dictionary& dict)
:
maxGlobalCells_(readLabel(dict.lookup("maxGlobalCells"))),
maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
keepPoints_(pointField(1, dict.lookup("locationInMesh")))
{
scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle")));
if (featAngle < 0 || featAngle > 180)
{
curvature_ = -GREAT;
}
else
{
curvature_ = Foam::cos(featAngle*mathematicalConstant::pi/180.0);
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh)
const
{
// Global calculation engine
globalIndex globalCells(mesh.nCells());
// Cell label per point
labelList cellLabels(keepPoints_.size());
forAll(keepPoints_, i)
{
const point& keepPoint = keepPoints_[i];
label localCellI = mesh.findCell(keepPoint);
label globalCellI = -1;
if (localCellI != -1)
{
Pout<< "Found point " << keepPoint << " in cell " << localCellI
<< " on processor " << Pstream::myProcNo() << endl;
globalCellI = globalCells.toGlobal(localCellI);
}
reduce(globalCellI, maxOp<label>());
if (globalCellI == -1)
{
FatalErrorIn
(
"refinementParameters::findCells(const polyMesh&) const"
) << "Point " << keepPoint
<< " is not inside the mesh or on a face or edge." << nl
<< "Bounding box of the mesh:" << mesh.bounds()
<< exit(FatalError);
}
if (globalCells.isLocal(globalCellI))
{
cellLabels[i] = localCellI;
}
else
{
cellLabels[i] = -1;
}
}
return cellLabels;
}
// ************************************************************************* //

154
src/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::refinementParameters
Description
Simple container to keep together refinement specific information.
SourceFiles
refinementParameters.C
\*---------------------------------------------------------------------------*/
#ifndef refinementParameters_H
#define refinementParameters_H
#include "dictionary.H"
#include "pointField.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
class polyMesh;
/*---------------------------------------------------------------------------*\
Class refinementParameters Declaration
\*---------------------------------------------------------------------------*/
class refinementParameters
{
// Private data
//- Total number of cells
const label maxGlobalCells_;
//- Per processor max number of cells
const label maxLocalCells_;
//- When to stop refining
const label minRefineCells_;
//- Curvature
scalar curvature_;
//- Number of layers between different refinement levels
const label nBufferLayers_;
//- Areas to keep
const pointField keepPoints_;
// Private Member Functions
//- Disallow default bitwise copy construct
refinementParameters(const refinementParameters&);
//- Disallow default bitwise assignment
void operator=(const refinementParameters&);
public:
// Constructors
//- Construct from dictionary - old syntax
refinementParameters(const dictionary& dict, const label dummy);
//- Construct from dictionary - new syntax
refinementParameters(const dictionary& dict);
// Member Functions
// Access
//- Total number of cells
label maxGlobalCells() const
{
return maxGlobalCells_;
}
//- Per processor max number of cells
label maxLocalCells() const
{
return maxLocalCells_;
}
//- When to stop refining
label minRefineCells() const
{
return minRefineCells_;
}
//- Curvature
scalar curvature() const
{
return curvature_;
}
//- Number of layers between different refinement levels
label nBufferLayers() const
{
return nBufferLayers_;
}
//- Areas to keep
const pointField& keepPoints() const
{
return keepPoints_;
}
// Other
//- Checks that cells are in mesh. Returns cells they are in.
labelList findCells(const polyMesh&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

54
src/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.C Normal file
View File

@ -0,0 +1,54 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "snapParameters.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::snapParameters::snapParameters(const dictionary& dict, const label dummy)
:
nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
snapTol_(readScalar(dict.lookup("snapTol"))),
nSmoothDispl_(readLabel(dict.lookup("nSmoothDispl"))),
nSnap_(readLabel(dict.lookup("nSnap")))
{}
// Construct from dictionary
Foam::snapParameters::snapParameters(const dictionary& dict)
:
nSmoothPatch_(readLabel(dict.lookup("nSmoothPatch"))),
snapTol_(readScalar(dict.lookup("tolerance"))),
nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
nSnap_(readLabel(dict.lookup("nRelaxIter")))
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// ************************************************************************* //

129
src/autoMesh/autoHexMesh/autoHexMeshDriver/snapParameters/snapParameters.H Normal file
View File

@ -0,0 +1,129 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::snapParameters
Description
Simple container to keep together snap specific information.
SourceFiles
snapParameters.C
\*---------------------------------------------------------------------------*/
#ifndef snapParameters_H
#define snapParameters_H
#include "dictionary.H"
#include "scalar.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Class forward declarations
/*---------------------------------------------------------------------------*\
Class snapParameters Declaration
\*---------------------------------------------------------------------------*/
class snapParameters
{
// Private data
const label nSmoothPatch_;
const scalar snapTol_;
const label nSmoothDispl_;
const label nSnap_;
// Private Member Functions
//- Disallow default bitwise copy construct
snapParameters(const snapParameters&);
//- Disallow default bitwise assignment
void operator=(const snapParameters&);
public:
// Constructors
//- Construct from dictionary - old syntax
snapParameters(const dictionary& dict, const label dummy);
//- Construct from dictionary - new syntax
snapParameters(const dictionary& dict);
// Member Functions
// Access
//- Number of patch smoothing iterations before finding
// correspondence to surface
label nSmoothPatch() const
{
return nSmoothPatch_;
}
//- Relative distance for points to be attracted by surface
// feature point
// or edge. True distance is this factor times local
// maximum edge length.
scalar snapTol() const
{
return snapTol_;
}
//- Number of mesh displacement smoothing iterations.
label nSmoothDispl() const
{
return nSmoothDispl_;
}
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
label nSnap() const
{
return nSnap_;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

711
src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C
View File

@ -24,9 +24,7 @@ License
\*----------------------------------------------------------------------------*/
#include "Pstream.H"
#include "meshRefinement.H"
#include "volMesh.H"
#include "volFields.H"
#include "surfaceMesh.H"
@ -34,7 +32,6 @@ License
#include "Time.H"
#include "refinementHistory.H"
#include "refinementSurfaces.H"
#include "cellSet.H"
#include "decompositionMethod.H"
#include "regionSplit.H"
#include "fvMeshDistribute.H"
@ -50,8 +47,9 @@ License
#include "globalPointPatchFields.H"
#include "calculatedPointPatchFields.H"
#include "processorPointPatch.H"
#include "featureEdgeMesh.H"
#include "globalIndex.H"
#include "meshTools.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -117,90 +115,6 @@ void Foam::meshRefinement::calcNeighbourData
}
//void Foam::meshRefinement::calcCanonicalBoundaryData
//(
// labelList& leftLevel,
// pointField& leftCc,
// labelList& rightLevel,
// pointField& rightCc
//) const
//{
// const labelList& cellLevel = meshCutter_.cellLevel();
// const pointField& cellCentres = mesh_.cellCentres();
// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
//
// // Get boundary face centre and level. Coupled aware.
// calcNeighbourData(rightLevel, rightCc);
//
// // Get owner face centre and level. Canonical sort for coupled faces.
// forAll(patches, patchI)
// {
// const polyPatch& pp = patches[patchI];
//
// const unallocLabelList& faceCells = pp.faceCells();
// const vectorField::subField faceCentres = pp.faceCentres();
//
// if (isA<processorPolyPatch>(pp))
// {
// const processorPolyPatch& procPp =
// refCast<const processorPolyPatch>(patches[patchI]);
//
// label bFaceI = pp.start()-mesh_.nInternalFaces();
// forAll(faceCells, i)
// {
// leftLevel[bFaceI] = cellLevel[faceCells[i]];
// leftCc[bFaceI] = cellCentres[faceCells[i]];
// bFaceI++;
// }
//
// // Swap if on neighbour so both sides have same values
// if (!procPp.owner())
// {
// forAll(leftLevel, i)
// {
// Swap(leftLevel[i], rightLevel[i]);
// Swap(leftCc[i], rightCc[i]);
// }
// }
// }
// else if (isA<cyclicPolyPatch>(pp))
// {
// const processorPolyPatch& cycPp =
// refCast<const processorPolyPatch>(patches[patchI]);
//
// label bFaceI = pp.start()-mesh_.nInternalFaces();
// forAll(faceCells, i)
// {
// leftLevel[bFaceI] = cellLevel[faceCells[i]];
// leftCc[bFaceI] = cellCentres[faceCells[i]];
// bFaceI++;
// }
//
// // First half has data in normal order: owner in leftLevel,leftCc
// // and coupled data in rightLevel, rightCc. Second half has it
// // swapped:
// bFaceI = pp.start()+cycPp.size()/2-mesh_.nInternalFaces();
// for (label i = 0; i < cycPp.size()/2; i++)
// {
// Swap(leftLevel[bFaceI], rightLevel[bFaceI]);
// Swap(leftCc[bFaceI], rightCc[bFaceI]);
// bFaceI++;
// }
// }
// else
// {
// label bFaceI = pp.start()-mesh_.nInternalFaces();
// forAll(faceCells, i)
// {
// leftLevel[bFaceI] = cellLevel[faceCells[i]];
// leftCc[bFaceI] = cellCentres[faceCells[i]];
// bFaceI++;
// }
// }
// }
//}
// Find intersections of edges (given by their two endpoints) with surfaces.
// Returns first intersection if there are more than one.
void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
@ -220,33 +134,46 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
// Collect segments we want to test for
pointField start(changedFaces.size());
pointField end(changedFaces.size());
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
label own = mesh_.faceOwner()[faceI];
pointIndexHit surfaceHitInfo;
start[i] = cellCentres[own];
if (mesh_.isInternalFace(faceI))
{
surfaceIndex_[faceI] = surfaces_.findAnyIntersection
(
cellCentres[own],
cellCentres[mesh_.faceNeighbour()[faceI]],
surfaceHitInfo
);
end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
}
else
{
surfaceIndex_[faceI] = surfaces_.findAnyIntersection
(
cellCentres[own],
neiCc[faceI-mesh_.nInternalFaces()],
surfaceHitInfo
);
end[i] = neiCc[faceI-mesh_.nInternalFaces()];
}
}
// Do tests in one go
labelList surfaceHit;
{
labelList surfaceLevel;
surfaces_.findHigherIntersection
(
start,
end,
labelList(start.size(), -1), // accept any intersection
surfaceHit,
surfaceLevel
);
}
// Keep just surface hit
forAll(surfaceHit, i)
{
surfaceIndex_[changedFaces[i]] = surfaceHit[i];
}
// Make sure both sides have same information. This should be
// case in general since same vectors but just to make sure.
syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>(), false);
@ -305,67 +232,72 @@ void Foam::meshRefinement::checkData()
}
// Check meshRefinement
{
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
const point& ownCc = mesh_.cellCentres()[own];
const point& neiCc = mesh_.cellCentres()[nei];
pointIndexHit surfaceHitInfo;
label surfI = surfaces_.findAnyIntersection
(
ownCc,
neiCc,
surfaceHitInfo
);
if (surfI != surfaceIndex_[faceI])
{
WarningIn("meshRefinement::checkData()")
<< "Internal face:" << faceI
<< " fc:" << mesh_.faceCentres()[faceI]
<< " cached surfaceIndex_:" << surfaceIndex_[faceI]
<< " current:" << surfI
<< " ownCc:" << ownCc << " neiCc:" << neiCc
<< endl;
}
}
// Get boundary face centre and level. Coupled aware.
labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
for
(
label faceI = mesh_.nInternalFaces();
faceI < mesh_.nFaces();
faceI++
)
// Collect segments we want to test for
pointField start(mesh_.nFaces());
pointField end(mesh_.nFaces());
forAll(start, faceI)
{
label own = mesh_.faceOwner()[faceI];
const point& ownCc = mesh_.cellCentres()[own];
start[faceI] = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
pointIndexHit surfaceHitInfo;
label surfI = surfaces_.findAnyIntersection
(
ownCc,
neiCc[faceI-mesh_.nInternalFaces()],
surfaceHitInfo
);
if (surfI != surfaceIndex_[faceI])
if (mesh_.isInternalFace(faceI))
{
WarningIn("meshRefinement::checkData()")
<< "Boundary face:" << faceI
<< " patch:" << mesh_.boundaryMesh().whichPatch(faceI)
<< " fc:" << mesh_.faceCentres()[faceI]
<< " cached surfaceIndex_:" << surfaceIndex_[faceI]
<< " current:" << surfI
<< " ownCc:" << ownCc
<< " neiCc:" << neiCc[faceI-mesh_.nInternalFaces()]
<< endl;
end[faceI] = mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]];
}
else
{
end[faceI] = neiCc[faceI-mesh_.nInternalFaces()];
}
}
// Do tests in one go
labelList surfaceHit;
{
labelList surfaceLevel;
surfaces_.findHigherIntersection
(
start,
end,
labelList(start.size(), -1), // accept any intersection
surfaceHit,
surfaceLevel
);
}
// Check
forAll(surfaceHit, faceI)
{
if (surfaceHit[faceI] != surfaceIndex_[faceI])
{
if (mesh_.isInternalFace(faceI))
{
WarningIn("meshRefinement::checkData()")
<< "Internal face:" << faceI
<< " fc:" << mesh_.faceCentres()[faceI]
<< " cached surfaceIndex_:" << surfaceIndex_[faceI]
<< " current:" << surfaceHit[faceI]
<< " ownCc:"
<< mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
<< " neiCc:"
<< mesh_.cellCentres()[mesh_.faceNeighbour()[faceI]]
<< endl;
}
else
{
WarningIn("meshRefinement::checkData()")
<< "Boundary face:" << faceI
<< " fc:" << mesh_.faceCentres()[faceI]
<< " cached surfaceIndex_:" << surfaceIndex_[faceI]
<< " current:" << surfaceHit[faceI]
<< " ownCc:"
<< mesh_.cellCentres()[mesh_.faceOwner()[faceI]]
<< endl;
}
}
}
}
@ -492,283 +424,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
}
// Get which are inside any closed surface. Note that all closed surfaces
// will have already been oriented to have keepPoint outside.
Foam::labelList Foam::meshRefinement::getInsideCells
(
const word& postfix
) const
{
const pointField& points = mesh_.points();
// This test can be done in various ways. The ultimate is the intersection
// of any edge of the mesh with any surface and intersection of any edge
// of any surface of the mesh with any mesh face.
// For now:
// - none of the edges of a cell intersects any surface
// - any point of the cell is inside.
// Is point inside any closed surface?
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PackedList<1> pointIsInside;
label nPointInside = surfaces_.markInsidePoints
(
points,
pointIsInside
);
if (debug)
{
Pout<< "getInsideCells : Points internal to closed surfaces :"
<< " local:" << nPointInside
<< " global:" << returnReduce(nPointInside, sumOp<label>())
<< endl;
}
// Find cells with all points inside closed surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
PackedList<1> cellIsInside(mesh_.nCells(), 1u);
label nCellInside = mesh_.nCells();
// Unmark cells with any point outside
forAll(points, pointI)
{
if (pointIsInside.get(pointI) == 0u)
{
// Outside point. Unmark cells.
const labelList& pCells = mesh_.pointCells()[pointI];
forAll(pCells, i)
{
label cellI = pCells[i];
if (cellIsInside.get(cellI) == 1u)
{
cellIsInside.set(cellI, 0u);
nCellInside--;
}
}
}
}
label totNCellInside = nCellInside;
reduce(totNCellInside, sumOp<label>());
if (debug)
{
Pout<< "getInsideCells : Cells using inside points only :"
<< " local:" << nCellInside << " global:" << totNCellInside
<< endl;
}
// Filter any cells with cc-cc edge intersection
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const cellList& cells = mesh_.cells();
cellSet insidePoints(mesh_, "insidePoints_"+postfix, mesh_.nCells()/100);
cellSet insidePointsButCut
(
mesh_,
"insidePointsButCut_"+postfix,
mesh_.nCells()/100
);
forAll(cells, cellI)
{
if (cellIsInside.get(cellI) == 1u)
{
insidePoints.insert(cellI);
const cell& cFaces = cells[cellI];
forAll(cFaces, i)
{
if (surfaceIndex_[cFaces[i]] != -1)
{
cellIsInside.set(cellI, 0u);
nCellInside--;
insidePointsButCut.insert(cellI);
break;
}
}
}
}
if (debug)
{
Pout<< "getInsideCells : cells with all points inside : "
<< insidePoints.objectPath() << endl;
insidePoints.write();
Pout<< "getInsideCells : cells with all points inside"
<< " but intersecting (cc) surface : "
<< insidePointsButCut.objectPath() << endl;
insidePointsButCut.write();
}
totNCellInside = nCellInside;
reduce(totNCellInside, sumOp<label>());
if (debug)
{
Pout<< "getInsideCells : After filtering cc-cc intersections :"
<< " local:" << nCellInside << " global:" << totNCellInside
<< endl;
}
// Filter cells with any normal edge intersection
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cellSet insidePointsButEdge
(
mesh_,
"insidePointsButEdge_"+postfix,
mesh_.nCells()/100
);
const labelListList& edgeCells = mesh_.edgeCells();
forAll(edgeCells, edgeI)
{
const labelList& eCells = edgeCells[edgeI];
// Does edge use any inside cells? Prefilter since edge intersection
// test expensive.
bool edgeHasInsideCells = false;
forAll(eCells, i)
{
if (cellIsInside.get(eCells[i]) == 1u)
{
edgeHasInsideCells = true;
break;
}
}
if (edgeHasInsideCells)
{
const edge& e = mesh_.edges()[edgeI];
// Find any pierces of surface by mesh edge.
pointIndexHit hit;
label surfI = surfaces().findHigherIntersection
(
points[e[0]],
points[e[1]],
-1, // minimum level.
hit
);
if (surfI != -1)
{
// This edge intersects some surface. Unmark all cells
// using this edge.
forAll(eCells, i)
{
if (cellIsInside.get(eCells[i]) == 1u)
{
cellIsInside.set(eCells[i], 0u);
nCellInside--;
insidePointsButEdge.insert(eCells[i]);
}
}
}
}
}
if (debug)
{
Pout<< "getInsideCells : cells with all points inside"
<< " but intersecting (edges) surface : "
<< insidePointsButEdge.objectPath() << endl;
insidePointsButEdge.write();
}
totNCellInside = nCellInside;
reduce(totNCellInside, sumOp<label>());
if (debug)
{
Pout<< "getInsideCells : After filtering edge interesections :"
<< " local:" << nCellInside << " global:" << totNCellInside
<< endl;
}
// Collect cells
// ~~~~~~~~~~~~~
DynamicList<label> cellsToRemove(nCellInside);
forAll(cellIsInside, cellI)
{
if (cellIsInside.get(cellI) == 1u)
{
cellsToRemove.append(cellI);
}
}
return cellsToRemove.shrink();
}
// Do all to remove inside cells
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeInsideCells
(
const string& msg,
const label exposedPatchI
)
{
labelList insideCells
(
getInsideCells
(
string::validate<word>(msg)
)
);
Info<< "Removing inside cells : "
<< returnReduce(insideCells.size(), sumOp<label>())
<< '.' << endl;
removeCells cellRemover(mesh_);
labelList exposedFaces(cellRemover.getExposedFaces(insideCells));
// Debug check
forAll(exposedFaces, i)
{
label faceI = exposedFaces[i];
if (surfaceIndex_[faceI] != -1)
{
FatalErrorIn
(
"removeInsideCells"
"(const label, const string&, const label)"
) << "Face:" << faceI
<< " fc:" << mesh_.faceCentres()[faceI]
<< " is actually cutting the surface and should never have been"
<< " included in the cells to remove."
<< abort(FatalError);
}
}
return doRemoveCells
(
insideCells,
exposedFaces,
labelList(exposedFaces.size(), exposedPatchI),
cellRemover
);
}
// Determine for multi-processor regions the lowest numbered cell on the lowest
// numbered processor.
void Foam::meshRefinement::getRegionMaster
@ -835,12 +490,14 @@ Foam::meshRefinement::meshRefinement
(
fvMesh& mesh,
const scalar tol,
const refinementSurfaces& surfaces
const refinementSurfaces& surfaces,
const shellSurfaces& shells
)
:
mesh_(mesh),
tol_(tol),
surfaces_(surfaces),
shells_(shells),
meshCutter_
(
mesh,
@ -926,109 +583,6 @@ Foam::label Foam::meshRefinement::countHits() const
}
//// Determine distribution to keep baffles together
//Foam::labelList Foam::meshRefinement::decomposePreserveBaffles
//(
// decompositionMethod& decomposer
//) const
//{
// // Find duplicate faces
// labelList duplicateFace
// (
// localPointRegion::findDuplicateFaces
// (
// mesh_,
// identity(mesh_.nFaces()-mesh_.nInternalFaces())
// + mesh_.nInternalFaces()
// )
// );
//
//
// // Agglomerate cells on both sides of duplicate face pair
//
// const pointField& cellCentres = mesh_.cellCentres();
// const labelList& faceOwner = mesh_.faceOwner();
//
// labelList toAgglom(mesh_.nCells(), -1);
// pointField agglomCellCentres(mesh_.nCells());
// labelList nAgglomCellCentres(mesh_.nCells(), 0);
//
// label agglomI = 0;
//
// // Do all baffle cells
// forAll(duplicateFace, i)
// {
// if (duplicateFace[i] != -1)
// {
// label own = faceOwner[i+mesh_.nInternalFaces()];
// label otherOwn =
// faceOwner[duplicateFace[i]+mesh_.nInternalFaces()];
//
// if (toAgglom[own] == -1 && toAgglom[otherOwn] == -1)
// {
// // Allocate new agglomeration
// agglomCellCentres[agglomI] =
// cellCentres[own]
// + cellCentres[otherOwn];
// nAgglomCellCentres[agglomI] = 2;
// toAgglom[own] = agglomI;
// toAgglom[otherOwn] = agglomI;
// agglomI++;
// }
// else if (toAgglom[own] != -1)
// {
// // Owner already part of agglomeration. Add otherOwn to it.
// label destAgglom = toAgglom[own];
// agglomCellCentres[destAgglom] += cellCentres[otherOwn];
// nAgglomCellCentres[destAgglom]++;
// toAgglom[otherOwn] = destAgglom;
// }
// else if (toAgglom[otherOwn] != -1)
// {
// // otherOwn already part of agglomeration. Add own to it.
// label destAgglom = toAgglom[otherOwn];
// agglomCellCentres[destAgglom] += cellCentres[own];
// nAgglomCellCentres[destAgglom]++;
// toAgglom[own] = destAgglom;
// }
// }
// }
//
// // Do all other cells
// forAll(toAgglom, cellI)
// {
// if (toAgglom[cellI] == -1)
// {
// agglomCellCentres[agglomI] = cellCentres[cellI];
// nAgglomCellCentres[agglomI] = 1;
// toAgglom[cellI] = agglomI;
// agglomI++;
// }
// }
//
// // Average
// agglomCellCentres.setSize(agglomI);
//
// forAll(agglomCellCentres, i)
// {
// agglomCellCentres[i] /= nAgglomCellCentres[i];
// }
//
// // Decompose based on agglomerated cell centres
// labelList agglomDistribution(decomposer.decompose(agglomCellCentres));
//
// // Rework back into decomposition for original mesh_
// labelList distribution(mesh_.nCells());
//
// forAll(toAgglom, cellI)
// {
// distribution[cellI] = agglomDistribution[toAgglom[cellI]];
// }
//
// return distribution;
//}
// Determine distribution to move connected regions onto one processor.
Foam::labelList Foam::meshRefinement::decomposeCombineRegions
(
@ -1213,9 +767,8 @@ Foam::labelList Foam::meshRefinement::intersectedPoints
forAll(f, fp)
{
if (isBoundaryPoint.get(f[fp]) == 0u)
if (isBoundaryPoint.set(f[fp], 1u))
{
isBoundaryPoint.set(f[fp], 1u);
nBoundaryPoints++;
}
}
@ -1239,9 +792,7 @@ Foam::labelList Foam::meshRefinement::intersectedPoints
//
// forAll(f, fp)
// {
// if (isBoundaryPoint.get(f[fp]) == 0u)
// {
// isBoundaryPoint.set(f[fp], 1u);
// if (isBoundaryPoint.set(f[fp], 1u))
// nBoundaryPoints++;
// }
// }
@ -1910,53 +1461,56 @@ void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
<< " Writing cellcentre-cellcentre intersections to file "
<< str.name() << endl;
// Internal faces
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
if (surfaceIndex_[faceI] != -1)
{
const point& ownCc = cellCentres[mesh_.faceOwner()[faceI]];
const point& neiCc = cellCentres[mesh_.faceNeighbour()[faceI]];
// Re-intersect to get position
pointIndexHit surfaceHitInfo;
surfaces_.findAnyIntersection(ownCc, neiCc, surfaceHitInfo);
meshTools::writeOBJ(str, ownCc);
vertI++;
meshTools::writeOBJ(str, surfaceHitInfo.hitPoint());
vertI++;
meshTools::writeOBJ(str, neiCc);
vertI++;
str << "l " << vertI-2 << ' ' << vertI-1 << nl
<< "l " << vertI-1 << ' ' << vertI << nl;
}
}
// Boundary faces
// Redo all intersections
// ~~~~~~~~~~~~~~~~~~~~~~
// Get boundary face centre and level. Coupled aware.
labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
forAll(neiCc, i)
labelList intersectionFaces(intersectedFaces());
// Collect segments we want to test for
pointField start(intersectionFaces.size());
pointField end(intersectionFaces.size());
forAll(intersectionFaces, i)
{
label faceI = i + mesh_.nInternalFaces();
label faceI = intersectionFaces[i];
start[i] = cellCentres[mesh_.faceOwner()[faceI]];
if (surfaceIndex_[faceI] != -1)
if (mesh_.isInternalFace(faceI))
{
const point& ownCc = cellCentres[mesh_.faceOwner()[faceI]];
end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
}
else
{
end[i] = neiCc[faceI-mesh_.nInternalFaces()];
}
}
// Re-intersect to get position
pointIndexHit surfaceHitInfo;
surfaces_.findAnyIntersection(ownCc, neiCc[i], surfaceHitInfo);
// Do tests in one go
labelList surfaceHit;
List<pointIndexHit> surfaceHitInfo;
surfaces_.findAnyIntersection
(
start,
end,
surfaceHit,
surfaceHitInfo
);
meshTools::writeOBJ(str, ownCc);
forAll(intersectionFaces, i)
{
if (surfaceHit[i] != -1)
{
meshTools::writeOBJ(str, start[i]);
vertI++;
meshTools::writeOBJ(str, surfaceHitInfo.hitPoint());
meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
vertI++;
meshTools::writeOBJ(str, neiCc[i]);
meshTools::writeOBJ(str, end[i]);
vertI++;
str << "l " << vertI-2 << ' ' << vertI-1 << nl
<< "l " << vertI-1 << ' ' << vertI << nl;
@ -1996,13 +1550,4 @@ void Foam::meshRefinement::write
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

40
src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
View File

@ -30,7 +30,7 @@ Description
(static) surfaces.
Maintains
- per face any intersections of this edge with any of the surfaces
- per face any intersections of the cc-cc segment with any of the surfaces
SourceFiles
meshRefinement.C
@ -62,6 +62,7 @@ class fvMesh;
class mapDistributePolyMesh;
class decompositionMethod;
class refinementSurfaces;
class shellSurfaces;
class removeCells;
class featureEdgeMesh;
class fvMeshDistribute;
@ -107,8 +108,12 @@ private:
//- tolerance used for sorting coordinates (used in 'less' routine)
const scalar tol_;
//- All surface-intersection interaction
const refinementSurfaces& surfaces_;
//- All shell-refinement interaction
const shellSurfaces& shells_;
//- refinement engine
hexRef8 meshCutter_;
@ -229,14 +234,18 @@ private:
//- Mark cells for refinement-shells based refinement.
label markInternalRefinement
(
const PtrList<searchableSurface>&,
const labelList& shellLevels,
const boolList& shellRefineInside,
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const;
//- Collect faces that are intersected and whose neighbours aren't
// yet marked for refinement.
labelList getRefineCandidateFaces
(
const labelList& refineCell
) const;
//- Mark cells for surface intersection based refinement.
label markSurfaceRefinement
(
@ -252,17 +261,13 @@ private:
// regions.
bool checkCurvature
(
const labelList& globalToPatch,
const scalar curvature,
const bool markDifferingRegions,
const label nAllowRefine,
const label newSurfI,
const label newTriI,
const label surfaceLevel,
const vector& surfaceNormal,
const label cellI,
label& maxCellSurfI,
label& maxCellTriI,
label& cellMaxLevel,
vector& cellMaxNormal,
labelList& refineCell,
label& nRefine
) const;
@ -273,9 +278,7 @@ private:
// (markDifferingRegions)
label markSurfaceCurvatureRefinement
(
const labelList& globalToPatch,
const scalar curvature,
const bool markDifferingRegions,
const label nAllowRefine,
const labelList& neiLevel,
const pointField& neiCc,
@ -381,7 +384,8 @@ public:
(
fvMesh& mesh,
const scalar tol,
const refinementSurfaces&
const refinementSurfaces&,
const shellSurfaces&
);
@ -482,22 +486,18 @@ public:
const labelList& adaptPatchIDs
);
// Refinement
//- Calculate list of cells to refine.
labelList refineCandidates
(
const point& keepPoint,
const labelList& globalToPatch,
const scalar curvature,
const PtrList<featureEdgeMesh>& featureMeshes,
const labelList& featureLevels,
const PtrList<searchableSurface>&,
const labelList& shellLevels,
const boolList& shellRefineInside,
const bool featureRefinement,
const bool internalRefinement,
const bool surfaceRefinement,

434
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
View File

@ -249,7 +249,7 @@ void Foam::meshRefinement::getBafflePatches
if (debug)
{
Pout<< "getBafflePatches : Not baffling surface "
<< surfaces_[surfI].searchableSurface::name() << endl;
<< surfaces_.names()[surfI] << endl;
}
}
else
@ -264,69 +264,95 @@ void Foam::meshRefinement::getBafflePatches
neiPatch.setSize(mesh_.nFaces());
neiPatch = -1;
forAll(surfaceIndex_, faceI)
// Collect candidate faces
// ~~~~~~~~~~~~~~~~~~~~~~~
labelList testFaces(intersectedFaces());
// Collect segments
// ~~~~~~~~~~~~~~~~
pointField start(testFaces.size());
pointField end(testFaces.size());
forAll(testFaces, i)
{
if (surfaceIndex_[faceI] != -1)
label faceI = testFaces[i];
label own = mesh_.faceOwner()[faceI];
if (mesh_.isInternalFace(faceI))
{
// Edge between owner and neighbour of face pierces a surface.
// Do closer inspection to find nearest intersection on both sides.
start[i] = cellCentres[own];
end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
}
else
{
start[i] = cellCentres[own];
end[i] = neiCc[faceI-mesh_.nInternalFaces()];
}
}
label own = mesh_.faceOwner()[faceI];
// Cc of neighbouring cell
point end
(
mesh_.isInternalFace(faceI)
? cellCentres[mesh_.faceNeighbour()[faceI]]
: neiCc[faceI-mesh_.nInternalFaces()]
);
// Do test for intersections
// ~~~~~~~~~~~~~~~~~~~~~~~~~
labelList surface1;
List<pointIndexHit> hit1;
labelList region1;
labelList surface2;
List<pointIndexHit> hit2;
labelList region2;
surfaces_.findNearestIntersection
(
surfacesToBaffle,
start,
end,
label surface1, surface2;
pointIndexHit hit1, hit2;
surfaces_.findNearestIntersection
(
surfacesToBaffle,
cellCentres[own],
end,
surface1,
hit1,
region1,
surface1,
hit1,
surface2,
hit2
);
surface2,
hit2,
region2
);
if (hit1.hit() && hit2.hit())
forAll(testFaces, i)
{
label faceI = testFaces[i];
if (hit1[i].hit() && hit2[i].hit())
{
if (str.valid())
{
if (str.valid())
{
meshTools::writeOBJ(str(), cellCentres[own]);
vertI++;
meshTools::writeOBJ(str(), hit1.rawPoint());
vertI++;
meshTools::writeOBJ(str(), hit2.rawPoint());
vertI++;
meshTools::writeOBJ(str(), end);
vertI++;
str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
meshTools::writeOBJ(str(), start[i]);
vertI++;
meshTools::writeOBJ(str(), hit1[i].rawPoint());
vertI++;
meshTools::writeOBJ(str(), hit2[i].rawPoint());
vertI++;
meshTools::writeOBJ(str(), end[i]);
vertI++;
str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
// Pick up the patches
ownPatch[faceI] = globalToPatch
[
surfaces_.globalRegion(surface1[i], region1[i])
];
neiPatch[faceI] = globalToPatch
[
surfaces_.globalRegion(surface2[i], region2[i])
];
// Pick up the patches
ownPatch[faceI] = globalToPatch
[
surfaces_.triangleRegion(surface1, hit1.index())
];
neiPatch[faceI] = globalToPatch
[
surfaces_.triangleRegion(surface2, hit2.index())
];
if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
{
FatalErrorIn("getBafflePatches()") << abort(FatalError);
}
if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
{
FatalErrorIn("getBafflePatches(..)")
<< "problem." << abort(FatalError);
}
}
}
@ -1203,15 +1229,20 @@ void Foam::meshRefinement::findCellZoneGeometric
const labelList& faceOwner = mesh_.faceOwner();
const labelList& faceNeighbour = mesh_.faceNeighbour();
forAll(cellToZone, cellI)
// Check if cell centre is inside
labelList insideSurfaces;
surfaces_.findInside
(
closedNamedSurfaces,
cellCentres,
insideSurfaces
);
forAll(insideSurfaces, cellI)
{
if (cellToZone[cellI] == -2)
{
label surfI = surfaces_.insideZone
(
closedNamedSurfaces,
cellCentres[cellI]
);
label surfI = insideSurfaces[cellI];
if (surfI != -1)
{
@ -1223,6 +1254,32 @@ void Foam::meshRefinement::findCellZoneGeometric
// Some cells with cell centres close to surface might have
// had been put into wrong surface. Recheck with perturbed cell centre.
// 1. Collect points
// Count points to test.
label nCandidates = 0;
forAll(namedSurfaceIndex, faceI)
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
{
if (mesh_.isInternalFace(faceI))
{
nCandidates += 2;
}
else
{
nCandidates += 1;
}
}
}
// Collect points.
pointField candidatePoints(nCandidates);
nCandidates = 0;
forAll(namedSurfaceIndex, faceI)
{
label surfI = namedSurfaceIndex[faceI];
@ -1236,55 +1293,65 @@ void Foam::meshRefinement::findCellZoneGeometric
{
label nei = faceNeighbour[faceI];
const point& neiCc = cellCentres[nei];
// Perturbed cc
const vector d = 1E-4*(neiCc - ownCc);
candidatePoints[nCandidates++] = ownCc-d;
candidatePoints[nCandidates++] = neiCc+d;
}
else
{
const point& neiFc = mesh_.faceCentres()[faceI];
// Perturbed cc
const vector d = 1E-4*(neiFc - ownCc);
candidatePoints[nCandidates++] = ownCc-d;
}
}
}
// Test perturbed owner
// 2. Test points for inside
surfaces_.findInside
(
closedNamedSurfaces,
candidatePoints,
insideSurfaces
);
// 3. Update zone information
nCandidates = 0;
forAll(namedSurfaceIndex, faceI)
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
{
label own = faceOwner[faceI];
if (mesh_.isInternalFace(faceI))
{
label ownSurfI = insideSurfaces[nCandidates++];
if (ownSurfI != -1)
{
label ownSurfI = surfaces_.insideZone
(
closedNamedSurfaces,
ownCc-d
);
if (ownSurfI != -1)
{
cellToZone[own] = surfaceToCellZone[ownSurfI];
}
cellToZone[own] = surfaceToCellZone[ownSurfI];
}
// Test perturbed neighbour
label neiSurfI = insideSurfaces[nCandidates++];
if (neiSurfI != -1)
{
label neiSurfI = surfaces_.insideZone
(
closedNamedSurfaces,
neiCc+d
);
label nei = faceNeighbour[faceI];
if (neiSurfI != -1)
{
cellToZone[nei] = surfaceToCellZone[neiSurfI];
}
cellToZone[nei] = surfaceToCellZone[neiSurfI];
}
}
else
{
const point& neiCc = mesh_.faceCentres()[faceI];
const vector d = 1E-4*(neiCc - ownCc);
// Test perturbed owner
label ownSurfI = insideSurfaces[nCandidates++];
if (ownSurfI != -1)
{
label ownSurfI = surfaces_.insideZone
(
closedNamedSurfaces,
ownCc-d
);
if (ownSurfI != -1)
{
cellToZone[own] = surfaceToCellZone[ownSurfI];
}
cellToZone[own] = surfaceToCellZone[ownSurfI];
}
}
}
@ -2291,7 +2358,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
isNamedSurface[surfI] = true;
Info<< "Surface : " << surfaces_[surfI].searchableSurface::name() << nl
Info<< "Surface : " << surfaces_.names()[surfI] << nl
<< " faceZone : " << faceZoneNames[surfI] << nl
<< " cellZone : " << cellZoneNames[surfI] << endl;
}
@ -2329,7 +2396,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (debug)
{
Pout<< "Faces on " << surfaces_[surfI].searchableSurface::name()
Pout<< "Faces on " << surfaces_.names()[surfI]
<< " will go into faceZone " << zoneI << endl;
}
surfaceToFaceZone[surfI] = zoneI;
@ -2387,8 +2454,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (debug)
{
Pout<< "Cells inside "
<< surfaces_[surfI].searchableSurface::name()
Pout<< "Cells inside " << surfaces_.names()[surfI]
<< " will go into cellZone " << zoneI << endl;
}
surfaceToCellZone[surfI] = zoneI;
@ -2444,78 +2510,83 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Note: for all internal faces? internal + coupled?
// Since zonify is run after baffling the surfaceIndex_ on baffles is
// not synchronised across both baffle faces. Fortunately we don't
// do zonify baffle faces anyway.
forAll(surfaceIndex_, faceI)
// do zonify baffle faces anyway (they are normal boundary faces).
// Collect candidate faces
// ~~~~~~~~~~~~~~~~~~~~~~~
labelList testFaces(intersectedFaces());
// Collect segments
// ~~~~~~~~~~~~~~~~
pointField start(testFaces.size());
pointField end(testFaces.size());
forAll(testFaces, i)
{
label surfI = surfaceIndex_[faceI];
label faceI = testFaces[i];
if
(
surfI != -1
&& (
mesh_.isInternalFace(faceI)
|| patches[patches.whichPatch(faceI)].coupled()
)
)
label own = mesh_.faceOwner()[faceI];
if (mesh_.isInternalFace(faceI))
{
if (isNamedSurface[surfI])
{
namedSurfaceIndex[faceI] = surfI;
nSurfFaces[surfI]++;
}
else
{
// Test more accurately for whether intersection is at
// zoneable surface
label surface1, surface2;
pointIndexHit hit1, hit2;
if (mesh_.isInternalFace(faceI))
{
surfaces_.findNearestIntersection
(
namedSurfaces,
cellCentres[faceOwner[faceI]],
cellCentres[faceNeighbour[faceI]],
surface1,
hit1,
surface2,
hit2
);
}
else
{
surfaces_.findNearestIntersection
(
namedSurfaces,
cellCentres[faceOwner[faceI]],
neiCc[faceI-mesh_.nInternalFaces()],
surface1,
hit1,
surface2,
hit2
);
}
if (hit1.hit())
{
// If both hit should probably choose nearest. For
// later.
namedSurfaceIndex[faceI] = surface1;
nSurfFaces[surface1]++;
}
else if (hit2.hit())
{
namedSurfaceIndex[faceI] = surface2;
nSurfFaces[surface2]++;
}
}
start[i] = cellCentres[own];
end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
}
else
{
start[i] = cellCentres[own];
end[i] = neiCc[faceI-mesh_.nInternalFaces()];
}
}
// Do test for intersections
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Note that we intersect all intersected faces again. Could reuse
// the information already in surfaceIndex_.
labelList surface1;
labelList surface2;
{
List<pointIndexHit> hit1;
labelList region1;
List<pointIndexHit> hit2;
labelList region2;
surfaces_.findNearestIntersection
(
namedSurfaces,
start,
end,
surface1,
hit1,
region1,
surface2,
hit2,
region2
);
}
forAll(testFaces, i)
{
label faceI = testFaces[i];
if (surface1[i] != -1)
{
// If both hit should probably choose nearest. For later.
namedSurfaceIndex[faceI] = surface1[i];
nSurfFaces[surface1[i]]++;
}
else if (surface2[i] != -1)
{
namedSurfaceIndex[faceI] = surface2[i];
nSurfFaces[surface2[i]]++;
}
}
// surfaceIndex migh have different surfaces on both sides if
// there happen to be a (obviously thin) surface with different
// regions between the cell centres. If one is on a named surface
@ -2534,7 +2605,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
forAll(nSurfFaces, surfI)
{
Pout<< "Surface:"
<< surfaces_[surfI].searchableSurface::name()
<< surfaces_.names()[surfI]
<< " nZoneFaces:" << nSurfFaces[surfI] << nl;
}
Pout<< endl;
@ -2591,7 +2662,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
(
mesh_.faces()[faceI], // modified face
faceI, // label of face
faceOwner[faceI], // owner
faceOwner[faceI], // owner
-1, // neighbour
false, // face flip
patchI, // patch for face
@ -2609,20 +2680,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// Put the cells into the correct zone
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
labelList closedNamedSurfaces(namedSurfaces.size());
label nClosed = 0;
forAll(namedSurfaces, i)
{
label surfI = namedSurfaces[i];
if (surfaces_.closed()[surfI])
{
closedNamedSurfaces[nClosed++] = surfI;
}
}
closedNamedSurfaces.setSize(nClosed);
// Closed surfaces with cellZone specified.
labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
// Zone per cell:
// -2 : unset
@ -2704,13 +2763,4 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

15
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementMerge.C
View File

@ -30,15 +30,9 @@ License
#include "removePoints.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Merge faces that are in-line.
@ -244,13 +238,4 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeEdges
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

842
src/autoMesh/autoHexMesh/meshRefinement/meshRefinementRefine.C
View File

File diff suppressed because it is too large Load Diff

989
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
View File

File diff suppressed because it is too large Load Diff

176
src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
View File

@ -38,34 +38,38 @@ SourceFiles
#ifndef refinementSurfaces_H
#define refinementSurfaces_H
#include "triSurfaceMeshes.H"
#include "PackedList.H"
#include "triSurfaceGeoMesh.H"
#include "triSurfaceFields.H"
#include "vectorList.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class searchableSurface;
class searchableSurfaces;
class shellSurfaces;
class triSurfaceMesh;
/*---------------------------------------------------------------------------*\
Class refinementSurfaces Declaration
\*---------------------------------------------------------------------------*/
class refinementSurfaces
:
public triSurfaceMeshes
{
// Private data
//- Reference to all geometry.
const searchableSurfaces& allGeometry_;
//- Indices of surfaces that are refinement ones
labelList surfaces_;
//- Surface name (word)
wordList names_;
//- Per surface whether is closed
boolList closed_;
//- Per 'interface' surface : name of faceZone to put faces into
wordList faceZoneNames_;
@ -86,24 +90,12 @@ class refinementSurfaces
//- From global region number to refinement level
labelList maxLevel_;
//- From global region number to added layers
labelList numLayers_;
//- From global region number to patch name
wordList patchName_;
//- From global region number to patch name
wordList patchType_;
//- Per surface refinement level adapted for shells.
PtrList<triSurfaceLabelField> minLevelFields_;
// Private Member Functions
static fileNameList extractFileNames(const PtrList<dictionary>&);
//- Disallow default bitwise copy construct
refinementSurfaces(const refinementSurfaces&);
@ -115,30 +107,36 @@ public:
// Constructors
//- Construct from directory (io.instance()) and dictionaries
//- Construct from surfaces and dictionaries
refinementSurfaces
(
const IOobject& io,
const searchableSurfaces& allGeometry,
const PtrList<dictionary>&
);
//- Construct from surfaces and dictionary
refinementSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary&
);
// Member Functions
// Access
const labelList& surfaces() const
{
return surfaces_;
}
//- Names of surfaces
const wordList& names() const
{
return names_;
}
//- Per surface whether is closed
const boolList& closed() const
{
return closed_;
}
//- Per 'interface' surface : name of faceZone to put faces into
const wordList& faceZoneNames() const
{
@ -151,15 +149,11 @@ public:
return cellZoneNames_;
}
//- Per 'interface' surface : if closed: zonify cells inside surface
const boolList& zoneInside() const
{
return zoneInside_;
}
//- Get indices of named surfaces (surfaces with cellZoneName)
labelList getNamedSurfaces() const;
//- Get indices of closed named surfaces
labelList getClosedNamedSurfaces() const;
//- From local region number to global region number
const labelList& regionOffset() const
@ -179,24 +173,6 @@ public:
return maxLevel_;
}
//- From global region number to added layers
const labelList& numLayers() const
{
return numLayers_;
}
//- From global region number to patch name. Patchnames can be empty
const wordList& patchName() const
{
return patchName_;
}
//- From global region number to patch name
const wordList& patchType() const
{
return patchType_;
}
// Helper
@ -206,14 +182,6 @@ public:
return regionOffset_[surfI]+regionI;
}
//- From triangle on surface to global region
label triangleRegion(const label surfI, const label triI) const
{
const triSurface& s = operator[](surfI);
return globalRegion(surfI, s[triI].region());
}
//- Min level for surface and region on surface
label minLevel(const label surfI, const label regionI) const
{
@ -231,23 +199,12 @@ public:
return minLevel_.size();
}
//- Minlevel updated for refinement shells
const triSurfaceLabelField& minLevelField(const label surfI) const
{
return minLevelFields_[surfI];
}
//- Calculate minLevelFields
void setMinLevelFields
(
const PtrList<searchableSurface>& shells,
const labelList& shellLevels,
const boolList& shellRefineInside
const shellSurfaces& shells
);
//- Helper: is surface closed?
static bool isSurfaceClosed(const triSurface&);
//- Helper: orient (closed only) surfaces so keepPoint is outside.
static void orientSurface
(
@ -263,36 +220,71 @@ public:
//- Find intersection of edge. Return -1 or first surface
// with higher (than currentLevel) minlevel.
// Return surface number and hit info.
label findHigherIntersection
// Return surface number and level.
void findHigherIntersection
(
const point& start,
const point& end,
const label currentLevel, // current cell refinement level
pointIndexHit&
const pointField& start,
const pointField& end,
const labelList& currentLevel, // current cell refinement level
labelList& surfaces,
labelList& surfaceLevel
) const;
//- Find intersection with max level. Return -1 or the surface
// with the highest maxLevel above currentLevel
label findHighestIntersection
//- Find all intersections of edge. Unsorted order.
void findAllHigherIntersections
(
const point& start,
const point& end,
const label currentLevel, // current cell refinement level
pointIndexHit&
const pointField& start,
const pointField& end,
const labelList& currentLevel, // current cell refinement level
List<vectorList>& surfaceNormal,
labelListList& surfaceLevel
) const;
//- Detect if a point is 'inside' (depending on zoneInside flag) a
// zoneable surface. Returns -1 if not, returns first surface it
// is.
label insideZone(const labelList& surfaces, const point& pt) const;
//- Find intersection nearest to the endpoints. surface1,2 are
// not indices into surfacesToTest but refinement surface indices.
void findNearestIntersection
(
const labelList& surfacesToTest,
const pointField& start,
const pointField& end,
labelList& surface1,
List<pointIndexHit>& hit1,
labelList& region1,
labelList& surface2,
List<pointIndexHit>& hit2,
labelList& region2
) const;
//- Mark for all points whether it is inside any closed surface
// Return number of inside points.
label markInsidePoints(const pointField&, PackedList<1>& isInside)
const;
//- Used for debugging only: find intersection of edge.
void findAnyIntersection
(
const pointField& start,
const pointField& end,
labelList& surfaces,
List<pointIndexHit>&
) const;
//- Find nearest point on surfaces.
void findNearest
(
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& surfaces,
List<pointIndexHit>&
) const;
//- Detect if a point is 'inside' (closed) surfaces.
// Returns -1 if not, returns first surface it is.
void findInside
(
const labelList& surfacesToTest,
const pointField& pt,
labelList& insideSurfaces
) const;
};

465
src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C Normal file
View File

@ -0,0 +1,465 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "searchableSurface.H"
#include "shellSurfaces.H"
#include "boundBox.H"
#include "triSurfaceMesh.H"
#include "refinementSurfaces.H"
#include "searchableSurfaces.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
template<>
const char*
NamedEnum<shellSurfaces::refineMode, 3>::
names[] =
{
"inside",
"outside",
"distance"
};
const NamedEnum<shellSurfaces::refineMode, 3> shellSurfaces::refineModeNames_;
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::shellSurfaces::setAndCheckLevels
(
const scalar shellI,
const List<Tuple2<scalar, label> >& distLevels
)
{
if (modes_[shellI] != DISTANCE && distLevels.size() != 1)
{
FatalErrorIn
(
"shellSurfaces::shellSurfaces"
"(const searchableSurfaces&, const dictionary&)"
) << "For refinement mode "
<< refineModeNames_[modes_[shellI]]
<< " specify only one distance+level."
<< " (its distance gets discarded)"
<< exit(FatalError);
}
// Extract information into separate distance and level
distances_[shellI].setSize(distLevels.size());
levels_[shellI].setSize(distLevels.size());
forAll(distLevels, j)
{
distances_[shellI][j] = distLevels[j].first();
levels_[shellI][j] = distLevels[j].second();
// Check in incremental order
if (j > 0)
{
if
(
(distances_[shellI][j] <= distances_[shellI][j-1])
|| (levels_[shellI][j] > levels_[shellI][j-1])
)
{
FatalErrorIn
(
"shellSurfaces::shellSurfaces"
"(const searchableSurfaces&, const dictionary&)"
) << "For refinement mode "
<< refineModeNames_[modes_[shellI]]
<< " : Refinement should be specified in order"
<< " of increasing distance"
<< " (and decreasing refinement level)." << endl
<< "Distance:" << distances_[shellI][j]
<< " refinementLevel:" << levels_[shellI][j]
<< exit(FatalError);
}
}
}
const searchableSurface& shell = allGeometry_[shells_[shellI]];
if (modes_[shellI] == DISTANCE)
{
Info<< "Refinement level according to distance to "
<< shell.name() << endl;
forAll(levels_[shellI], j)
{
Info<< " level " << levels_[shellI][j]
<< " for all cells within " << distances_[shellI][j]
<< " meter." << endl;
}
}
else
{
if (!allGeometry_[shells_[shellI]].hasVolumeType())
{
FatalErrorIn
(
"shellSurfaces::shellSurfaces"
"(const searchableSurfaces&"
", const PtrList<dictionary>&)"
) << "Shell " << shell.name()
<< " does not support testing for "
<< refineModeNames_[modes_[shellI]] << endl
<< "Probably it is not closed."
<< exit(FatalError);
}
if (modes_[shellI] == INSIDE)
{
Info<< "Refinement level " << levels_[shellI][0]
<< " for all cells inside " << shell.name() << endl;
}
else
{
Info<< "Refinement level " << levels_[shellI][0]
<< " for all cells outside " << shell.name() << endl;
}
}
}
// Specifically orient triSurfaces using a calculated point outside.
// Done since quite often triSurfaces not of consistent orientation which
// is (currently) necessary for sideness calculation
void Foam::shellSurfaces::orient()
{
// Determine outside point.
boundBox overallBb
(
point(GREAT, GREAT, GREAT),
point(-GREAT, -GREAT, -GREAT)
);
bool hasSurface = false;
forAll(shells_, shellI)
{
const searchableSurface& s = allGeometry_[shells_[shellI]];
if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
{
const triSurfaceMesh& shell = refCast<const triSurfaceMesh>(s);
if (shell.triSurface::size() > 0)
{
const pointField& points = shell.points();
hasSurface = true;
boundBox shellBb(points[0], points[0]);
// Assume surface is compact!
for (label i = 0; i < points.size(); i++)
{
const point& pt = points[i];
shellBb.min() = min(shellBb.min(), pt);
shellBb.max() = max(shellBb.max(), pt);
}
overallBb.min() = min(overallBb.min(), shellBb.min());
overallBb.max() = max(overallBb.max(), shellBb.max());
}
}
}
if (hasSurface)
{
const point outsidePt(2*overallBb.max() - overallBb.min());
//Info<< "Using point " << outsidePt << " to orient shells" << endl;
forAll(shells_, shellI)
{
const searchableSurface& s = allGeometry_[shells_[shellI]];
if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
{
triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
(
refCast<const triSurfaceMesh>(s)
);
refinementSurfaces::orientSurface(outsidePt, shell);
}
}
}
}
// Find maximum level of a shell.
void Foam::shellSurfaces::findHigherLevel
(
const pointField& pt,
const label shellI,
labelList& maxLevel
) const
{
const labelList& levels = levels_[shellI];
if (modes_[shellI] == DISTANCE)
{
// Distance mode.
const scalarField& distances = distances_[shellI];
// Collect all those points that have a current maxLevel less than
// (any of) the shell. Also collect the furthest distance allowable
// to any shell with a higher level.
pointField candidates(pt.size());
labelList candidateMap(pt.size());
scalarField candidateDistSqr(pt.size());
label candidateI = 0;
forAll(maxLevel, pointI)
{
forAllReverse(levels, levelI)
{
if (levels[levelI] > maxLevel[pointI])
{
candidates[candidateI] = pt[pointI];
candidateMap[candidateI] = pointI;
candidateDistSqr[candidateI] = sqr(distances[levelI]);
candidateI++;
break;
}
}
}
candidates.setSize(candidateI);
candidateMap.setSize(candidateI);
candidateDistSqr.setSize(candidateI);
// Do the expensive nearest test only for the candidate points.
List<pointIndexHit> nearInfo;
allGeometry_[shells_[shellI]].findNearest
(
candidates,
candidateDistSqr,
nearInfo
);
// Update maxLevel
forAll(nearInfo, candidateI)
{
if (nearInfo[candidateI].hit())
{
// Check which level it actually is in.
label minDistI = findLower
(
distances,
mag(nearInfo[candidateI].hitPoint()-candidates[candidateI])
);
label pointI = candidateMap[candidateI];
// pt is inbetween shell[minDistI] and shell[minDistI+1]
maxLevel[pointI] = levels[minDistI+1];
}
}
}
else
{
// Inside/outside mode
// Collect all those points that have a current maxLevel less than the
// shell.
pointField candidates(pt.size());
labelList candidateMap(pt.size());
label candidateI = 0;
forAll(maxLevel, pointI)
{
if (levels[0] > maxLevel[pointI])
{
candidates[candidateI] = pt[pointI];
candidateMap[candidateI] = pointI;
candidateI++;
}
}
candidates.setSize(candidateI);
candidateMap.setSize(candidateI);
// Do the expensive nearest test only for the candidate points.
List<searchableSurface::volumeType> volType;
allGeometry_[shells_[shellI]].getVolumeType(candidates, volType);
forAll(volType, i)
{
label pointI = candidateMap[i];
if
(
(
modes_[shellI] == INSIDE
&& volType[i] == searchableSurface::INSIDE
)
|| (
modes_[shellI] == OUTSIDE
&& volType[i] == searchableSurface::OUTSIDE
)
)
{
maxLevel[pointI] = levels[0];
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::shellSurfaces::shellSurfaces
(
const searchableSurfaces& allGeometry,
const PtrList<dictionary>& shellDicts
)
:
allGeometry_(allGeometry)
{
shells_.setSize(shellDicts.size());
modes_.setSize(shellDicts.size());
distances_.setSize(shellDicts.size());
levels_.setSize(shellDicts.size());
forAll(shellDicts, shellI)
{
const dictionary& dict = shellDicts[shellI];
const word name = dict.lookup("name");
const word type = dict.lookup("type");
shells_[shellI] = allGeometry_.findSurfaceID(name);
if (shells_[shellI] == -1)
{
FatalErrorIn
(
"shellSurfaces::shellSurfaces"
"(const searchableSurfaces&, const PtrList<dictionary>&)"
) << "No surface called " << name << endl
<< "Valid surfaces are " << allGeometry_.names()
<< exit(FatalError);
}
modes_[shellI] = refineModeNames_.read(dict.lookup("refineMode"));
// Read pairs of distance+level
setAndCheckLevels(shellI, dict.lookup("levels"));
}
// Orient shell surfaces before any searching is done. Note that this
// only needs to be done for inside or outside. Orienting surfaces
// constructs lots of addressing which we want to avoid.
orient();
}
Foam::shellSurfaces::shellSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& shellsDict
)
:
allGeometry_(allGeometry)
{
shells_.setSize(shellsDict.size());
modes_.setSize(shellsDict.size());
distances_.setSize(shellsDict.size());
levels_.setSize(shellsDict.size());
label shellI = 0;
forAllConstIter(dictionary, shellsDict, iter)
{
shells_[shellI] = allGeometry_.findSurfaceID(iter().keyword());
if (shells_[shellI] == -1)
{
FatalErrorIn
(
"shellSurfaces::shellSurfaces"
"(const searchableSurfaces&, const dictionary>&"
) << "No surface called " << iter().keyword() << endl
<< "Valid surfaces are " << allGeometry_.names()
<< exit(FatalError);
}
const dictionary& dict = shellsDict.subDict(iter().keyword());
modes_[shellI] = refineModeNames_.read(dict.lookup("mode"));
// Read pairs of distance+level
setAndCheckLevels(shellI, dict.lookup("levels"));
shellI++;
}
// Orient shell surfaces before any searching is done. Note that this
// only needs to be done for inside or outside. Orienting surfaces
// constructs lots of addressing which we want to avoid.
orient();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Highest shell level
Foam::label Foam::shellSurfaces::maxLevel() const
{
label overallMax = 0;
forAll(levels_, shellI)
{
overallMax = max(overallMax, max(levels_[shellI]));
}
return overallMax;
}
void Foam::shellSurfaces::findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const
{
// Maximum level of any shell. Start off with level of point.
maxLevel = ptLevel;
forAll(shells_, shellI)
{
findHigherLevel(pt, shellI, maxLevel);
}
}
// ************************************************************************* //

182
src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.H Normal file
View File

@ -0,0 +1,182 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
shellSurfaces
Description
Encapsulates queries for volume refinement ('refine all cells within
shell').
SourceFiles
shellSurfaces.C
\*---------------------------------------------------------------------------*/
#ifndef shellSurfaces_H
#define shellSurfaces_H
#include "searchableSurface.H"
#include "Tuple2.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class searchableSurfaces;
/*---------------------------------------------------------------------------*\
Class shellSurfaces Declaration
\*---------------------------------------------------------------------------*/
class shellSurfaces
{
public:
// Public data types
//- Volume refinement controls
enum refineMode
{
INSIDE, // Refine all inside shell
OUTSIDE, // ,, outside
DISTANCE // Refine based on distance to shell
};
private:
// Private data
//- Reference to all geometry.
const searchableSurfaces& allGeometry_;
//- Indices of surfaces that are shells
labelList shells_;
//- Per shell whether to refine inside or outside
List<refineMode> modes_;
//- Per shell the list of ranges
List<scalarField> distances_;
//- Per shell per distance the refinement level
labelListList levels_;
// Private data
//- refineMode names
static const NamedEnum<refineMode, 3> refineModeNames_;
// Private Member Functions
//- Helper function for initialisation.
void setAndCheckLevels
(
const scalar shellI,
const List<Tuple2<scalar, label> >&
);
void orient();
void findHigherLevel
(
const pointField& pt,
const label shellI,
labelList& maxLevel
) const;
public:
// Constructors
//- Construct from components
shellSurfaces
(
const searchableSurfaces& allGeometry,
const labelList& shells,
const List<refineMode>& modes,
const List<scalarField>& distances,
const labelListList& levels
);
//- Construct from geometry and dictionaries
shellSurfaces
(
const searchableSurfaces& allGeometry,
const PtrList<dictionary>& shellDicts
);
//- Construct from geometry and dictionary
shellSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& shellsDict
);
// Member Functions
// Access
//const List<scalarField>& distances() const
//{
// return distances_;
//}
//
////- Per shell per distance the refinement level
//const labelListList& levels() const
//{
// return levels_;
//}
// Query
//- Highest shell level
label maxLevel() const;
//- Find shell level higher than ptLevel
void findHigherLevel
(
const pointField& pt,
const labelList& ptLevel,
labelList& maxLevel
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

2
src/autoMesh/autoHexMesh/trackedParticle/ExactParticle.C
View File

@ -256,7 +256,7 @@ Foam::Ostream& Foam::operator<<
const ExactParticle<ParticleType>& p
)
{
return operator<<(os, static_cast<Particle<ParticleType> >(p));
return operator<<(os, static_cast<const Particle<ParticleType>&>(p));
}

4
src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C
View File

@ -553,7 +553,7 @@ Foam::labelListList Foam::addPatchCellLayer::addedCells() const
void Foam::addPatchCellLayer::setRefinement
(
const scalar expansionRatio,
const scalarField& expansionRatio,
const indirectPrimitivePatch& pp,
const labelList& nFaceLayers,
const labelList& nPointLayers,
@ -885,7 +885,7 @@ void Foam::addPatchCellLayer::setRefinement
addedPoints_[patchPointI][i] = addedVertI;
disp *= expansionRatio;
disp *= expansionRatio[patchPointI];
}
}
}

4
src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H
View File

@ -305,7 +305,7 @@ public:
// (instead of e.g. from patch faces)
void setRefinement
(
const scalar expansionRatio,
const scalarField& expansionRatio,
const indirectPrimitivePatch& pp,
const labelList& nFaceLayers,
const labelList& nPointLayers,
@ -325,7 +325,7 @@ public:
{
setRefinement
(
1.0, // expansion ration
scalarField(pp.nPoints(), 1.0), // expansion ration
pp,
labelList(pp.size(), nLayers),
labelList(pp.nPoints(), nLayers),

18
src/meshTools/Make/files
View File

@ -53,6 +53,14 @@ indexedOctree/treeDataFace.C
indexedOctree/treeDataPoint.C
indexedOctree/treeDataTriSurface.C
searchableSurface = searchableSurface
$(searchableSurface)/searchableBox.C
$(searchableSurface)/searchableSphere.C
$(searchableSurface)/searchableSurface.C
$(searchableSurface)/searchableSurfaces.C
$(searchableSurface)/searchableSurfacesQueries.C
$(searchableSurface)/triSurfaceMesh.C
topoSets = sets/topoSets
$(topoSets)/cellSet.C
$(topoSets)/topoSet.C
@ -117,20 +125,10 @@ $(intersectedSurface)/intersectedSurface.C
$(intersectedSurface)/edgeSurface.C
triSurface/triSurfaceSearch/triSurfaceSearch.C
triSurface/octreeData/octreeDataTriSurface.C
triSurface/octreeData/octreeDataTriSurfaceTreeLeaf.C
triSurface/triangleFuncs/triangleFuncs.C
triSurface/searchableSurface/searchableSurface.C
triSurface/searchableSurface/triSurfaceMesh.C
triSurface/searchableSurface/searchableBox.C
triSurface/surfaceFeatures/surfaceFeatures.C
triSurface/triSurfaceMeshes/triSurfaceMeshes.C
triSurface/triSurfaceTools/triSurfaceTools.C
triSurface/triSurfaceTools/geompack/geompack.C

541
src/meshTools/searchableSurface/searchableBox.C Normal file
View File

@ -0,0 +1,541 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "searchableBox.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableBox, 0);
addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::searchableBox::projectOntoCoordPlane
(
const direction dir,
const point& planePt,
pointIndexHit& info
) const
{
// Set point
info.rawPoint()[dir] = planePt[dir];
// Set face
if (planePt[dir] == min()[dir])
{
info.setIndex(dir*2);
}
else if (planePt[dir] == max()[dir])
{
info.setIndex(dir*2+1);
}
else
{
FatalErrorIn("searchableBox::projectOntoCoordPlane(..)")
<< "Point on plane " << planePt
<< " is not on coordinate " << min()[dir]
<< " nor " << max()[dir] << abort(FatalError);
}
}
// Returns miss or hit with face (0..5) and region(always 0)
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const point& bbMid,
const point& sample,
const scalar nearestDistSqr
) const
{
// Point can be inside or outside. For every component direction can be
// left of min, right of max or inbetween.
// - outside points: project first one x plane (either min().x()
// or max().x()), then onto y plane and finally z. You should be left
// with intersection point
// - inside point: find nearest side (compare to mid point). Project onto
// that.
// The face is set to the last projected face.
// Outside point projected onto cube. Assume faces 0..5.
pointIndexHit info(true, sample, -1);
bool outside = false;
// (for internal points) per direction what nearest cube side is
point near;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (info.rawPoint()[dir] < min()[dir])
{
projectOntoCoordPlane(dir, min(), info);
outside = true;
}
else if (info.rawPoint()[dir] > max()[dir])
{
projectOntoCoordPlane(dir, max(), info);
outside = true;
}
else if (info.rawPoint()[dir] > bbMid[dir])
{
near[dir] = max()[dir];
}
else
{
near[dir] = min()[dir];
}
}
// For outside points the info will be correct now. Handle inside points
// using the three near distances. Project onto the nearest plane.
if (!outside)
{
vector dist(cmptMag(info.rawPoint() - near));
if (dist.x() < dist.y())
{
if (dist.x() < dist.z())
{
// Project onto x plane
projectOntoCoordPlane(vector::X, near, info);
}
else
{
projectOntoCoordPlane(vector::Z, near, info);
}
}
else
{
if (dist.y() < dist.z())
{
projectOntoCoordPlane(vector::Y, near, info);
}
else
{
projectOntoCoordPlane(vector::Z, near, info);
}
}
}
// Check if outside. Optimisation: could do some checks on distance already
// on components above
if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
{
info.setMiss();
info.setIndex(-1);
}
return info;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableBox::searchableBox
(
const IOobject& io,
const treeBoundBox& bb
)
:
searchableSurface(io),
treeBoundBox(bb)
{}
Foam::searchableBox::searchableBox
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
treeBoundBox(dict.lookup("min"), dict.lookup("max"))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableBox::~searchableBox()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchableBox::regions() const
{
if (regions_.size() == 0)
{
regions_.setSize(1);
regions_[0] = "region0";
}
return regions_;
}
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
return findNearest(mid(), sample, nearestDistSqr);
}
Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const
{
const point bbMid(mid());
// Outside point projected onto cube. Assume faces 0..5.
pointIndexHit info(true, sample, -1);
bool outside = false;
// (for internal points) per direction what nearest cube side is
point near;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (info.rawPoint()[dir] < min()[dir])
{
projectOntoCoordPlane(dir, min(), info);
outside = true;
}
else if (info.rawPoint()[dir] > max()[dir])
{
projectOntoCoordPlane(dir, max(), info);
outside = true;
}
else if (info.rawPoint()[dir] > bbMid[dir])
{
near[dir] = max()[dir];
}
else
{
near[dir] = min()[dir];
}
}
// For outside points the info will be correct now. Handle inside points
// using the three near distances. Project onto the nearest two planes.
if (!outside)
{
// Get the per-component distance to nearest wall
vector dist(cmptMag(info.rawPoint() - near));
SortableList<scalar> sortedDist(3);
sortedDist[0] = dist[0];
sortedDist[1] = dist[1];
sortedDist[2] = dist[2];
sortedDist.sort();
// Project onto nearest
projectOntoCoordPlane(sortedDist.indices()[0], near, info);
// Project onto second nearest
projectOntoCoordPlane(sortedDist.indices()[1], near, info);
}
// Check if outside. Optimisation: could do some checks on distance already
// on components above
if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
{
info.setMiss();
info.setIndex(-1);
}
return info;
}
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const
{
notImplemented
(
"searchableBox::findNearest"
"(const linePointRef&, treeBoundBox&, point&)"
);
return pointIndexHit();
}
Foam::pointIndexHit Foam::searchableBox::findLine
(
const point& start,
const point& end
) const
{
pointIndexHit info(false, start, -1);
bool foundInter;
if (posBits(start) == 0)
{
if (posBits(end) == 0)
{
// Both start and end inside.
foundInter = false;
}
else
{
// end is outside. Clip to bounding box.
foundInter = intersects(end, start, info.rawPoint());
}
}
else
{
// start is outside. Clip to bounding box.
foundInter = intersects(start, end, info.rawPoint());
}
// Classify point
if (foundInter)
{
info.setHit();
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (info.rawPoint()[dir] == min()[dir])
{
info.setIndex(2*dir);
break;
}
else if (info.rawPoint()[dir] == max()[dir])
{
info.setIndex(2*dir+1);
break;
}
}
if (info.index() == -1)
{
FatalErrorIn("searchableBox::findLine(const point&, const point&)")
<< "point " << info.rawPoint()
<< " on segment " << start << end
<< " should be on face of " << *this
<< " but it isn't." << abort(FatalError);
}
}
return info;
}
Foam::pointIndexHit Foam::searchableBox::findLineAny
(
const point& start,
const point& end
) const
{
return findLine(start, end);
}
void Foam::searchableBox::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
info.setSize(samples.size());
const point bbMid(mid());
forAll(samples, i)
{
info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
}
}
void Foam::searchableBox::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
forAll(start, i)
{
info[i] = findLine(start[i], end[i]);
}
}
void Foam::searchableBox::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
forAll(start, i)
{
info[i] = findLineAny(start[i], end[i]);
}
}
void Foam::searchableBox::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
info.setSize(start.size());
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
// Tolerances
const vectorField dirVec(end-start);
const scalarField magSqrDirVec(magSqr(dirVec));
const vectorField smallVec
(
Foam::sqrt(SMALL)*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
forAll(start, pointI)
{
hits.clear();
// Current starting point of ray.
point pt = start[pointI];
while (true)
{
// See if any intersection between pt and end
pointIndexHit inter = findLine(pt, end[pointI]);
if (!inter.hit())
{
break;
}
hits.append(inter);
pt = inter.hitPoint() + smallVec[pointI];
if (((pt-start[pointI])&dirVec[pointI]) > magSqrDirVec[pointI])
{
// Adding smallVec has taken us beyond end
break;
}
}
hits.shrink();
info[pointI].transfer(hits);
hits.clear();
}
}
void Foam::searchableBox::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
region = 0;
}
void Foam::searchableBox::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
normal.setSize(info.size());
normal = vector::zero;
forAll(info, i)
{
if (info[i].hit())
{
normal[i] = treeBoundBox::faceNormals[info[i].index()];
}
else
{
// Set to what?
}
}
}
void Foam::searchableBox::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
volType.setSize(points.size());
volType = INSIDE;
forAll(points, pointI)
{
const point& pt = points[pointI];
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
{
volType[pointI] = OUTSIDE;
break;
}
}
}
}
// ************************************************************************* //

249
src/meshTools/searchableSurface/searchableBox.H Normal file
View File

@ -0,0 +1,249 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableBox
Description
Searching on bounding box
SourceFiles
searchableBox.C
\*---------------------------------------------------------------------------*/
#ifndef searchableBox_H
#define searchableBox_H
#include "searchableSurface.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableBox Declaration
\*---------------------------------------------------------------------------*/
class searchableBox
:
public searchableSurface,
public treeBoundBox
{
private:
// Private Member Data
mutable wordList regions_;
// Private Member Functions
//- Project onto component dir of planePt and update index() (=face)
void projectOntoCoordPlane
(
const direction dir,
const point& planePt,
pointIndexHit& info
) const;
//- Returns miss or hit with face (0..5)
pointIndexHit findNearest
(
const point& bbMid,
const point& sample,
const scalar nearestDistSqr
) const;
//- Disallow default bitwise copy construct
searchableBox(const searchableBox&);
//- Disallow default bitwise assignment
void operator=(const searchableBox&);
public:
//- Runtime type information
TypeName("searchableBox");
// Constructors
//- Construct from components
searchableBox(const IOobject& io, const treeBoundBox& bb);
//- Construct from dictionary (used by searchableSurface)
searchableBox
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableBox();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return true;
}
// Single point queries.
//- Calculate nearest point on surface. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface (=face 0..5)
// - point: actual nearest point found
pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Calculate nearest point on edge. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface(=?)
// - point: actual nearest point found
pointIndexHit findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Find nearest to segment. Returns
// - bool : any point found?
// - label: relevant index in shapes (=face 0..5)
// - point: actual nearest point found
// sets:
// - tightest : bounding box
// - linePoint : corresponding nearest point on line
pointIndexHit findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const;
//- Find nearest intersection of line between start and end.
pointIndexHit findLine
(
const point& start,
const point& end
) const;
//- Find any intersection of line between start and end.
pointIndexHit findLineAny
(
const point& start,
const point& end
) const;
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("searchableBox::writeData(Ostream&) const");
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

328
src/meshTools/searchableSurface/searchableSphere.C Normal file
View File

@ -0,0 +1,328 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "searchableSphere.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableSphere, 0);
addToRunTimeSelectionTable(searchableSurface, searchableSphere, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointIndexHit Foam::searchableSphere::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
pointIndexHit info(false, sample, -1);
const vector n(sample-centre_);
scalar magN = mag(n);
if (nearestDistSqr > sqr(magN-radius_))
{
info.rawPoint() = centre_ + n/magN*radius_;
info.setHit();
info.setIndex(0);
}
return info;
}
// From Graphics Gems - intersection of sphere with ray
void Foam::searchableSphere::findLineAll
(
const point& start,
const point& end,
pointIndexHit& near,
pointIndexHit& far
) const
{
near.setMiss();
far.setMiss();
vector dir(end-start);
scalar magSqrDir = magSqr(dir);
if (magSqrDir > ROOTVSMALL)
{
const vector toCentre(centre_-start);
scalar magSqrToCentre = magSqr(toCentre);
dir /= Foam::sqrt(magSqrDir);
scalar v = (toCentre & dir);
scalar disc = sqr(radius_) - (magSqrToCentre - sqr(v));
if (disc >= 0)
{
scalar d = Foam::sqrt(disc);
scalar nearParam = v-d;
if (nearParam >= 0 && sqr(nearParam) <= magSqrDir)
{
near.setHit();
near.setPoint(start + nearParam*dir);
near.setIndex(0);
}
scalar farParam = v+d;
if (farParam >= 0 && sqr(farParam) <= magSqrDir)
{
far.setHit();
far.setPoint(start + farParam*dir);
far.setIndex(0);
}
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSphere::searchableSphere
(
const IOobject& io,
const point& centre,
const scalar radius
)
:
searchableSurface(io),
centre_(centre),
radius_(radius)
{}
Foam::searchableSphere::searchableSphere
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
centre_(dict.lookup("centre")),
radius_(readScalar(dict.lookup("radius")))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableSphere::~searchableSphere()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchableSphere::regions() const
{
if (regions_.size() == 0)
{
regions_.setSize(1);
regions_[0] = "region0";
}
return regions_;
}
void Foam::searchableSphere::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
info.setSize(samples.size());
forAll(samples, i)
{
info[i] = findNearest(samples[i], nearestDistSqr[i]);
}
}
void Foam::searchableSphere::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Pick nearest intersection. If none intersected take second one.
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{
info[i] = b;
}
}
}
void Foam::searchableSphere::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Discard far intersection
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{
info[i] = b;
}
}
}
void Foam::searchableSphere::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
info.setSize(start.size());
pointIndexHit near, far;
forAll(start, i)
{
findLineAll(start[i], end[i], near, far);
if (near.hit())
{
if (far.hit())
{
info[i].setSize(2);
info[i][0] = near;
info[i][1] = far;
}
else
{
info[i].setSize(1);
info[i][0] = near;
}
}
else
{
if (far.hit())
{
info[i].setSize(1);
info[i][0] = far;
}
}
}
}
void Foam::searchableSphere::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
region = 0;
}
void Foam::searchableSphere::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
normal.setSize(info.size());
normal = vector::zero;
forAll(info, i)
{
if (info[i].hit())
{
normal[i] = info[i].hitPoint() - centre_;
normal[i] /= mag(normal[i]);
}
else
{
// Set to what?
}
}
}
void Foam::searchableSphere::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
volType.setSize(points.size());
volType = INSIDE;
forAll(points, pointI)
{
const point& pt = points[pointI];
if (magSqr(pt - centre_) <= sqr(radius_))
{
volType[pointI] = INSIDE;
}
else
{
volType[pointI] = OUTSIDE;
}
}
}
// ************************************************************************* //

204
src/meshTools/searchableSurface/searchableSphere.H Normal file
View File

@ -0,0 +1,204 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableSphere
Description
Searching on sphere
SourceFiles
searchableSphere.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSphere_H
#define searchableSphere_H
#include "searchableSurface.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableSphere Declaration
\*---------------------------------------------------------------------------*/
class searchableSphere
:
public searchableSurface
{
private:
// Private Member Data
//- Centre point
const point centre_;
//- Radius squared
const scalar radius_;
//- Names of regions
mutable wordList regions_;
// Private Member Functions
//- Find nearest point on sphere.
pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Find intersection with sphere
void findLineAll
(
const point& start,
const point& end,
pointIndexHit& near,
pointIndexHit& far
) const;
//- Disallow default bitwise copy construct
searchableSphere(const searchableSphere&);
//- Disallow default bitwise assignment
void operator=(const searchableSphere&);
public:
//- Runtime type information
TypeName("searchableSphere");
// Constructors
//- Construct from components
searchableSphere(const IOobject& io, const point&, const scalar radius);
//- Construct from dictionary (used by searchableSurface)
searchableSphere
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableSphere();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return true;
}
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("searchableSphere::writeData(Ostream&) const");
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

43
src/meshTools/triSurface/searchableSurface/searchableSurface.C → src/meshTools/searchableSurface/searchableSurface.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,15 +34,13 @@ namespace Foam
defineTypeNameAndDebug(searchableSurface, 0);
defineRunTimeSelectionTable(searchableSurface, dict);
//defineRunTimeSelectionTable(searchableSurface, istream);
// Construct named object from dictionary
autoPtr<searchableSurface> searchableSurface::New
(
const word& searchableSurfaceType,
const word& name,
const objectRegistry& obj,
const IOobject& io,
const dictionary& dict
)
{
@ -55,7 +53,7 @@ autoPtr<searchableSurface> searchableSurface::New
FatalErrorIn
(
"searchableSurface::New(const word&, const word&"
", const objectRegistry&, const dictionary&)"
", const IOobject&, const dictionary&)"
) << "Unknown searchableSurface type " << searchableSurfaceType
<< endl << endl
<< "Valid searchableSurface types : " << endl
@ -63,44 +61,15 @@ autoPtr<searchableSurface> searchableSurface::New
<< exit(FatalError);
}
return autoPtr<searchableSurface>(cstrIter()(name, obj, dict));
return autoPtr<searchableSurface>(cstrIter()(io, dict));
}
//// Construct named object from Istream
//autoPtr<searchableSurface> searchableSurface::New
//(
// const word& searchableSurfaceType,
// const objectRegistry& obj,
// Istream& is
//)
//{
// istreamConstructorTable::iterator cstrIter =
// istreamConstructorTablePtr_
// ->find(searchableSurfaceType);
//
// if (cstrIter == istreamConstructorTablePtr_->end())
// {
// FatalErrorIn
// (
// "searchableSurface::New(const word&, const objectRegistry&"
// ", Istream&)"
// ) << "Unknown searchableSurface type " << searchableSurfaceType
// << endl << endl
// << "Valid searchableSurface types : " << endl
// << istreamConstructorTablePtr_->toc()
// << exit(FatalError);
// }
//
// return autoPtr<searchableSurface>(cstrIter()(obj, is));
//}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSurface::searchableSurface(const word& name)
Foam::searchableSurface::searchableSurface(const IOobject& io)
:
name_(name)
regIOobject(io)
{}

322
src/meshTools/searchableSurface/searchableSurface.H Normal file
View File

@ -0,0 +1,322 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableSurface
Description
Base class of (analytical or triangulated) surface.
Encapsulates all the search routines. WIP.
Information returned is usually a pointIndexHit:
- bool : was intersection/nearest found?
- point : intersection point or nearest point
- index : unique index on surface (e.g. triangle for triSurfaceMesh)
SourceFiles
searchableSurface.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurface_H
#define searchableSurface_H
#include "pointField.H"
#include "typeInfo.H"
#include "runTimeSelectionTables.H"
#include "pointIndexHit.H"
#include "linePointRef.H"
#include "objectRegistry.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class objectRegistry;
class treeBoundBox;
/*---------------------------------------------------------------------------*\
Class searchableSurface Declaration
\*---------------------------------------------------------------------------*/
class searchableSurface
:
public regIOobject
{
public:
// Data types
//- volume types
enum volumeType
{
UNKNOWN = 0,
MIXED = 1, // not used. only here to maintain consistency with
// indexedOctree volumeType.
INSIDE = 2,
OUTSIDE = 3
};
private:
// Private data
const word name_;
// Private Member Functions
//- Disallow default bitwise copy construct
searchableSurface(const searchableSurface&);
//- Disallow default bitwise assignment
void operator=(const searchableSurface&);
public:
//- Runtime type information
TypeName("searchableSurface");
// Declare run-time constructor selection table
// For the dictionary constructor
declareRunTimeSelectionTable
(
autoPtr,
searchableSurface,
dict,
(
const IOobject& io,
const dictionary& dict
),
(io, dict)
);
//- Class used for the read-construction of
// PtrLists of searchableSurface.
class iNew
{
IOobject& io_;
public:
iNew(IOobject& io)
:
io_(io)
{}
autoPtr<searchableSurface> operator()(Istream& is) const
{
word surfaceType(is);
word readName(is);
dictionary dict(is);
autoPtr<IOobject> namedIO(io_.clone());
namedIO().rename(readName);
return searchableSurface::New(surfaceType, namedIO(), dict);
}
};
// Constructors
searchableSurface(const IOobject& io);
//- Clone
virtual autoPtr<searchableSurface> clone() const
{
notImplemented("autoPtr<searchableSurface> clone() const");
return autoPtr<searchableSurface>(NULL);
}
// Selectors
//- Return a reference to the selected searchableSurface
static autoPtr<searchableSurface> New
(
const word& surfaceType,
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableSurface();
// Member Functions
//- Names of regions.
virtual const wordList& regions() const = 0;
//- Whether supports volume type below.
virtual bool hasVolumeType() const = 0;
// Single point queries.
////- Calculate nearest point on surface. Returns
//// - bool : any point found nearer than nearestDistSqr
//// - label: relevant index in surface
//// - label: region in surface
//// - point: actual nearest point found
//virtual pointIndexHit findNearest
//(
// const point& sample,
// const scalar nearestDistSqr
//) const = 0;
//
////- Calculate nearest point on edge. Returns
//// - bool : any point found nearer than nearestDistSqr
//// - label: relevant index in surface
//// - label: region in surface
//// - point: actual nearest point found
//virtual pointIndexHit findNearestOnEdge
//(
// const point& sample,
// const scalar nearestDistSqr
//) const = 0;
//
////- Find nearest to segment. Returns
//// - bool : any point found?
//// - label: relevant index in shapes
//// - label: region in surface
//// - point: actual nearest point found
//// sets:
//// - tightest : bounding box
//// - linePoint : corresponding nearest point on line
//virtual pointIndexHit findNearest
//(
// const linePointRef& ln,
// treeBoundBox& tightest,
// point& linePoint
//) const = 0;
//
////- Find nearest intersection of line between start and end.
//virtual pointIndexHit findLine
//(
// const point& start,
// const point& end
//) const = 0;
//
////- Find any intersection of line between start and end.
//virtual pointIndexHit findLineAny
//(
// const point& start,
// const point& end
//) const = 0;
// Multiple point queries. When surface is distributed the index
// should be a global index. Not done yet.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const = 0;
//- Find first intersection on segment from start to end.
// Note: searchableSurfacesQueries expects no
// intersection to be found if start==end. Is problem?
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const = 0;
//- Return any intersection on segment from start to end.
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const = 0;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const = 0;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const = 0;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const = 0;
//- Determine type (inside/outside) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const = 0;
// Other
////- Get bounding box.
//const boundBox& bounds() const = 0;
////- Set bounding box.
//void setBounds
//(
// const boundBox&,
// autoPtr<mapDistribute>& faceMap,
// autoPtr<mapDistribute>& pointMap
//) = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

355
src/meshTools/searchableSurface/searchableSurfaces.C Normal file
View File

@ -0,0 +1,355 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "searchableSurfaces.H"
#include "searchableSurfacesQueries.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableSurfaces, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct with length.
Foam::searchableSurfaces::searchableSurfaces(const label size)
:
PtrList<searchableSurface>(size),
regionNames_(size),
allSurfaces_(identity(size))
{}
//Foam::searchableSurfaces::searchableSurfaces
//(
// const IOobject& io,
// const PtrList<dictionary>& dicts
//)
//:
// PtrList<searchableSurface>(dicts.size()),
// regionNames_(dicts.size()),
// allSurfaces_(identity(dicts.size()))
//{
// forAll(dicts, surfI)
// {
// const dictionary& dict = dicts[surfI];
//
// // Make IOobject with correct name
// autoPtr<IOobject> namedIO(io.clone());
// namedIO().rename(dict.lookup("name"));
//
// // Create and hook surface
// set
// (
// surfI,
// searchableSurface::New
// (
// dict.lookup("type"),
// namedIO(),
// dict
// )
// );
// const searchableSurface& s = operator[](surfI);
//
// // Construct default region names by prepending surface name
// // to region name.
// const wordList& localNames = s.regions();
//
// wordList globalNames(localNames.size());
// forAll(localNames, regionI)
// {
// globalNames[regionI] = s.name() + '_' + localNames[regionI];
// }
//
// // See if dictionary provides any global region names.
// if (dict.found("regions"))
// {
// const dictionary& regionsDict = dict.subDict("regions");
//
// forAllConstIter(dictionary, regionsDict, iter)
// {
// const word& key = iter().keyword();
//
// if (regionsDict.isDict(key))
// {
// // Get the dictionary for region iter.key()
// const dictionary& regionDict = regionsDict.subDict(key);
//
// label index = findIndex(localNames, key);
//
// if (index == -1)
// {
// FatalErrorIn
// (
// "searchableSurfaces::searchableSurfaces"
// "( const IOobject&, const dictionary&)"
// ) << "Unknown region name " << key
// << " for surface " << s.name() << endl
// << "Valid region names are " << localNames
// << exit(FatalError);
// }
//
// globalNames[index] = word(regionDict.lookup("name"));
// }
// }
// }
//
// // Now globalNames contains the names of the regions.
// Info<< "Surface:" << s.name() << " has regions:"
// << endl;
// forAll(globalNames, regionI)
// {
// Info<< " " << globalNames[regionI] << endl;
// }
//
// // Create reverse lookup
// forAll(globalNames, regionI)
// {
// regionNames_.insert
// (
// globalNames[regionI],
// labelPair(surfI, regionI)
// );
// }
// }
//}
Foam::searchableSurfaces::searchableSurfaces
(
const IOobject& io,
const dictionary& topDict
)
:
PtrList<searchableSurface>(topDict.size()),
names_(topDict.size()),
regionNames_(topDict.size()),
allSurfaces_(identity(topDict.size()))
{
label surfI = 0;
forAllConstIter(dictionary, topDict, iter)
{
const word& key = iter().keyword();
if (!topDict.isDict(key))
{
FatalErrorIn
(
"searchableSurfaces::searchableSurfaces"
"( const IOobject&, const dictionary&)"
) << "Found non-dictionary entry " << iter()
<< " in top-level dictionary " << topDict
<< exit(FatalError);
}
const dictionary& dict = topDict.subDict(key);
names_[surfI] = key;
if (dict.found("name"))
{
dict.lookup("name") >> names_[surfI];
}
// Make IOobject with correct name
autoPtr<IOobject> namedIO(io.clone());
// Note: we would like to e.g. register triSurface 'sphere.stl' as
// 'sphere'. Unfortunately
// no support for having object read from different location than
// their object name. Maybe have stlTriSurfaceMesh which appends .stl
// when reading/writing?
namedIO().rename(key); // names_[surfI]
// Create and hook surface
set
(
surfI,
searchableSurface::New
(
dict.lookup("type"),
namedIO(),
dict
)
);
const searchableSurface& s = operator[](surfI);
// Construct default region names by prepending surface name
// to region name.
const wordList& localNames = s.regions();
wordList& rNames = regionNames_[surfI];
rNames.setSize(localNames.size());
forAll(localNames, regionI)
{
rNames[regionI] = names_[surfI] + '_' + localNames[regionI];
}
// See if dictionary provides any global region names.
if (dict.found("regions"))
{
const dictionary& regionsDict = dict.subDict("regions");
forAllConstIter(dictionary, regionsDict, iter)
{
const word& key = iter().keyword();
if (regionsDict.isDict(key))
{
// Get the dictionary for region iter.keyword()
const dictionary& regionDict = regionsDict.subDict(key);
label index = findIndex(localNames, key);
if (index == -1)
{
FatalErrorIn
(
"searchableSurfaces::searchableSurfaces"
"( const IOobject&, const dictionary&)"
) << "Unknown region name " << key
<< " for surface " << s.name() << endl
<< "Valid region names are " << localNames
<< exit(FatalError);
}
rNames[index] = word(regionDict.lookup("name"));
}
}
}
surfI++;
}
// Trim (not really necessary since we don't allow non-dictionary entries)
PtrList<searchableSurface>::setSize(surfI);
names_.setSize(surfI);
regionNames_.setSize(surfI);
allSurfaces_.setSize(surfI);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::searchableSurfaces::findSurfaceID(const word& wantedName)
const
{
return findIndex(names_, wantedName);
}
// Find any intersection
void Foam::searchableSurfaces::findAnyIntersection
(
const pointField& start,
const pointField& end,
labelList& hitSurfaces,
List<pointIndexHit>& hitInfo
) const
{
searchableSurfacesQueries::findAnyIntersection
(
*this,
allSurfaces_,
start,
end,
hitSurfaces,
hitInfo
);
}
// Find intersections of edge nearest to both endpoints.
void Foam::searchableSurfaces::findAllIntersections
(
const pointField& start,
const pointField& end,
labelListList& hitSurfaces,
List<List<pointIndexHit> >& hitInfo
) const
{
searchableSurfacesQueries::findAllIntersections
(
*this,
allSurfaces_,
start,
end,
hitSurfaces,
hitInfo
);
}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfaces::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
) const
{
return searchableSurfacesQueries::findNearest
(
*this,
allSurfaces_,
samples,
nearestDistSqr,
nearestSurfaces,
nearestInfo
);
}
//- Calculate point which is on a set of surfaces.
Foam::pointIndexHit Foam::searchableSurfaces::facesIntersection
(
const scalar initDistSqr,
const scalar convergenceDistSqr,
const point& start
) const
{
return searchableSurfacesQueries::facesIntersection
(
*this,
allSurfaces_,
initDistSqr,
convergenceDistSqr,
start
);
}
// ************************************************************************* //

187
src/meshTools/searchableSurface/searchableSurfaces.H Normal file
View File

@ -0,0 +1,187 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableSurfaces
Description
Container for searchableSurfaces.
SourceFiles
searchableSurfaces.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurfaces_H
#define searchableSurfaces_H
#include "searchableSurface.H"
#include "labelPair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableSurfaces Declaration
\*---------------------------------------------------------------------------*/
class searchableSurfaces
:
public PtrList<searchableSurface>
{
// Private data
//- Surface names
wordList names_;
//- Region names per surface
List<wordList> regionNames_;
////- From global region name to surface and region on surface
//HashTable<labelPair> regionNames_;
//- Indices of all surfaces. Precalculated and stored.
labelList allSurfaces_;
//- Disallow default bitwise copy construct
searchableSurfaces(const searchableSurfaces&);
//- Disallow default bitwise assignment
void operator=(const searchableSurfaces&);
public:
ClassName("searchableSurfaces");
// Constructors
//- Construct with length specified. Fill later.
explicit searchableSurfaces(const label);
////- Construct from list of dictionaries
//searchableSurfaces(const IOobject&, const PtrList<dictionary>&);
//- Construct from dictionary
searchableSurfaces(const IOobject&, const dictionary&);
// Member Functions
const wordList& names() const
{
return names_;
}
wordList& names()
{
return names_;
}
const List<wordList>& regionNames() const
{
return regionNames_;
}
List<wordList>& regionNames()
{
return regionNames_;
}
////- If adding surfaces 'by hand'
//HashTable<labelPair>& regionNames()
//{
// return regionNames_;
//}
////- Get surface and region for a name
//const labelPair& surfaceRegion(const word& globalRegion) const
//{
// return regionNames_[globalRegion];
//}
//- Find index of surface. Return -1 if not found.
label findSurfaceID(const word& name) const;
// Multiple point queries.
//- Find any intersection. Return hit point information and
// surface number. If multiple surfaces hit the first surface
// is returned, not necessarily the nearest (to start).
void findAnyIntersection
(
const pointField& start,
const pointField& end,
labelList& surfaces,
List<pointIndexHit>&
) const;
//- Find all intersections in order from start to end. Returns for
// every hit the surface and the hit info.
void findAllIntersections
(
const pointField& start,
const pointField& end,
labelListList& surfaces,
List<List<pointIndexHit> >& surfaceHits
) const;
//- Find nearest. Return -1 (and a miss()) or surface and nearest
// point.
void findNearest
(
const pointField&,
const scalarField& nearestDistSqr,
labelList& surfaces,
List<pointIndexHit>&
) const;
// Single point queries
//- Calculate point which is on a set of surfaces.
pointIndexHit facesIntersection
(
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const point& start
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

822
src/meshTools/searchableSurface/searchableSurfacesQueries.C Normal file
View File

@ -0,0 +1,822 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "searchableSurfacesQueries.H"
#include "SortableList.H"
#include "OFstream.H"
#include "meshTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableSurfacesQueries, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointIndexHit Foam::searchableSurfacesQueries::tempFindNearest
(
const searchableSurface& surf,
const point& pt,
const scalar initDistSqr
)
{
pointField onePoint(1, pt);
scalarField oneDist(1, initDistSqr);
List<pointIndexHit> oneHit(1);
surf.findNearest(onePoint, oneDist, oneHit);
return oneHit[0];
}
// Calculate sum of distance to surfaces.
Foam::scalar Foam::searchableSurfacesQueries::sumDistSqr
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const scalar initDistSqr,
const point& pt
)
{
scalar sum = 0;
forAll(surfacesToTest, testI)
{
label surfI = surfacesToTest[testI];
pointIndexHit hit
(
tempFindNearest(allSurfaces[surfI], pt, initDistSqr)
);
// Note: make it fall over if not hit.
sum += magSqr(hit.hitPoint()-pt);
}
return sum;
}
// Reflects the point furthest away around the triangle centre by a factor fac.
// (triangle centre is the average of all points but the ihi. pSum is running
// sum of all points)
Foam::scalar Foam::searchableSurfacesQueries::tryMorphTet
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const scalar initDistSqr,
List<vector>& p,
List<scalar>& y,
vector& pSum,
const label ihi,
const scalar fac
)
{
scalar fac1 = (1.0-fac)/vector::nComponents;
scalar fac2 = fac1-fac;
vector ptry = pSum*fac1-p[ihi]*fac2;
scalar ytry = sumDistSqr(allSurfaces, surfacesToTest, initDistSqr, ptry);
if (ytry < y[ihi])
{
y[ihi] = ytry;
pSum += ptry - p[ihi];
p[ihi] = ptry;
}
return ytry;
}
bool Foam::searchableSurfacesQueries::morphTet
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const scalar initDistSqr,
const scalar convergenceDistSqr,
const label maxIter,
List<vector>& p,
List<scalar>& y
)
{
vector pSum = sum(p);
autoPtr<OFstream> str;
label vertI = 0;
if (debug)
{
wordList names(surfacesToTest.size());
forAll(surfacesToTest, i)
{
names[i] = allSurfaces[surfacesToTest[i]].name();
}
Pout<< "searchableSurfacesQueries::morphTet : intersection of "
<< names << " starting from points:" << p << endl;
str.reset(new OFstream("track.obj"));
meshTools::writeOBJ(str(), p[0]);
vertI++;
}
for (label iter = 0; iter < maxIter; iter++)
{
// Get the indices of highest, second-highest and lowest values.
label ihi, inhi, ilo;
{
SortableList<scalar> sortedY(y);
ilo = sortedY.indices()[0];
ihi = sortedY.indices()[sortedY.size()-1];
inhi = sortedY.indices()[sortedY.size()-2];
}
if (debug)
{
Pout<< "Iteration:" << iter
<< " lowest:" << y[ilo] << " highest:" << y[ihi]
<< " points:" << p << endl;
meshTools::writeOBJ(str(), p[ilo]);
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
if (y[ihi] < convergenceDistSqr)
{
// Get point on 0th surface.
Swap(p[0], p[ilo]);
Swap(y[0], y[ilo]);
return true;
}
// Reflection: point furthest away gets reflected.
scalar ytry = tryMorphTet
(
allSurfaces,
surfacesToTest,
10*y[ihi], // search box.
p,
y,
pSum,
ihi,
-1.0
);
if (ytry <= y[ilo])
{
// If in right direction (y lower) expand by two.
ytry = tryMorphTet
(
allSurfaces,
surfacesToTest,
10*y[ihi],
p,
y,
pSum,
ihi,
2.0
);
}
else if (ytry >= y[inhi])
{
// If inside tet try contraction.
scalar ysave = y[ihi];
ytry = tryMorphTet
(
allSurfaces,
surfacesToTest,
10*y[ihi],
p,
y,
pSum,
ihi,
0.5
);
if (ytry >= ysave)
{
// Contract around lowest point.
forAll(p, i)
{
if (i != ilo)
{
p[i] = 0.5*(p[i] + p[ilo]);
y[i] = sumDistSqr
(
allSurfaces,
surfacesToTest,
y[ihi],
p[i]
);
}
}
pSum = sum(p);
}
}
}
if (debug)
{
meshTools::writeOBJ(str(), p[0]);
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
// Failure to converge. Return best guess so far.
label ilo = findMin(y);
Swap(p[0], p[ilo]);
Swap(y[0], y[ilo]);
return false;
}
//// Get all intersections (in order) for single surface.
//void Foam::searchableSurfacesQueries::findAllIntersections
//(
// const searchableSurface& s,
// const pointField& start,
// const pointField& end,
// const vectorField& smallVec,
// List<List<pointIndexHit> >& surfaceHitInfo
//)
//{
// surfaceHitInfo.setSize(start.size());
//
// // Current start point of vector
// pointField p0(start);
//
// List<pointIndexHit> intersectInfo(start.size());
//
// // For test whether finished doing vector.
// const vectorField dirVec(end-start);
// const scalarField magSqrDirVec(magSqr(dirVec));
//
// while (true)
// {
// // Find first intersection. Synced.
// s.findLine(p0, end, intersectInfo);
//
// label nHits = 0;
//
// forAll(intersectInfo, i)
// {
// if (intersectInfo[i].hit())
// {
// nHits++;
//
// label sz = surfaceHitInfo[i].size();
// surfaceHitInfo[i].setSize(sz+1);
// surfaceHitInfo[i][sz] = intersectInfo[i];
//
// p0[i] = intersectInfo[i].hitPoint() + smallVec[i];
//
// // If beyond endpoint set to endpoint so as not to pick up
// // any intersections. Could instead just filter out hits.
// if (((p0[i]-start[i])&dirVec[i]) > magSqrDirVec[i])
// {
// p0[i] = end[i];
// }
// }
// else
// {
// // Set to endpoint to stop intersection test. See above.
// p0[i] = end[i];
// }
// }
//
// // returnReduce(nHits) ?
// if (nHits == 0)
// {
// break;
// }
// }
//}
// Given current set of hits (allSurfaces, allInfo) merge in those coming from
// surface surfI.
void Foam::searchableSurfacesQueries::mergeHits
(
const point& start,
const scalar mergeDist,
const label testI, // index of surface
const List<pointIndexHit>& surfHits, // hits on surface
labelList& allSurfaces,
List<pointIndexHit>& allInfo,
scalarList& allDistSqr
)
{
// Precalculate distances
scalarList surfDistSqr(surfHits.size());
forAll(surfHits, i)
{
surfDistSqr[i] = magSqr(surfHits[i].hitPoint()-start);
}
forAll(surfDistSqr, i)
{
label index = findLower(allDistSqr, surfDistSqr[i]);
// Check if equal to lower.
if
(
index >= 0
&& (mag(allDistSqr[index]-surfDistSqr[i]) < mergeDist)
)
{
// Same. Do not count.
//Pout<< "point:" << surfHits[i].hitPoint()
// << " considered same as:" << allInfo[index].hitPoint()
// << " within tol:" << mergeDist
// << endl;
}
else
{
// Check if equal to higher
label next = index+1;
if
(
next < allDistSqr.size()
&& (mag(allDistSqr[next]-surfDistSqr[i]) < mergeDist)
)
{
//Pout<< "point:" << surfHits[i].hitPoint()
// << " considered same as:" << allInfo[next].hitPoint()
// << " within tol:" << mergeDist
// << endl;
}
else
{
// Insert after index
label sz = allSurfaces.size();
allSurfaces.setSize(sz+1);
allInfo.setSize(allSurfaces.size());
allDistSqr.setSize(allSurfaces.size());
// Make space.
for (label j = sz-1; j > index; --j)
{
allSurfaces[j+1] = allSurfaces[j];
allInfo[j+1] = allInfo[j];
allDistSqr[j+1] = allDistSqr[j];
}
// Insert new value
allSurfaces[index+1] = testI;
allInfo[index+1] = surfHits[i];
allDistSqr[index+1] = surfDistSqr[i];
}
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Find any intersection
void Foam::searchableSurfacesQueries::findAnyIntersection
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& start,
const pointField& end,
labelList& hitSurfaces,
List<pointIndexHit>& hitInfo
)
{
hitSurfaces.setSize(start.size());
hitSurfaces = -1;
hitInfo.setSize(start.size());
// Work arrays
labelList hitMap(identity(start.size()));
pointField p0(start);
pointField p1(end);
List<pointIndexHit> intersectInfo(start.size());
forAll(surfacesToTest, testI)
{
// Do synchronised call to all surfaces.
allSurfaces[surfacesToTest[testI]].findLineAny(p0, p1, intersectInfo);
// Copy all hits into arguments, continue with misses
label newI = 0;
forAll(intersectInfo, i)
{
if (intersectInfo[i].hit())
{
hitInfo[hitMap[i]] = intersectInfo[i];
hitSurfaces[hitMap[i]] = testI;
}
else
{
if (i != newI)
{
hitMap[newI] = hitMap[i];
p0[newI] = p0[i];
p1[newI] = p1[i];
}
newI++;
}
}
// All done? Note that this decision should be synchronised
if (newI == 0)
{
break;
}
// Trim and continue
hitMap.setSize(newI);
p0.setSize(newI);
p1.setSize(newI);
intersectInfo.setSize(newI);
}
}
void Foam::searchableSurfacesQueries::findAllIntersections
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& start,
const pointField& end,
labelListList& hitSurfaces,
List<List<pointIndexHit> >& hitInfo
)
{
// Note: maybe move the single-surface all intersections test into
// searchable surface? Some of the tolerance issues might be
// lessened.
// 2. Currently calling searchableSurface::findLine with start==end
// is expected to find no intersection. Problem if it does.
hitSurfaces.setSize(start.size());
hitInfo.setSize(start.size());
if (surfacesToTest.size() == 0)
{
return;
}
// Test first surface
allSurfaces[surfacesToTest[0]].findLineAll(start, end, hitInfo);
// Set hitSurfaces and distance
List<scalarList> hitDistSqr(hitInfo.size());
forAll(hitInfo, pointI)
{
const List<pointIndexHit>& pHits = hitInfo[pointI];
labelList& pSurfaces = hitSurfaces[pointI];
pSurfaces.setSize(pHits.size());
pSurfaces = 0;
scalarList& pDistSqr = hitDistSqr[pointI];
pDistSqr.setSize(pHits.size());
forAll(pHits, i)
{
pDistSqr[i] = magSqr(pHits[i].hitPoint() - start[pointI]);
}
}
if (surfacesToTest.size() > 1)
{
// Small vector to increment start vector by to find next intersection
// along line. Constant factor added to make sure that start==end still
// ends iteration in findAllIntersections. Also SMALL is just slightly
// too small.
const vectorField smallVec
(
1E2*SMALL*(end-start)
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
// Tolerance used to check whether points are equal. Note: used to
// compare distance^2. Note that we use the maximum possible tolerance
// (reached at intersections close to the end point)
const scalarField mergeDist(2*mag(smallVec)*mag(end-start));
// Test the other surfaces and merge (according to distance from start).
for (label testI = 1; testI < surfacesToTest.size(); testI++)
{
List<List<pointIndexHit> > surfHits;
allSurfaces[surfacesToTest[testI]].findLineAll
(
start,
end,
surfHits
);
forAll(surfHits, pointI)
{
mergeHits
(
start[pointI], // Current segment
mergeDist[pointI],
testI, // Surface and its hits
surfHits[pointI],
hitSurfaces[pointI], // Merge into overall hit info
hitInfo[pointI],
hitDistSqr[pointI]
);
}
}
}
}
//// Find intersections of edge nearest to both endpoints.
//void Foam::searchableSurfacesQueries::findNearestIntersection
//(
// const PtrList<searchableSurface>& allSurfaces,
// const labelList& surfacesToTest,
// const pointField& start,
// const pointField& end,
//
// labelList& surface1,
// List<pointIndexHit>& hit1,
// labelList& surface2,
// List<pointIndexHit>& hit2
//)
//{
// // 1. intersection from start to end
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// // Initialize arguments
// surface1.setSize(start.size());
// surface1 = -1;
// hit1.setSize(start.size());
//
// // Current end of segment to test.
// pointField nearest(end);
// // Work array
// List<pointIndexHit> nearestInfo(start.size());
//
// forAll(surfacesToTest, testI)
// {
// // See if any intersection between start and current nearest
// allSurfaces[surfacesToTest[testI]].findLine
// (
// start,
// nearest,
// nearestInfo
// );
//
// forAll(nearestInfo, pointI)
// {
// if (nearestInfo[pointI].hit())
// {
// hit1[pointI] = nearestInfo[pointI];
// surface1[pointI] = testI;
// nearest[pointI] = hit1[pointI].hitPoint();
// }
// }
// }
//
//
// // 2. intersection from end to last intersection
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// // Find the nearest intersection from end to start. Note that we
// // initialize to the first intersection (if any).
// surface2 = surface1;
// hit2 = hit1;
//
// // Set current end of segment to test.
// forAll(nearest, pointI)
// {
// if (hit1[pointI].hit())
// {
// nearest[pointI] = hit1[pointI].hitPoint();
// }
// else
// {
// // Disable testing by setting to end.
// nearest[pointI] = end[pointI];
// }
// }
//
// forAll(surfacesToTest, testI)
// {
// // See if any intersection between end and current nearest
// allSurfaces[surfacesToTest[i]].findLine(end, nearest, nearestInfo);
//
// forAll(nearestInfo, pointI)
// {
// if (nearestInfo[pointI].hit())
// {
// hit2[pointI] = nearestInfo[pointI];
// surface2[pointI] = testI;
// nearest[pointI] = hit2[pointI].hitPoint();
// }
// }
// }
//}
// Find nearest. Return -1 or nearest point
void Foam::searchableSurfacesQueries::findNearest
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const pointField& samples,
const scalarField& nearestDistSqr,
labelList& nearestSurfaces,
List<pointIndexHit>& nearestInfo
)
{
// Initialise
nearestSurfaces.setSize(samples.size());
nearestSurfaces = -1;
nearestInfo.setSize(samples.size());
// Work arrays
scalarField minDistSqr(nearestDistSqr);
List<pointIndexHit> hitInfo(samples.size());
forAll(surfacesToTest, testI)
{
allSurfaces[surfacesToTest[testI]].findNearest
(
samples,
minDistSqr,
hitInfo
);
// Update minDistSqr and arguments
forAll(hitInfo, pointI)
{
if (hitInfo[pointI].hit())
{
minDistSqr[pointI] = magSqr
(
hitInfo[pointI].hitPoint()
- samples[pointI]
);
nearestInfo[pointI] = hitInfo[pointI];
nearestSurfaces[pointI] = testI;
}
}
}
}
//- Calculate point which is on a set of surfaces.
Foam::pointIndexHit Foam::searchableSurfacesQueries::facesIntersection
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const scalar initDistSqr,
const scalar convergenceDistSqr,
const point& start
)
{
// Get four starting points. Take these as the projection of the
// starting point onto the surfaces and the mid point
List<point> nearest(surfacesToTest.size()+1);
point sumNearest = vector::zero;
forAll(surfacesToTest, i)
{
pointIndexHit hit
(
tempFindNearest(allSurfaces[surfacesToTest[i]], start, initDistSqr)
);
if (hit.hit())
{
nearest[i] = hit.hitPoint();
sumNearest += nearest[i];
}
else
{
FatalErrorIn
(
"searchableSurfacesQueries::facesIntersection"
"(const labelList&, const scalar, const scalar, const point&)"
) << "Did not find point within distance "
<< initDistSqr << " of starting point " << start
<< " on surface "
<< allSurfaces[surfacesToTest[i]].IOobject::name()
<< abort(FatalError);
}
}
nearest[nearest.size()-1] = sumNearest / surfacesToTest.size();
// Get the sum of distances (initial evaluation)
List<scalar> nearestDist(nearest.size());
forAll(nearestDist, i)
{
nearestDist[i] = sumDistSqr
(
allSurfaces,
surfacesToTest,
initDistSqr,
nearest[i]
);
}
//- Downhill Simplex method
bool converged = morphTet
(
allSurfaces,
surfacesToTest,
initDistSqr,
convergenceDistSqr,
2000,
nearest,
nearestDist
);
pointIndexHit intersection;
if (converged)
{
// Project nearest onto 0th surface.
intersection = tempFindNearest
(
allSurfaces[surfacesToTest[0]],
nearest[0],
nearestDist[0]
);
}
//if (!intersection.hit())
//{
// // Restart
// scalar smallDist = Foam::sqr(convergenceDistSqr);
// nearest[0] = intersection.hitPoint();
// nearest[1] = nearest[0];
// nearest[1].x() += smallDist;
// nearest[2] = nearest[0];
// nearest[2].y() += smallDist;
// nearest[3] = nearest[0];
// nearest[3].z() += smallDist;
//
// forAll(nearestDist, i)
// {
// nearestDist[i] = sumDistSqr
// (
// surfacesToTest,
// initDistSqr,
// nearest[i]
// );
// }
//
// intersection = morphTet
// (
// allSurfaces,
// surfacesToTest,
// initDistSqr,
// convergenceDistSqr,
// 1000,
// nearest,
// nearestDist
// );
//}
return intersection;
}
// ************************************************************************* //

198
src/meshTools/searchableSurface/searchableSurfacesQueries.H Normal file
View File

@ -0,0 +1,198 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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::searchableSurfacesQueries
Description
A collection of tools for searchableSurfaces.
SourceFiles
searchableSurfacesQueries.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurfacesQueries_H
#define searchableSurfacesQueries_H
#include "searchableSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class plane;
/*---------------------------------------------------------------------------*\
Class searchableSurfacesQueries Declaration
\*---------------------------------------------------------------------------*/
class searchableSurfacesQueries
{
// Private data
// Private Member Functions
//- Temporary wrapper around findNearest. Used in facesIntersection only
static pointIndexHit tempFindNearest
(
const searchableSurface&,
const point& pt,
const scalar initDistSqr
);
//- Calculate sum of distances to nearest point on surfaces. Is used
// in minimisation to find intersection. Returns sum of (square of)
// distances to the surfaces.
static scalar sumDistSqr
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const scalar initialDistSqr, // search box
const point& pt
);
//- Takes the tet (points p) and reflects the point with the
// highest value around the centre (pSum). Checks if it gets closer
// and updates p, y if so.
static scalar tryMorphTet
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const scalar initialDistSqr,
List<vector>& p,
List<scalar>& y,
vector& pSum,
const label ihi,
const scalar fac
);
//- Downhill simplex method: find the point with min cumulative
// distance to all surfaces. Does so by morphing a tet (points p).
// Returns the point on the 0th surface or hit if not reached within
// maxIters iterations.
static bool morphTet
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const label maxIter,
List<vector>& p,
List<scalar>& y
);
//static void findAllIntersections
//(
// const searchableSurface& s,
// const pointField& start,
// const pointField& end,
// const vectorField& smallVec,
// List<List<pointIndexHit> >&
//);
static void mergeHits
(
const point& start,
const scalar mergeDist,
const label surfI,
const List<pointIndexHit>& surfHits,
labelList& allSurfaces,
List<pointIndexHit>& allInfo,
scalarList& allDistSqr
);
public:
// Declare name of the class and its debug switch
ClassName("searchableSurfacesQueries");
// Multiple point queries.
//- Find any intersection. Return hit point information and
// index in surfacesToTest. If multiple surfaces hit the first
// surface is returned, not necessarily the nearest (to start).
static void findAnyIntersection
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const pointField& start,
const pointField& end,
labelList& surfaces,
List<pointIndexHit>&
);
//- Find all intersections in order from start to end. Returns for
// every hit the index in surfacesToTest and the hit info.
static void findAllIntersections
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const pointField& start,
const pointField& end,
labelListList& surfaces,
List<List<pointIndexHit> >& surfaceHits
);
//- Find nearest. Return -1 (and a miss()) or surface and nearest
// point.
static void findNearest
(
const PtrList<searchableSurface>&,
const labelList& surfacesToTest,
const pointField&,
const scalarField& nearestDistSqr,
labelList& surfaces,
List<pointIndexHit>&
);
// Single point queries
//- Calculate point which is on a set of surfaces. WIP.
static pointIndexHit facesIntersection
(
const PtrList<searchableSurface>& allSurfaces,
const labelList& surfacesToTest,
const scalar initDistSqr,
const scalar convergenceDistSqr,
const point& start
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

552
src/meshTools/searchableSurface/triSurfaceMesh.C Normal file
View File

@ -0,0 +1,552 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 "triSurfaceMesh.H"
#include "Random.H"
#include "addToRunTimeSelectionTable.H"
#include "EdgeMap.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(triSurfaceMesh, 0);
addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Check file existence
const Foam::fileName& Foam::triSurfaceMesh::checkFile
(
const fileName& fName,
const fileName& objectName
)
{
if (fName == fileName::null)
{
FatalErrorIn
(
"triSurfaceMesh::checkFile(const fileName&, const fileName&)"
) << "Cannot find triSurfaceMesh starting from "
<< objectName << exit(FatalError);
}
return fName;
}
bool Foam::triSurfaceMesh::isSurfaceClosed() const
{
// Construct pointFaces. Let's hope surface has compact point
// numbering ...
labelListList pointFaces;
invertManyToMany(points().size(), *this, pointFaces);
// Loop over all faces surrounding point. Count edges emanating from point.
// Every edge should be used by two faces exactly.
// To prevent doing work twice per edge only look at edges to higher
// point
EdgeMap<label> facesPerEdge(100);
forAll(pointFaces, pointI)
{
const labelList& pFaces = pointFaces[pointI];
facesPerEdge.clear();
forAll(pFaces, i)
{
const labelledTri& f = triSurface::operator[](pFaces[i]);
label fp = findIndex(f, pointI);
// Forward edge
{
label p1 = f[f.fcIndex(fp)];
if (p1 > pointI)
{
const edge e(pointI, p1);
EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
if (eFnd != facesPerEdge.end())
{
if (eFnd() == 2)
{
return false;
}
eFnd()++;
}
else
{
facesPerEdge.insert(e, 1);
}
}
}
// Reverse edge
{
label p1 = f[f.rcIndex(fp)];
if (p1 > pointI)
{
const edge e(pointI, p1);
EdgeMap<label>::iterator eFnd = facesPerEdge.find(e);
if (eFnd != facesPerEdge.end())
{
if (eFnd() == 2)
{
return false;
}
eFnd()++;
}
else
{
facesPerEdge.insert(e, 1);
}
}
}
}
// Check for any edges used only once.
forAllConstIter(EdgeMap<label>, facesPerEdge, iter)
{
if (iter() != 2)
{
return false;
}
}
}
return true;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
:
searchableSurface(io),
objectRegistry(io),
triSurface(s),
surfaceClosed_(-1)
{}
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
:
searchableSurface(io),
objectRegistry(io),
triSurface
(
checkFile
(
searchableSurface::filePath(),
searchableSurface::objectPath()
)
),
surfaceClosed_(-1)
{}
Foam::triSurfaceMesh::triSurfaceMesh
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
objectRegistry(io),
triSurface
(
checkFile
(
searchableSurface::filePath(),
searchableSurface::objectPath()
)
),
surfaceClosed_(-1)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::triSurfaceMesh::~triSurfaceMesh()
{
clearOut();
}
void Foam::triSurfaceMesh::clearOut()
{
tree_.clear();
edgeTree_.clear();
triSurface::clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
{
tree_.clear();
edgeTree_.clear();
triSurface::movePoints(newPoints);
}
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceMesh::tree() const
{
if (!tree_.valid())
{
treeBoundBox bb(points(), meshPoints());
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
tree_.reset
(
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface(*this),
bb.extend(rndGen, 1E-3), // slightly randomize bb
10, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
return tree_();
}
const Foam::indexedOctree<Foam::treeDataEdge>&
Foam::triSurfaceMesh::edgeTree() const
{
if (!edgeTree_.valid())
{
treeBoundBox bb(localPoints());
// Boundary edges
labelList bEdges
(
identity
(
nEdges()
-nInternalEdges()
)
+ nInternalEdges()
);
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
edgeTree_.reset
(
new indexedOctree<treeDataEdge>
(
treeDataEdge
(
false, // cachebb
edges(), // edges
localPoints(), // points
bEdges // selected edges
),
bb.extend(rndGen, 1E-3), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
return edgeTree_();
}
const Foam::wordList& Foam::triSurfaceMesh::regions() const
{
if (regions_.size() == 0)
{
regions_.setSize(patches().size());
forAll(regions_, regionI)
{
regions_[regionI] = patches()[regionI].name();
}
}
return regions_;
}
//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
//(
// const point& sample,
// const scalar nearestDistSqr
//) const
//{
// return tree().findNearest(sample, nearestDistSqr);
//}
//
//
//Foam::pointIndexHit Foam::triSurfaceMesh::findNearestOnEdge
//(
// const point& sample,
// const scalar nearestDistSqr
//) const
//{
// return = edgeTree().findNearest(sample, nearestDistSqr);
//}
//
//
//Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
//(
// const linePointRef& ln,
// treeBoundBox& tightest,
// point& linePoint
//) const
//{
// return tree().findNearest(ln, tightest, linePoint);
//}
//
//
//Foam::pointIndexHit Foam::triSurfaceMesh::findLine
//(
// const point& start,
// const point& end
//) const
//{
// return tree().findLine(start, end);
//}
//
//
//Foam::pointIndexHit Foam::triSurfaceMesh::findLineAny
//(
// const point& start,
// const point& end
//) const
//{
// return tree().findLineAny(start, end);
//}
// Find out if surface is closed.
bool Foam::triSurfaceMesh::hasVolumeType() const
{
if (surfaceClosed_ == -1)
{
if (isSurfaceClosed())
{
surfaceClosed_ = 1;
}
else
{
surfaceClosed_ = 0;
}
}
return surfaceClosed_ == 1;
}
void Foam::triSurfaceMesh::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(samples.size());
forAll(samples, i)
{
static_cast<pointIndexHit&>(info[i]) =
octree.findNearest(samples[i], nearestDistSqr[i]);
}
}
void Foam::triSurfaceMesh::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) = octree.findLine
(
start[i],
end[i]
);
}
}
void Foam::triSurfaceMesh::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
forAll(start, i)
{
static_cast<pointIndexHit&>(info[i]) =
octree.findLineAny(start[i], end[i]);
}
}
void Foam::triSurfaceMesh::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
const indexedOctree<treeDataTriSurface>& octree = tree();
info.setSize(start.size());
// Work array
DynamicList<pointIndexHit, 1, 1> hits;
// Tolerances
const vectorField dirVec(end-start);
const scalarField magSqrDirVec(magSqr(dirVec));
const vectorField smallVec
(
Foam::sqrt(SMALL)*dirVec
+ vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
);
forAll(start, pointI)
{
hits.clear();
// Current starting point of ray.
point pt = start[pointI];
while (true)
{
// See if any intersection between pt and end
pointIndexHit inter = octree.findLine(pt, end[pointI]);
if (!inter.hit())
{
break;
}
hits.append(inter);
pt = inter.hitPoint() + smallVec[pointI];
if (((pt-start[pointI])&dirVec[pointI]) > magSqrDirVec[pointI])
{
// Adding smallVec has taken us beyond end
break;
}
}
hits.shrink();
info[pointI].transfer(hits);
hits.clear();
}
}
void Foam::triSurfaceMesh::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
region[i] = triSurface::operator[](info[i].index()).region();
}
else
{
region[i] = -1;
}
}
}
void Foam::triSurfaceMesh::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
normal.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
normal[i] = faceNormals()[info[i].index()];
}
else
{
// Set to what?
normal[i] = vector::zero;
}
}
}
void Foam::triSurfaceMesh::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
volType.setSize(points.size());
forAll(points, pointI)
{
const point& pt = points[pointI];
// - use cached volume type per each tree node
// - cheat conversion since same values
volType[pointI] = static_cast<volumeType>(tree().getVolumeType(pt));
}
}
// ************************************************************************* //

103
src/meshTools/triSurface/searchableSurface/triSurfaceMesh.H → src/meshTools/searchableSurface/triSurfaceMesh.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -54,7 +54,7 @@ namespace Foam
class triSurfaceMesh
:
public searchableSurface,
public objectRegistry,
public objectRegistry, // so we can store fields
public triSurface
{
private:
@ -67,6 +67,12 @@ private:
//- Search tree for boundary edges.
mutable autoPtr<indexedOctree<treeDataEdge> > edgeTree_;
//- Names of regions
mutable wordList regions_;
//- Is surface closed
mutable label surfaceClosed_;
// Private Member Functions
//- Check file existence
@ -76,6 +82,10 @@ private:
const fileName& objectName
);
//- Check whether surface is closed without calculating any permanent
// addressing.
bool isSurfaceClosed() const;
//- Disallow default bitwise copy construct
triSurfaceMesh(const triSurfaceMesh&);
@ -97,11 +107,10 @@ public:
//- Construct read
triSurfaceMesh(const IOobject& io);
//- Construct as searchableSurface
//- Construct from dictionary (used by searchableSurface)
triSurfaceMesh
(
const word& name,
const objectRegistry& obj,
const IOobject& io,
const dictionary& dict
);
@ -128,57 +137,73 @@ public:
// searchableSurface implementation
//- Calculate nearest point on surface. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearest
virtual const wordList& regions() const;
//- Whether supports volume type below. I.e. whether is closed.
virtual bool hasVolumeType() const;
// Multiple point queries.
virtual void findNearest
(
const point& sample,
const scalar nearestDistSqr
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
//- Calculate nearest point on edge. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearestOnEdge
virtual void findLine
(
const point& sample,
const scalar nearestDistSqr
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Find nearest to line. Returns
// - bool : any point found?
// - label: relevant index in shapes
// - point: actual nearest point found
// sets:
// - tightest : bounding box
// - linePoint : corresponding nearest point on line
virtual pointIndexHit findNearest
virtual void findLineAny
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Find nearest intersection of line between start and end.
virtual pointIndexHit findLine
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const point& start,
const point& end
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- Find any intersection of line between start and end.
virtual pointIndexHit findLineAny
//- From a set of points and indices get the region
virtual void getRegion
(
const point& start,
const point& end
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual volumeType getVolumeType(const point&) const;
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("triSurfaceMesh::writeData(Ostream&) const");
return false;
}
};

247
src/meshTools/triSurface/searchableSurface/searchableBox.C
View File

@ -1,247 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "searchableBox.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableBox, 0);
addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableBox::searchableBox
(
const word& name,
const treeBoundBox& bb
)
:
searchableSurface(name),
treeBoundBox(bb)
{}
Foam::searchableBox::searchableBox
(
const word& name,
const objectRegistry& obj,
const dictionary& dict
)
:
searchableSurface(name),
treeBoundBox(dict.lookup("min"), dict.lookup("max"))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableBox::~searchableBox()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
// Point can be inside or outside. For every component direction can be
// left of min, right of max or inbetween.
// - outside points: project first one x plane (either min().x()
// or max().x()), then onto y plane and finally z. You should be left
// with intersection point
// - inside point: find nearest side (compare to mid point). Pick any
// one of three points.
const point bbMid(mid());
// Outside point projected onto cube
point interPt(sample);
bool outside = false;
// (for internal points) Per direction what nearest cube side is
point near;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (interPt[dir] < min()[dir])
{
interPt[dir] = min()[dir];
outside = true;
}
else if (interPt[dir] > max()[dir])
{
interPt[dir] = max()[dir];
outside = true;
}
else if (interPt[dir] > bbMid[dir])
{
near[dir] = max()[dir];
}
else
{
near[dir] = min()[dir];
}
}
// For outside points the interPt will be correct now. Handle inside points
// using the three near distances. Project onto the nearest plane.
if (!outside)
{
vector dist(cmptMag(interPt - near));
if (dist.x() < dist.y())
{
if (dist.x() < dist.z())
{
interPt.x() = near.x();
}
else
{
interPt.z() = near.z();
}
}
else
{
if (dist.y() < dist.z())
{
interPt.y() = near.y();
}
else
{
interPt.z() = near.z();
}
}
}
return pointIndexHit(true, interPt, 0);
}
Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const
{
const point bbMid(mid());
point interPt(sample);
for (direction dir = 0; dir < vector::nComponents; dir++)
{
// Project onto left or right depending on mid
if (interPt[dir] > bbMid[dir])
{
interPt[dir] = max()[dir];
}
else
{
interPt[dir] = min()[dir];
}
}
return pointIndexHit(true, interPt, 0);
}
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const
{
notImplemented
(
"searchableBox::findNearest"
"(const linePointRef&, treeBoundBox&, point&)"
);
return pointIndexHit();
}
Foam::pointIndexHit Foam::searchableBox::findLine
(
const point& start,
const point& end
) const
{
point intPt;
bool foundInter = intersects(start, end, intPt);
return pointIndexHit(foundInter, intPt, 0);
}
Foam::pointIndexHit Foam::searchableBox::findLineAny
(
const point& start,
const point& end
) const
{
return findLine(start, end);
}
Foam::searchableSurface::volumeType Foam::searchableBox::getVolumeType
(
const point& pt
) const
{
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
{
return OUTSIDE;
}
}
return INSIDE;
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

161
src/meshTools/triSurface/searchableSurface/searchableBox.H
View File

@ -1,161 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::searchableBox
Description
Searching on bounding box
SourceFiles
searchableBox.C
\*---------------------------------------------------------------------------*/
#ifndef searchableBox_H
#define searchableBox_H
#include "searchableSurface.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableBox Declaration
\*---------------------------------------------------------------------------*/
class searchableBox
:
public searchableSurface,
public treeBoundBox
{
private:
// Private member data
// Private Member Functions
//- Disallow default bitwise copy construct
searchableBox(const searchableBox&);
//- Disallow default bitwise assignment
void operator=(const searchableBox&);
public:
//- Runtime type information
TypeName("searchableBox");
// Constructors
//- Construct from components
searchableBox(const word& name, const treeBoundBox& bb);
searchableBox
(
const word& name,
const objectRegistry& obj,
const dictionary& dict
);
// Destructor
virtual ~searchableBox();
// Member Functions
//- Calculate nearest point on surface. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Calculate nearest point on edge. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const;
//- Find nearest to line. Returns
// - bool : any point found?
// - label: relevant index in shapes
// - point: actual nearest point found
// sets:
// - tightest : bounding box
// - linePoint : corresponding nearest point on line
virtual pointIndexHit findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const;
//- Find nearest intersection of line between start and end.
virtual pointIndexHit findLine
(
const point& start,
const point& end
) const;
//- Find any intersection of line between start and end.
virtual pointIndexHit findLineAny
(
const point& start,
const point& end
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual volumeType getVolumeType(const point&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

257
src/meshTools/triSurface/searchableSurface/searchableSurface.H
View File

@ -1,257 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::searchableSurface
Description
Base class of (analytical or triangulated) surface.
Encapsulates all the search routines.
SourceFiles
searchableSurface.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurface_H
#define searchableSurface_H
#include "pointField.H"
#include "typeInfo.H"
#include "runTimeSelectionTables.H"
#include "pointIndexHit.H"
#include "linePointRef.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class objectRegistry;
class treeBoundBox;
/*---------------------------------------------------------------------------*\
Class searchableSurface Declaration
\*---------------------------------------------------------------------------*/
class searchableSurface
{
public:
// Data types
//- volume types
enum volumeType
{
UNKNOWN = 0,
MIXED = 1,
INSIDE = 2,
OUTSIDE = 3
};
private:
// Private data
const word name_;
// Private Member Functions
//- Disallow default bitwise copy construct
searchableSurface(const searchableSurface&);
//- Disallow default bitwise assignment
void operator=(const searchableSurface&);
public:
//- Runtime type information
TypeName("searchableSurface");
// Declare run-time constructor selection table
// For the dictionary constructor
declareRunTimeSelectionTable
(
autoPtr,
searchableSurface,
dict,
(
const word& name,
const objectRegistry& obj,
const dictionary& dict
),
(name, obj, dict)
);
//// For the Istream constructor
//declareRunTimeSelectionTable
//(
// autoPtr,
// searchableSurface,
// istream,
// (
// const objectRegistry& obj,
// Istream& is
// ),
// (obj, is)
//);
//- Class used for the read-construction of
// PtrLists of searchableSurface
class iNew
{
const objectRegistry& obj_;
public:
iNew(const objectRegistry& obj)
:
obj_(obj)
{}
autoPtr<searchableSurface> operator()(Istream& is) const
{
word surfaceType(is);
word name(is);
dictionary dict(is);
return searchableSurface::New(surfaceType, name, obj_, dict);
}
};
// Constructors
searchableSurface(const word& name);
//- Clone
virtual autoPtr<searchableSurface> clone() const
{
notImplemented("autoPtr<searchableSurface> clone() const");
return autoPtr<searchableSurface>(NULL);
}
// Selectors
//- Return a reference to the selected searchableSurface
static autoPtr<searchableSurface> New
(
const word& surfaceType,
const word& name,
const objectRegistry& obj,
const dictionary& dict
);
////- Return a reference to the selected searchableSurface
//static autoPtr<searchableSurface> New
//(
// const word& surfaceType,
// const objectRegistry& obj,
// Istream& is
//);
// Destructor
virtual ~searchableSurface();
// Member Functions
//- Return name
const word& name() const
{
return name_;
}
//- Calculate nearest point on surface. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const = 0;
//- Calculate nearest point on edge. Returns
// - bool : any point found nearer than nearestDistSqr
// - label: relevant index in surface
// - point: actual nearest point found
virtual pointIndexHit findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const = 0;
//- Find nearest to segment. Returns
// - bool : any point found?
// - label: relevant index in shapes
// - point: actual nearest point found
// sets:
// - tightest : bounding box
// - linePoint : corresponding nearest point on line
virtual pointIndexHit findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const = 0;
//- Find nearest intersection of line between start and end.
virtual pointIndexHit findLine
(
const point& start,
const point& end
) const = 0;
//- Find any intersection of line between start and end.
virtual pointIndexHit findLineAny
(
const point& start,
const point& end
) const = 0;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual volumeType getVolumeType(const point&) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

271
src/meshTools/triSurface/searchableSurface/triSurfaceMesh.C
View File

@ -1,271 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "triSurfaceMesh.H"
#include "Random.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(triSurfaceMesh, 0);
addToRunTimeSelectionTable(searchableSurface, triSurfaceMesh, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
//- Check file existence
const Foam::fileName& Foam::triSurfaceMesh::checkFile
(
const fileName& fName,
const fileName& objectName
)
{
if (fName == fileName::null)
{
FatalErrorIn
(
"triSurfaceMesh::checkFile(const fileName&, const fileName&)"
) << "Cannot find triSurfaceMesh starting from "
<< objectName << exit(FatalError);
}
return fName;
}
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceMesh::tree() const
{
if (!tree_.valid())
{
treeBoundBox bb(points(), meshPoints());
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
tree_.reset
(
new indexedOctree<treeDataTriSurface>
(
treeDataTriSurface(*this),
bb.extend(rndGen, 1E-3), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
return tree_();
}
const Foam::indexedOctree<Foam::treeDataEdge>&
Foam::triSurfaceMesh::edgeTree() const
{
if (!edgeTree_.valid())
{
treeBoundBox bb(localPoints());
// Boundary edges
labelList bEdges
(
identity
(
nEdges()
-nInternalEdges()
)
+ nInternalEdges()
);
// Random number generator. Bit dodgy since not exactly random ;-)
Random rndGen(65431);
edgeTree_.reset
(
new indexedOctree<treeDataEdge>
(
treeDataEdge
(
false, // cachebb
edges(), // edges
localPoints(), // points
bEdges // selected edges
),
bb.extend(rndGen, 1E-3), // slightly randomize bb
8, // maxLevel
10, // leafsize
3.0 // duplicity
)
);
}
return edgeTree_();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
//- Construct from triangles, patches, points.
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io, const triSurface& s)
:
searchableSurface(io.name()),
objectRegistry(io),
triSurface(s)
{}
Foam::triSurfaceMesh::triSurfaceMesh(const IOobject& io)
:
searchableSurface(io.name()),
objectRegistry(io),
triSurface(checkFile(filePath(), objectPath()))
{}
Foam::triSurfaceMesh::triSurfaceMesh
(
const word& name,
const objectRegistry& obj,
const dictionary& dict
)
:
searchableSurface(name),
objectRegistry
(
IOobject
(
name,
"constant",
"triSurface",
obj,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
triSurface(checkFile(filePath(), objectPath()))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::triSurfaceMesh::~triSurfaceMesh()
{}
void Foam::triSurfaceMesh::clearOut()
{
tree_.clear();
edgeTree_.clear();
triSurface::clearOut();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
{
tree_.clear();
edgeTree_.clear();
triSurface::movePoints(newPoints);
}
Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
return tree().findNearest(sample, nearestDistSqr);
}
Foam::pointIndexHit Foam::triSurfaceMesh::findNearestOnEdge
(
const point& sample,
const scalar nearestDistSqr
) const
{
return edgeTree().findNearest(sample, nearestDistSqr);
}
Foam::pointIndexHit Foam::triSurfaceMesh::findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint
) const
{
return tree().findNearest(ln, tightest, linePoint);
}
Foam::pointIndexHit Foam::triSurfaceMesh::findLine
(
const point& start,
const point& end
) const
{
return tree().findLine(start, end);
}
Foam::pointIndexHit Foam::triSurfaceMesh::findLineAny
(
const point& start,
const point& end
) const
{
return tree().findLineAny(start, end);
}
Foam::searchableSurface::volumeType
Foam::triSurfaceMesh::getVolumeType
(
const point& pt
) const
{
// - use cached volume type per each tree node
// - cheat conversion since same values
return static_cast<volumeType>(tree().getVolumeType(pt));
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
// ************************************************************************* //

916
src/meshTools/triSurface/triSurfaceMeshes/triSurfaceMeshes.C
View File

@ -1,916 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "triSurfaceMeshes.H"
#include "Random.H"
#include "Time.H"
#include "SortableList.H"
#include "IOmanip.H"
#include "plane.H"
#include "SortableList.H"
#include "triSurfaceTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(triSurfaceMeshes, 0);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Calculate sum of distance to surfaces.
Foam::scalar Foam::triSurfaceMeshes::sumDistSqr
(
const labelList& surfacesToTest,
const scalar initDistSqr,
const point& pt
) const
{
scalar sum = 0;
forAll(surfacesToTest, i)
{
label surfI = surfacesToTest[i];
pointIndexHit hit(operator[](surfI).findNearest(pt, initDistSqr));
// Note: make it fall over if not hit.
sum += magSqr(hit.hitPoint()-pt);
}
return sum;
}
// Reflects the point furthest away around the triangle centre by a factor fac.
// (triangle centre is the average of all points but the ihi. pSum is running
// sum of all points)
Foam::scalar Foam::triSurfaceMeshes::tryMorphTet
(
const labelList& surfacesToTest,
const scalar initDistSqr,
List<vector>& p,
List<scalar>& y,
vector& pSum,
const label ihi,
const scalar fac
) const
{
scalar fac1 = (1.0-fac)/vector::nComponents;
scalar fac2 = fac1-fac;
vector ptry = pSum*fac1-p[ihi]*fac2;
scalar ytry = sumDistSqr(surfacesToTest, initDistSqr, ptry);
if (ytry < y[ihi])
{
y[ihi] = ytry;
pSum += ptry - p[ihi];
p[ihi] = ptry;
}
return ytry;
}
bool Foam::triSurfaceMeshes::morphTet
(
const labelList& surfacesToTest,
const scalar initDistSqr,
const scalar convergenceDistSqr,
const label maxIter,
List<vector>& p,
List<scalar>& y
) const
{
vector pSum = sum(p);
autoPtr<OFstream> str;
label vertI = 0;
if (debug)
{
Pout<< "triSurfaceMeshes::morphTet : intersection of "
<< IndirectList<fileName>(names(), surfacesToTest)()
<< " starting from points:" << p << endl;
str.reset(new OFstream("track.obj"));
meshTools::writeOBJ(str(), p[0]);
vertI++;
}
for (label iter = 0; iter < maxIter; iter++)
{
// Get the indices of highest, second-highest and lowest values.
label ihi, inhi, ilo;
{
SortableList<scalar> sortedY(y);
ilo = sortedY.indices()[0];
ihi = sortedY.indices()[sortedY.size()-1];
inhi = sortedY.indices()[sortedY.size()-2];
}
if (debug)
{
Pout<< "Iteration:" << iter
<< " lowest:" << y[ilo] << " highest:" << y[ihi]
<< " points:" << p << endl;
meshTools::writeOBJ(str(), p[ilo]);
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
if (y[ihi] < convergenceDistSqr)
{
// Get point on 0th surface.
Swap(p[0], p[ilo]);
Swap(y[0], y[ilo]);
return true;
}
// Reflection: point furthest away gets reflected.
scalar ytry = tryMorphTet
(
surfacesToTest,
10*y[ihi], // search box.
p,
y,
pSum,
ihi,
-1.0
);
if (ytry <= y[ilo])
{
// If in right direction (y lower) expand by two.
ytry = tryMorphTet(surfacesToTest, 10*y[ihi], p, y, pSum, ihi, 2.0);
}
else if (ytry >= y[inhi])
{
// If inside tet try contraction.
scalar ysave = y[ihi];
ytry = tryMorphTet(surfacesToTest, 10*y[ihi], p, y, pSum, ihi, 0.5);
if (ytry >= ysave)
{
// Contract around lowest point.
forAll(p, i)
{
if (i != ilo)
{
p[i] = 0.5*(p[i] + p[ilo]);
y[i] = sumDistSqr(surfacesToTest, y[ihi], p[i]);
}
}
pSum = sum(p);
}
}
}
if (debug)
{
meshTools::writeOBJ(str(), p[0]);
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
// Failure to converge. Return best guess so far.
label ilo = findMin(y);
Swap(p[0], p[ilo]);
Swap(y[0], y[ilo]);
return false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
Foam::triSurfaceMeshes::triSurfaceMeshes
(
const IOobject& io,
const fileNameList& names
)
:
PtrList<triSurfaceMesh>(names.size()),
allSurfaces_(identity(names.size()))
{
forAll(names, i)
{
autoPtr<IOobject> surfaceIO = io.clone();
surfaceIO().rename(names[i]);
Info<< "Loading surface " << surfaceIO().name() << endl;
//fileName fullPath = surfaceIO().filePath();
//
//if (fullPath.size() == 0)
//{
// FatalErrorIn
// (
// "triSurfaceMeshes::triSurfaceMeshes"
// "(const IOobject&, const fileNameList&)"
// ) << "Cannot load surface " << surfaceIO().name()
// << " starting from path " << surfaceIO().path()
// << exit(FatalError);
//}
set(i, new triSurfaceMesh(surfaceIO()));
if (Pstream::master())
{
string oldPrefix(Pout.prefix());
Pout.prefix() += " ";
operator[](i).writeStats(Pout);
Pout.prefix() = oldPrefix;
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::fileNameList Foam::triSurfaceMeshes::allNames(const IOobject& io)
{
return readDir(io.path(), fileName::FILE);
}
Foam::fileNameList Foam::triSurfaceMeshes::names() const
{
fileNameList surfNames(size());
forAll(surfNames, surfI)
{
surfNames[surfI] = operator[](surfI).IOobject::name();
}
return surfNames;
}
// Find any intersection
Foam::label Foam::triSurfaceMeshes::findAnyIntersection
(
const labelList& surfaces,
const point& start,
const point& end,
pointIndexHit& hitInfo
) const
{
forAll(surfaces, i)
{
label surfI = surfaces[i];
hitInfo = operator[](surfI).findLineAny(start, end);
if (hitInfo.hit())
{
return surfI;
}
}
return -1;
}
Foam::label Foam::triSurfaceMeshes::findAnyIntersection
(
const point& start,
const point& end,
pointIndexHit& hitInfo
) const
{
return findAnyIntersection
(
allSurfaces_,
start,
end,
hitInfo
);
}
// Find intersections of edge nearest to both endpoints.
void Foam::triSurfaceMeshes::findAllIntersections
(
const labelList& surfaces,
const point& start,
const point& end,
labelList& surfacesIndex,
List<pointIndexHit>& surfaceHitInfo
) const
{
DynamicList<label> hitSurfaces(surfaces.size());
DynamicList<pointIndexHit> hitInfos(surfaces.size());
DynamicList<scalar> hitDistSqr(surfaces.size());
const vector dirVec(end-start);
const scalar magSqrDirVec(magSqr(dirVec));
const vector smallVec(Foam::sqrt(SMALL)*dirVec);
forAll(surfaces, i)
{
label surfI = surfaces[i];
// Current starting point of ray.
point pt = start;
while (true)
{
// See if any intersection between pt and end
pointIndexHit inter = operator[](surfI).findLine(pt, end);
if (!inter.hit())
{
break;
}
hitSurfaces.append(surfI);
hitInfos.append(inter);
hitDistSqr.append(magSqr(inter.hitPoint() - start));
pt = inter.hitPoint() + smallVec;
if (((pt-start)&dirVec) > magSqrDirVec)
{
// Adding smallVec has taken us beyond end
break;
}
}
}
surfacesIndex.setSize(hitSurfaces.size());
surfaceHitInfo.setSize(hitSurfaces.size());
if (hitSurfaces.size() > 0)
{
// Sort and transfer to arguments
hitSurfaces.shrink();
hitInfos.shrink();
hitDistSqr.shrink();
// Sort from start to end.
SortableList<scalar> sortedDist(hitDistSqr);
forAll(sortedDist.indices(), newPos)
{
label oldPos = sortedDist.indices()[newPos];
surfacesIndex[newPos] = hitSurfaces[oldPos];
surfaceHitInfo[newPos] = hitInfos[oldPos];
}
}
}
void Foam::triSurfaceMeshes::findAllIntersections
(
const point& start,
const point& end,
labelList& surfacesIndex,
List<pointIndexHit>& surfaceHitInfo
) const
{
findAllIntersections
(
allSurfaces_,
start,
end,
surfacesIndex,
surfaceHitInfo
);
}
// Find intersections of edge nearest to both endpoints.
void Foam::triSurfaceMeshes::findNearestIntersection
(
const labelList& surfaces,
const point& start,
const point& end,
label& surface1,
pointIndexHit& hit1,
label& surface2,
pointIndexHit& hit2
) const
{
surface1 = -1;
// Initialize to endpoint
hit1 = pointIndexHit(false, end, -1);
forAll(surfaces, i)
{
label surfI = surfaces[i];
if (hit1.rawPoint() == start)
{
break;
}
// See if any intersection between start and current nearest
pointIndexHit inter = operator[](surfI).findLine
(
start,
hit1.rawPoint()
);
if (inter.hit())
{
hit1 = inter;
surface1 = surfI;
}
}
// Find the nearest intersection from end to start. Note that we initialize
// to the first intersection (if any).
surface2 = surface1;
hit2 = pointIndexHit(hit1);
if (hit1.hit())
{
// Test from the end side.
forAll(surfaces, i)
{
label surfI = surfaces[i];
if (hit2.rawPoint() == end)
{
break;
}
// See if any intersection between end and current nearest
pointIndexHit inter = operator[](surfI).findLine
(
end,
hit2.rawPoint()
);
if (inter.hit())
{
hit2 = inter;
surface2 = surfI;
}
}
}
}
void Foam::triSurfaceMeshes::findNearestIntersection
(
const point& start,
const point& end,
label& surface1,
pointIndexHit& hit1,
label& surface2,
pointIndexHit& hit2
) const
{
findNearestIntersection
(
allSurfaces_,
start,
end,
surface1,
hit1,
surface2,
hit2
);
}
// Find nearest. Return -1 or nearest point
Foam::label Foam::triSurfaceMeshes::findNearest
(
const labelList& surfaces,
const point& pt,
const scalar nearestDistSqr,
pointIndexHit& nearestHit
) const
{
// nearest surface
label minSurface = -1;
scalar minDistSqr = Foam::sqr(GREAT);
forAll(surfaces, i)
{
label surfI = surfaces[i];
pointIndexHit hit
(
operator[](surfI).findNearest(pt, nearestDistSqr)
);
if (hit.hit())
{
scalar distSqr = magSqr(hit.hitPoint()-pt);
if (distSqr < minDistSqr)
{
minDistSqr = distSqr;
minSurface = surfI;
nearestHit = hit;
}
}
}
if (minSurface == -1)
{
// maxLevel unchanged. No interesting surface hit.
nearestHit.setMiss();
}
return minSurface;
}
// Find nearest. Return -1 or nearest point
Foam::label Foam::triSurfaceMeshes::findNearest
(
const point& pt,
const scalar nearestDistSqr,
pointIndexHit& nearestHit
) const
{
return findNearest
(
allSurfaces_,
pt,
nearestDistSqr,
nearestHit
);
}
Foam::label Foam::triSurfaceMeshes::findNearestAndClassify
(
const labelList& surfacesToTest,
const point& pt,
const scalar nearestDistSqr,
surfaceLocation& nearest
) const
{
pointIndexHit nearestInfo;
label surfI = findNearest
(
surfacesToTest,
pt,
nearestDistSqr,
nearestInfo
);
if (surfI != -1)
{
nearest = triSurfaceTools::classify
(
operator[](surfI),
nearestInfo.index(), // triangle
nearestInfo.rawPoint() // point on edge/inside triangle
);
}
return surfI;
}
//- Calculate point which is on a set of surfaces.
Foam::surfaceLocation Foam::triSurfaceMeshes::facesIntersectionTrack
(
const labelList& surfacesToTest,
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const point& start
) const
{
autoPtr<OFstream> str;
label vertI = 0;
if (debug)
{
str.reset(new OFstream("track.obj"));
}
surfaceLocation current
(
pointIndexHit
(
false,
start,
-1
),
triPointRef::NONE,
-1
);
// Dump start point
if (str.valid())
{
Pout<< "Starting at " << current.info()
<< " written as point " << vertI
<< endl;
meshTools::writeOBJ(str(), current.rawPoint());
vertI++;
}
// Now slide current along all faces until it settles
label iter = 0;
const label maxIter = 10;
for (; iter < maxIter; iter++)
{
// - store old position
// - slide it along all faces
// - if it hasn't changed position exit loop
if (debug)
{
Pout<< endl
<< "Iteration " << iter << endl
<< "-----------" << endl;
}
point oldCurrent = current.rawPoint();
forAll(surfacesToTest, i)
{
label surfI = surfacesToTest[i];
// Project pt onto surf
// ~~~~~~~~~~~~~~~~~~~~
point copy(current.rawPoint()); // need to work on copy
if
(
findNearestAndClassify
(
labelList(1, surfI),
copy,
initialDistSqr, // initialDistSqr
current
)
== -1
)
{
FatalErrorIn("triSurfaceMeshes::facesIntersectionTrack(..)")
<< "Did not find a point on surface " << surfI
<< " which is within sqrt(" << initialDistSqr
<< ") of " << copy
<< abort(FatalError);
}
if (debug)
{
Pout<< "Nearest onto surface " << surfI
<< ' ' << operator[](surfI).IOobject::name() << " : "
<< current.info() << " written as point " << vertI
<< endl;
}
// Dump current
if (str.valid())
{
meshTools::writeOBJ(str(), current.rawPoint());
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
// Find corresponding point on next surface
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
label nextSurfI = surfacesToTest[surfacesToTest.fcIndex(i)];
surfaceLocation next;
findNearestAndClassify
(
labelList(1, nextSurfI),
point(current.rawPoint()),
initialDistSqr, // initialDistSqr
next
);
// Since next is on a different surface we cannot compare
// indices (like in snapToEnd) so make sure this does not
// happen.
next.elementType() = triPointRef::NONE;
next.setIndex(-123);
if (debug)
{
Pout<< "Nearest onto next surface " << nextSurfI
<< ' ' << operator[](nextSurfI).IOobject::name() << " : "
<< next.info() << endl;
}
if
(
magSqr(current.rawPoint() - next.rawPoint())
> convergenceDistSqr
)
{
// Track from current to next
// ~~~~~~~~~~~~~~~~~~~~~~~~~~
//vector n = current.normal(surfaces[surfI]);
current = triSurfaceTools::trackToEdge
(
operator[](surfI),
current,
next,
plane(start, current.rawPoint(), next.rawPoint())
);
if (debug)
{
Pout<< "Sliding along surface "
<< surfI << ' ' << operator[](surfI).IOobject::name()
<< " in direction of "
<< nextSurfI << ' '
<< operator[](nextSurfI).IOobject::name()
<< " stopped at " << current.info()
<< " written as point " << vertI << endl;
}
if (str.valid())
{
meshTools::writeOBJ(str(), current.rawPoint());
vertI++;
str()<< "l " << vertI-1 << ' ' << vertI << nl;
}
}
}
scalar distSqr = magSqr(oldCurrent - current.rawPoint());
if (debug)
{
Pout<< "distSqr:" << distSqr
<< " convergenceDistSqr:" << convergenceDistSqr << endl;
}
if (distSqr < convergenceDistSqr)
{
break;
}
}
if (iter == maxIter)
{
FatalErrorIn("triSurfaceMeshes::facesIntersectionTrack(..)")
<< "Did not converge in " << iter
<< " iterations to find a point which is on surfaces "
<< IndirectList<fileName>(names(), surfacesToTest)() << nl
<< "Start point:" << start << nl
<< "Current nearest:" << current.info() << nl
<< "Please check that your surfaces actually intersect and"
<< " that the starting point is close to the intersection point."
<< abort(FatalError);
}
return current;
}
//- Calculate point which is on a set of surfaces.
Foam::pointIndexHit Foam::triSurfaceMeshes::facesIntersection
(
const labelList& surfacesToTest,
const scalar initDistSqr,
const scalar convergenceDistSqr,
const point& start
) const
{
// Get four starting points. Take these as the projection of the
// starting point onto the surfaces and the mid point
List<point> nearest(surfacesToTest.size()+1);
point sumNearest = vector::zero;
forAll(surfacesToTest, i)
{
label surfI = surfacesToTest[i];
pointIndexHit hit
(
operator[](surfI).findNearest(start, initDistSqr)
);
if (hit.hit())
{
nearest[i] = hit.hitPoint();
sumNearest += nearest[i];
}
else
{
FatalErrorIn
(
"triSurfaceMeshes::facesIntersection"
"(const labelList&, const scalar, const scalar, const point&)"
) << "Did not find point within distance "
<< initDistSqr << " of starting point " << start
<< " on surface " << operator[](surfI).IOobject::name()
<< abort(FatalError);
}
}
nearest[nearest.size()-1] = sumNearest / surfacesToTest.size();
// Get the sum of distances (initial evaluation)
List<scalar> nearestDist(nearest.size());
forAll(nearestDist, i)
{
nearestDist[i] = sumDistSqr(surfacesToTest, initDistSqr, nearest[i]);
}
//- Downhill Simplex method
bool converged = morphTet
(
surfacesToTest,
initDistSqr,
convergenceDistSqr,
2000,
nearest,
nearestDist
);
pointIndexHit intersection;
if (converged)
{
// Project nearest onto 0th surface.
intersection = operator[](surfacesToTest[0]).findNearest
(
nearest[0],
nearestDist[0]
);
}
//if (!intersection.hit())
//{
// // Restart
// scalar smallDist = Foam::sqr(convergenceDistSqr);
// nearest[0] = intersection.hitPoint();
// nearest[1] = nearest[0];
// nearest[1].x() += smallDist;
// nearest[2] = nearest[0];
// nearest[2].y() += smallDist;
// nearest[3] = nearest[0];
// nearest[3].z() += smallDist;
//
// forAll(nearestDist, i)
// {
// nearestDist[i] = sumDistSqr
// (
// surfacesToTest,
// initDistSqr,
// nearest[i]
// );
// }
//
// intersection = morphTet
// (
// surfacesToTest,
// initDistSqr,
// convergenceDistSqr,
// 1000,
// nearest,
// nearestDist
// );
//}
return intersection;
}
// ************************************************************************* //

261
src/meshTools/triSurface/triSurfaceMeshes/triSurfaceMeshes.H
View File

@ -1,261 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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::triSurfaceMeshes
Description
Container for triSurfaces read from files.
SourceFiles
triSurfaceMeshes.C
\*---------------------------------------------------------------------------*/
#ifndef triSurfaceMeshes_H
#define triSurfaceMeshes_H
#include "triSurfaceMesh.H"
#include "fileNameList.H"
#include "treeDataTriSurface.H"
#include "indexedOctree.H"
#include "surfaceLocation.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
class plane;
/*---------------------------------------------------------------------------*\
Class triSurfaceMeshes Declaration
\*---------------------------------------------------------------------------*/
class triSurfaceMeshes
:
public PtrList<triSurfaceMesh>
{
// Private data
//- Indices of all surfaces. Precalculated and stored.
const labelList allSurfaces_;
// Private Member Functions
//- Calculate sum of distances to nearest point on surfaces. Is used
// in minimisation to find intersection. Returns sum of (square of)
// distances to the surfaces.
scalar sumDistSqr
(
const labelList& surfacesToTest,
const scalar initialDistSqr, // search box
const point& pt
) const;
//- Takes the tet (points p) and reflects the point with the
// highest value around the centre (pSum). Checks if it gets closer
// and updates p, y if so.
scalar tryMorphTet
(
const labelList& surfacesToTest,
const scalar initialDistSqr,
List<vector>& p,
List<scalar>& y,
vector& pSum,
const label ihi,
const scalar fac
) const;
//- Downhill simplex method: find the point with min cumulative
// distance to all surfaces. Does so by morphing a tet (points p).
// Returns the point on the 0th surface or hit if not reached within
// maxIters iterations.
bool morphTet
(
const labelList& surfacesToTest,
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const label maxIter,
List<vector>& p,
List<scalar>& y
) const;
//- Disallow default bitwise copy construct
triSurfaceMeshes(const triSurfaceMeshes&);
//- Disallow default bitwise assignment
void operator=(const triSurfaceMeshes&);
public:
ClassName("triSurfaceMeshes");
// Constructors
//- Construct from directory (io.instance()) and list of local filenames
triSurfaceMeshes(const IOobject& io, const fileNameList&);
// Member Functions
//- Get all surfaces in directory
static fileNameList allNames(const IOobject&);
//- Get names of surfaces
fileNameList names() const;
//- Find any intersection. Return hit point information and surface
// number
label findAnyIntersection
(
const labelList& surfacesToTest,
const point& start,
const point& end,
pointIndexHit&
) const;
label findAnyIntersection
(
const point& start,
const point& end,
pointIndexHit&
) const;
//- Find all intersections in order from start to end. Returns for
// every hit the surface and the hit info.
void findAllIntersections
(
const labelList& surfacesToTest,
const point& start,
const point& end,
labelList& surfaces,
List<pointIndexHit>& surfaceHits
) const;
void findAllIntersections
(
const point& start,
const point& end,
labelList& surfaces,
List<pointIndexHit>& surfaceHits
) const;
//- Find intersections of edge nearest to both endpoints.
void findNearestIntersection
(
const labelList& surfacesToTest,
const point& start,
const point& end,
label& surface1, // surface index
pointIndexHit& hit1, // hit point information
label& surface2,
pointIndexHit& hit2
) const;
void findNearestIntersection
(
const point& start,
const point& end,
label& surface1, // surface index
pointIndexHit& hit1, // hit point information
label& surface2,
pointIndexHit& hit2
) const;
//- Find nearest. Return -1 (and a miss()) or surface and nearest point.
label findNearest
(
const labelList& surfacesToTest,
const point&,
const scalar nearestDistSqr,
pointIndexHit&
) const;
label findNearest
(
const point&,
const scalar nearestDistSqr,
pointIndexHit&
) const;
//- Find nearest point (like findNearest) but also return whether
// point is on edge or on point. Returns -1 if not found
// (within nearestDistSqr) on any of the surfaces.
// Returns surfaceLocation
// - hit : nearest point is inside triangle
// - elementType : none, edge or point
// - index : triI, edgeI or localPointI
label findNearestAndClassify
(
const labelList& surfacesToTest,
const point& pt,
const scalar nearestDistSqr,
surfaceLocation& nearest
) const;
//- Calculate point which is on a set of surfaces. Takes
// - initialDistSqr : search dimensions to find point on surface
// - convergenceDistSqr : when point is converged
// Will bomb out if no stable point reached after certain number
// of iterations.
surfaceLocation facesIntersectionTrack
(
const labelList& surfacesToTest,
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const point& start
) const;
//- Calculate point which is on a set of surfaces.
pointIndexHit facesIntersection
(
const labelList& surfacesToTest,
const scalar initialDistSqr,
const scalar convergenceDistSqr,
const point& start
) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //