sampledSet/midPoint, midPointAndFace: Improved robustness of the mid-point cell seaching and selecting

This commit is contained in:
Henry Weller
2016-03-19 21:22:09 +00:00
parent 649128313b
commit a7e410396a
4 changed files with 119 additions and 155 deletions

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-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,6 +25,7 @@ License
#include "midPointSet.H" #include "midPointSet.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "meshSearch.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,58 +48,46 @@ void Foam::midPointSet::genSamples()
labelList midSegments(2*size()); labelList midSegments(2*size());
scalarList midCurveDist(2*size()); scalarList midCurveDist(2*size());
label midI = 0; label mSamplei = 0;
label samplei = 0;
label sampleI = 0; while (size() > 0)
while(true && size()>0)
{ {
// calculate midpoint between sampleI and sampleI+1 (if in same segment) // Calculate midpoint between samplei and samplei+1 (if in same segment)
while while
( (
(sampleI < size() - 1) (samplei < size() - 1)
&& (segments_[sampleI] == segments_[sampleI+1]) && (segments_[samplei] == segments_[samplei+1])
) )
{ {
midPoints[midI] = point midPoint(0.5*(operator[](samplei) + operator[](samplei+1)));
0.5*(operator[](sampleI) + operator[](sampleI+1)); label cellm = pointInCell(midPoint, samplei);
label cell1 = getCell(faces_[sampleI], midPoints[midI]); if (cellm != -1)
label cell2 = getCell(faces_[sampleI+1], midPoints[midI]);
if (cell1 != cell2)
{ {
FatalErrorInFunction midPoints[mSamplei] = midPoint;
<< " midI:" << midI midCells[mSamplei] = cellm;
<< " sampleI:" << sampleI midSegments[mSamplei] = segments_[samplei];
<< " pts[sampleI]:" << operator[](sampleI) midCurveDist[mSamplei] = mag(midPoints[mSamplei] - start());
<< " face[sampleI]:" << faces_[sampleI] mSamplei++;
<< " pts[sampleI+1]:" << operator[](sampleI+1)
<< " face[sampleI+1]:" << faces_[sampleI+1]
<< " cell1:" << cell1
<< " cell2:" << cell2
<< abort(FatalError);
} }
midCells[midI] = cell1; samplei++;
midSegments[midI] = segments_[sampleI];
midCurveDist[midI] = mag(midPoints[midI] - start());
midI++;
sampleI++;
} }
if (sampleI == size() - 1) if (samplei == size() - 1)
{ {
break; break;
} }
sampleI++;
samplei++;
} }
midPoints.setSize(midI); midPoints.setSize(mSamplei);
midCells.setSize(midI); midCells.setSize(mSamplei);
midSegments.setSize(midI); midSegments.setSize(mSamplei);
midCurveDist.setSize(midI); midCurveDist.setSize(mSamplei);
setSamples setSamples
( (
midPoints, midPoints,

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-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -39,101 +39,83 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Rework faceOnlySet samples.
// Take two consecutive samples
void Foam::midPointAndFaceSet::genSamples() void Foam::midPointAndFaceSet::genSamples()
{ {
// Generate midpoints and add to face points // Generate midpoints and add to face points
List<point> newSamplePoints(3*size()); List<point> mpfSamplePoints(3*size());
labelList newSampleCells(3*size()); labelList mpfSampleCells(3*size());
labelList newSampleFaces(3*size()); labelList mpfSampleFaces(3*size());
labelList newSampleSegments(3*size()); labelList mpfSampleSegments(3*size());
scalarList newSampleCurveDist(3*size()); scalarList mpfSampleCurveDist(3*size());
label newSampleI = 0; label mpfSamplei = 0;
label samplei = 0;
label sampleI = 0; while (size() > 0)
while(true && size()>0)
{ {
// sampleI is start of segment // Add first face
mpfSamplePoints[mpfSamplei] = operator[](samplei);
// Add sampleI mpfSampleCells[mpfSamplei] = cells_[samplei];
newSamplePoints[newSampleI] = operator[](sampleI); mpfSampleFaces[mpfSamplei] = faces_[samplei];
newSampleCells[newSampleI] = cells_[sampleI]; mpfSampleSegments[mpfSamplei] = segments_[samplei];
newSampleFaces[newSampleI] = faces_[sampleI]; mpfSampleCurveDist[mpfSamplei] = curveDist_[samplei];
newSampleSegments[newSampleI] = segments_[sampleI]; mpfSamplei++;
newSampleCurveDist[newSampleI] = curveDist_[sampleI];
newSampleI++;
while while
( (
(sampleI < size() - 1) (samplei < size() - 1)
&& (segments_[sampleI] == segments_[sampleI+1]) && (segments_[samplei] == segments_[samplei+1])
) )
{ {
// Add mid point point midPoint(0.5*(operator[](samplei) + operator[](samplei+1)));
const point mid = 0.5*(operator[](sampleI) + operator[](sampleI+1)); label cellm = pointInCell(midPoint, samplei);
label cell1 = getCell(faces_[sampleI], mid); if (cellm != -1)
label cell2 = getCell(faces_[sampleI+1], mid);
if (cell1 != cell2)
{ {
FatalErrorInFunction mpfSamplePoints[mpfSamplei] = midPoint;
<< " newSampleI:" << newSampleI mpfSampleCells[mpfSamplei] = cellm;
<< " pts[sampleI]:" << operator[](sampleI) mpfSampleFaces[mpfSamplei] = -1;
<< " face[sampleI]:" << faces_[sampleI] mpfSampleSegments[mpfSamplei] = segments_[samplei];
<< " pts[sampleI+1]:" << operator[](sampleI+1) mpfSampleCurveDist[mpfSamplei] =
<< " face[sampleI+1]:" << faces_[sampleI+1] mag(mpfSamplePoints[mpfSamplei] - start());
<< " cell1:" << cell1
<< " cell2:" << cell2 mpfSamplei++;
<< abort(FatalError);
} }
newSamplePoints[newSampleI] = mid; // Add second face
newSampleCells[newSampleI] = cell1; mpfSamplePoints[mpfSamplei] = operator[](samplei+1);
newSampleFaces[newSampleI] = -1; mpfSampleCells[mpfSamplei] = cells_[samplei+1];
newSampleSegments[newSampleI] = segments_[sampleI]; mpfSampleFaces[mpfSamplei] = faces_[samplei+1];
newSampleCurveDist[newSampleI] = mpfSampleSegments[mpfSamplei] = segments_[samplei+1];
mag(newSamplePoints[newSampleI] - start()); mpfSampleCurveDist[mpfSamplei] =
mag(mpfSamplePoints[mpfSamplei] - start());
newSampleI++; mpfSamplei++;
// Add sampleI+1 samplei++;
newSamplePoints[newSampleI] = operator[](sampleI+1);
newSampleCells[newSampleI] = cells_[sampleI+1];
newSampleFaces[newSampleI] = faces_[sampleI+1];
newSampleSegments[newSampleI] = segments_[sampleI+1];
newSampleCurveDist[newSampleI] =
mag(newSamplePoints[newSampleI] - start());
newSampleI++;
sampleI++;
} }
if (sampleI == size() - 1) if (samplei == size() - 1)
{ {
break; break;
} }
sampleI++; samplei++;
} }
newSamplePoints.setSize(newSampleI); mpfSamplePoints.setSize(mpfSamplei);
newSampleCells.setSize(newSampleI); mpfSampleCells.setSize(mpfSamplei);
newSampleFaces.setSize(newSampleI); mpfSampleFaces.setSize(mpfSamplei);
newSampleSegments.setSize(newSampleI); mpfSampleSegments.setSize(mpfSamplei);
newSampleCurveDist.setSize(newSampleI); mpfSampleCurveDist.setSize(mpfSamplei);
setSamples setSamples
( (
newSamplePoints, mpfSamplePoints,
newSampleCells, mpfSampleCells,
newSampleFaces, mpfSampleFaces,
newSampleSegments, mpfSampleSegments,
newSampleCurveDist mpfSampleCurveDist
); );
} }

View File

@ -47,61 +47,57 @@ Foam::label Foam::sampledSet::getBoundaryCell(const label faceI) const
} }
Foam::label Foam::sampledSet::getCell Foam::label Foam::sampledSet::pointInCell
( (
const label faceI, const point& p,
const point& sample const label samplei
) const ) const
{ {
if (faceI == -1) const label cellio = mesh().faceOwner()[faces_[samplei]];
const label cellin = mesh().faceNeighbour()[faces_[samplei]];
const label celljo = mesh().faceOwner()[faces_[samplei+1]];
const label celljn = mesh().faceNeighbour()[faces_[samplei+1]];
// If mid-point is in the cell including the sampled faces
// include in list otherwise ignore
label cellm =
(cellio == celljo || cellio == celljn) ? cellio
: (cellin == celljo || cellin == celljn) ? cellin
: -1;
if (cellm != -1)
{ {
FatalErrorInFunction // Check the mid-point is in the cell otherwise ignore
<< "Illegal face label " << faceI if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
<< abort(FatalError);
}
if (faceI >= mesh().nInternalFaces())
{
label cellI = getBoundaryCell(faceI);
if (!mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{ {
FatalErrorInFunction cellm = -1;
<< "Found cell " << cellI << " using face " << faceI
<< ". But cell does not contain point " << sample
<< abort(FatalError);
}
return cellI;
}
else
{
// Try owner and neighbour to see which one contains sample
label cellI = mesh().faceOwner()[faceI]; if (debug)
if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{
return cellI;
}
else
{
cellI = mesh().faceNeighbour()[faceI];
if (mesh().pointInCell(sample, cellI, searchEngine_.decompMode()))
{ {
return cellI; WarningInFunction
} << "Could not find mid-point " << p
else << " cell " << cellm << endl;
{
FatalErrorInFunction
<< "None of the neighbours of face "
<< faceI << " contains point " << sample
<< abort(FatalError);
return -1;
} }
} }
} }
else if (debug)
{
WarningInFunction
<< "Could not find cell for mid-point" << nl
<< " samplei: " << samplei
<< " pts[samplei]: " << operator[](samplei)
<< " face[samplei]: " << faces_[samplei]
<< " pts[samplei+1]: " << operator[](samplei+1)
<< " face[samplei+1]: " << faces_[samplei+1]
<< " cellio: " << cellio
<< " cellin: " << cellin
<< " celljo: " << celljo
<< " celljn: " << celljn
<< endl;
}
return cellm;
} }

View File

@ -92,12 +92,9 @@ protected:
//- Returns cell next to boundary face //- Returns cell next to boundary face
label getBoundaryCell(const label) const; label getBoundaryCell(const label) const;
//- Returns cell using face and containing sample //- Return the cell in which the point on the sample line
label getCell // resides if found otherwise return -1
( label pointInCell(const point& p, const label samplei) const;
const label faceI,
const point& sample
) const;
//- Calculates inproduct of face normal and vector sample-face centre //- Calculates inproduct of face normal and vector sample-face centre
// <0 if sample inside. // <0 if sample inside.