exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Immersed boundary robustness improvements

Browse source
This commit is contained in:
Hrvoje Jasak 2019-03-13 12:19:27 +00:00
parent 5fce429505
commit 6fe4ba3aa3
7 changed files with 335 additions and 210 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

198
src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.C
View file

@ -53,14 +53,6 @@ namespace Foam
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::debug::tolerancesSwitch
Foam::immersedBoundaryPolyPatch::liveFactor_
(
"immersedBoundaryLiveFactor",
1e-6
);
const Foam::debug::tolerancesSwitch
Foam::immersedBoundaryPolyPatch::spanFactor_
(
@ -261,7 +253,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
// If the nearest triangle cannot be found within span than this is
// most probably a tri surface search error. Mark unknown and check
// later. (IG 22/Nov/2018)
if(tss.nearest(C[cellI], span/spanFactor_()).index() == -1)
if (tss.nearest(C[cellI], span/spanFactor_()).index() == -1)
{
intersectedCell[cellI] = immersedPoly::UNKNOWN;
}
@ -329,7 +321,7 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
{
// There are either no wet or dry negbours or there are both.
// This should not be possible. NOTE: the check is not
// parallelised and this can theoretically lead to failiures in
// parallelised and this can theoretically lead to failures in
// strange arrangaments.
// Issue a warning, mark CUT and hope for the best.
// (IG 22/Nov/2018)
@ -433,60 +425,50 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
}
else
{
// Real intesection. Check cut
const scalar cellFactor = cutCell.wetVolume()/V[cellI];
// True intersection. Cut the cell and store all
// derived data
if (cellFactor < liveFactor_())
// Note: volumetric check is not allowed because true
// intersection guarantees that the faces of the cell
// have been cut. Therefore, the cell MUST be an IB cell.
// If the cut is invalid, Marooney Maneouvre shall correct
// the error in sum(Sf). HJ, 12/Mar/2019
// Store ibFace with local points. Points merge will
// take place later
const face& cutFace = cutCell.faces()[0];
const pointField& cutPoints = cutCell.points();
// Collect the renumbered face, using the point labels
// from the unmergedPoints list
face renumberedFace(cutFace.size());
// Insert points and renumber the face
forAll (cutFace, cpI)
{
// Thin cut: cell is dry
intersectedCell[cellI] = immersedPoly::DRY;
unmergedPoints.append(cutPoints[cutFace[cpI]]);
renumberedFace[cpI] = nIbPoints;
nIbPoints++;
}
else if (cellFactor > (1 - liveFactor_()))
{
// Thick cut: cell is wet
intersectedCell[cellI] = immersedPoly::WET;
}
else
{
// True intersection. Cut the cell and store all
// derived data
// Store ibFace with local points. Points merge will
// take place later
const face& cutFace = cutCell.faces()[0];
// Record the face
unmergedFaces[nIbCells] = renumberedFace;
const pointField& cutPoints = cutCell.points();
// Collect cut cell index
ibCells[nIbCells] = cellI;
// Collect the renumbered face, using the point labels
// from the unmergedPoints list
face renumberedFace(cutFace.size());
// Record the live centre
ibCellCentres[nIbCells] = cutCell.wetVolumeCentre();
// Insert points and renumber the face
forAll (cutFace, cpI)
{
unmergedPoints.append(cutPoints[cutFace[cpI]]);
renumberedFace[cpI] = nIbPoints;
nIbPoints++;
}
// Record the live volume
ibCellVolumes[nIbCells] = cutCell.wetVolume();
// Record the face
unmergedFaces[nIbCells] = renumberedFace;
// Record the nearest triangle to the face centre
nearestTri[nIbCells] =
tss.nearest(cutFace.centre(cutPoints), span).index();
// Collect cut cell index
ibCells[nIbCells] = cellI;
// Record the live centre
ibCellCentres[nIbCells] = cutCell.wetVolumeCentre();
// Record the live volume
ibCellVolumes[nIbCells] = cutCell.wetVolume();
// Record the nearest triangle to the face centre
nearestTri[nIbCells] =
tss.nearest(cutFace.centre(cutPoints), span).index();
nIbCells++;
}
nIbCells++;
}
}
}
@ -1079,36 +1061,24 @@ void Foam::immersedBoundaryPolyPatch::calcImmersedBoundary() const
}
else
{
// Real intesection. Check cut
// Real intesection. Check cut. Rejection on thin cut is
// performed by ImmersedFace. HJ, 13/Mar/2019
const scalar faceFactor =
cutFace.wetAreaMag()/mag(S[faceI]);
if (faceFactor < liveFactor_())
{
// Thin cut: face is dry
intersectedFace[faceI] = immersedPoly::DRY;
}
else if (faceFactor > (1 - liveFactor_()))
{
// Thick cut: face is wet
intersectedFace[faceI] = immersedPoly::WET;
}
else
{
// True intersection. Collect data
intersectedFace[faceI] = immersedPoly::CUT;
// True intersection. Collect data
intersectedFace[faceI] = immersedPoly::CUT;
// Get intersected face index
ibFaces[nIbFaces] = faceI;
// Get intersected face index
ibFaces[nIbFaces] = faceI;
// Get wet centre
ibFaceCentres[nIbFaces] = cutFace.wetAreaCentre();
// Get wet centre
ibFaceCentres[nIbFaces] = cutFace.wetAreaCentre();
// Get wet area, preserving original normal direction
ibFaceAreas[nIbFaces] = faceFactor*S[faceI];
// Get wet area, preserving original normal direction
ibFaceAreas[nIbFaces] = faceFactor*S[faceI];
nIbFaces++;
}
nIbFaces++;
}
}
}
@ -1260,41 +1230,55 @@ void Foam::immersedBoundaryPolyPatch::calcCorrectedGeometry() const
// by the Marooney Maneouvre, where the face sum imbalance is compensated
// in the cut face. HJ, 11/Dec/2017
vectorField sumSf(mesh.nCells(), vector::zero);
// Get face addressing
const labelList& owner = mesh.faceOwner();
const labelList& neighbour = mesh.faceNeighbour();
forAll (owner, faceI)
{
sumSf[owner[faceI]] += Sf[faceI];
}
forAll (neighbour, faceI)
{
sumSf[neighbour[faceI]] -= Sf[faceI];
}
forAll (cutCells, cutCellI)
{
const label ccc = cutCells[cutCellI];
if
(
mag(sumSf[ccc]) > 1e-12
// The below criteria is still not firmly established, potential
// place for cutting errors (IG 19/Dec/2018)
&& mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI] > 1e-6
)
{
// Info<< "Marooney Maneouvre for cell " << ccc
// << " error: " << mag(sumSf[ccc] + ibSf[cutCellI]) << " "
// << mag(sumSf[ccc] + ibSf[cutCellI])/cutCellVolumes[cutCellI]
// << " " << sumSf[ccc] << " "
// << " V: " << cutCellVolumes[cutCellI]
// << " Sf: " << ibSf[cutCellI] << endl;
ibSf[cutCellI] = -sumSf[ccc];
// Calculate sum Sf and sumMagSf for the cell
const cell& curCell = mesh.cells()[ccc];
vector curSumSf = vector::zero;
scalar curSumMagSf = 0;
// Collect from regular faces
forAll (curCell, cfI)
{
const vector& curSf = Sf[curCell[cfI]];
// Check owner/neighbour
if (owner[curCell[cfI]] == ccc)
{
curSumSf += curSf;
}
else
{
curSumSf -= curSf;
}
curSumMagSf += mag(curSf);
}
// Add cut face only into mag. The second part is handled in the
// if-statement
curSumMagSf += mag(ibSf[cutCellI]);
// Adjustment is peformed when the openness is greater than a certain
// fraction of surface area. Criterion by IG, 13/Mar/2019
// Switched to using absolute check from primitiveMeshCheck.
// HJ, 13/Mar/2019
// if (mag(curSumSf + ibSf[cutCellI]) > 1e-6*curSumMagSf)
if (mag(curSumSf + ibSf[cutCellI]) > primitiveMesh::closedThreshold_)
{
Info<< "Marooney Maneouvre for cell " << ccc
<< " error: " << mag(curSumSf + ibSf[cutCellI]) << " "
<< " S: " << curSumMagSf
<< " V: " << cutCellVolumes[cutCellI]
<< " Sf: " << ibSf[cutCellI] << endl;
// Create IB face to ideally close the cell
ibSf[cutCellI] = -curSumSf;
}
}
}

5
src/immersedBoundary/immersedBoundary/immersedBoundaryPolyPatch/immersedBoundaryPolyPatch.H
View file

@ -86,11 +86,6 @@ class immersedBoundaryPolyPatch
// Static data
//- Live volume factor
// Fraction of face/cell size within the IB for which the face/cell is
// considered to be live
static const debug::tolerancesSwitch liveFactor_;
//- Search span factor
// Factor used to multiply the span of the cell to get the tri-seach
// span

197
src/immersedBoundary/immersedBoundary/immersedPoly/ImmersedCell.C
View file

@ -74,7 +74,7 @@ void Foam::ImmersedCell<Distance>::getBase
template<class Distance>
Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
void Foam::ImmersedCell<Distance>::insertIntersectionPoints()
{
// Get list of edges
const edgeList& edges = this->edges();
@ -99,16 +99,19 @@ Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
const scalar edgeLength = mag(points_[end] - points_[start]);
// Check if there is a legitimate cut to be found
// Note: synced tolerances in ImmersedCell and ImmersedFace
// HJ, 13/Mar/2019
if
(
depth_[start]*depth_[end] < 0
&& edgeLength > SMALL
&& mag(depth_[start]) > absTol_
&& mag(depth_[end]) > absTol_
&& mag(depth_[start]) > edgeLength*immersedPoly::tolerance_()
&& mag(depth_[end]) > edgeLength*immersedPoly::tolerance_()
)
{
// Prepare a new point to insert
point cutPoint;
scalar depthAtCut = 0;
// Intersection is along the edge length (pf[end] - pf[start])
// times the ratio of the depth at start and the difference
@ -131,14 +134,14 @@ Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
scalar d1 = depth_[end];
// Convergence criterion is the depth at newP
scalar depthAtCut = dist_.distance(cutPoint);
depthAtCut = dist_.distance(cutPoint);
// initialize loop counter
label iters = 0;
while
(
(mag(depthAtCut) > absTol_)
(mag(depthAtCut) > immersedPoly::tolerance_())
&& (iters < immersedPoly::nIter_())
)
{
@ -184,7 +187,7 @@ Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
// Count points added to new face
label nfp = 0;
// Loop throuhg old face. If this edge is found, add the
// Loop through old face. If this edge is found, add the
// cut point label into the edge
forAll (oldFace, fpI)
{
@ -210,7 +213,8 @@ Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
// Debug: check if point insertion was successful
if (nfp < newFace.size())
{
FatalErrorIn("badInsertion")
FatalErrorInFunction
<< "badInsertion"
<< abort(FatalError);
}
@ -228,41 +232,42 @@ Foam::label Foam::ImmersedCell<Distance>::insertIntersectionPoints()
// Extra depths are all zero
depth_.setSize(nPoints);
// For all cut points set depth to exactly zero
for (label i = oldSize; i < depth_.size(); i++)
{
depth_[i] = 0;
}
// Check for successful intersection: more than 3 added points
// can form an internal face
return nPoints - oldSize;
}
template<class Distance>
Foam::face Foam::ImmersedCell<Distance>::createInternalFace
(
const label nIntersections
) const
Foam::face Foam::ImmersedCell<Distance>::createInternalFace() const
{
// Sanity check: Do we have at least 3 intersection points?
if (nIntersections < 3)
// Declare internal face with mixed-up point ordering
face unorderedInternalFace(points_.size());
// Collect all points with zero distance to surface
label nPif = 0;
forAll (depth_, pointI)
{
FatalErrorIn
(
"ImmersedCell::createInternalFace(const label nIntersections)"
) << "Less than 3 intersection points between cell and free surface"
<< abort(FatalError);
if (mag(depth_[pointI]) < absTol_)
{
// Found point on zero plane
unorderedInternalFace[nPif] = pointI;
nPif++;
}
}
// Declare internal face with mixed-up point ordering
face unorderedInternalFace(nIntersections);
unorderedInternalFace.setSize(nPif);
forAll (unorderedInternalFace, i)
// Sanity check: Do we have at least 3 points at zero distance?
if (nPif < 3)
{
label pointID = points_.size() - nIntersections + i;
unorderedInternalFace[i] = pointID;
FatalErrorInFunction
<< "Less than 3 intersection points in cell on free surface." << nl
<< "depth: " << depth_
<< abort(FatalError);
}
// Order points, so that they form a polygon
@ -308,7 +313,8 @@ Foam::face Foam::ImmersedCell<Distance>::createInternalFace
if (pointID == -1)
{
// All intersection points are colinear
FatalErrorIn("Colinear points in cut")
FatalErrorInFunction
<< "Colinear points in cut"
<< abort(FatalError);
}
@ -373,12 +379,10 @@ Foam::face Foam::ImmersedCell<Distance>::createInternalFace
}
}
if(nUndecided == depth_.size())
if (nUndecided == depth_.size())
{
FatalErrorIn
(
"ImmersedCell::createInternalFace(const label nIntersections)"
) << "All points lay on the tri surface, zero volume cell?"
FatalErrorInFunction
<< "All points lay on the tri surface, zero volume cell?"
<< nl << "Points: " << points_
<< abort(FatalError);
}
@ -427,6 +431,7 @@ Foam::ImmersedCell<Distance>::ImmersedCell
absTol_(0),
isAllWet_(false),
isAllDry_(false),
isBadCut_(false),
// Initialize points_ with points from cell
points_(mesh_.cells()[cellID_].points(mesh_.faces(), mesh_.points())),
faces_(),
@ -435,14 +440,14 @@ Foam::ImmersedCell<Distance>::ImmersedCell
faceNeighbour_(),
depth_(dist.distance(points_))
{
const cell& origCell = mesh.cells()[cellID];
const cell& origCell = mesh_.cells()[cellID];
// Build a valid 1-cell mesh in local addressing
// Create hash table that maps points on global mesh to local point list
HashTable<label, label, Hash<label> > pointMapTable(points_.size());
labelList origCellPointLabels = origCell.labels(mesh.faces());
labelList origCellPointLabels = origCell.labels(mesh_.faces());
forAll (points_, pointI)
{
@ -458,11 +463,11 @@ Foam::ImmersedCell<Distance>::ImmersedCell
forAll (origCell, faceI)
{
// Get old point list of faceI
face origFace(mesh.faces()[origCell[faceI]]);
face origFace(mesh_.faces()[origCell[faceI]]);
// Make sure that all faces point outward,
// since they are going to be outside cells
if (!(mesh.faceOwner()[origCell[faceI]] == cellID))
if (!(mesh_.faceOwner()[origCell[faceI]] == cellID))
{
// Cell is not owner of face, revert face orientation
// for the use in a 1-cell mesh
@ -536,9 +541,11 @@ Foam::ImmersedCell<Distance>::ImmersedCell
// Created expanded point and face lists
// Insert intersection points and adjust depth for intersections
// This will add further points into the intersected face
// Depth at intersectin will be zero. HJ, 5/Dec/2017
label nIntersections = insertIntersectionPoints();
// This will add further points into the intersected face if needed
// Depth at intersection will be zero. HJ, 5/Dec/2017
// Note that it is possible to have the cut face even if no new points
// have been introduced. HJ, 13/Mar/2019
insertIntersectionPoints();
// Update primitiveMesh parameters
this->reset
@ -558,7 +565,18 @@ Foam::ImmersedCell<Distance>::ImmersedCell
// At this point, a 1-cell mesh with faces enriched for intersections
// is valid. HJ, 5/Dec/2017
// Recheck, if there has been a successful cut at all
// Check if there has been a successful cut at all
// For a good cut there should be at least 3 points at zero level
label nIntersections = 0;
forAll (depth_, pointI)
{
if (mag(depth_[pointI]) < absTol_)
{
nIntersections++;
}
}
if (nIntersections < 3)
{
// Check if cell centre is wet or dry, depending on greatest distance
@ -581,7 +599,6 @@ Foam::ImmersedCell<Distance>::ImmersedCell
return;
}
}
// From here on, there exists a valid intersection
// Resize the face list. Each face can be split into two, with one
@ -599,7 +616,7 @@ Foam::ImmersedCell<Distance>::ImmersedCell
// with one point or edge, insert internal face that
// connects all intersection points
// create internal face, which gets inserted at front of faces_ list
faces_[0] = createInternalFace(nIntersections);
faces_[0] = createInternalFace();
// Internal face points out of the wet cell. Make the wet cell its owner
faceOwner_[0] = WET;
@ -611,9 +628,12 @@ Foam::ImmersedCell<Distance>::ImmersedCell
// For all faces with inserted points, do face splitting
forAll (enrichedFaces, oldFaceI)
{
const face& oldFace = mesh.faces()[origCell[oldFaceI]];
const face& oldFace = mesh_.faces()[origCell[oldFaceI]];
const face& newFace = enrichedFaces[oldFaceI];
// Calculate old face area locally to avoid triggering polyMesh
const scalar oldFaceArea = mag(mesh_.faceAreas()[origCell[oldFaceI]]);
// If a face has been modified, it will have extra points
if (newFace.size() != oldFace.size())
{
@ -637,13 +657,13 @@ Foam::ImmersedCell<Distance>::ImmersedCell
dryFace[nDry] = newFace[pointI];
nDry++;
}
else if (depth_[newFace[pointI]] < 0)
else if (depth_[newFace[pointI]] < -absTol_)
{
// Point is submerged, add to wetFace
wetFace[nWet] = newFace[pointI];
nWet++;
}
else
else // depth_[newFace[pointI]] > absTol_
{
// Otherwise point must be dry, add to dryFace
dryFace[nDry] = newFace[pointI];
@ -660,6 +680,23 @@ Foam::ImmersedCell<Distance>::ImmersedCell
faceOwner_[nFaces] = WET;
nFaces++;
// Check for bad wet face cut
if
(
wetFace.mag(points_)
> (1 + immersedPoly::badCutFactor_())*oldFaceArea
)
{
// Wet face area is greater than original face area
// This is a bad cut
Info<< "Bad cell face cut: wet = ("
<< wetFace.mag(points_) << " "
<< oldFaceArea
<< ")" << endl;
isBadCut_ = true;
}
}
if (nDry >= 3)
@ -670,6 +707,23 @@ Foam::ImmersedCell<Distance>::ImmersedCell
faceOwner_[nFaces] = DRY;
nFaces++;
// Check for bad dry face cut
if
(
dryFace.mag(points_)
> (1 + immersedPoly::badCutFactor_())*oldFaceArea
)
{
// Dry face area is greater than original face area
// This is a bad cut
Info<< "Bad cell face cut: dry = ("
<< dryFace.mag(points_) << " "
<< oldFaceArea
<< ")" << endl;
isBadCut_ = true;
}
}
}
else
@ -685,7 +739,13 @@ Foam::ImmersedCell<Distance>::ImmersedCell
// Create face depth distance as a subset
scalarField faceDepth(depth_, newFace);
if (mag(min(faceDepth)) > mag(max(faceDepth)))
// Since the face has not been cut, all faceDepth should have the
// same sign. Otherwise, the face should straddle the immersed
// surface. Check on minimum.
// Note: this is a very precise check on purpose: there is no cut
// and the face belongs either to a wet cell or a dry cell
// HJ, 12/Mar/2019
if (min(faceDepth) < scalar(0))
{
// Negative distance: wet face
faceOwner_[nFaces] = WET;
@ -721,20 +781,45 @@ Foam::ImmersedCell<Distance>::ImmersedCell
<< "depth: " << depth_ << endl;
# endif
// Recheck, if the cut cell has significant volume. If not, reset it
scalar wetCut = cellVolumes()[WET]/mesh_.cellVolumes()[cellID_];
const scalar oldCellVolume = mesh_.cellVolumes()[cellID_];
if (wetCut > 1 - immersedPoly::tolerance_())
// Note: is it legal to cut a zero volume cell? HJ, 11/Mar/2019
scalar wetCut = cellVolumes()[WET]/oldCellVolume;
scalar dryCut = cellVolumes()[DRY]/oldCellVolume;
// Check for bad cell cut based on volume
if
(
wetCut < -immersedPoly::badCutFactor_()
|| wetCut > (1 + immersedPoly::badCutFactor_())
|| dryCut < -immersedPoly::badCutFactor_()
|| dryCut > (1 + immersedPoly::badCutFactor_())
)
{
// The cell is practically wet
isAllWet_ = true;
isBadCut_ = true;
}
if (wetCut < immersedPoly::tolerance_())
if (isBadCut_)
{
// The cell is practically dry
isAllDry_ = true;
Info<< "Bad cell cut: volume = (" << wetCut << " " << dryCut
<< ") = " << wetCut + dryCut << nl
// << "Points: " << nl << this->points() << nl
// << "Faces: " << nl << this->faces() << nl
// << "Owner: " << nl << this->faceOwner() << nl
// << "Neighbour: " << nl << this->faceNeighbour() << nl
// << "Cut (wet dry) = (" << isAllWet_ << " " << isAllDry_ << ")"
<< endl;
}
// Correction on cutting is not allowed, as it results in an open cell
// if faces are cut and the cell is not.
// Previous check confirmed more than 3 valid cut points in the cell,
// which means that some of the faces were cut.
// Cutting tolerances for the cell and face have been adjusted to make sure
// identical cut has been produced.
// HJ, 11/Mar/2019
}

47
src/immersedBoundary/immersedBoundary/immersedPoly/ImmersedCell.H
View file

@ -35,6 +35,21 @@ Implementation
The class builds a primitiveMesh with the local point, face and cell list
which is then manipulated as needed.
Note
The new immersed boundary algorithm requires that a valid intersection
is calculated for the cell even if no such intersection exists, eg there
are only 1 or 2 faces protruding through the surface.
The reason for this is thta due to the Marooney Maneouvre correction,
the cell cut status is determined based on the face cut status:
if there are properly cut faces in the cell, a valid cell intersection
needs to be provided.
It is expected that such errors will happen only under extreme
circumstances.
Further, the ImmersedCell class is enhanced by the capacity to detect
bad intersections, eg. two detached corners of a single cells are cut off.
SourceFiles
ImmersedCell.C
@ -79,7 +94,7 @@ class ImmersedCell
//- Point field with original cell points and intersection points
pointField cellPointsAndIntersections_;
//- Minimum edge length in cell
//- Edge length tolerance based on minimum edge length in cell
scalar absTol_;
//- All wet
@ -88,6 +103,9 @@ class ImmersedCell
//- All dry
bool isAllDry_;
//- Bad cut
bool isBadCut_;
// Primitive mesh data for two cell mesh
@ -113,13 +131,12 @@ class ImmersedCell
//- Copied from class geomCellLooper
void getBase(const vector& n, vector& e0, vector& e1) const;
//- Create points where the water surface intersects face edges
// Return number of inserted points
// and if so, insert additional points into points_ and adjust depth
label insertIntersectionPoints();
//- Create points where the surface intersects face edges.
// Insert additional points into points_ and adjust depth
void insertIntersectionPoints();
//- Make an internal face out of the intersection points
face createInternalFace(const label nIntersections) const;
face createInternalFace() const;
public:
@ -186,6 +203,12 @@ public:
return isAllDry_;
}
//- Does the cell have a bad cut
inline bool isBadCut() const
{
return isBadCut_;
}
//- Return center of the cell's wet part
inline point wetVolumeCentre() const
{
@ -212,7 +235,17 @@ public:
}
else
{
return this->cellVolumes()[WET];
// Handle possible bad cut: cut volume must be between zero
// and original cell volume
return Foam::max
(
scalar(0),
Foam::min
(
mesh_.cellVolumes()[cellID_],
this->cellVolumes()[WET]
)
);
}
}
};

62
src/immersedBoundary/immersedBoundary/immersedPoly/ImmersedFace.C
View file

@ -65,12 +65,14 @@ void Foam::ImmersedFace<Distance>::createSubfaces
const scalar edgeLength = curEdge.mag(localPoints);
// Check if there is a legitimate cut to be found
// Note: synced tolerances in ImmersedCell and ImmersedFace
// HJ, 13/Mar/2019
if
(
depth[start]*depth[end] < 0
&& edgeLength > SMALL
&& mag(depth[start]) > edgeLength*tolerance_()
&& mag(depth[end]) > edgeLength*tolerance_()
&& mag(depth[start]) > edgeLength*immersedPoly::tolerance_()
&& mag(depth[end]) > edgeLength*immersedPoly::tolerance_()
)
{
// Prepare a new point to insert and determine its location
@ -177,8 +179,8 @@ void Foam::ImmersedFace<Distance>::createSubfaces
// Analyse new depth
// For each point, determine if it is submerged(=-1), dry(=1) or
// on the surface (=0)
// For each point, determine if it is submerged( = -1), dry( = 1) or
// on the surface ( = 0)
labelField isSubmerged(newDepth.size());
forAll (newDepth, pointI)
@ -270,19 +272,19 @@ void Foam::ImmersedFace<Distance>::createSubfaces
{
drySubface_.setSize(nDry);
// Check area: if it is very small, reset the face
if
(
drySubface_.mag(facePointsAndIntersections_)
< immersedPoly::tolerance_()*
localFace.mag(facePointsAndIntersections_)
)
{
// The face is practically wet
isAllWet_ = true;
// // Check area: if it is very small, reset the face
// if
// (
// drySubface_.mag(facePointsAndIntersections_)
// < immersedPoly::liveFactor_()*
// localFace.mag(facePointsAndIntersections_)
// )
// {
// // The face is practically wet
// isAllWet_ = true;
drySubface_.clear();
}
// drySubface_.clear();
// }
}
if (nWet < 3)
@ -297,21 +299,21 @@ void Foam::ImmersedFace<Distance>::createSubfaces
{
wetSubface_.setSize(nWet);
// Check area: if it is very small, reset the face
// Note: must use relative tolerance, ie divide with original face
// area magnitude. HJ, 30/Aug/2018
if
(
wetSubface_.mag(facePointsAndIntersections_)
< immersedPoly::tolerance_()*
localFace.mag(facePointsAndIntersections_)
)
{
// The face is practically dry
isAllDry_ = true;
// // Check area: if it is very small, reset the face
// // Note: must use relative tolerance, ie divide with original face
// // area magnitude. HJ, 30/Aug/2018
// if
// (
// wetSubface_.mag(facePointsAndIntersections_)
// < immersedPoly::liveFactor_()*
// localFace.mag(facePointsAndIntersections_)
// )
// {
// // The face is practically dry
// isAllDry_ = true;
wetSubface_.clear();
}
// wetSubface_.clear();
// }
}
}
}

20
src/immersedBoundary/immersedBoundary/immersedPoly/immersedPoly.C
View file

@ -31,7 +31,7 @@ const Foam::debug::optimisationSwitch
Foam::immersedPoly::nIter_
(
"immersedPolyNIter",
10
3
);
@ -39,7 +39,23 @@ const Foam::debug::tolerancesSwitch
Foam::immersedPoly::tolerance_
(
"immersedPolyTolerance",
1e-7
1e-4
);
const Foam::debug::tolerancesSwitch
Foam::immersedPoly::liveFactor_
(
"immersedPolyLiveFactor",
1e-6
);
const Foam::debug::tolerancesSwitch
Foam::immersedPoly::badCutFactor_
(
"immersedPolyBadCutFactor",
0.01
);

16
src/immersedBoundary/immersedBoundary/immersedPoly/immersedPoly.H
View file

@ -26,7 +26,8 @@ Class
Description
immersedPoly holds basic data for immersed face/cell objects
Wet side is "inside", indicated by negative distance to surface
WET side is "inside", indicated by negative distance to surface.
Positive distance to surface is DRY.
SourceFiles
immersedPoly.C
@ -67,9 +68,18 @@ public:
//- Number of iterations in the iterative intersection
static const debug::optimisationSwitch nIter_;
//- Tolerance for considering a point on surface. Used as a fraction
// of edge length
//- Tolerance for considering a point on surface.
// Used as a fraction of edge length
static const debug::tolerancesSwitch tolerance_;
//- Live factor
// Fraction of face/cell size within the IB for which the face/cell is
// considered to be live
static const debug::tolerancesSwitch liveFactor_;
//- Bad cut factor
// Fraction of face/cell size error that declares a bad cut
static const debug::tolerancesSwitch badCutFactor_;
};