ENH: snappyHexMesh: implicit feature edge snapping

This commit is contained in:
mattijs
2012-09-10 17:36:29 +01:00
parent d996f4218a
commit 5fb4fb855c
8 changed files with 1448 additions and 529 deletions

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -663,18 +663,11 @@ void Foam::autoRefineDriver::mergePatchFaces
const dictionary& motionDict
)
{
const fvMesh& mesh = meshRefiner_.mesh();
Info<< nl
<< "Merge refined boundary faces" << nl
<< "----------------------------" << nl
<< endl;
if (debug)
{
const_cast<Time&>(mesh.time())++;
}
meshRefiner_.mergePatchFacesUndo
(
Foam::cos(degToRad(45.0)),

View File

@ -1465,6 +1465,7 @@ void Foam::autoSnapDriver::doSnap
{
disp = calcNearestSurfaceFeature
(
snapParams,
iter,
featureCos,
scalar(iter+1)/nFeatIter,

View File

@ -111,17 +111,43 @@ class autoSnapDriver
// Feature line snapping
//- Is point on two feature edges that make a largish angle?
bool isFeaturePoint
(
const scalar featureCos,
const indirectPrimitivePatch& pp,
const PackedBoolList& isFeatureEdge,
const label pointI
) const;
void smoothAndConstrain
(
const indirectPrimitivePatch& pp,
const List<pointConstraint>& constraints,
vectorField& disp
) const;
//void calcNearest
void smoothAndConstrain2
(
const bool applyConstraints,
const indirectPrimitivePatch& pp,
const List<pointConstraint>& constraints,
vectorField& disp
) const;
void calcNearest
(
const label iter,
const indirectPrimitivePatch& pp,
vectorField& pointDisp,
vectorField& pointSurfaceNormal,
vectorField& pointRotation
) const;
//void calcNearestFace
//(
// const pointField& points,
// vectorField& disp,
// vectorField& surfaceNormal
// const label iter,
// const indirectPrimitivePatch& pp,
// vectorField& faceDisp,
// vectorField& faceSurfaceNormal,
// vectorField& faceRotation
//) const;
void calcNearestFace
(
@ -129,19 +155,18 @@ class autoSnapDriver
const indirectPrimitivePatch& pp,
vectorField& faceDisp,
vectorField& faceSurfaceNormal,
labelList& faceSurfaceRegion,
vectorField& faceRotation
) const;
void interpolateFaceToPoint
(
const label iter,
const indirectPrimitivePatch& pp,
const vectorField& faceSurfaceNormal,
const vectorField& faceDisp,
const vectorField& faceRotation,
vectorField& patchDisp,
vectorField& patchRotationDisp
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceRotation,
const List<List<point> >& pointFaceCentres,
vectorField& patchDisp
//vectorField& patchRotationDisp
) const;
void correctAttraction
(
@ -152,6 +177,14 @@ class autoSnapDriver
const point& pt,
vector& edgeOffset // offset from pt to point on edge
) const;
//- Return hit if on multiple points
pointIndexHit findMultiPatchPoint
(
const point& pt,
const labelList& patchIDs,
const List<point>& faceCentres
) const;
void binFeatureFace
(
const label iter,
@ -186,6 +219,25 @@ class autoSnapDriver
DynamicList<label>& surfaceCount
) const;
void featureAttractionUsingReconstruction
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const List<List<point> >& pointFaceSurfNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
const labelListList& pointFacePatchID,
vector& patchAttraction,
pointConstraint& patchConstraint
) const;
void featureAttractionUsingReconstruction
(
const label iter,
@ -196,31 +248,12 @@ class autoSnapDriver
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
const labelListList& pointFacePatchID,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void determineAllFeatures
(
const label iter,
const scalar featureCos,
const indirectPrimitivePatch&,
const scalarField&,
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
List<labelList>& pointAttractor,
List<List<pointConstraint> >& pointConstraints,
// Feature-edge to pp point
List<List<DynamicList<point> > >& edgeAttractors,
List<List<DynamicList<pointConstraint> > >& edgeConstraints,
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void determineFeatures
(
const label iter,
@ -232,6 +265,7 @@ class autoSnapDriver
const List<List<point> >& pointFaceNormals,
const List<List<point> >& pointFaceDisp,
const List<List<point> >& pointFaceCentres,
const labelListList& pointFacePatchID,
List<labelList>& pointAttractor,
List<List<pointConstraint> >& pointConstraints,
@ -242,6 +276,47 @@ class autoSnapDriver
List<pointConstraint>& patchConstraints
) const;
//- Find point on nearest feature edge (within searchDist).
// Return point and feature
// and store feature-edge to mesh-point and vice versa
pointIndexHit findNearFeatureEdge
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const point& estimatedPt,
label& featI,
List<List<DynamicList<point> > >&,
List<List<DynamicList<pointConstraint> > >&,
vectorField&,
List<pointConstraint>&
) const;
//- Find nearest feature point (within searchDist).
// Return feature point
// and store feature-point to mesh-point and vice versa.
// If another mesh point already referring to this feature
// point and further away, reset that one to a near feature
// edge (using findNearFeatureEdge above)
labelPair findNearFeaturePoint
(
const indirectPrimitivePatch& pp,
const scalarField& snapDist,
const label pointI,
const point& estimatedPt,
// Feature-point to pp point
List<labelList>& pointAttractor,
List<List<pointConstraint> >& pointConstraints,
// Feature-edge to pp point
List<List<DynamicList<point> > >& edgeAttractors,
List<List<DynamicList<pointConstraint> > >& edgeConstraints,
// pp point to nearest feature
vectorField& patchAttraction,
List<pointConstraint>& patchConstraints
) const;
void featureAttractionUsingFeatureEdges
(
const label iter,
@ -270,6 +345,7 @@ class autoSnapDriver
vectorField calcNearestSurfaceFeature
(
const snapParameters& snapParams,
const label iter,
const scalar featureCos,
const scalar featureAttract,
@ -342,7 +418,6 @@ public:
motionSmoother& meshMover
) const;
//- Smooth the displacement field to the internal.
void smoothDisplacement
(

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,7 +34,9 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
snapTol_(readScalar(dict.lookup("tolerance"))),
nSmoothDispl_(readLabel(dict.lookup("nSolveIter"))),
nSnap_(readLabel(dict.lookup("nRelaxIter"))),
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1))
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false))
{}

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -37,6 +37,7 @@ SourceFiles
#include "dictionary.H"
#include "scalar.H"
#include "Switch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,6 +64,10 @@ class snapParameters
const label nFeatureSnap_;
const Switch explicitFeatureSnap_;
const Switch implicitFeatureSnap_;
// Private Member Functions
@ -119,6 +124,16 @@ public:
return nFeatureSnap_;
}
Switch explicitFeatureSnap() const
{
return explicitFeatureSnap_;
}
Switch implicitFeatureSnap() const
{
return implicitFeatureSnap_;
}
};

View File

@ -326,6 +326,16 @@ private:
const labelList& globalToPatch
) const;
//- Geometric test on see whether face needs to be baffled:
// is face boundary face and perpendicular to surface normal?
bool validBaffleTopology
(
const label faceI,
const vector& n1,
const vector& n2,
const vector& testDir
) const;
//- Determine patches for baffles
void getBafflePatches
(
@ -414,7 +424,7 @@ private:
//- Extract those baffles (duplicate) faces that are on the edge
// of a baffle region. These are candidates for merging.
List<labelPair> filterDuplicateFaces(const List<labelPair>&) const;
List<labelPair> freeStandingBaffles(const List<labelPair>&) const;
// Zone handling

View File

@ -43,6 +43,7 @@ License
#include "OFstream.H"
#include "regionSplit.H"
#include "removeCells.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -216,6 +217,45 @@ Foam::label Foam::meshRefinement::getBafflePatch
}
// Check if we are a boundary face and normal of surface does
// not align with test vector. In this case there'd probably be
// a freestanding 'baffle' so we might as well not create it.
// Note that since it is not a proper baffle we cannot detect it
// afterwards so this code cannot be merged with the
// filterDuplicateFaces code.
bool Foam::meshRefinement::validBaffleTopology
(
const label faceI,
const vector& n1,
const vector& n2,
const vector& testDir
) const
{
label patchI = mesh_.boundaryMesh().whichPatch(faceI);
if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
{
return true;
}
else if (mag(n1&n2) > cos(degToRad(30)))
{
// Both normals aligned. Check that test vector perpendicularish to
// surface normal
scalar magTestDir = mag(testDir);
if (magTestDir > VSMALL)
{
if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
{
//Pout<< "** disabling baffling face "
// << mesh_.faceCentres()[faceI] << endl;
return false;
}
}
}
return true;
}
// Determine patches for baffles on all intersected unnamed faces
void Foam::meshRefinement::getBafflePatches
(
@ -745,7 +785,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
// Done by counting the number of baffles faces per mesh edge. If edge
// has 2 boundary faces and both are baffle faces it is the edge of a baffle
// region.
Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
(
const List<labelPair>& couples
) const
@ -852,12 +892,112 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
}
filteredCouples.setSize(filterI);
//Info<< "filterDuplicateFaces : from "
//Info<< "freeStandingBaffles : from "
// << returnReduce(couples.size(), sumOp<label>())
// << " down to "
// << returnReduce(filteredCouples.size(), sumOp<label>())
// << " baffles." << nl << endl;
//XXXXXX
// {
// // Collect segments
// // ~~~~~~~~~~~~~~~~
//
// pointField start(filteredCouples.size());
// pointField end(filteredCouples.size());
//
// const pointField& cellCentres = mesh_.cellCentres();
//
// forAll(filteredCouples, i)
// {
// const labelPair& couple = couples[i];
// start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
// end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
// }
//
// // Extend segments a bit
// {
// const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
// start -= smallVec;
// end += smallVec;
// }
//
//
// // Do test for intersections
// // ~~~~~~~~~~~~~~~~~~~~~~~~~
// labelList surface1;
// List<pointIndexHit> hit1;
// labelList region1;
// vectorField normal1;
//
// labelList surface2;
// List<pointIndexHit> hit2;
// labelList region2;
// vectorField normal2;
//
// surfaces_.findNearestIntersection
// (
// surfacesToBaffle,
// start,
// end,
//
// surface1,
// hit1,
// region1,
// normal1,
//
// surface2,
// hit2,
// region2,
// normal2
// );
//
// forAll(testFaces, i)
// {
// if (hit1[i].hit() && hit2[i].hit())
// {
// bool createBaffle = true;
//
// label faceI = couples[i].first();
// label patchI = mesh_.boundaryMesh().whichPatch(faceI);
// if (patchI != -1 && !mesh_.boundaryMesh()[patchI].coupled())
// {
// // Check if we are a boundary face and normal of surface
// // does
// // not align with test vector. In this case there'd
// // probably be
// // a freestanding 'baffle' so we might as well not
// // create it.
// // Note that since it is not a proper baffle we cannot
// // detect it
// // afterwards so this code cannot be merged with the
// // filterDuplicateFaces code.
// if (mag(normal1[i]&normal2[i]) > cos(degToRad(30)))
// {
// // Both normals aligned
// vector n = end[i]-start[i];
// scalar magN = mag(n);
// if (magN > VSMALL)
// {
// n /= magN;
//
// if (mag(normal1[i]&n) < cos(degToRad(45)))
// {
// Pout<< "** disabling baffling face "
// << mesh_.faceCentres()[faceI] << endl;
// createBaffle = false;
// }
// }
// }
// }
//
//
// }
//XXXXXX
return filteredCouples;
}
@ -1717,11 +1857,6 @@ void Foam::meshRefinement::baffleAndSplitMesh
neiPatch
);
if (debug)
{
runTime++;
}
createBaffles(ownPatch, neiPatch);
if (debug)
@ -1874,15 +2009,11 @@ void Foam::meshRefinement::baffleAndSplitMesh
<< "---------------------------" << nl
<< endl;
if (debug)
{
runTime++;
}
// List of pairs of freestanding baffle faces.
List<labelPair> couples
(
filterDuplicateFaces // filter out freestanding baffles
freeStandingBaffles // filter out freestanding baffles
(
getDuplicateFaces // get all baffles
(
@ -2484,11 +2615,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
// the information already in surfaceIndex_.
labelList surface1;
List<pointIndexHit> hit1;
vectorField normal1;
labelList surface2;
List<pointIndexHit> hit2;
vectorField normal2;
{
List<pointIndexHit> hit1;
labelList region1;
List<pointIndexHit> hit2;
labelList region2;
surfaces_.findNearestIntersection
(
@ -2511,9 +2644,36 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (surface1[i] != -1)
{
// If both hit should probably choose nearest. For later.
namedSurfaceIndex[faceI] = surface1[i];
nSurfFaces[surface1[i]]++;
//- Not allowed not to create baffle - is vital for regioning.
// Have logic instead at erosion!
//bool createBaffle = validBaffleTopology
//(
// faceI,
// normal1[i],
// normal2[i],
// end[i]-start[i]
//);
//
// If both hit should probably choose 'nearest'
if
(
surface2[i] != -1
&& (
magSqr(hit2[i].hitPoint())
< magSqr(hit1[i].hitPoint())
)
)
{
namedSurfaceIndex[faceI] = surface2[i];
nSurfFaces[surface2[i]]++;
}
else
{
namedSurfaceIndex[faceI] = surface1[i];
nSurfFaces[surface1[i]]++;
}
}
else if (surface2[i] != -1)
{