diff --git a/cellFinders/CellFinder/CellFinder.C b/cellFinders/CellFinder/CellFinder.C index 8324807..323b3a3 100644 --- a/cellFinders/CellFinder/CellFinder.C +++ b/cellFinders/CellFinder/CellFinder.C @@ -26,7 +26,6 @@ License #include "objectRegistry.H" #include "IOField.H" #include "codeRules.H" -#include "patchDistMethod.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -83,87 +82,9 @@ Foam::autoPtr Foam::CellFinder::New void Foam::CellFinder::createFields() const { - } -Foam::tmp Foam::CellFinder::distanceField() const -{ - - // Grab h for the current patch - const volScalarField & h = mesh_.lookupObject ("h"); - if (debug) - { - Info<< "CellFinder: Creating dist field" << nl; - } - - tmp dist - ( - new volScalarField - ( - IOobject - ( - "dist", - mesh_.time().timeName(), - mesh_, - IOobject::READ_IF_PRESENT, - IOobject::NO_WRITE - ), - mesh_, - dimensionedScalar("dist", dimLength, 0), - h.boundaryField().types() - ) - ); - - bool precomputedDist = - #ifdef FOAM_NEW_GEOMFIELD_RULES - mag(max(dist().primitiveField())) > VSMALL; - #else - mag(max(dist().internalField())) > VSMALL; - #endif - - if (debug) - { - Info<<"CellFinder: using precumputed distance field" << nl; - } - - if (!precomputedDist) - { - labelHashSet patchIDs(1); - patchIDs.insert(patch().index()); - - dictionary methodDict = dictionary(); - methodDict.lookupOrAddDefault(word("method"), word("meshWave")); - - if (debug) - { - Info<< "Initializing patchDistanceMethod" << nl; - } - - autoPtr pdm - ( - patchDistMethod::New - ( - methodDict, - mesh_, - patchIDs - ) - ); - - if (debug) - { - Info<< "CellFinder: Computing dist field" << nl; - } -#ifdef FOAM_NEW_TMP_RULES - pdm->correct(dist.ref()); -#else - pdm->correct(dist()); -#endif - } - - return dist; -} - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // diff --git a/cellFinders/CellFinder/CellFinder.H b/cellFinders/CellFinder/CellFinder.H index 89110f2..4a55528 100644 --- a/cellFinders/CellFinder/CellFinder.H +++ b/cellFinders/CellFinder/CellFinder.H @@ -63,9 +63,6 @@ protected: // Protected Member Functions - //- Compute distance field - tmp distanceField() const; - void createFields() const; public: diff --git a/cellFinders/TreeCellFinder/TreeCellFinder.C b/cellFinders/TreeCellFinder/TreeCellFinder.C index 4f63f96..57db9a6 100644 --- a/cellFinders/TreeCellFinder/TreeCellFinder.C +++ b/cellFinders/TreeCellFinder/TreeCellFinder.C @@ -189,17 +189,31 @@ void Foam::TreeCellFinder::findCellIndices point = faceCentres[i] - faceNormals[i]*h[i]; // If h is zero, or the point is outside the domain, - // set it to distance to adjacent cell's centre - // Set the cellIndexList component accordingly + // grab the wall-adjacent cell bool inside = #ifdef FOAM_VOLUMETYPE_NOT_CAPITAL boundaryTreePtr->getVolumeType(point) == volumeType::inside; #else boundaryTreePtr->getVolumeType(point) == volumeType::INSIDE; #endif - if ((h[i] == 0) || (!inside) || (treePtr->nodes().empty())) + if ((h[i] <= 0) || (!inside) || (treePtr->nodes().empty())) { - //h_[i] = mag(cellCentres[i] - faceCentres[i]); + if (h[i] < 0) + { + Warning + << "TreeCellFinder: " << h[i] + << " is negative and thus not a valid distance. " + << "Will fall back to wall-adjacent cell for face " + << i << " on patch " << patch().name() << nl; + } + else if (!inside) + { + Warning + << "TreeCellFinder: the point " << h[i] + << " away from the wall is outside the domain. " + << "Will fall back to wall-adjacent cell for face " + << i << " on patch " << patch().name() << nl; + } indexList[i] = faceCells[i]; } else @@ -207,7 +221,6 @@ void Foam::TreeCellFinder::findCellIndices pih = treePtr->findNearest(point, treePtr->bb().mag()); indexList[i] = searchCellLabels[pih.index()]; - //h_[i] = mag(C[indexList_[i]] - faceCentres[i]); } } if (debug) @@ -477,4 +490,82 @@ Foam::TreeCellFinder::findCandidateCellLabels return tCandidates; } +Foam::tmp Foam::TreeCellFinder::distanceField() const +{ + + // Grab h for the current patch + const volScalarField & h = mesh_.lookupObject ("h"); + if (debug) + { + Info<< "CellFinder: Creating dist field" << nl; + } + + tmp dist + ( + new volScalarField + ( + IOobject + ( + "dist", + mesh_.time().timeName(), + mesh_, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh_, + dimensionedScalar("dist", dimLength, 0), + h.boundaryField().types() + ) + ); + + bool precomputedDist = + #ifdef FOAM_NEW_GEOMFIELD_RULES + mag(max(dist().primitiveField())) > VSMALL; + #else + mag(max(dist().internalField())) > VSMALL; + #endif + + if (debug) + { + Info<<"CellFinder: using precumputed distance field" << nl; + } + + if (!precomputedDist) + { + labelHashSet patchIDs(1); + patchIDs.insert(patch().index()); + + dictionary methodDict = dictionary(); + methodDict.lookupOrAddDefault(word("method"), word("meshWave")); + + if (debug) + { + Info<< "Initializing patchDistanceMethod" << nl; + } + + autoPtr pdm + ( + patchDistMethod::New + ( + methodDict, + mesh_, + patchIDs + ) + ); + + if (debug) + { + Info<< "CellFinder: Computing dist field" << nl; + } +#ifdef FOAM_NEW_TMP_RULES + pdm->correct(dist.ref()); +#else + pdm->correct(dist()); +#endif + } + + return dist; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/cellFinders/TreeCellFinder/TreeCellFinder.H b/cellFinders/TreeCellFinder/TreeCellFinder.H index 4cd1893..281311c 100644 --- a/cellFinders/TreeCellFinder/TreeCellFinder.H +++ b/cellFinders/TreeCellFinder/TreeCellFinder.H @@ -39,6 +39,7 @@ SourceFiles #include "runTimeSelectionTables.H" #include "addToRunTimeSelectionTable.H" #include "CellFinder.H" +#include "patchDistMethod.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -88,6 +89,9 @@ public: // Member functions + //- Compute distance field + tmp distanceField() const; + //- Find the sampling cell indices for a single cell sampler void findCellIndices ( diff --git a/tests/cellFinders/TreeCellFinder/testTreeCellFinder.C b/tests/cellFinders/TreeCellFinder/testTreeCellFinder.C index f42616d..b57e050 100644 --- a/tests/cellFinders/TreeCellFinder/testTreeCellFinder.C +++ b/tests/cellFinders/TreeCellFinder/testTreeCellFinder.C @@ -26,13 +26,352 @@ TEST_F(TreeCellFinderTest, FullConstructor) const fvPatch & patch = mesh.boundary()["bottomWall"]; h.boundaryFieldRef()[patch.index()] == 0.25; TreeCellFinder finder(patch); +} + + +TEST_F(TreeCellFinderTest, FindDistanceBasedBottomWall) +{ + extern argList * mainArgs; + const argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + volScalarField & h = const_cast + ( + mesh.thisDb().lookupObject("h") + ); + const volVectorField & C = mesh.C(); + + const fvPatch & patch = mesh.boundary()["bottomWall"]; + TreeCellFinder finder(patch); labelList indexList(patch.size(), 0); + + h.boundaryFieldRef()[patch.index()] == 0.1; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.1); + } + + h.boundaryFieldRef()[patch.index()] == 0.32; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.3); + } + + h.boundaryFieldRef()[patch.index()] == 1.9; finder.findCellIndices ( indexList, h.boundaryFieldRef()[patch.index()] ); - Info << indexList; + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 1.9); + } } + + +TEST_F(TreeCellFinderTest, FindDistanceBasedTopWall) +{ + extern argList * mainArgs; + const argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + volScalarField & h = const_cast + ( + mesh.thisDb().lookupObject("h") + ); + const volVectorField & C = mesh.C(); + + const fvPatch & patch = mesh.boundary()["topWall"]; + TreeCellFinder finder(patch); + labelList indexList(patch.size(), 0); + + h.boundaryFieldRef()[patch.index()] == 0.1; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 1.9); + } + + h.boundaryFieldRef()[patch.index()] == 0.32; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 1.7); + } + + h.boundaryFieldRef()[patch.index()] == 1.9; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.1); + } +} + + +TEST_F(TreeCellFinderTest, FindDistanceInvalidDistance) +{ + extern argList * mainArgs; + const argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + volScalarField & h = const_cast + ( + mesh.thisDb().lookupObject("h") + ); + const volVectorField & C = mesh.C(); + + const fvPatch & patch = mesh.boundary()["bottomWall"]; + TreeCellFinder finder(patch); + labelList indexList(patch.size(), 0); + + h.boundaryFieldRef()[patch.index()] == -0.1; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.1); + } +} + +TEST_F(TreeCellFinderTest, FindDistanceZeroDistance) +{ + extern argList * mainArgs; + const argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + volScalarField & h = const_cast + ( + mesh.thisDb().lookupObject("h") + ); + const volVectorField & C = mesh.C(); + + const fvPatch & patch = mesh.boundary()["bottomWall"]; + TreeCellFinder finder(patch); + labelList indexList(patch.size(), 0); + + h.boundaryFieldRef()[patch.index()] == 0.0; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.1); + } +} + + +TEST_F(TreeCellFinderTest, FindDistanceTooLargeDistance) +{ + extern argList * mainArgs; + const argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + volScalarField & h = const_cast + ( + mesh.thisDb().lookupObject("h") + ); + const volVectorField & C = mesh.C(); + + const fvPatch & patch = mesh.boundary()["bottomWall"]; + TreeCellFinder finder(patch); + labelList indexList(patch.size(), 0); + + h.boundaryFieldRef()[patch.index()] == 10.0; + finder.findCellIndices + ( + indexList, + h.boundaryFieldRef()[patch.index()] + ); + + forAll(indexList, i) + { + ASSERT_FLOAT_EQ(C[indexList[i]][1], 0.1); + } +} + + +TEST_F(TreeCellFinderTest, DistanceFieldBottom) +{ + extern argList * mainArgs; + argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + + const fvPatch & patch = mesh.boundary()["bottomWall"]; + TreeCellFinder finder(patch); + + tmp tDist = finder.distanceField(); + const volScalarField & dist = tDist(); + + const vectorField C = mesh.C(); + + forAll (C, i) + { + ASSERT_NEAR(dist[i], C[i][1], 1e-8); + } +} + +TEST_F(TreeCellFinderTest, DistanceFieldTop) +{ + extern argList * mainArgs; + argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + + const fvPatch & patch = mesh.boundary()["topWall"]; + TreeCellFinder finder(patch); + + tmp tDist = finder.distanceField(); + const volScalarField & dist = tDist(); + + const vectorField C = mesh.C(); + + forAll (C, i) + { + ASSERT_NEAR(dist[i], 2 - C[i][1], 1e-8); + } +} + + +TEST_F(TreeCellFinderTest, FindCandidateCellLabelsBottom) +{ + extern argList * mainArgs; + argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + const fvPatch & patchBottom = mesh.boundary()["topWall"]; + TreeCellFinder finderBottom(patchBottom); + + scalarField h(finderBottom.patch().size()); + + h = 0.0; + tmp distField(finderBottom.distanceField()); + + const scalarField & dist = +#ifdef FOAM_NEW_GEOMFIELD_RULES + distField().primitiveField(); +#else + distField().internalField(); +#endif + tmp tCellLabelsBottom = + finderBottom.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsBottom().size(), 0); + + h = 0.1; + tCellLabelsBottom = finderBottom.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsBottom().size(), 9); + + h = 0.21; + tCellLabelsBottom = finderBottom.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsBottom().size(), 18); + + h = 1; + tCellLabelsBottom = finderBottom.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsBottom().size(), 90); + +} + + +TEST_F(TreeCellFinderTest, FindCandidateCellLabelsTop) +{ + extern argList * mainArgs; + argList & args = *mainArgs; + Time runTime(Foam::Time::controlDictName, args); + autoPtr meshPtr = createMesh(runTime); + const fvMesh & mesh = meshPtr(); + createSamplingHeightField(mesh); + + const fvPatch & patchTop = mesh.boundary()["topWall"]; + TreeCellFinder finderTop(patchTop); + + scalarField h(finderTop.patch().size()); + + tmp distField = finderTop.distanceField(); + const scalarField & dist = +#ifdef FOAM_NEW_GEOMFIELD_RULES + distField().primitiveField(); +#else + distField().internalField(); +#endif + + h = 0.0; + tmp tCellLabelsTop = finderTop.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsTop().size(), 0) << "h = 0"; + + h = 0.1; + tCellLabelsTop = finderTop.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsTop().size(), 9) << "h = 0.1"; + + h = 0.21; + tCellLabelsTop = finderTop.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsTop().size(), 18) << "h = 0.21"; + + h = 1; + tCellLabelsTop = finderTop.findCandidateCellLabels(dist, h); + ASSERT_EQ(tCellLabelsTop().size(), 90) << "h = 1"; +} \ No newline at end of file diff --git a/tests/samplers/Sampler/testSampler.C b/tests/samplers/Sampler/testSampler.C index 906f20c..978ba7e 100644 --- a/tests/samplers/Sampler/testSampler.C +++ b/tests/samplers/Sampler/testSampler.C @@ -279,103 +279,3 @@ TEST_F(SamplerTest, AddField) ASSERT_EQ(sampler.Sampler::nSampledFields(), 1); } - -// TEST_F(SamplerTest, DistanceFieldBottom) -// { -// extern argList * mainArgs; -// argList & args = *mainArgs; -// Time runTime(Foam::Time::controlDictName, args); -// autoPtr meshPtr = createMesh(runTime); -// const fvMesh & mesh = meshPtr(); -// createSamplingHeightField(mesh); - - -// const fvPatch & patch = mesh.boundary()["bottomWall"]; -// DummySampler sampler("SingleCellSampler", patch, 3.0); - -// tmp tDist = sampler.wrapDistanceField(); -// const volScalarField & dist = tDist(); - -// const vectorField C = mesh.C(); - -// forAll (C, i) -// { -// ASSERT_NEAR(dist[i], C[i][1], 1e-8); -// } -// } - -// TEST_F(SamplerTest, DistanceFieldTop) -// { -// extern argList * mainArgs; -// argList & args = *mainArgs; -// Time runTime(Foam::Time::controlDictName, args); -// autoPtr meshPtr = createMesh(runTime); -// const fvMesh & mesh = meshPtr(); -// createSamplingHeightField(mesh); - - -// const fvPatch & patch = mesh.boundary()["topWall"]; -// DummySampler sampler("SingleCellSampler", patch, 3.0); - -// tmp tDist = sampler.wrapDistanceField(); -// const volScalarField & dist = tDist(); - -// const vectorField C = mesh.C(); - -// forAll (C, i) -// { -// ASSERT_NEAR(dist[i], 2 - C[i][1], 1e-8); -// } -// } - - -// TEST_F(SamplerTest, FindSearchCellLabels) -// { -// extern argList * mainArgs; -// argList & args = *mainArgs; -// Time runTime(Foam::Time::controlDictName, args); -// autoPtr meshPtr = createMesh(runTime); -// const fvMesh & mesh = meshPtr(); -// createSamplingHeightField(mesh); -// volScalarField & h = -// const_cast(mesh.lookupObject ("h")); - -// const fvPatch & patchBottom = mesh.boundary()["topWall"]; - -// h.boundaryFieldRef()[patchBottom.index()] == 0; -// DummySampler samplerBottom("SingleCellSampler", patchBottom, 3.0); -// tmp tCellLabelsBottom = samplerBottom.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsBottom().size(), 0); - -// h.boundaryFieldRef()[patchBottom.index()] == 0.1; -// tCellLabelsBottom = samplerBottom.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsBottom().size(), 9); - -// h.boundaryFieldRef()[patchBottom.index()] == 0.21; -// tCellLabelsBottom = samplerBottom.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsBottom().size(), 18); - -// h.boundaryFieldRef()[patchBottom.index()] == 1; -// tCellLabelsBottom = samplerBottom.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsBottom().size(), 90); - -// const fvPatch & patchTop = mesh.boundary()["topWall"]; - -// h.boundaryFieldRef()[patchTop.index()] == 0; -// DummySampler samplerTop("SingleCellSampler", patchTop, 3.0); -// tmp tCellLabelsTop = samplerTop.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsTop().size(), 0) << "h = 0"; - -// h.boundaryFieldRef()[patchTop.index()] == 0.1; -// tCellLabelsTop = samplerTop.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsTop().size(), 9) << "h = 0.1"; - -// h.boundaryFieldRef()[patchTop.index()] == 0.21; -// tCellLabelsTop = samplerTop.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsTop().size(), 18) << "h = 0.21"; - -// h.boundaryFieldRef()[patchTop.index()] == 1; -// tCellLabelsTop = samplerTop.wrapFindSearchCellLabels(); -// ASSERT_EQ(tCellLabelsTop().size(), 90) << "h = 1"; - -// }