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

Moving to a cell approach for cells that are concave. WIP - still

Browse Source 
problems.
This commit is contained in:
graham
2009-10-06 16:59:46 +01:00
parent 9ef96ef0a8
commit 8fd827aeb3
5 changed files with 307 additions and 450 deletions
Download Patch File Download Diff File Expand all files Collapse all files

201
src/lagrangian/basic/Cloud/Cloud.C
View File

@ -42,7 +42,7 @@ template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::trackingRescueTolerance = 1e-3;
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-6;
const Foam::scalar Foam::Cloud<ParticleType>::intersectionTolerance = 1e-9;
template<class ParticleType>
const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
@ -52,18 +52,27 @@ const Foam::scalar Foam::Cloud<ParticleType>::planarCosAngle = (1 - 1e-6);
// Calculate using face properties only
template<class ParticleType>
Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
(
const label cellI,
const label f0I,
const label f1I
) const
bool Foam::Cloud<ParticleType>::isConcaveCell(const label cellI) const
{
const cellList& cells = pMesh().cells();
const labelList& fOwner = pMesh().faceOwner();
const vectorField& fAreas = pMesh().faceAreas();
const pointField& fCentres = pMesh().faceCentres();
const labelList& fOwner = pMesh().faceOwner();
const pointField& points = pMesh().points();
const cell& cFaces = cells[cellI];
forAll(cFaces, i)
{
label f0I = cFaces[i];
forAll(cFaces, j)
{
if (j != i)
{
label f1I = cFaces[j];
// Get f0 centre
const point& fc0 = fCentres[f0I];
@ -95,17 +104,22 @@ Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
const point& pt = points[ptI];
// Normal distance from fc0 to pt
// If the cell is concave, the point on f1 will be
// on the outside of the plane of f0, defined by
// its centre and normal, and the angle between
// (fc0 -pt) and fn0 will be greater than 90
// degrees, so the dot product will be negative.
scalar d = ((fc0 - pt) & fn0);
if (d < 0)
{
// Proper concave: f1 on wrong side of f0
return 2;
return true;
}
}
// Could be a convex cell, but check angle for planar faces.
// Check for co-planar faces, which are also treated
// as concave, as they are not strictly convex.
vector fn1 = fAreas[f1I];
fn1 /= (mag(fn1) + VSMALL);
@ -120,159 +134,62 @@ Foam::label Foam::Cloud<ParticleType>::isIntrudingFace
if ((fn0 & fn1) > planarCosAngle)
{
// Planar face
return 1;
return true;
}
}
}
}
return 0;
return false;
}
template<class ParticleType>
void Foam::Cloud<ParticleType>::calcIntrudingFaces() const
void Foam::Cloud<ParticleType>::calcConcaveCells() const
{
const labelList& fOwner = pMesh().faceOwner();
const cellList& cells = pMesh().cells();
intrudesIntoOwnerPtr_.reset(new PackedBoolList(pMesh().nFaces()));
PackedBoolList& intrudesIntoOwner = intrudesIntoOwnerPtr_();
intrudesIntoNeighbourPtr_.reset(new PackedBoolList(pMesh().nFaces()));
PackedBoolList& intrudesIntoNeighbour = intrudesIntoNeighbourPtr_();
PackedBoolList cellsWithProblemFaces(cells.size());
DynamicList<label> planarFaces;
DynamicList<label> intrudingFaces;
label nIntruded = 0;
label nPlanar = 0;
concaveCellPtr_.reset(new PackedBoolList(pMesh().nCells()));
PackedBoolList& concaveCell = concaveCellPtr_();
forAll(cells, cellI)
{
const cell& cFaces = cells[cellI];
forAll(cFaces, i)
if (isConcaveCell(cellI))
{
label f0 = cFaces[i];
bool own0 = (fOwner[f0] == cellI);
// Now check that all other faces of the cell are on the 'inside'
// of the face.
bool intrudes = false;
forAll(cFaces, j)
{
if (j != i)
{
label state = isIntrudingFace
(
cellI,
f0,
cFaces[j]
);
if (state == 1)
{
// planar
planarFaces.append(f0);
intrudes = true;
nPlanar++;
break;
}
else if (state == 2)
{
// concave
intrudingFaces.append(f0);
intrudes = true;
nIntruded++;
break;
}
}
concaveCell[cellI] = 1;
}
if (intrudes)
{
cellsWithProblemFaces[cellI] = 1;
if (own0)
{
intrudesIntoOwner[f0] = 1;
}
else
{
intrudesIntoNeighbour[f0] = 1;
}
// Force all cells to be decomposed
concaveCell[cellI] = 1;
}
// intrudesIntoOwner[f0] = 1;
// intrudesIntoNeighbour[f0] = 1;
}
}
// if (debug)
{
Pout<< "Cloud<ParticleType>::calcIntrudingFaces() :"
<< " overall faces in mesh : "
<< pMesh().nFaces() << nl
<< " of these planar : "
<< nPlanar << nl
<< " of these intruding into their cell (concave) : "
<< nIntruded << endl;
// Write cells that are a problem to file
// Info<< "Hard coded picking up every face." << endl;
// Write the faces and cells that are a problem to file to be
// compared and made into sets
{
DynamicList<label> cellsWithIntrudingFaces;
DynamicList<label> tmpConcaveCells;
forAll(cells, cellI)
{
if (cellsWithProblemFaces[cellI])
if (concaveCell[cellI])
{
cellsWithIntrudingFaces.append(cellI);
tmpConcaveCells.append(cellI);
}
}
fileName fName = pMesh().time().path()/"cellsWithProblemFaces";
Pout<< "Cloud<ParticleType>::calcConcaveCells() :"
<< " overall cells in mesh : " << pMesh().nCells() << nl
<< " of these concave : " << tmpConcaveCells.size() << endl;
fileName fName = pMesh().time().path()/"concaveCells";
Pout<< " Writing " << fName.name() << endl;
OFstream file(fName);
file << cellsWithIntrudingFaces;
file << tmpConcaveCells;
file.flush();
}
{
fileName fName = pMesh().time().path()/"planarFaces";
Pout<< " Writing " << fName.name() << endl;
OFstream file(fName);
file << planarFaces;
file.flush();
}
{
fileName fName = pMesh().time().path()/"intrudingFaces";
Pout<< " Writing " << fName.name() << endl;
OFstream file(fName);
file << intrudingFaces;
file.flush();
}
}
}
@ -318,8 +235,7 @@ Foam::Cloud<ParticleType>::Cloud
template<class ParticleType>
void Foam::Cloud<ParticleType>::clearOut()
{
intrudesIntoOwnerPtr_.clear();
intrudesIntoNeighbourPtr_.clear();
concaveCellPtr_.clear();
labels_.clearStorage();
}
@ -327,26 +243,15 @@ void Foam::Cloud<ParticleType>::clearOut()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class ParticleType>
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoOwner()
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::concaveCell()
const
{
if (!intrudesIntoOwnerPtr_.valid())
if (!concaveCellPtr_.valid())
{
calcIntrudingFaces();
}
return intrudesIntoOwnerPtr_();
calcConcaveCells();
}
template<class ParticleType>
const Foam::PackedBoolList& Foam::Cloud<ParticleType>::intrudesIntoNeighbour()
const
{
if (!intrudesIntoNeighbourPtr_.valid())
{
calcIntrudingFaces();
}
return intrudesIntoNeighbourPtr_();
return concaveCellPtr_();
}

29
src/lagrangian/basic/Cloud/Cloud.H
View File

@ -78,10 +78,8 @@ class Cloud
const bool concaveCheck_;
//- Does face intrude into owner cell
mutable autoPtr<PackedBoolList> intrudesIntoOwnerPtr_;
//- Does face intrude into neighbour cell
mutable autoPtr<PackedBoolList> intrudesIntoNeighbourPtr_;
//- Is the cell concave
mutable autoPtr<PackedBoolList> concaveCellPtr_;
//- Overall count of particles ever created. Never decreases.
mutable label particleCount_;
@ -92,20 +90,11 @@ class Cloud
// Private member functions
//- Helper: is face f1 on wrong side of f0 (both faces are on same cell)
// 0 : convex
// 1 : planar to within planarCosAngle
// 2 : concave
label isIntrudingFace
(
const label cellI,
const label f0,
const label f1
) const;
//- Determine if a particular cell is concave
bool isConcaveCell(const label cellI) const;
//- Analyse mesh for faces not being convex.
// Sets intrudesIntoOwner_,intrudesIntoNeighbour_.
void calcIntrudingFaces() const;
void calcConcaveCells() const;
//- Initialise cloud on IO constructor
void initCloud(const bool checkClass);
@ -229,12 +218,8 @@ public:
return IDLList<ParticleType>::size();
};
//- Per face whether it intrudes into the owner cell.(demand driven)
const PackedBoolList& intrudesIntoOwner() const;
//- Per face whether it intrudes into the neighbour cell. (demand
// driven.
const PackedBoolList& intrudesIntoNeighbour() const;
//- Whether each cell is concave (demand driven data)
const PackedBoolList& concaveCell() const;
// Iterators

5
src/lagrangian/basic/Cloud/CloudIO.C
View File

@ -90,11 +90,6 @@ Foam::Cloud<ParticleType>::Cloud
particleCount_(0)
{
initCloud(checkClass);
// if (concaveCheck_)
// {
// intrudesIntoOwner();
// }
}

405
src/lagrangian/basic/Particle/Particle.C
View File

@ -48,67 +48,10 @@ void Foam::Particle<ParticleType>::findFaces
faceList.clear();
// forAll(faces, i)
// {
// label facei = faces[i];
// scalar lam = lambda(C, position, facei);
// if ((lam > 0) && (lam < 1.0))
// {
// faceList.append(facei);
// }
// }
vector deltaPosition = position - C;
scalar deltaPositionMag = mag(deltaPosition);
scalar lam;
forAll(faces, i)
{
label facei = faces[i];
lam = GREAT;
if (cloud_.concaveCheck_)
{
// Check for intruding face
// ~~~~~~~~~~~~~~~~~~~~~~~~
label isOwner = (celli_ == mesh.faceOwner()[facei]);
if
(
(isOwner && cloud_.intrudesIntoOwner()[facei])
|| (!isOwner && cloud_.intrudesIntoNeighbour()[facei])
)
{
// Concave face. Use exact intersection.
pointHit inter = mesh.faces()[facei].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[facei],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
if (inter.hit())
{
lam = inter.distance()/deltaPositionMag;
}
}
else
{
// Convex face. Can use lambda.
lam = lambda(C, position, facei);
}
}
else
{
lam = lambda(C, position, facei);
}
scalar lam = lambda(C, position, facei);
if ((lam > 0) && (lam < 1.0))
{
@ -145,6 +88,163 @@ void Foam::Particle<ParticleType>::findFaces
}
template<class ParticleType>
void Foam::Particle<ParticleType>::trackToFaceByDecomposition
(
scalar& trackFraction,
const vector& endPosition
)
{
const polyMesh& mesh = cloud_.polyMesh_;
const labelList& faces = mesh.cells()[celli_];
vector deltaPosition = endPosition - position_;
scalar lambdaMin = GREAT;
label workingFace = -1;
scalar hitLambda = GREAT;
forAll(faces, i)
{
label facei = faces[i];
pointHit inter = mesh.faces()[facei].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[facei],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
if (inter.hit())
{
hitLambda = inter.distance();
// There is a problem if the particle is essentially on
// the face, it will flip.
if (hitLambda <= 1.0 && hitLambda < lambdaMin)
{
Info<< origProc_ << " "
<< origId_ << " "
<< position_ << " "
<< endPosition << " "
<< trackFraction << " "
<< stepFraction_ << " "
<< hitLambda << " "
<< inter.hitPoint() << " "
<< celli_ << " "
<< facei << endl;
lambdaMin = hitLambda;
workingFace = facei;
}
}
}
if (workingFace == -1)
{
// inside cell
facei_ = workingFace;
trackFraction = 1.0;
position_ = endPosition;
}
else
{
if (lambdaMin < 0.0)
{
FatalErrorIn("Particle::trackToFaceByDecomposition")
<< "Decomposed tracking failure, lambdaMin = " << lambdaMin
<< nl << abort(FatalError);
}
trackFraction = lambdaMin;
position_ += trackFraction*deltaPosition;
// What about softImpact?
}
}
template<class ParticleType>
void Foam::Particle<ParticleType>::trackToFaceByPlanes
(
scalar& trackFraction,
const vector& endPosition
)
{
facei_ = -1;
DynamicList<label>& faces = cloud_.labels_;
findFaces(endPosition, faces);
if (faces.empty())
{
// inside cell
trackFraction = 1.0;
position_ = endPosition;
}
else
{
// hit face
scalar lambdaMin = GREAT;
if (faces.size() == 1)
{
lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
facei_ = faces[0];
}
else
{
forAll(faces, i)
{
scalar lam =
lambda(position_, endPosition, faces[i], stepFraction_);
if (lam < lambdaMin)
{
lambdaMin = lam;
facei_ = faces[i];
}
}
}
// For warped faces the particle can be 'outside' the cell.
// This will yield a lambda larger than 1, or smaller than 0
// For values < 0, the particle travels away from the cell
// and we don't move the particle, only change cell.
// For values larger than 1, we move the particle to endPosition only.
if (lambdaMin > 0.0)
{
if (lambdaMin <= 1.0)
{
trackFraction = lambdaMin;
position_ += trackFraction*(endPosition - position_);
}
else
{
trackFraction = 1.0;
position_ = endPosition;
}
}
else if (static_cast<ParticleType&>(*this).softImpact())
{
// Soft-sphere particles can travel outside the domain
// but we don't use lambda since this the particle
// is going away from face
trackFraction = 1.0;
position_ = endPosition;
}
}
}
template<class ParticleType>
template<class TrackData>
void Foam::Particle<ParticleType>::prepareForParallelTransfer
@ -275,7 +375,6 @@ Foam::label Foam::Particle<ParticleType>::track
}
template<class ParticleType>
Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
{
@ -283,6 +382,7 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
return track(endPosition, dummyTd);
}
template<class ParticleType>
template<class TrackData>
Foam::scalar Foam::Particle<ParticleType>::trackToFace
@ -293,154 +393,25 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
{
const polyMesh& mesh = cloud_.polyMesh_;
DynamicList<label>& faces = cloud_.labels_;
findFaces(endPosition, faces);
facei_ = -1;
scalar trackFraction = 0.0;
scalar lambdaMin = GREAT;
// Temporary addition to report on infinite loops ->
List<scalar> lambdas(faces.size());
forAll(faces, i)
if (cloud_.concaveCell()[celli_])
{
lambdas[i] = lambda(position_, endPosition, faces[i], stepFraction_);
}
// <- Temporary addition to report on infinite loops
if (faces.empty()) // inside cell
{
trackFraction = 1.0;
position_ = endPosition;
}
else // hit face
{
// scalar lambdaMin = GREAT;
if (faces.size() == 1)
{
lambdaMin = lambda(position_, endPosition, faces[0], stepFraction_);
facei_ = faces[0];
trackToFaceByDecomposition(trackFraction, endPosition);
}
else
{
// // If the particle has to cross more than one cell to reach the
// // endPosition, we check which way to go.
// // If one of the faces is a boundary face and the particle is
// // outside, we choose the boundary face.
// The particle is outside if one of the lambda's is > 1 or < 0
scalar fallBackLambdaMin = GREAT;
label fallBackFacei = -1;
vector deltaPosition = endPosition - position_;
scalar deltaPositionMag = mag(deltaPosition);
scalar hitLambda = GREAT;
forAll(faces, i)
{
scalar lam =
lambda(position_, endPosition, faces[i], stepFraction_);
// Calculate fallback
if (lam < fallBackLambdaMin)
{
fallBackLambdaMin = lam;
fallBackFacei = faces[i];
trackToFaceByPlanes(trackFraction, endPosition);
}
if (cloud_.concaveCheck_)
if (facei_ > -1)
{
// Check for intruding face
// ~~~~~~~~~~~~~~~~~~~~~~~~
// A face was hit
label isOwner = (celli_ == mesh.faceOwner()[faces[i]]);
if (cloud_.internalFace(facei_))
{
// Internal face, change cell
if
(
(isOwner && cloud_.intrudesIntoOwner()[faces[i]])
|| (!isOwner && cloud_.intrudesIntoNeighbour()[faces[i]])
)
{
// Concave face. Use exact intersection.
pointHit inter = mesh.faces()[faces[i]].intersection
(
position_,
deltaPosition,
mesh.faceCentres()[faces[i]],
mesh.points(),
intersection::HALF_RAY,
Cloud<ParticleType>::intersectionTolerance
);
if (inter.hit())
{
hitLambda = inter.distance()/deltaPositionMag;
if (hitLambda < lambdaMin)
{
lambdaMin = hitLambda;
facei_ = faces[i];
}
}
}
else
{
// Convex face. Can use lambda.
if (lam < lambdaMin)
{
lambdaMin = lam;
facei_ = faces[i];
}
}
}
}
// Fall back to the original algorithm if no
// intersection was found
if (facei_ == -1)
{
lambdaMin = fallBackLambdaMin;
facei_ = fallBackFacei;
}
}
bool internalFace = cloud_.internalFace(facei_);
// For warped faces the particle can be 'outside' the cell.
// This will yield a lambda larger than 1, or smaller than 0
// For values < 0, the particle travels away from the cell
// and we don't move the particle, only change cell.
// For values larger than 1, we move the particle to endPosition only.
if (lambdaMin > 0.0)
{
if (lambdaMin <= 1.0)
{
trackFraction = lambdaMin;
position_ += trackFraction*(endPosition - position_);
}
else
{
trackFraction = 1.0;
position_ = endPosition;
}
}
else if (static_cast<ParticleType&>(*this).softImpact())
{
// Soft-sphere particles can travel outside the domain
// but we don't use lambda since this the particle
// is going away from face
trackFraction = 1.0;
position_ = endPosition;
}
// change cell
if (internalFace) // Internal face
{
if (celli_ == mesh.faceOwner()[facei_])
{
celli_ = mesh.faceNeighbour()[facei_];
@ -451,10 +422,8 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
}
else
{
FatalErrorIn
(
"Particle::trackToFace(const vector&, TrackData&)"
)<< "addressing failure" << nl
FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
<< "addressing failure" << nl
<< abort(FatalError);
}
}
@ -526,31 +495,17 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
if (trackFraction < Cloud<ParticleType>::minValidTrackFraction)
{
Pout<< "tracking error "
<< origId_ << " "
<< origProc_ << " "
<< position_ << " "
<< endPosition << " "
<< position_ + 1000*(endPosition - position_) << " "
<< position_ - 1000*(endPosition - position_) << " "
<< lambdaMin << " "
<< trackFraction << " "
<< stepFraction_ << " "
<< celli_ << " "
<< facei_ << " "
<< cloud_.intrudesIntoOwner()[facei_] << " "
<< cloud_.intrudesIntoNeighbour()[facei_] << " "
<< faces << " "
<< lambdas
<< endl;
// if (stepFraction_ < SMALL)
// {
// Pout<< "Intervene to mark face" << endl;
// cloud_.intrudesIntoOwnerPtr_()[facei_] = 1;
// cloud_.intrudesIntoNeighbourPtr_()[facei_] = 1;
// }
// Pout<< "tracking error "
// << origId_ << " "
// << origProc_ << " "
// << position_ << " "
// << endPosition << " "
// << trackFraction << " "
// << stepFraction_ << " "
// << celli_ << " "
// << facei_ << " "
// << cloud_.concaveCell()[celli_]
// << endl;
const point& cc = mesh.cellCentres()[celli_];
position_ +=

17
src/lagrangian/basic/Particle/Particle.H
View File

@ -173,6 +173,23 @@ protected:
DynamicList<label>& faceList
) const;
//- Track to cell face using face decomposition, used for
// concave cells.
void trackToFaceByDecomposition
(
scalar& trackFraction,
const vector& endPosition
);
//- Track to cell face using the infinite planes of the faces.
// The cell *must* be convex.
void trackToFaceByPlanes
(
scalar& trackFraction,
const vector& endPosition
);
// Patch interactions