ENH: snappyHexMesh: enable free-standing baffle merging

This commit is contained in:
mattijs
2013-05-16 12:49:16 +01:00
parent 5f68e587e6
commit 25a373bd39
7 changed files with 445 additions and 260 deletions

View File

@ -180,6 +180,12 @@ castellatedMeshControls
// - used if feature snapping (see snapControls below) is used // - used if feature snapping (see snapControls below) is used
resolveFeatureAngle 30; resolveFeatureAngle 30;
// Planar angle:
// - used to determine if neighbouring surface normals
// are roughly the same so e.g. free-standing baffles can be merged
// If not specified same as resolveFeatureAngle
//planarAngle 30;
// Region-wise refinement // Region-wise refinement
// ~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~

View File

@ -522,6 +522,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
false, // perpendicular edge connected cells false, // perpendicular edge connected cells
scalarField(0), // per region perpendicular angle scalarField(0), // per region perpendicular angle
!handleSnapProblems, // merge free standing baffles? !handleSnapProblems, // merge free standing baffles?
refineParams.planarAngle(),
motionDict, motionDict,
const_cast<Time&>(mesh.time()), const_cast<Time&>(mesh.time()),
globalToMasterPatch_, globalToMasterPatch_,
@ -606,8 +607,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
handleSnapProblems, handleSnapProblems,
handleSnapProblems, // remove perp edge connected cells handleSnapProblems, // remove perp edge connected cells
perpAngle, // perp angle perpAngle, // perp angle
false, // merge free standing baffles? true, // merge free standing baffles?
//true, // merge free standing baffles? refineParams.planarAngle(), // planar angle
motionDict, motionDict,
const_cast<Time&>(mesh.time()), const_cast<Time&>(mesh.time()),
globalToMasterPatch_, globalToMasterPatch_,

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -41,6 +41,7 @@ Foam::refinementParameters::refinementParameters
maxLocalCells_(readLabel(dict.lookup("procCellLimit"))), maxLocalCells_(readLabel(dict.lookup("procCellLimit"))),
minRefineCells_(readLabel(dict.lookup("minimumRefine"))), minRefineCells_(readLabel(dict.lookup("minimumRefine"))),
curvature_(readScalar(dict.lookup("curvature"))), curvature_(readScalar(dict.lookup("curvature"))),
planarAngle_(dict.lookupOrDefault("planarAngle", curvature_)),
nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))), nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))),
keepPoints_(dict.lookup("keepPoints")), keepPoints_(dict.lookup("keepPoints")),
allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
@ -53,6 +54,14 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
maxGlobalCells_(readLabel(dict.lookup("maxGlobalCells"))), maxGlobalCells_(readLabel(dict.lookup("maxGlobalCells"))),
maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))), maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))),
minRefineCells_(readLabel(dict.lookup("minRefinementCells"))), minRefineCells_(readLabel(dict.lookup("minRefinementCells"))),
planarAngle_
(
dict.lookupOrDefault
(
"planarAngle",
readScalar(dict.lookup("resolveFeatureAngle"))
)
),
nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))), nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))),
keepPoints_(pointField(1, dict.lookup("locationInMesh"))), keepPoints_(pointField(1, dict.lookup("locationInMesh"))),
allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),

View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -67,6 +67,9 @@ class refinementParameters
//- Curvature //- Curvature
scalar curvature_; scalar curvature_;
//- Planarity criterion
scalar planarAngle_;
//- Number of layers between different refinement levels //- Number of layers between different refinement levels
const label nBufferLayers_; const label nBufferLayers_;
@ -129,6 +132,12 @@ public:
return curvature_; return curvature_;
} }
//- Angle when two intersections are considered to be planar
scalar planarAngle() const
{
return planarAngle_;
}
//- Number of layers between different refinement levels //- Number of layers between different refinement levels
label nBufferLayers() const label nBufferLayers() const
{ {

View File

@ -253,6 +253,14 @@ private:
label& nRefine label& nRefine
); );
//- Mark every cell with level of feature passing through it
// (or -1 if not passed through). Uses tracking.
void markFeatureCellLevel
(
const point& keepPoint,
labelList& maxFeatureLevel
) const;
//- Calculate list of cells to refine based on intersection of //- Calculate list of cells to refine based on intersection of
// features. // features.
label markFeatureRefinement label markFeatureRefinement
@ -337,16 +345,6 @@ private:
const labelList& globalToSlavePatch const labelList& globalToSlavePatch
) const; ) 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 //- Determine patches for baffles
void getBafflePatches void getBafflePatches
( (
@ -435,7 +433,11 @@ private:
//- Extract those baffles (duplicate) faces that are on the edge //- Extract those baffles (duplicate) faces that are on the edge
// of a baffle region. These are candidates for merging. // of a baffle region. These are candidates for merging.
List<labelPair> freeStandingBaffles(const List<labelPair>&) const; List<labelPair> freeStandingBaffles
(
const List<labelPair>&,
const scalar freeStandingAngle
) const;
// Zone handling // Zone handling
@ -489,6 +491,16 @@ private:
labelList& namedSurfaceIndex labelList& namedSurfaceIndex
) const; ) const;
//- Remove any loose standing cells
void handleSnapProblems
(
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch
);
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
meshRefinement(const meshRefinement&); meshRefinement(const meshRefinement&);
@ -707,6 +719,7 @@ public:
const bool removeEdgeConnectedCells, const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle, const scalarField& perpendicularAngle,
const bool mergeFreeStanding, const bool mergeFreeStanding,
const scalar freeStandingAngle,
const dictionary& motionDict, const dictionary& motionDict,
Time& runTime, Time& runTime,
const labelList& globalToMasterPatch, const labelList& globalToMasterPatch,

View File

@ -44,6 +44,7 @@ License
#include "regionSplit.H" #include "regionSplit.H"
#include "removeCells.H" #include "removeCells.H"
#include "unitConversion.H" #include "unitConversion.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -217,43 +218,43 @@ Foam::label Foam::meshRefinement::getBafflePatch
} }
// Check if we are a boundary face and normal of surface does //// Check if we are a boundary face and normal of surface does
// not align with test vector. In this case there'd probably be //// not align with test vector. In this case there'd probably be
// a freestanding 'baffle' so we might as well not create it. //// a freestanding 'baffle' so we might as well not create it.
// Note that since it is not a proper baffle we cannot detect it //// Note that since it is not a proper baffle we cannot detect it
// afterwards so this code cannot be merged with the //// afterwards so this code cannot be merged with the
// filterDuplicateFaces code. //// filterDuplicateFaces code.
bool Foam::meshRefinement::validBaffleTopology //bool Foam::meshRefinement::validBaffleTopology
( //(
const label faceI, // const label faceI,
const vector& n1, // const vector& n1,
const vector& n2, // const vector& n2,
const vector& testDir // const vector& testDir
) const //) const
{ //{
//
label patchI = mesh_.boundaryMesh().whichPatch(faceI); // label patchI = mesh_.boundaryMesh().whichPatch(faceI);
if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled()) // if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
{ // {
return true; // return true;
} // }
else if (mag(n1&n2) > cos(degToRad(30))) // else if (mag(n1&n2) > cos(degToRad(30)))
{ // {
// Both normals aligned. Check that test vector perpendicularish to // // Both normals aligned. Check that test vector perpendicularish to
// surface normal // // surface normal
scalar magTestDir = mag(testDir); // scalar magTestDir = mag(testDir);
if (magTestDir > VSMALL) // if (magTestDir > VSMALL)
{ // {
if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45))) // if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
{ // {
//Pout<< "** disabling baffling face " // //Pout<< "** disabling baffling face "
// << mesh_.faceCentres()[faceI] << endl; // // << mesh_.faceCentres()[faceI] << endl;
return false; // return false;
} // }
} // }
} // }
return true; // return true;
} //}
// Determine patches for baffles on all intersected unnamed faces // Determine patches for baffles on all intersected unnamed faces
@ -800,7 +801,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
// region. // region.
Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
( (
const List<labelPair>& couples const List<labelPair>& couples,
const scalar planarAngle
) const ) const
{ {
// All duplicate faces on edge of the patch are to be merged. // All duplicate faces on edge of the patch are to be merged.
@ -809,6 +811,22 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
labelList nBafflesPerEdge(mesh_.nEdges(), 0); labelList nBafflesPerEdge(mesh_.nEdges(), 0);
// This algorithm is quite tricky. We don't want to use edgeFaces and
// also want it to run in parallel so it is now an algorithm over
// all (boundary) faces instead.
// We want to pick up any edges that are only used by the baffle
// or internal faces but not by any other boundary faces. So
// - increment count on an edge by 1 if it is used by any (uncoupled)
// boundary face.
// - increment count on an edge by 1000000 if it is used by a baffle face
// - sum in parallel
//
// So now any edge that is used by baffle faces only will have the
// value 2*1000000+2*1.
const label baffleValue = 1000000;
// Count number of boundary faces per edge // Count number of boundary faces per edge
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@ -847,18 +865,22 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
forAll(couples, i) forAll(couples, i)
{ {
const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0); {
label f0 = couples[i].first();
const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
forAll(fEdges0, fEdgeI) forAll(fEdges0, fEdgeI)
{ {
nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000; nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
}
} }
const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1); {
label f1 = couples[i].second();
const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
forAll(fEdges1, fEdgeI) forAll(fEdges1, fEdgeI)
{ {
nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000; nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
}
} }
} }
@ -873,8 +895,8 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
// Baffles which are not next to other boundaries and baffles will have // Baffles which are not next to other boundaries and baffles will have
// nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which
// both baffle faces) // are both baffle faces)
List<labelPair> filteredCouples(couples.size()); List<labelPair> filteredCouples(couples.size());
label filterI = 0; label filterI = 0;
@ -889,15 +911,15 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
== patches.whichPatch(couple.second()) == patches.whichPatch(couple.second())
) )
{ {
const labelList& fEdges = mesh_.faceEdges(couples[i].first()); const labelList& fEdges = mesh_.faceEdges(couple.first());
forAll(fEdges, fEdgeI) forAll(fEdges, fEdgeI)
{ {
label edgeI = fEdges[fEdgeI]; label edgeI = fEdges[fEdgeI];
if (nBafflesPerEdge[edgeI] == 2*1000000+2*1) if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
{ {
filteredCouples[filterI++] = couples[i]; filteredCouples[filterI++] = couple;
break; break;
} }
} }
@ -905,111 +927,127 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
} }
filteredCouples.setSize(filterI); filteredCouples.setSize(filterI);
//Info<< "freeStandingBaffles : from "
// << returnReduce(couples.size(), sumOp<label>()) label nFiltered = returnReduce(filteredCouples.size(), sumOp<label>());
// << " down to "
// << returnReduce(filteredCouples.size(), sumOp<label>()) Info<< "freeStandingBaffles : detected "
// << " baffles." << nl << endl; << nFiltered
<< " free-standing baffles out of "
<< returnReduce(couples.size(), sumOp<label>())
<< nl << endl;
if (nFiltered > 0)
{
// Collect segments
// ~~~~~~~~~~~~~~~~
//XXXXXX pointField start(filteredCouples.size());
// { pointField end(filteredCouples.size());
// // Collect segments
// // ~~~~~~~~~~~~~~~~ const pointField& cellCentres = mesh_.cellCentres();
//
// pointField start(filteredCouples.size()); forAll(filteredCouples, i)
// pointField end(filteredCouples.size()); {
// const labelPair& couple = filteredCouples[i];
// const pointField& cellCentres = mesh_.cellCentres(); start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
// end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
// forAll(filteredCouples, i) }
// {
// const labelPair& couple = couples[i]; // Extend segments a bit
// start[i] = cellCentres[mesh_.faceOwner()[couple.first()]]; {
// end[i] = cellCentres[mesh_.faceOwner()[couple.second()]]; const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
// } start -= smallVec;
// end += smallVec;
// // Extend segments a bit }
// {
// const vectorField smallVec(Foam::sqrt(SMALL)*(end-start));
// start -= smallVec; // Do test for intersections
// end += smallVec; // ~~~~~~~~~~~~~~~~~~~~~~~~~
// } labelList surface1;
// List<pointIndexHit> hit1;
// labelList region1;
// // Do test for intersections vectorField normal1;
// // ~~~~~~~~~~~~~~~~~~~~~~~~~
// labelList surface1; labelList surface2;
// List<pointIndexHit> hit1; List<pointIndexHit> hit2;
// labelList region1; labelList region2;
// vectorField normal1; vectorField normal2;
//
// labelList surface2; surfaces_.findNearestIntersection
// List<pointIndexHit> hit2; (
// labelList region2; identity(surfaces_.surfaces().size()),
// vectorField normal2; start,
// end,
// surfaces_.findNearestIntersection
surface1,
hit1,
region1,
normal1,
surface2,
hit2,
region2,
normal2
);
//mkDir(mesh_.time().path()/timeName());
//OBJstream str
//( //(
// surfacesToBaffle, // mesh_.time().path()/timeName()/"flatBaffles.obj"
// 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
const scalar planarAngleCos = Foam::cos(degToRad(planarAngle));
label filterI = 0;
forAll(filteredCouples, i)
{
const labelPair& couple = filteredCouples[i];
if
(
hit1[i].hit()
&& hit2[i].hit()
&& (
surface1[i] != surface2[i]
|| hit1[i].index() != hit2[i].index()
)
)
{
// Two different hits. Check angle.
//str.write
//(
// linePointRef(hit1[i].hitPoint(), hit2[i].hitPoint()),
// normal1[i],
// normal2[i]
//);
if ((normal1[i]&normal2[i]) > planarAngleCos)
{
// Both normals aligned
vector n = end[i]-start[i];
scalar magN = mag(n);
if (magN > VSMALL)
{
filteredCouples[filterI++] = couple;
}
}
}
else if (hit1[i].hit() || hit2[i].hit())
{
// Single hit. Do not include in freestanding baffles.
}
}
filteredCouples.setSize(filterI);
Info<< "freeStandingBaffles : detected "
<< returnReduce(filterI, sumOp<label>())
<< " planar (within " << planarAngle
<< " degrees) free-standing baffles out of "
<< nFiltered
<< nl << endl;
}
return filteredCouples; return filteredCouples;
} }
@ -1832,77 +1870,15 @@ void Foam::meshRefinement::makeConsistentFaceIndex
} }
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // void Foam::meshRefinement::handleSnapProblems
// Split off unreachable areas of mesh.
void Foam::meshRefinement::baffleAndSplitMesh
( (
const bool handleSnapProblems,
const bool removeEdgeConnectedCells, const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle, const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
const dictionary& motionDict, const dictionary& motionDict,
Time& runTime, Time& runTime,
const labelList& globalToMasterPatch, const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch, const labelList& globalToSlavePatch
const point& keepPoint
) )
{
// Introduce baffles
// ~~~~~~~~~~~~~~~~~
// Split the mesh along internal faces wherever there is a pierce between
// two cell centres.
Info<< "Introducing baffles for "
<< returnReduce(countHits(), sumOp<label>())
<< " faces that are intersected by the surface." << nl << endl;
// Swap neighbouring cell centres and cell level
labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
labelList ownPatch, neiPatch;
getBafflePatches
(
globalToMasterPatch,
neiLevel,
neiCc,
ownPatch,
neiPatch
);
createBaffles(ownPatch, neiPatch);
if (debug)
{
// Debug:test all is still synced across proc patches
checkData();
}
Info<< "Created baffles in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
printMeshInfo(debug, "After introducing baffles");
if (debug&meshRefinement::MESH)
{
Pout<< "Writing baffled mesh to time " << timeName()
<< endl;
write(debug, runTime.path()/"baffles");
Pout<< "Dumped debug data in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
}
// Introduce baffles to delete problem cells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Create some additional baffles where we want surface cells removed.
if (handleSnapProblems)
{ {
Info<< nl Info<< nl
<< "Introducing baffles to block off problem cells" << nl << "Introducing baffles to block off problem cells" << nl
@ -1918,7 +1894,6 @@ void Foam::meshRefinement::baffleAndSplitMesh
perpendicularAngle, perpendicularAngle,
globalToMasterPatch globalToMasterPatch
) )
//markFacesOnProblemCellsGeometric(motionDict)
); );
Info<< "Analyzed problem cells in = " Info<< "Analyzed problem cells in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl; << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
@ -1981,6 +1956,166 @@ void Foam::meshRefinement::baffleAndSplitMesh
} }
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Split off unreachable areas of mesh.
void Foam::meshRefinement::baffleAndSplitMesh
(
const bool doHandleSnapProblems,
const bool removeEdgeConnectedCells,
const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
const scalar freeStandingAngle,
const dictionary& motionDict,
Time& runTime,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch,
const point& keepPoint
)
{
// Introduce baffles
// ~~~~~~~~~~~~~~~~~
// Split the mesh along internal faces wherever there is a pierce between
// two cell centres.
Info<< "Introducing baffles for "
<< returnReduce(countHits(), sumOp<label>())
<< " faces that are intersected by the surface." << nl << endl;
// Swap neighbouring cell centres and cell level
labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
calcNeighbourData(neiLevel, neiCc);
labelList ownPatch, neiPatch;
getBafflePatches
(
globalToMasterPatch,
neiLevel,
neiCc,
ownPatch,
neiPatch
);
createBaffles(ownPatch, neiPatch);
if (debug)
{
// Debug:test all is still synced across proc patches
checkData();
}
Info<< "Created baffles in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
printMeshInfo(debug, "After introducing baffles");
if (debug&meshRefinement::MESH)
{
Pout<< "Writing baffled mesh to time " << timeName()
<< endl;
write(debug, runTime.path()/"baffles");
Pout<< "Dumped debug data in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
}
// Introduce baffles to delete problem cells
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Create some additional baffles where we want surface cells removed.
if (doHandleSnapProblems)
{
//Info<< nl
// << "Introducing baffles to block off problem cells" << nl
// << "----------------------------------------------" << nl
// << endl;
//
//labelList facePatch
//(
// markFacesOnProblemCells
// (
// motionDict,
// removeEdgeConnectedCells,
// perpendicularAngle,
// globalToMasterPatch
// )
// //markFacesOnProblemCellsGeometric(motionDict)
//);
//Info<< "Analyzed problem cells in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//
//if (debug&meshRefinement::MESH)
//{
// faceSet problemTopo(mesh_, "problemFacesTopo", 100);
//
// const labelList facePatchTopo
// (
// markFacesOnProblemCells
// (
// motionDict,
// removeEdgeConnectedCells,
// perpendicularAngle,
// globalToMasterPatch
// )
// );
// forAll(facePatchTopo, faceI)
// {
// if (facePatchTopo[faceI] != -1)
// {
// problemTopo.insert(faceI);
// }
// }
// problemTopo.instance() = timeName();
// Pout<< "Dumping " << problemTopo.size()
// << " problem faces to " << problemTopo.objectPath() << endl;
// problemTopo.write();
//}
//
//Info<< "Introducing baffles to delete problem cells." << nl << endl;
//
//if (debug)
//{
// runTime++;
//}
//
//// Create baffles with same owner and neighbour for now.
//createBaffles(facePatch, facePatch);
//
//if (debug)
//{
// // Debug:test all is still synced across proc patches
// checkData();
//}
//Info<< "Created baffles in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//
//printMeshInfo(debug, "After introducing baffles");
//
//if (debug&meshRefinement::MESH)
//{
// Pout<< "Writing extra baffled mesh to time "
// << timeName() << endl;
// write(debug, runTime.path()/"extraBaffles");
// Pout<< "Dumped debug data in = "
// << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
//}
handleSnapProblems
(
removeEdgeConnectedCells,
perpendicularAngle,
motionDict,
runTime,
globalToMasterPatch,
globalToSlavePatch
);
}
// Select part of mesh // Select part of mesh
// ~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~
@ -2036,7 +2171,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
( (
identity(mesh_.nFaces()-mesh_.nInternalFaces()) identity(mesh_.nFaces()-mesh_.nInternalFaces())
+mesh_.nInternalFaces() +mesh_.nInternalFaces()
) ),
freeStandingAngle
) )
); );
@ -2052,6 +2188,18 @@ void Foam::meshRefinement::baffleAndSplitMesh
// from them. // from them.
mergeBaffles(couples); mergeBaffles(couples);
// Detect any problem cells resulting from merging of baffles
// and delete them
handleSnapProblems
(
removeEdgeConnectedCells,
perpendicularAngle,
motionDict,
runTime,
globalToMasterPatch,
globalToSlavePatch
);
if (debug) if (debug)
{ {
// Debug:test all is still synced across proc patches // Debug:test all is still synced across proc patches
@ -2663,18 +2811,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
if (surface1[i] != -1) if (surface1[i] != -1)
{ {
//- 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 both hit should probably choose 'nearest'
if if
( (

View File

@ -244,14 +244,10 @@ bool Foam::meshRefinement::markForRefine
} }
// Calculates list of cells to refine based on intersection with feature edge. void Foam::meshRefinement::markFeatureCellLevel
Foam::label Foam::meshRefinement::markFeatureRefinement
( (
const point& keepPoint, const point& keepPoint,
const label nAllowRefine, labelList& maxFeatureLevel
labelList& refineCell,
label& nRefine
) const ) const
{ {
// We want to refine all cells containing a feature edge. // We want to refine all cells containing a feature edge.
@ -384,7 +380,7 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
// Largest refinement level of any feature passed through // Largest refinement level of any feature passed through
labelList maxFeatureLevel(mesh_.nCells(), -1); maxFeatureLevel = labelList(mesh_.nCells(), -1);
// Database to pass into trackedParticle::move // Database to pass into trackedParticle::move
trackedParticle::trackingData td(cloud, maxFeatureLevel); trackedParticle::trackingData td(cloud, maxFeatureLevel);
@ -467,8 +463,23 @@ Foam::label Foam::meshRefinement::markFeatureRefinement
// Track all particles to their end position. // Track all particles to their end position.
cloud.move(td, GREAT); cloud.move(td, GREAT);
} }
}
// Calculates list of cells to refine based on intersection with feature edge.
Foam::label Foam::meshRefinement::markFeatureRefinement
(
const point& keepPoint,
const label nAllowRefine,
labelList& refineCell,
label& nRefine
) const
{
// Largest refinement level of any feature passed through
labelList maxFeatureLevel;
markFeatureCellLevel(keepPoint, maxFeatureLevel);
// See which cells to refine. maxFeatureLevel will hold highest level // See which cells to refine. maxFeatureLevel will hold highest level
// of any feature edge that passed through. // of any feature edge that passed through.