49 for (label celli = 0; celli < mesh_.
nCells(); ++celli)
51 if (ignoreCells.
test(celli))
59 ignoreCells.
set(celli);
70 const bitSet& ignoreCells
75 bitSet blockedFaces(mesh_.nFaces());
77 const labelList& faceOwn = mesh_.faceOwner();
78 const labelList& faceNei = mesh_.faceNeighbour();
81 for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
88 ignoreCells.test(faceOwn[facei])
89 != ignoreCells.test(faceNei[facei])
92 blockedFaces.
set(facei);
96 for (
const polyPatch&
patch : mesh_.boundaryMesh())
105 const label facei =
patch.start() + patchFacei;
106 if (ignoreCells.test(faceOwn[facei]))
108 blockedFaces.set(facei);
125 bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
129 blockedFaces.clearStorage();
137 forAll(regionColour, celli)
139 if (!ignoreCells.
test(celli))
141 ++nCutsPerRegion[regionColour[celli]];
150 boolList keepRegion(rs.nRegions(),
false);
154 const label largest =
findMax(nCutsPerRegion);
163 keepRegion[largest] =
true;
168 Info<<
"Had " <<
sum(nCutsPerRegion) <<
" cuts, in " 169 << nCutsPerRegion.size() <<
" regions, largest is " 176 forAll(regionColour, celli)
178 if (!keepRegion.test(regionColour[celli]))
191 if (nearestPoints_.empty())
194 <<
"Ignoring nearestPoints - no points provided" <<
nl 200 bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
203 regionSplit rs(mesh_, blockedFaces);
204 blockedFaces.clearStorage();
212 typedef Tuple2<scalar, label> nearInfo;
213 List<nearInfo> nearest(nearPts.size(), nearInfo(GREAT, -1));
221 if (ignoreCells.test(celli))
226 const point& pt = cc[celli];
227 const label regioni = regionColour[celli];
229 ++nCutsPerRegion[regioni];
232 for (nearInfo& near : nearest)
234 const scalar distSqr =
magSqr(nearPts[pointi] - pt);
237 if (distSqr < near.first())
239 near.first() = distSqr;
240 near.second() = regioni;
254 boolList keepRegion(rs.nRegions(),
false);
258 const label largest =
findMax(nCutsPerRegion);
260 for (
const nearInfo& near : nearest)
262 const scalar distSqr = near.first();
263 const label regioni = near.second();
265 if (regioni != -1 && distSqr < maxDistanceSqr_)
267 keepRegion[regioni] =
true;
273 Info<<
"Had " <<
sum(nCutsPerRegion) <<
" cuts, in " 274 << nCutsPerRegion.size() <<
" regions, largest is " 277 Info<<
"nearestPoints (max distance = " 278 <<
sqrt(maxDistanceSqr_) <<
")" <<
nl;
282 const scalar dist =
sqrt(nearest[pointi].first());
283 const label regioni = nearest[pointi].second();
285 Info<<
" " << nearPts[pointi] <<
" region " 286 << regioni <<
" distance " 289 if (!keepRegion.test(regioni))
300 forAll(regionColour, celli)
302 if (!keepRegion.test(regionColour[celli]))
304 ignoreCells.set(celli);
315 const searchableSurface& geom = geometryPtr_();
318 const pointField& fc = surface_.faceCentres();
321 bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
324 regionSplit rs(mesh_, blockedFaces);
325 blockedFaces.clearStorage();
333 List<pointIndexHit> nearest;
341 calcAbsoluteDistance(faceDistance, fc, nearest);
350 const label celli = meshCells_[facei];
351 const label regioni = regionColour[celli];
353 const scalar faceArea = surface_[facei].mag(surface_.points());
354 distRegion[regioni] += (faceDistance[facei] * faceArea);
355 areaRegion[regioni] += (faceArea);
363 forAll(distRegion, regioni)
365 distRegion[regioni] /= (areaRegion[regioni] + VSMALL);
374 bitSet acceptFaces(fc.size());
379 const label celli = meshCells_[facei];
380 const label regioni = regionColour[celli];
386 if (absProximity_ < distRegion[regioni])
392 acceptFaces.set(facei);
399 const fileName
outputName(surfaceName() +
"-region-proximity-filter");
403 surfaceWriters::vtkWriter
writer 410 writer.write(
"absolute-distance", faceDistance);
415 ListOps::create<label>
418 [&](
const label celli){
return regionColour[celli]; }
421 writer.write(
"face-region", faceRegion);
426 ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
428 writer.write(
"filter-state", faceFilterState);
440 forAll(faceDistance, facei)
442 if (absProximity_ < faceDistance[facei])
448 acceptFaces.set(facei);
458 surface_.subsetMesh(acceptFaces, pointMap,
faceMap)
469 const searchableSurface& geom = geometryPtr_();
472 const pointField& fc = surface_.faceCentres();
478 List<pointIndexHit> nearest;
486 calcAbsoluteDistance(faceDistance, fc, nearest);
492 geom.getNormal(nearest, norms);
494 faceNormalDistance.
resize(fc.size());
513 bitSet acceptFaces(fc.size());
517 forAll(faceDistance, facei)
519 if (absProximity_ < faceDistance[facei])
524 Pout<<
"trim reject: " 525 << faceDistance[facei] <<
nl;
530 acceptFaces.set(facei);
538 const fileName
outputName(surfaceName() +
"-face-proximity-filter");
542 surfaceWriters::vtkWriter
writer 549 writer.write(
"absolute-distance", faceDistance);
550 writer.write(
"normal-distance", faceNormalDistance);
555 ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
557 writer.write(
"filter-state", faceFilterState);
566 surface_.subsetMesh(acceptFaces, pointMap,
faceMap)
570 meshCells_ = UIndirectList<label>(meshCells_,
faceMap)();
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
label findMax(const ListType &input, label start=0)
Linear search for the index of the max element, similar to std::max_element but for lists and returns...
scalar diff(const triad &A, const triad &B)
Return a quantity of the difference between two triads.
cutType getCellCutType(const label celli) const
Cell cut for an individual cell, with special handling for TETCUT and SPHERE cuts.
Field< label > labelField
Specialisation of Field<T> for label.
bitSet filterPrepareRegionSplit(const bitSet &ignoreCells) const
Prepare blockedFaces for region split.
void set(const bitSet &bitset)
Set specified bits from another bitset.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void filterKeepNearestRegions(bitSet &ignoreCells) const
Keep region(s) closest to the nearest points.
void resize(const label len)
Adjust allocated size of list.
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
constexpr char nl
The newline '\n' character (0x0a)
bool refineBlockedCells(bitSet &ignoreCells, const isoSurfaceBase &isoContext) const
Re-filter the blocked cells based on the anticipated cuts.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
static void broadcast(Type &value, const label comm=UPstream::worldComm)
Broadcast content (contiguous or non-contiguous) to all processes in communicator.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
Low-level components common to various iso-surface algorithms.
#define forAll(list, i)
Loop across all elements in list.
bool returnReduceAnd(const bool value, const label comm=UPstream::worldComm)
Perform logical (and) MPI Allreduce on a copy. Uses UPstream::reduceAnd.
word outputName("finiteArea-edges.obj")
vectorField pointField
pointField is a vectorField.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
MeshedSurface< face > meshedSurface
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
int debug
Static debugging option.
vtk::internalMeshWriter writer(topoMesh, topoCells, vtk::formatType::INLINE_ASCII, runTime.path()/"blockTopology")
void filterKeepLargestRegion(bitSet &ignoreCells) const
Keep region with the most cuts (after region split)
vector point
Point is a vector.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
void filterFaceProximity()
Adjust extracted iso-surface to remove far faces.
A List with indirect addressing. Like IndirectList but does not store addressing. ...
static bool master(const label communicator=worldComm)
Am I the master rank.
const std::string patch
OpenFOAM patch number as a std::string.
void transfer(pointField &pointLst, List< Face > &faceLst)
Transfer the components.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static void listCombineGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
List< label > labelList
A List of labels.
List< bool > boolList
A List of bools.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
FlatOutput::OutputAdaptor< Container, Delimiters > flatOutput(const Container &obj, Delimiters delim)
Global flatOutput() function with specified output delimiters.
void filterRegionProximity(bitSet &ignoreCells) const
Remove region(s) that have far faces.
static constexpr const zero Zero
Global zero (0)