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
227 isEdgeIntersecting =
false;
231 const edge&
e = edges[edgeI];
245 isEdgeIntersecting.
set(edgeI);
251 label detectIntersectionPoints
253 const scalar tolerance,
256 const scalar extendFactor,
260 const bool checkSelfIntersect,
261 const bitSet& initialIsEdgeIntersecting,
267 const pointField initialLocalPoints(initialPoints,
s.meshPoints());
272 isPointOnHitEdge.
setSize(
s.nPoints());
273 isPointOnHitEdge =
false;
278 scalar tol =
max(tolerance, 10*
s.tolerance());
284 pointField start(initialLocalPoints+tol*displacement);
285 pointField end(initialLocalPoints+extendFactor*displacement);
288 s.findLineAny(start,
end, hits);
295 && !localFaces[hits[pointI].index()].
found(pointI)
298 scale[pointI] =
max(0.0, scale[pointI]-0.2);
300 isPointOnHitEdge.
set(pointI);
308 if (checkSelfIntersect)
310 bitSet isEdgeIntersecting;
311 detectSelfIntersections(
s, isEdgeIntersecting);
318 const edge&
e = edges[edgeI];
320 if (isEdgeIntersecting[edgeI] && !initialIsEdgeIntersecting[edgeI])
322 if (isPointOnHitEdge.
set(
e[0]))
324 label start =
s.meshPoints()[
e[0]];
327 Pout<<
"Additional self intersection at " 331 scale[
e[0]] =
max(0.0, scale[
e[0]]-0.2);
334 if (isPointOnHitEdge.
set(
e[1]))
336 label
end =
s.meshPoints()[
e[1]];
339 Pout<<
"Additional self intersection at " 343 scale[
e[1]] =
max(0.0, scale[
e[1]]-0.2);
370 const edge&
e = edges[edgeI];
371 const scalar w = edgeWeights[edgeI];
373 res[
e[0]] += w*
fld[
e[1]];
374 sumWeight[
e[0]] += w;
376 res[
e[1]] += w*
fld[
e[0]];
377 sumWeight[
e[1]] += w;
385 if (
mag(sumWeight[pointI]) < VSMALL)
388 res[pointI] =
fld[pointI];
392 res[pointI] /= sumWeight[pointI];
403 const bitSet& isAffectedPoint,
416 const edge&
e = edges[edgeI];
419 edgeWeights[edgeI] = 1.0/(
max(w, SMALL));
428 if (isAffectedPoint.
test(pointI))
433 0.5*
fld[pointI] + 0.5*avgFld[pointI]
445 const bitSet& isFeaturePoint,
447 const bitSet& isAffectedPoint
454 bitSet isSmoothPoint(isAffectedPoint);
457 bitSet newIsSmoothPoint(isSmoothPoint);
460 const edge&
e = edges[edgeI];
461 if (isSmoothPoint.test(
e[0]))
463 newIsSmoothPoint.
set(
e[1]);
465 if (isSmoothPoint.test(
e[1]))
467 newIsSmoothPoint.set(
e[0]);
470 Info<<
"From points-to-smooth " << isSmoothPoint.count()
471 <<
" to " << newIsSmoothPoint.count()
473 isSmoothPoint.transfer(newIsSmoothPoint);
477 for (label i = 0; i < nSmooth; i++)
483 forAll(isSmoothPoint, pointI)
485 if (isSmoothPoint[pointI] && !isFeaturePoint[pointI])
492 avg += faceCentres[
pFaces[pFaceI]];
495 newPoints[meshPoints[pointI]] = avg;
511 const edge&
e = edges[edgeI];
515 pointSum[
e[0]] += eMid;
517 pointSum[
e[1]] += eMid;
526 isSmoothPoint[pointI]
527 && isFeaturePoint[pointI]
528 && nPointSum[pointI] >= 2
531 newPoints[meshPoints[pointI]] =
532 pointSum[pointI]/nPointSum[pointI];
537 s.movePoints(newPoints);
542 bitSet newIsSmoothPoint(isSmoothPoint);
545 const edge&
e = edges[edgeI];
546 if (isSmoothPoint[
e[0]])
548 newIsSmoothPoint.
set(
e[1]);
550 if (isSmoothPoint[
e[1]])
552 newIsSmoothPoint.set(
e[0]);
555 Info<<
"From points-to-smooth " << isSmoothPoint.count()
556 <<
" to " << newIsSmoothPoint.count()
558 isSmoothPoint.transfer(newIsSmoothPoint);
567 int main(
int argc,
char *argv[])
571 "Inflates surface according to point normals." 577 "Creates inflated version of surface using point normals. " 578 "Takes surface, distance to inflate and additional safety factor" 582 "checkSelfIntersection",
583 "Also check for self-intersection" 589 "Number of smoothing iterations (default 20)" 600 "Switch on additional debug information" 614 const auto extendFactor =
args.
get<scalar>(3);
615 const bool checkSelfIntersect =
args.
found(
"checkSelfIntersection");
621 Info<<
"Inflating surface " << inputName
622 <<
" according to point normals" <<
nl 624 <<
" safety factor : " << extendFactor <<
nl 625 <<
" self intersection test : " << checkSelfIntersect <<
nl 626 <<
" smoothing iterations : " << nSmooth <<
nl 627 <<
" feature angle : " << featureAngle <<
nl 631 if (extendFactor < 1 || extendFactor > 10)
634 <<
"Illegal safety factor " << extendFactor
635 <<
". It is usually 1..2" 658 Info<<
"Detected features:" <<
nl 659 <<
" feature points : " << features.featurePoints().size()
660 <<
" out of " <<
s.nPoints() <<
nl 661 <<
" feature edges : " << features.featureEdges().size()
662 <<
" out of " <<
s.nEdges() <<
nl 665 bitSet isFeaturePoint(
s.nPoints(), features.featurePoints());
674 Info<<
"Constructing field scale\n" <<
endl;
693 Info<<
"Calculating vertex normals\n" <<
endl;
706 Info<<
"Constructing field pointDisplacement\n" <<
endl;
720 (
distance*scale.field()*pointNormals)
728 bitSet isScaledPoint(
s.nPoints());
732 bitSet initialIsEdgeIntersecting;
733 if (checkSelfIntersect)
735 detectSelfIntersections(
s, initialIsEdgeIntersecting);
741 initialIsEdgeIntersecting.setSize(
s.nEdges(),
true);
754 newPoints[meshPoints[i]] += pointDisplacement[i];
756 s.movePoints(newPoints);
757 Info<<
"Bounding box : " <<
s.bounds() <<
endl;
763 if (
s.pointFaces()[pointI].size())
765 scale[pointI] =
mag(pointDisplacement[pointI])/
distance;
781 label nIntersections = detectIntersectionPoints
790 initialIsEdgeIntersecting,
795 Info<<
"Detected " << nIntersections <<
" intersections" <<
nl 798 if (nIntersections == 0)
806 forAll(isAffectedPoint, pointI)
808 if (isAffectedPoint.
test(pointI))
810 isScaledPoint.set(pointI);
816 for (label i = 0; i < nSmooth; i++)
819 oldScale.rename(
"oldScale");
832 pointDisplacement.normalise();
833 pointDisplacement.field() *=
distance*scale.field();
837 lloydsSmoothing(nSmooth,
s, isFeaturePoint, edgeStat, isAffectedPoint);
844 label meshPointI = meshPoints[i];
845 pointDisplacement[i] =
pts[meshPointI]-initialPoints[meshPointI];
856 Info<<
"Detected overall intersecting points " << isScaledPoint.count()
857 <<
" out of " <<
s.nPoints() <<
nl <<
endl;
864 for (
const label pointi : isScaledPoint)
866 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.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
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)
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)