68 scalar calcVertexNormalWeight
76 label index =
f.
find(pI);
96 auto& pointNormals = tpointNormals.ref();
108 const label faceI =
pFaces[fI];
113 scalar weight = calcVertexNormalWeight
121 pointNormals[pI] += weight * areaNorm;
124 pointNormals[pI].normalise();
127 return tpointNormals;
134 const bitSet& isFeaturePoint,
152 const edge&
e =
s.edges()[edgeI];
155 if (!isFeaturePoint.
test(
e[i]))
157 pointNormals[
e[i]] =
Zero;
167 const labelList& eFaces = edgeFaces[edgeI];
171 for (
const label facei : eFaces)
173 n +=
s.faceNormals()[facei];
179 const edge&
e =
s.edges()[edgeI];
182 if (!isFeaturePoint.
test(
e[i]))
184 pointNormals[
e[i]] +=
n;
193 if (nNormals[pointI] > 0)
201 forAll(pointNormals, pointI)
203 if (
mag(
mag(pointNormals[pointI])-1) > SMALL)
211 return tpointNormals;
215 void detectSelfIntersections
218 bitSet& isEdgeIntersecting
228 isEdgeIntersecting =
false;
232 const edge&
e = edges[edgeI];
246 isEdgeIntersecting.
set(edgeI);
252 label detectIntersectionPoints
254 const scalar tolerance,
257 const scalar extendFactor,
261 const bool checkSelfIntersect,
262 const bitSet& initialIsEdgeIntersecting,
268 const pointField initialLocalPoints(initialPoints,
s.meshPoints());
273 isPointOnHitEdge.
setSize(
s.nPoints());
274 isPointOnHitEdge =
false;
279 scalar tol =
max(tolerance, 10*
s.tolerance());
285 pointField start(initialLocalPoints+tol*displacement);
286 pointField end(initialLocalPoints+extendFactor*displacement);
289 s.findLineAny(start,
end, hits);
296 && !localFaces[hits[pointI].index()].
found(pointI)
299 scale[pointI] =
max(0.0, scale[pointI]-0.2);
301 isPointOnHitEdge.
set(pointI);
309 if (checkSelfIntersect)
311 bitSet isEdgeIntersecting;
312 detectSelfIntersections(
s, isEdgeIntersecting);
320 const edge&
e = edges[edgeI];
322 if (isEdgeIntersecting[edgeI] && !initialIsEdgeIntersecting[edgeI])
324 if (isPointOnHitEdge.
set(
e[0]))
326 label start =
s.meshPoints()[
e[0]];
329 Pout<<
"Additional self intersection at " 333 scale[
e[0]] =
max(0.0, scale[
e[0]]-0.2);
336 if (isPointOnHitEdge.
set(
e[1]))
338 label
end =
s.meshPoints()[
e[1]];
341 Pout<<
"Additional self intersection at " 345 scale[
e[1]] =
max(0.0, scale[
e[1]]-0.2);
364 auto& res = tres.
ref();
372 const edge&
e = edges[edgeI];
373 const scalar w = edgeWeights[edgeI];
375 res[
e[0]] += w*
fld[
e[1]];
376 sumWeight[
e[0]] += w;
378 res[
e[1]] += w*
fld[
e[0]];
379 sumWeight[
e[1]] += w;
387 if (
mag(sumWeight[pointI]) < VSMALL)
390 res[pointI] =
fld[pointI];
394 res[pointI] /= sumWeight[pointI];
405 const bitSet& isAffectedPoint,
418 const edge&
e = edges[edgeI];
421 edgeWeights[edgeI] = 1.0/(
max(w, SMALL));
430 if (isAffectedPoint.
test(pointI))
435 0.5*
fld[pointI] + 0.5*avgFld[pointI]
447 const bitSet& isFeaturePoint,
449 const bitSet& isAffectedPoint
456 bitSet isSmoothPoint(isAffectedPoint);
459 bitSet newIsSmoothPoint(isSmoothPoint);
462 const edge&
e = edges[edgeI];
463 if (isSmoothPoint.test(
e[0]))
465 newIsSmoothPoint.
set(
e[1]);
467 if (isSmoothPoint.test(
e[1]))
469 newIsSmoothPoint.set(
e[0]);
472 Info<<
"From points-to-smooth " << isSmoothPoint.count()
473 <<
" to " << newIsSmoothPoint.count()
475 isSmoothPoint.transfer(newIsSmoothPoint);
479 for (label i = 0; i < nSmooth; i++)
485 forAll(isSmoothPoint, pointI)
487 if (isSmoothPoint[pointI] && !isFeaturePoint[pointI])
494 avg += faceCentres[
pFaces[pFaceI]];
497 newPoints[meshPoints[pointI]] = avg;
513 const edge&
e = edges[edgeI];
517 pointSum[
e[0]] += eMid;
519 pointSum[
e[1]] += eMid;
528 isSmoothPoint[pointI]
529 && isFeaturePoint[pointI]
530 && nPointSum[pointI] >= 2
533 newPoints[meshPoints[pointI]] =
534 pointSum[pointI]/nPointSum[pointI];
539 s.movePoints(newPoints);
544 bitSet newIsSmoothPoint(isSmoothPoint);
547 const edge&
e = edges[edgeI];
548 if (isSmoothPoint[
e[0]])
550 newIsSmoothPoint.
set(
e[1]);
552 if (isSmoothPoint[
e[1]])
554 newIsSmoothPoint.set(
e[0]);
557 Info<<
"From points-to-smooth " << isSmoothPoint.count()
558 <<
" to " << newIsSmoothPoint.count()
560 isSmoothPoint.transfer(newIsSmoothPoint);
569 int main(
int argc,
char *argv[])
573 "Inflates surface according to point normals." 579 "Creates inflated version of surface using point normals. " 580 "Takes surface, distance to inflate and additional safety factor" 584 "checkSelfIntersection",
585 "Also check for self-intersection" 591 "Number of smoothing iterations (default 20)" 602 "Switch on additional debug information" 616 const auto extendFactor =
args.
get<scalar>(3);
617 const bool checkSelfIntersect =
args.
found(
"checkSelfIntersection");
623 Info<<
"Inflating surface " << inputName
624 <<
" according to point normals" <<
nl 626 <<
" safety factor : " << extendFactor <<
nl 627 <<
" self intersection test : " << checkSelfIntersect <<
nl 628 <<
" smoothing iterations : " << nSmooth <<
nl 629 <<
" feature angle : " << featureAngle <<
nl 633 if (extendFactor < 1 || extendFactor > 10)
636 <<
"Illegal safety factor " << extendFactor
637 <<
". It is usually 1..2" 660 Info<<
"Detected features:" <<
nl 661 <<
" feature points : " << features.featurePoints().size()
662 <<
" out of " <<
s.nPoints() <<
nl 663 <<
" feature edges : " << features.featureEdges().size()
664 <<
" out of " <<
s.nEdges() <<
nl 667 bitSet isFeaturePoint(
s.nPoints(), features.featurePoints());
676 Info<<
"Constructing field scale\n" <<
endl;
695 Info<<
"Calculating vertex normals\n" <<
endl;
708 Info<<
"Constructing field pointDisplacement\n" <<
endl;
722 (
distance*scale.field()*pointNormals)
730 bitSet isScaledPoint(
s.nPoints());
734 bitSet initialIsEdgeIntersecting;
735 if (checkSelfIntersect)
737 detectSelfIntersections(
s, initialIsEdgeIntersecting);
743 initialIsEdgeIntersecting.setSize(
s.nEdges(),
true);
756 newPoints[meshPoints[i]] += pointDisplacement[i];
758 s.movePoints(newPoints);
759 Info<<
"Bounding box : " <<
s.bounds() <<
endl;
765 if (
s.pointFaces()[pointI].size())
767 scale[pointI] =
mag(pointDisplacement[pointI])/
distance;
783 label nIntersections = detectIntersectionPoints
792 initialIsEdgeIntersecting,
797 Info<<
"Detected " << nIntersections <<
" intersections" <<
nl 800 if (nIntersections == 0)
808 forAll(isAffectedPoint, pointI)
810 if (isAffectedPoint.
test(pointI))
812 isScaledPoint.set(pointI);
818 for (label i = 0; i < nSmooth; i++)
821 oldScale.rename(
"oldScale");
834 pointDisplacement.normalise();
835 pointDisplacement.field() *=
distance*scale.field();
839 lloydsSmoothing(nSmooth,
s, isFeaturePoint, edgeStat, isAffectedPoint);
847 label meshPointI = meshPoints[i];
848 pointDisplacement[i] =
pts[meshPointI]-initialPoints[meshPointI];
859 Info<<
"Detected overall intersecting points " << isScaledPoint.count()
860 <<
" out of " <<
s.nPoints() <<
nl <<
endl;
867 for (
const label pointi : isScaledPoint)
869 str.
write(initialPoints[meshPoints[pointi]]);
label nPoints() const
Number of points supporting patch faces.
static void noFunctionObjects(bool addWithOption=false)
Remove '-noFunctionObjects' option and ignore any occurrences.
static void addNote(const string ¬e)
Add extra notes for the usage information.
void size(const label n)
Older name for setAddressableSize.
fileName path() const
Return path = rootPath/caseName. Same as TimePaths::path()
void set(const bitSet &bitset)
Set specified bits from another bitset.
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
virtual Ostream & write(const char c) override
Write character.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
T & ref() const
Return non-const reference to the contents of a non-null managed pointer.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
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)
void setSize(const label n, unsigned int val=0u)
Alias for resize()
virtual bool loop()
Return true if run should continue and if so increment time.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static void addBoolOption(const word &optName, const string &usage="", bool advanced=false)
Add a bool option to validOptions with usage information.
static void noParallel()
Remove the parallel options.
scalar distance(const vector &p1, const vector &p2)
This class describes the interaction of an object (often a face) and a point. It carries the info of ...
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
#define forAll(list, i)
Loop across all elements in list.
IOoject and searching on triSurface.
label fcIndex(const label i) const noexcept
The forward circular index. The next index in the list which returns to the first at the end of the l...
Not a classified feature edge.
const dimensionedScalar e
Elementary charge.
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
A class for handling words, derived from Foam::string.
A triangular face using a FixedList of labels corresponding to mesh vertices.
Tree tree(triangles.begin(), triangles.end())
Reading is optional [identical to LAZY_READ].
virtual bool write(const bool writeOnProc=true) const
Write using setting from DB.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
static void addOption(const word &optName, const string ¶m="", const string &usage="", bool advanced=false)
Add an option to validOptions with usage information.
const Field< point_type > & points() const noexcept
Return reference to global points.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
errorManip< error > abort(error &err)
bool test(const label pos) const
Test for True value at specified position, never auto-vivify entries.
label find(const T &val) const
Find index of the first occurrence of the value.
const labelListList & pointFaces() const
Return point-face addressing.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
An OFstream that keeps track of vertices and provides convenience output methods for OBJ files...
static word timeName(const scalar t, const int precision=precision_)
Return a time name for the given scalar time value formatted with the given precision.
const word & constant() const noexcept
Return constant name.
int debug
Static debugging option.
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Ostream & printExecutionTime(OSstream &os) const
Print the elapsed ExecutionTime (cpu-time), ClockTime.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=cellModel::ref(cellModel::HEX);labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells].reset(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< SMALL) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &oldCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
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;for(const word &name :lagrangianScalarNames){ IOField< scalar > fld(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
T get(const label index) const
Get a value from the argument at index.
Non-pointer based hierarchical recursive searching.
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Type gAverage(const FieldField< Field, Type > &f)
void resize_nocopy(const label numElem)
Currently identical to resize. Subject to future change (Oct-2021)
Automatically write from objectRegistry::writeObject()
static void addArgument(const string &argName, const string &usage="")
Append a (mandatory) argument to validArgs.
messageStream Info
Information stream (stdout output on master, null elsewhere)
A class for managing temporary objects.
void normalise()
Inplace normalise this field. Generally a no-op except for vector fields.
Triangulated surface description with patch information.
label rcIndex(const label i) const noexcept
The reverse circular index. The previous index in the list which returns to the last at the beginning...
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))
Foam::argList args(argc, argv)
Defines the attributes of an object for which implicit objectRegistry management is supported...
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
const dimensionedScalar mp
Proton mass.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
bool found(const word &optName) const
Return true if the named option is found.
Holds feature edges/points of surface.
static constexpr const zero Zero
Global zero (0)