52 const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1
e-3;
56 Foam::scalar Foam::geomCellLooper::snapTol_ = 0.1;
61 Foam::scalar Foam::geomCellLooper::minEdgeLen(
const label vertI)
const 63 scalar minLen = GREAT;
69 const edge&
e =
mesh().
edges()[pEdges[pEdgeI]];
79 bool Foam::geomCellLooper::cutEdge
81 const plane& cutPlane,
90 scalar
s = cutPlane.normalIntersect(
pts[
e.start()],
e.vec(
pts));
92 if ((
s > -pointEqualTol_) && (
s < 1 + pointEqualTol_))
109 Foam::label Foam::geomCellLooper::snapToVert
122 else if (weight > (1-tol))
137 scalar nComp =
n & base;
139 if (
mag(nComp) > 0.8)
148 if (
mag(nComp) > 0.8)
173 bool Foam::geomCellLooper::edgeEndsCut
179 label edgeI =
getEdge(loop[index]);
183 const label prevCut = loop[loop.rcIndex(index)];
184 const label nextCut = loop[loop.fcIndex(index)];
186 if (!isEdge(prevCut) && !isEdge(nextCut))
190 label v0 = getVertex(prevCut);
191 label v1 = getVertex(nextCut);
195 (v0 ==
e[0] && v1 ==
e[1])
196 || (v0 ==
e[1] && v1 ==
e[0])
208 Foam::geomCellLooper::geomCellLooper(
const polyMesh&
mesh)
231 plane(
mesh().cellCentres()[celli], refDir),
244 const plane& cutPlane,
281 label edgeI = cellEdges[i];
283 const edge&
e = edges[edgeI];
285 bool useStart =
false;
293 if (checkedPoints.insert(
e.start()))
295 const scalar typLen = pointEqualTol_ * minEdgeLen(
e.start());
301 localLoop.
append(vertToEVert(
e.start()));
302 localLoopWeights.append(-GREAT);
308 if (checkedPoints.insert(
e.end()))
310 const scalar typLen = pointEqualTol_ * minEdgeLen(
e.end());
316 localLoop.append(vertToEVert(
e.end()));
317 localLoopWeights.append(-GREAT);
327 if (!useEnd && !useStart)
333 if (cutEdge(cutPlane, edgeI, cutWeight))
336 label cutVertI = snapToVert(snapTol_, edgeI, cutWeight);
341 localLoop.append(edgeToEVert(edgeI));
342 localLoopWeights.append(cutWeight);
349 label cut = vertToEVert(cutVertI);
351 if (!localLoop.found(cut))
353 localLoop.append(vertToEVert(cutVertI));
354 localLoopWeights.append(-GREAT);
361 if (localLoop.size() <= 2)
367 localLoopWeights.shrink();
376 loopPoints[i] = coord(localLoop[i], localLoopWeights[i]);
377 ctr += loopPoints[i];
379 ctr /= localLoop.size();
384 getBase(cutPlane.
normal(), e0, e1);
388 SortableList<scalar> sortedAngles(localLoop.size());
398 loop.
setSize(localLoop.size());
404 const labelList& indices = sortedAngles.indices();
408 loop[i] = localLoop[indices[i]];
409 loopWeights[i] = localLoopWeights[indices[i]];
414 bool filterLoop =
false;
420 if (isEdge(cut) && edgeEndsCut(loop, i))
439 if (isEdge(cut) && edgeEndsCut(loop, i))
445 filteredLoop[filterI] = loop[i];
446 filteredLoopWeights[filterI] = loopWeights[i];
451 filteredLoopWeights.setSize(filterI);
454 loopWeights.
transfer(filteredLoopWeights);
462 Pout<<
"At angle:" << sortedAngles[i] <<
endl 465 writeCut(
Pout, loop[i], loopWeights[i]);
467 Pout<<
" coord:" << coord(loop[i], loopWeights[i]);
const labelListList & cellEdges() const
void size(const label n)
Older name for setAddressableSize.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
void append(const T &val)
Append an element at the end of the list.
const labelListList & pointEdges() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
const cellList & cells() const
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Macros for easy insertion into run-time selection tables.
Various functions to operate on Lists.
virtual const pointField & points() const
Return raw points.
#define forAll(list, i)
Loop across all elements in list.
static label getEdge(List< DynamicList< label >> &pe, DynamicList< edge > &es, const label pointi, const label nextPointi)
vectorField pointField
pointField is a vectorField.
const dimensionedScalar e
Elementary charge.
void setSize(const label n)
Alias for resize()
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
scalar distance(const point &p) const
Return distance (magnitude) from the given point to the plane.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
const polyMesh & mesh() const
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
int debug
Static debugging option.
defineTypeNameAndDebug(combustionModel, 0)
virtual bool cut(const vector &refDir, const label celli, const boolList &vertIsCut, const boolList &edgeIsCut, const scalarField &edgeWeight, labelList &loop, scalarField &loopWeights) const
Create cut along circumference of celli. Gets current mesh cuts.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
vector point
Point is a vector.
List< label > labelList
A List of labels.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const vector & normal() const noexcept
The plane unit normal.
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
static constexpr const zero Zero
Global zero (0)