surfaceInflate.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2015 OpenFOAM Foundation
9 Copyright (C) 2020-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 surfaceInflate
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Inflates surface. WIP. Checks for overlaps and locally lowers
35 inflation distance
36
37Usage
38 - surfaceInflate [OPTION]
39
40 \param -checkSelfIntersection \n
41 Includes checks for self-intersection.
42
43 \param -nSmooth
44 Specify number of smoothing iterations
45
46 \param -featureAngle
47 Specify a feature angle
48
49
50 E.g. inflate surface by 20mm with 1.5 safety factor:
51 surfaceInflate DTC-scaled.obj 0.02 1.5 -featureAngle 45 -nSmooth 2
52
53\*---------------------------------------------------------------------------*/
54
55#include "argList.H"
56#include "Time.H"
57#include "triSurfaceFields.H"
58#include "triSurfaceMesh.H"
59#include "triSurfaceGeoMesh.H"
60#include "bitSet.H"
61#include "OBJstream.H"
62#include "surfaceFeatures.H"
63
64using namespace Foam;
65
66// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
67
68scalar calcVertexNormalWeight
69(
70 const triFace& f,
71 const label pI,
72 const vector& fN,
73 const pointField& points
74)
75{
76 label index = f.find(pI);
77
78 if (index == -1)
79 {
81 << "Point not in face" << abort(FatalError);
82 }
83
84 const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
85 const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
86
87 return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
88}
89
90
91tmp<vectorField> calcVertexNormals(const triSurface& surf)
92{
93 // Weighted average of normals of faces attached to the vertex
94 // Weight = fA / (mag(e0)^2 * mag(e1)^2);
95 auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero);
96 auto& pointNormals = tpointNormals.ref();
97
98 const pointField& points = surf.points();
99 const labelListList& pointFaces = surf.pointFaces();
100 const labelList& meshPoints = surf.meshPoints();
101
102 forAll(pointFaces, pI)
103 {
104 const labelList& pFaces = pointFaces[pI];
105
106 forAll(pFaces, fI)
107 {
108 const label faceI = pFaces[fI];
109 const triFace& f = surf[faceI];
110
111 vector areaNorm = f.areaNormal(points);
112
113 scalar weight = calcVertexNormalWeight
114 (
115 f,
116 meshPoints[pI],
117 areaNorm,
118 points
119 );
120
121 pointNormals[pI] += weight * areaNorm;
122 }
123
124 pointNormals[pI].normalise();
125 }
126
127 return tpointNormals;
128}
129
130
131tmp<vectorField> calcPointNormals
132(
133 const triSurface& s,
134 const bitSet& isFeaturePoint,
136)
137{
138 //const pointField pointNormals(s.pointNormals());
139 tmp<vectorField> tpointNormals(calcVertexNormals(s));
140 vectorField& pointNormals = tpointNormals.ref();
141
142
143 // feature edges: create edge normals from edgeFaces only.
144 {
145 const labelListList& edgeFaces = s.edgeFaces();
146
147 labelList nNormals(s.nPoints(), Zero);
148 forAll(edgeStat, edgeI)
149 {
150 if (edgeStat[edgeI] != surfaceFeatures::NONE)
151 {
152 const edge& e = s.edges()[edgeI];
153 forAll(e, i)
154 {
155 if (!isFeaturePoint.test(e[i]))
156 {
157 pointNormals[e[i]] = Zero;
158 }
159 }
160 }
161 }
162
163 forAll(edgeStat, edgeI)
164 {
165 if (edgeStat[edgeI] != surfaceFeatures::NONE)
166 {
167 const labelList& eFaces = edgeFaces[edgeI];
168
169 // Get average edge normal
170 vector n = Zero;
171 for (const label facei : eFaces)
172 {
173 n += s.faceNormals()[facei];
174 }
175 n /= eFaces.size();
176
177
178 // Sum to points
179 const edge& e = s.edges()[edgeI];
180 forAll(e, i)
181 {
182 if (!isFeaturePoint.test(e[i]))
183 {
184 pointNormals[e[i]] += n;
185 nNormals[e[i]]++;
186 }
187 }
188 }
189 }
190
191 forAll(nNormals, pointI)
192 {
193 if (nNormals[pointI] > 0)
194 {
195 pointNormals[pointI].normalise();
196 }
197 }
198 }
199
200
201 forAll(pointNormals, pointI)
202 {
203 if (mag(mag(pointNormals[pointI])-1) > SMALL)
204 {
206 << "unitialised"
207 << exit(FatalError);
208 }
209 }
210
211 return tpointNormals;
212}
213
214
215void detectSelfIntersections
216(
217 const triSurfaceMesh& s,
218 bitSet& isEdgeIntersecting
219)
220{
221 const edgeList& edges = s.edges();
222 const indexedOctree<treeDataTriSurface>& tree = s.tree();
223 const labelList& meshPoints = s.meshPoints();
224 const pointField& points = s.points();
225
226 isEdgeIntersecting.setSize(edges.size());
227 isEdgeIntersecting = false;
228
229 forAll(edges, edgeI)
230 {
231 const edge& e = edges[edgeI];
232
233 pointIndexHit hitInfo
234 (
235 tree.findLine
236 (
237 points[meshPoints[e[0]]],
238 points[meshPoints[e[1]]],
240 )
241 );
242
243 if (hitInfo.hit())
244 {
245 isEdgeIntersecting.set(edgeI);
246 }
247 }
248}
249
250
251label detectIntersectionPoints
252(
253 const scalar tolerance,
254 const triSurfaceMesh& s,
255
256 const scalar extendFactor,
257 const pointField& initialPoints,
258 const vectorField& displacement,
259
260 const bool checkSelfIntersect,
261 const bitSet& initialIsEdgeIntersecting,
262
263 bitSet& isPointOnHitEdge,
264 scalarField& scale
265)
266{
267 const pointField initialLocalPoints(initialPoints, s.meshPoints());
268 const List<labelledTri>& localFaces = s.localFaces();
269
270
271 label nHits = 0;
272 isPointOnHitEdge.setSize(s.nPoints());
273 isPointOnHitEdge = false;
274
275
276 // 1. Extrusion offset vectors intersecting new surface location
277 {
278 scalar tol = max(tolerance, 10*s.tolerance());
279
280 // Collect all the edge vectors. Slightly shorten the edges to prevent
281 // finding lots of intersections. The fast triangle intersection routine
282 // has problems with rays passing through a point of the
283 // triangle.
284 pointField start(initialLocalPoints+tol*displacement);
285 pointField end(initialLocalPoints+extendFactor*displacement);
286
288 s.findLineAny(start, end, hits);
289
290 forAll(hits, pointI)
291 {
292 if
293 (
294 hits[pointI].hit()
295 && !localFaces[hits[pointI].index()].found(pointI)
296 )
297 {
298 scale[pointI] = max(0.0, scale[pointI]-0.2);
299
300 isPointOnHitEdge.set(pointI);
301 nHits++;
302 }
303 }
304 }
305
306
307 // 2. (new) surface self intersections
308 if (checkSelfIntersect)
309 {
310 bitSet isEdgeIntersecting;
311 detectSelfIntersections(s, isEdgeIntersecting);
312
313 const edgeList& edges = s.edges();
314 const pointField& points = s.points();
315
316 forAll(edges, edgeI)
317 {
318 const edge& e = edges[edgeI];
319
320 if (isEdgeIntersecting[edgeI] && !initialIsEdgeIntersecting[edgeI])
321 {
322 if (isPointOnHitEdge.set(e[0]))
323 {
324 label start = s.meshPoints()[e[0]];
325 const point& pt = points[start];
326
327 Pout<< "Additional self intersection at "
328 << pt
329 << endl;
330
331 scale[e[0]] = max(0.0, scale[e[0]]-0.2);
332 nHits++;
333 }
334 if (isPointOnHitEdge.set(e[1]))
335 {
336 label end = s.meshPoints()[e[1]];
337 const point& pt = points[end];
338
339 Pout<< "Additional self intersection at "
340 << pt
341 << endl;
342
343 scale[e[1]] = max(0.0, scale[e[1]]-0.2);
344 nHits++;
345 }
346 }
347 }
348 }
349
350 return nHits;
351}
352
353
355(
356 const triSurface& s,
357 const scalarField& fld,
358 const scalarField& edgeWeights
359)
360{
361 tmp<scalarField> tres(new scalarField(s.nPoints(), Zero));
362 scalarField& res = tres.ref();
363
364 scalarField sumWeight(s.nPoints(), Zero);
365
366 const edgeList& edges = s.edges();
367
368 forAll(edges, edgeI)
369 {
370 const edge& e = edges[edgeI];
371 const scalar w = edgeWeights[edgeI];
372
373 res[e[0]] += w*fld[e[1]];
374 sumWeight[e[0]] += w;
375
376 res[e[1]] += w*fld[e[0]];
377 sumWeight[e[1]] += w;
378 }
379
380 // Average
381 // ~~~~~~~
382
383 forAll(res, pointI)
384 {
385 if (mag(sumWeight[pointI]) < VSMALL)
386 {
387 // Unconnected point. Take over original value
388 res[pointI] = fld[pointI];
389 }
390 else
391 {
392 res[pointI] /= sumWeight[pointI];
393 }
394 }
395
396 return tres;
397}
398
399
400void minSmooth
401(
402 const triSurface& s,
403 const bitSet& isAffectedPoint,
404 const scalarField& fld,
405 scalarField& newFld
406)
407{
408
409 const edgeList& edges = s.edges();
410 const pointField& points = s.points();
411 const labelList& mp = s.meshPoints();
412
413 scalarField edgeWeights(edges.size());
414 forAll(edges, edgeI)
415 {
416 const edge& e = edges[edgeI];
417 scalar w = mag(points[mp[e[0]]]-points[mp[e[1]]]);
418
419 edgeWeights[edgeI] = 1.0/(max(w, SMALL));
420 }
421
422 tmp<scalarField> tavgFld = avg(s, fld, edgeWeights);
423
424 const scalarField& avgFld = tavgFld();
425
426 forAll(fld, pointI)
427 {
428 if (isAffectedPoint.test(pointI))
429 {
430 newFld[pointI] = min
431 (
432 fld[pointI],
433 0.5*fld[pointI] + 0.5*avgFld[pointI]
434 //avgFld[pointI]
435 );
436 }
437 }
438}
439
440
441void lloydsSmoothing
442(
443 const label nSmooth,
444 triSurface& s,
445 const bitSet& isFeaturePoint,
446 const List<surfaceFeatures::edgeStatus>& edgeStat,
447 const bitSet& isAffectedPoint
448)
449{
450 const labelList& meshPoints = s.meshPoints();
451 const edgeList& edges = s.edges();
452
453
454 bitSet isSmoothPoint(isAffectedPoint);
455 // Extend isSmoothPoint
456 {
457 bitSet newIsSmoothPoint(isSmoothPoint);
458 forAll(edges, edgeI)
459 {
460 const edge& e = edges[edgeI];
461 if (isSmoothPoint.test(e[0]))
462 {
463 newIsSmoothPoint.set(e[1]);
464 }
465 if (isSmoothPoint.test(e[1]))
466 {
467 newIsSmoothPoint.set(e[0]);
468 }
469 }
470 Info<< "From points-to-smooth " << isSmoothPoint.count()
471 << " to " << newIsSmoothPoint.count()
472 << endl;
473 isSmoothPoint.transfer(newIsSmoothPoint);
474 }
475
476 // Do some smoothing (Lloyds algorithm) around problematic points
477 for (label i = 0; i < nSmooth; i++)
478 {
479 const labelListList& pointFaces = s.pointFaces();
480 const pointField& faceCentres = s.faceCentres();
481
482 pointField newPoints(s.points());
483 forAll(isSmoothPoint, pointI)
484 {
485 if (isSmoothPoint[pointI] && !isFeaturePoint[pointI])
486 {
487 const labelList& pFaces = pointFaces[pointI];
488
489 point avg(Zero);
490 forAll(pFaces, pFaceI)
491 {
492 avg += faceCentres[pFaces[pFaceI]];
493 }
494 avg /= pFaces.size();
495 newPoints[meshPoints[pointI]] = avg;
496 }
497 }
498
499
500 // Move points on feature edges only according to feature edges.
501
502 const pointField& points = s.points();
503
504 vectorField pointSum(s.nPoints(), Zero);
505 labelList nPointSum(s.nPoints(), Zero);
506
507 forAll(edges, edgeI)
508 {
509 if (edgeStat[edgeI] != surfaceFeatures::NONE)
510 {
511 const edge& e = edges[edgeI];
512 const point& start = points[meshPoints[e[0]]];
513 const point& end = points[meshPoints[e[1]]];
514 point eMid = 0.5*(start+end);
515 pointSum[e[0]] += eMid;
516 nPointSum[e[0]]++;
517 pointSum[e[1]] += eMid;
518 nPointSum[e[1]]++;
519 }
520 }
521
522 forAll(pointSum, pointI)
523 {
524 if
525 (
526 isSmoothPoint[pointI]
527 && isFeaturePoint[pointI]
528 && nPointSum[pointI] >= 2
529 )
530 {
531 newPoints[meshPoints[pointI]] =
532 pointSum[pointI]/nPointSum[pointI];
533 }
534 }
535
536
537 s.movePoints(newPoints);
538
539
540 // Extend isSmoothPoint
541 {
542 bitSet newIsSmoothPoint(isSmoothPoint);
543 forAll(edges, edgeI)
544 {
545 const edge& e = edges[edgeI];
546 if (isSmoothPoint[e[0]])
547 {
548 newIsSmoothPoint.set(e[1]);
549 }
550 if (isSmoothPoint[e[1]])
551 {
552 newIsSmoothPoint.set(e[0]);
553 }
554 }
555 Info<< "From points-to-smooth " << isSmoothPoint.count()
556 << " to " << newIsSmoothPoint.count()
557 << endl;
558 isSmoothPoint.transfer(newIsSmoothPoint);
559 }
560 }
561}
562
563
564
565// Main program:
566
567int main(int argc, char *argv[])
568{
569 argList::addNote
570 (
571 "Inflates surface according to point normals."
572 );
573
574 argList::noParallel();
575 argList::addNote
576 (
577 "Creates inflated version of surface using point normals. "
578 "Takes surface, distance to inflate and additional safety factor"
579 );
580 argList::addBoolOption
581 (
582 "checkSelfIntersection",
583 "Also check for self-intersection"
584 );
585 argList::addOption
586 (
587 "nSmooth",
588 "integer",
589 "Number of smoothing iterations (default 20)"
590 );
591 argList::addOption
592 (
593 "featureAngle",
594 "scalar",
595 "Feature angle"
596 );
597 argList::addBoolOption
598 (
599 "debug",
600 "Switch on additional debug information"
601 );
602
603 argList::addArgument("input", "The input surface file");
604 argList::addArgument("distance", "The inflate distance");
605 argList::addArgument("factor", "The extend safety factor [1,10]");
606
607 argList::noFunctionObjects(); // Never use function objects
608
609 #include "setRootCase.H"
610 #include "createTime.H"
611
612 const auto inputName = args.get<word>(1);
613 const auto distance = args.get<scalar>(2);
614 const auto extendFactor = args.get<scalar>(3);
615 const bool checkSelfIntersect = args.found("checkSelfIntersection");
616 const auto nSmooth = args.getOrDefault<label>("nSmooth", 10);
617 const auto featureAngle = args.getOrDefault<scalar>("featureAngle", 180);
618 const bool debug = args.found("debug");
619
620
621 Info<< "Inflating surface " << inputName
622 << " according to point normals" << nl
623 << " distance : " << distance << nl
624 << " safety factor : " << extendFactor << nl
625 << " self intersection test : " << checkSelfIntersect << nl
626 << " smoothing iterations : " << nSmooth << nl
627 << " feature angle : " << featureAngle << nl
628 << endl;
629
630
631 if (extendFactor < 1 || extendFactor > 10)
632 {
634 << "Illegal safety factor " << extendFactor
635 << ". It is usually 1..2"
636 << exit(FatalError);
637 }
638
639
640
641 // Load triSurfaceMesh
643 (
645 (
646 inputName, // name
647 runTime.constant(), // instance
648 "triSurface", // local
649 runTime, // registry
650 IOobject::MUST_READ,
651 IOobject::AUTO_WRITE
652 )
653 );
654
655 // Mark features
656 const surfaceFeatures features(s, 180.0-featureAngle);
657
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
663 << endl;
664
665 bitSet isFeaturePoint(s.nPoints(), features.featurePoints());
666
667 const List<surfaceFeatures::edgeStatus> edgeStat(features.toStatus());
668
669
670 const pointField initialPoints(s.points());
671
672
673 // Construct scale
674 Info<< "Constructing field scale\n" << endl;
676 (
678 (
679 "scale", // name
680 runTime.timeName(), // instance
681 s, // registry
682 IOobject::READ_IF_PRESENT,
683 IOobject::AUTO_WRITE
684 ),
685 s,
686 dimensionedScalar("scale", dimLength, 1.0)
687 );
688
689
690 // Construct unit normals
691
692 Info<< "Calculating vertex normals\n" << endl;
693 const pointField pointNormals
694 (
695 calcPointNormals
696 (
697 s,
698 isFeaturePoint,
699 edgeStat
700 )
701 );
702
703
704 // Construct pointDisplacement
705 Info<< "Constructing field pointDisplacement\n" << endl;
706 triSurfacePointVectorField pointDisplacement
707 (
709 (
710 "pointDisplacement", // name
711 runTime.timeName(), // instance
712 s, // registry
713 IOobject::NO_READ,
714 IOobject::AUTO_WRITE
715 ),
716 s,
717 dimLength,
718 distance*scale*pointNormals
719 );
720
721
722 const labelList& meshPoints = s.meshPoints();
723
724
725 // Any point on any intersected edge in any of the iterations
726 bitSet isScaledPoint(s.nPoints());
727
728
729 // Detect any self intersections on initial mesh
730 bitSet initialIsEdgeIntersecting;
731 if (checkSelfIntersect)
732 {
733 detectSelfIntersections(s, initialIsEdgeIntersecting);
734 }
735 else
736 {
737 // Mark all edges as already self intersecting so avoid detecting any
738 // new ones
739 initialIsEdgeIntersecting.setSize(s.nEdges(), true);
740 }
741
742
743 // Inflate
744 while (runTime.loop())
745 {
746 Info<< "Time = " << runTime.timeName() << nl << endl;
747
748 // Move to new location
749 pointField newPoints(initialPoints);
750 forAll(meshPoints, i)
751 {
752 newPoints[meshPoints[i]] += pointDisplacement[i];
753 }
754 s.movePoints(newPoints);
755 Info<< "Bounding box : " << s.bounds() << endl;
756
757
758 // Work out scale from pointDisplacement
759 forAll(scale, pointI)
760 {
761 if (s.pointFaces()[pointI].size())
762 {
763 scale[pointI] = mag(pointDisplacement[pointI])/distance;
764 }
765 else
766 {
767 scale[pointI] = 1.0;
768 }
769 }
770
771
772 Info<< "Scale : " << gAverage(scale) << endl;
773 Info<< "Displacement : " << gAverage(pointDisplacement) << endl;
774
775
776
777 // Detect any intersections and scale back
778 bitSet isAffectedPoint;
779 label nIntersections = detectIntersectionPoints
780 (
781 1e-9, // intersection tolerance
782 s,
783 extendFactor,
784 initialPoints,
785 pointDisplacement,
786
787 checkSelfIntersect,
788 initialIsEdgeIntersecting,
789
790 isAffectedPoint,
791 scale
792 );
793 Info<< "Detected " << nIntersections << " intersections" << nl
794 << endl;
795
796 if (nIntersections == 0)
797 {
798 runTime.write();
799 break;
800 }
801
802
803 // Accumulate all affected points
804 forAll(isAffectedPoint, pointI)
805 {
806 if (isAffectedPoint.test(pointI))
807 {
808 isScaledPoint.set(pointI);
809 }
810 }
811
812 // Smear out lowering of scale so any edges not found are
813 // still included
814 for (label i = 0; i < nSmooth; i++)
815 {
816 triSurfacePointScalarField oldScale(scale);
817 oldScale.rename("oldScale");
818 minSmooth
819 (
820 s,
821 bitSet(s.nPoints(), true),
822 oldScale,
823 scale
824 );
825 }
826
827
828 // From scale update the pointDisplacement
829 pointDisplacement *= distance*scale/mag(pointDisplacement);
830
831
832 // Do some smoothing (Lloyds algorithm)
833 lloydsSmoothing(nSmooth, s, isFeaturePoint, edgeStat, isAffectedPoint);
834
835
836 // Update pointDisplacement
837 const pointField& pts = s.points();
838 forAll(meshPoints, i)
839 {
840 label meshPointI = meshPoints[i];
841 pointDisplacement[i] = pts[meshPointI]-initialPoints[meshPointI];
842 }
843
844
845 // Write
846 runTime.write();
847
848 runTime.printExecutionTime(Info);
849 }
850
851
852 Info<< "Detected overall intersecting points " << isScaledPoint.count()
853 << " out of " << s.nPoints() << nl << endl;
854
855
856 if (debug)
857 {
858 OBJstream str(runTime.path()/"isScaledPoint.obj");
859 forAll(isScaledPoint, pointI)
860 {
861 if (isScaledPoint[pointI])
862 {
863 str.write(initialPoints[meshPoints[pointI]]);
864 }
865 }
866 }
867
868
869 Info<< "End\n" << endl;
870
871 return 0;
872}
873
874
875// ************************************************************************* //
bool found
label n
Y[inertIndex] max(0.0)
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
void normalise()
Definition: Field.C:538
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
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.
Definition: List.H:330
OFstream that keeps track of vertices.
Definition: OBJstream.H:61
void setSize(const label n, unsigned int val=0u)
Alias for resize()
Definition: PackedList.H:557
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label nPoints() const
Number of points supporting patch faces.
const labelList & meshPoints() const
Return labelList of mesh points in patch.
const Field< point_type > & points() const noexcept
Return reference to global points.
const labelListList & pointFaces() const
Return point-face addressing.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:66
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:590
bool test(const label pos) const
Test value at specified position, never auto-vivify entries.
Definition: bitSetI.H:521
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
Holds feature edges/points of surface.
A class for managing temporary objects.
Definition: tmp.H:65
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:72
IOoject and searching on triSurface.
Triangulated surface description with patch information.
Definition: triSurface.H:79
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
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))
const dimensionedScalar mp
Proton mass.
int debug
Static debugging option.
Namespace for OpenFOAM.
scalar distance(const vector &p1, const vector &p2)
Definition: curveTools.C:12
messageStream Info
Information stream (stdout output on master, null elsewhere)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Type gAverage(const FieldField< Field, Type > &f)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
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]
Definition: readKivaGrid.H:235
labelList f(nPoints)
Foam::argList args(argc, argv)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Fields for triSurface.