conformalVoronoiMesh.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) 2012-2016 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
27\*---------------------------------------------------------------------------*/
28
30#include "initialPointsMethod.H"
31#include "relaxationModel.H"
32#include "faceAreaWeightModel.H"
33#include "meshSearch.H"
34#include "vectorTools.H"
35#include "IOmanip.H"
36#include "indexedCellChecks.H"
39#include "OBJstream.H"
40#include "indexedVertexOps.H"
41#include "DelaunayMeshTools.H"
42
43// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44
45namespace Foam
46{
47 defineTypeNameAndDebug(conformalVoronoiMesh, 0);
48}
49
50const Foam::Enum
51<
53>
55({
56 { dualMeshPointType::internal, "internal" },
57 { dualMeshPointType::surface, "surface" },
58 { dualMeshPointType::featureEdge, "featureEdge" },
59 { dualMeshPointType::featurePoint, "featurePoint" },
60 { dualMeshPointType::constrained, "constrained" },
61});
62
63
64// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
65
66void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground() const
67{
68 const cellShapeControlMesh& cellSizeMesh =
69 cellShapeControl_.shapeControlMesh();
70
71 DynamicList<Foam::point> pts(number_of_vertices());
72
73 for
74 (
75 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
76 vit != finite_vertices_end();
77 ++vit
78 )
79 {
80 if (vit->internalOrBoundaryPoint() && !vit->referred())
81 {
82 pts.append(topoint(vit->point()));
83 }
84 }
85
86 boundBox bb(pts);
87
88 boundBox cellSizeMeshBb = cellSizeMesh.bounds();
89
90 bool fullyContained = true;
91
92 if (!cellSizeMeshBb.contains(bb))
93 {
94 Pout<< "Triangulation not fully contained in cell size mesh."
95 << endl;
96
97 Pout<< "Cell Size Mesh Bounds = " << cellSizeMesh.bounds() << endl;
98 Pout<< "foamyHexMesh Bounds = " << bb << endl;
99
100 fullyContained = false;
101 }
102
103 reduce(fullyContained, andOp<unsigned int>());
104
105 Info<< "Triangulation is "
106 << (fullyContained ? "fully" : "not fully")
107 << " contained in the cell size mesh"
108 << endl;
109}
110
111
112void Foam::conformalVoronoiMesh::insertInternalPoints
113(
114 List<Point>& points,
115 bool distribute
116)
117{
118 label nPoints = points.size();
119
120 if (Pstream::parRun())
121 {
122 reduce(nPoints, sumOp<label>());
123 }
124
125 Info<< " " << nPoints << " points to insert..." << endl;
126
127 if (Pstream::parRun() && distribute)
128 {
129 List<Foam::point> transferPoints(points.size());
130
131 forAll(points, pI)
132 {
133 transferPoints[pI] = topoint(points[pI]);
134 }
135
136 // Send the points that are not on this processor to the appropriate
137 // place
139 (
140 decomposition_().distributePoints(transferPoints)
141 );
142
143 transferPoints.clear();
144
145 map().distribute(points);
146 }
147
148 label nVert = number_of_vertices();
149
151
152 label nInserted(number_of_vertices() - nVert);
153
154 if (Pstream::parRun())
155 {
156 reduce(nInserted, sumOp<label>());
157 }
158
159 Info<< " " << nInserted << " points inserted"
160 << ", failed to insert " << nPoints - nInserted
161 << " ("
162 << 100.0*(nPoints - nInserted)/(nInserted + SMALL)
163 << " %)"<< endl;
164
165 for
166 (
167 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
168 vit != finite_vertices_end();
169 ++vit
170 )
171 {
173 {
174 vit->index() = getNewVertexIndex();
175 vit->type() = Vb::vtInternal;
176 }
177 }
178}
179
180
181Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
182(
183 List<Vb>& vertices,
184 bool distribute,
185 bool reIndex
186)
187{
188 if (Pstream::parRun() && distribute)
189 {
190 autoPtr<mapDistribute> mapDist =
191 decomposition_().distributePoints(vertices);
192
193 // Re-index the point pairs if one or both have been distributed.
194 // If both, remove
195
196 // If added a point, then need to know its point pair
197 // If one moved, then update procIndex locally
198
199 forAll(vertices, vI)
200 {
201 vertices[vI].procIndex() = Pstream::myProcNo();
202 }
203 }
204
205 label preReinsertionSize(number_of_vertices());
206
207 Map<label> oldToNewIndices =
208 this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
209
210 const label nReinserted = returnReduce
211 (
212 label(number_of_vertices()) - preReinsertionSize,
213 sumOp<label>()
214 );
215
216 Info<< " Reinserted " << nReinserted << " vertices out of "
217 << returnReduce(vertices.size(), sumOp<label>())
218 << endl;
219
220 return oldToNewIndices;
221}
222
223
225(
226 const pointIndexHitAndFeatureList& surfaceHits,
227 const fileName fName,
228 DynamicList<Vb>& pts
229)
230{
231 forAll(surfaceHits, i)
232 {
233 vectorField norm(1);
234
235 const pointIndexHit surfaceHit = surfaceHits[i].first();
236 const label featureIndex = surfaceHits[i].second();
237
238 allGeometry_[featureIndex].getNormal
239 (
240 List<pointIndexHit>(1, surfaceHit),
241 norm
242 );
243
244 const vector& normal = norm[0];
245
246 const Foam::point& surfacePt(surfaceHit.hitPoint());
247
249 geometryToConformTo_.meshableSide(featureIndex, surfaceHit);
250
251 if (meshableSide == extendedFeatureEdgeMesh::BOTH)
252 {
253 createBafflePointPair
254 (
255 pointPairDistance(surfacePt),
256 surfacePt,
257 normal,
258 true,
259 pts
260 );
261 }
262 else if (meshableSide == extendedFeatureEdgeMesh::INSIDE)
263 {
264 createPointPair
265 (
266 pointPairDistance(surfacePt),
267 surfacePt,
268 normal,
269 true,
270 pts
271 );
272 }
273 else if (meshableSide == extendedFeatureEdgeMesh::OUTSIDE)
274 {
275 createPointPair
276 (
277 pointPairDistance(surfacePt),
278 surfacePt,
279 -normal,
280 true,
281 pts
282 );
283 }
284 else
285 {
287 << meshableSide << ", bad"
288 << endl;
289 }
290 }
291
292 if (foamyHexMeshControls().objOutput() && !fName.empty())
293 {
294 DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
295 }
296}
297
298
299void Foam::conformalVoronoiMesh::insertEdgePointGroups
300(
301 const pointIndexHitAndFeatureList& edgeHits,
302 const fileName fName,
303 DynamicList<Vb>& pts
304)
305{
306 forAll(edgeHits, i)
307 {
308 if (edgeHits[i].first().hit())
309 {
310 const extendedFeatureEdgeMesh& feMesh
311 (
312 geometryToConformTo_.features()[edgeHits[i].second()]
313 );
314
315// const bool isBaffle =
316// geometryToConformTo_.isBaffleFeature(edgeHits[i].second());
317
318 createEdgePointGroup
319 (
320 feMesh,
321 edgeHits[i].first(),
322 pts
323 );
324 }
325 }
326
327 if (foamyHexMeshControls().objOutput() && !fName.empty())
328 {
329 DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
330 }
331}
332
333
334bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
335{
336 scalar exclusionRangeSqr = featurePointExclusionDistanceSqr(pt);
337
338 pointIndexHit info;
339 label featureHit;
340
341 geometryToConformTo_.findFeaturePointNearest
342 (
343 pt,
344 exclusionRangeSqr,
345 info,
346 featureHit
347 );
348
349 return info.hit();
350}
351
352
353bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
354(
355 const Foam::point& pt
356) const
357{
358 scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
359
360 pointIndexHit info;
361 label featureHit;
362
363 geometryToConformTo_.findEdgeNearest
364 (
365 pt,
366 exclusionRangeSqr,
367 info,
368 featureHit
369 );
370
371 return info.hit();
372}
373
374
375void Foam::conformalVoronoiMesh::insertInitialPoints()
376{
377 Info<< nl << "Inserting initial points" << endl;
378
379 timeCheck("Before initial points call");
380
381 List<Point> initPts = initialPointsMethod_->initialPoints();
382
383 timeCheck("After initial points call");
384
385 // Assume that the initial points method made the correct decision for
386 // which processor each point should be on, so give distribute = false
387 insertInternalPoints(initPts);
388
389 if (initialPointsMethod_->fixInitialPoints())
390 {
391 for
392 (
393 Finite_vertices_iterator vit = finite_vertices_begin();
394 vit != finite_vertices_end();
395 ++vit
396 )
397 {
398 vit->fixed() = true;
399 }
400 }
401
402 if (foamyHexMeshControls().objOutput())
403 {
405 (
406 time().path()/"initialPoints.obj",
407 *this,
409 );
410 }
411}
412
413
414void Foam::conformalVoronoiMesh::distribute()
415{
416 if (!Pstream::parRun())
417 {
418 return;
419 }
420
421 DynamicList<Foam::point> points(number_of_vertices());
422 DynamicList<Foam::indexedVertexEnum::vertexType> types
423 (
424 number_of_vertices()
425 );
426 DynamicList<scalar> sizes(number_of_vertices());
427 DynamicList<tensor> alignments(number_of_vertices());
428
429 for
430 (
431 Finite_vertices_iterator vit = finite_vertices_begin();
432 vit != finite_vertices_end();
433 ++vit
434 )
435 {
436 if (vit->real())
437 {
438 points.append(topoint(vit->point()));
439 types.append(vit->type());
440 sizes.append(vit->targetCellSize());
441 alignments.append(vit->alignment());
442 }
443 }
444
445 autoPtr<mapDistribute> mapDist =
446 DistributedDelaunayMesh<Delaunay>::distribute(decomposition_(), points);
447
448 mapDist().distribute(types);
449 mapDist().distribute(sizes);
450 mapDist().distribute(alignments);
451
452 // Reset the entire tessellation
454
455 Info<< nl << " Inserting distributed tessellation" << endl;
456
457 // Internal points have to be inserted first
458
459 DynamicList<Vb> verticesToInsert(points.size());
460
461 forAll(points, pI)
462 {
463 verticesToInsert.append
464 (
465 Vb
466 (
467 toPoint(points[pI]),
468 -1,
469 types[pI],
471 )
472 );
473
474 verticesToInsert.last().targetCellSize() = sizes[pI];
475 verticesToInsert.last().alignment() = alignments[pI];
476 }
477
478 this->rangeInsertWithInfo
479 (
480 verticesToInsert.begin(),
481 verticesToInsert.end(),
482 true
483 );
484
485 Info<< " Total number of vertices after redistribution "
486 << returnReduce
487 (
488 label(number_of_vertices()), sumOp<label>()
489 )
490 << endl;
491}
492
493
494void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
495{
496 controlMeshRefinement meshRefinement
497 (
498 cellShapeControl_
499 );
500
501 smoothAlignmentSolver meshAlignmentSmoother
502 (
503 cellShapeControl_.shapeControlMesh()
504 );
505
506 meshRefinement.initialMeshPopulation(decomposition_);
507
508 cellShapeControlMesh& cellSizeMesh = cellShapeControl_.shapeControlMesh();
509
510 if (Pstream::parRun())
511 {
512 if (!distributeBackground(cellSizeMesh))
513 {
514 // Synchronise the cell size mesh if it has not been distributed
515 cellSizeMesh.distribute(decomposition_());
516 }
517 }
518
519 const dictionary& motionControlDict
520 = foamyHexMeshControls().foamyHexMeshDict().subDict("motionControl");
521
522 const label nMaxIter =
523 motionControlDict.get<label>("maxRefinementIterations");
524
525 Info<< "Maximum number of refinement iterations : " << nMaxIter << endl;
526
527 for (label i = 0; i < nMaxIter; ++i)
528 {
529 label nAdded = meshRefinement.refineMesh(decomposition_);
530 //cellShapeControl_.refineMesh(decomposition_);
531 reduce(nAdded, sumOp<label>());
532
533 if (Pstream::parRun())
534 {
535 cellSizeMesh.distribute(decomposition_());
536 }
537
538 Info<< " Iteration " << i
539 << " Added = " << nAdded << " points"
540 << endl;
541
542 if (nAdded == 0)
543 {
544 break;
545 }
546 }
547
548 if (Pstream::parRun())
549 {
550 // Need to distribute the cell size mesh to cover the background mesh
551 if (!distributeBackground(cellSizeMesh))
552 {
553 cellSizeMesh.distribute(decomposition_());
554 }
555 }
556
557 meshAlignmentSmoother.smoothAlignments
558 (
559 motionControlDict.get<label>("maxSmoothingIterations")
560 );
561
562 Info<< "Background cell size and alignment mesh:" << endl;
563 cellSizeMesh.printInfo(Info);
564
565 Info<< "Triangulation is "
566 << (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
567 << endl;
568
569 if (foamyHexMeshControls().writeCellShapeControlMesh())
570 {
571 //cellSizeMesh.writeTriangulation();
572 cellSizeMesh.write();
573 }
574
575 if (foamyHexMeshControls().printVertexInfo())
576 {
577 cellSizeMesh.printVertexInfo(Info);
578 }
579
580// Info<< "Estimated number of cells in final mesh = "
581// << returnReduce
582// (
583// cellSizeMesh.estimateCellCount(decomposition_),
584// sumOp<label>()
585// )
586// << endl;
587}
588
589
590void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
591{
592 Info<< nl << "Calculating target cell alignment and size" << endl;
593
594 for
595 (
596 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
597 vit != finite_vertices_end();
598 vit++
599 )
600 {
601 if (vit->internalOrBoundaryPoint())
602 {
603 pointFromPoint pt = topoint(vit->point());
604
605 cellShapeControls().cellSizeAndAlignment
606 (
607 pt,
608 vit->targetCellSize(),
609 vit->alignment()
610 );
611 }
612 }
613}
614
615
616Foam::face Foam::conformalVoronoiMesh::buildDualFace
617(
618 const Delaunay::Finite_edges_iterator& eit
619) const
620{
621 Cell_circulator ccStart = incident_cells(*eit);
622 Cell_circulator cc1 = ccStart;
623 Cell_circulator cc2 = cc1;
624
625 // Advance the second circulator so that it always stays on the next
626 // cell around the edge;
627 cc2++;
628
629 DynamicList<label> verticesOnFace;
630
631 label nUniqueVertices = 0;
632
633 do
634 {
635 if
636 (
637 cc1->hasFarPoint() || cc2->hasFarPoint()
638 || is_infinite(cc1) || is_infinite(cc2)
639 )
640 {
641 Cell_handle c = eit->first;
642 Vertex_handle vA = c->vertex(eit->second);
643 Vertex_handle vB = c->vertex(eit->third);
644
645// DelaunayMeshTools::drawDelaunayCell(Pout, cc1);
646// DelaunayMeshTools::drawDelaunayCell(Pout, cc2);
647
649 << "Dual face uses circumcenter defined by a "
650 << "Delaunay tetrahedron with no internal "
651 << "or boundary points. Defining Delaunay edge ends: "
652 << vA->info() << " "
653 << vB->info() << nl
654 <<endl;//<< exit(FatalError);
655 }
656 else
657 {
658 label cc1I = cc1->cellIndex();
659 label cc2I = cc2->cellIndex();
660
661 if (cc1I != cc2I)
662 {
663 if (!verticesOnFace.found(cc1I))
664 {
665 nUniqueVertices++;
666 }
667
668 verticesOnFace.append(cc1I);
669 }
670 }
671
672 cc1++;
673 cc2++;
674
675 } while (cc1 != ccStart);
676
677 verticesOnFace.shrink();
678
679 if (verticesOnFace.size() >= 3 && nUniqueVertices < 3)
680 {
681 // There are not enough unique vertices on this face to
682 // justify its size, it may have a form like:
683
684 // Vertices:
685 // A B
686 // A B
687
688 // Face:
689 // ABAB
690
691 // Setting the size to be below 3, so that it will not be
692 // created
693
694 verticesOnFace.setSize(nUniqueVertices);
695 }
696
697 return face(verticesOnFace);
698}
699
700
701Foam::label Foam::conformalVoronoiMesh::maxFilterCount
702(
703 const Delaunay::Finite_edges_iterator& eit
704) const
705{
706 Cell_circulator ccStart = incident_cells(*eit);
707 Cell_circulator cc = ccStart;
708
709 label maxFC = 0;
710
711 do
712 {
713 if (cc->hasFarPoint())
714 {
715 Cell_handle c = eit->first;
716 Vertex_handle vA = c->vertex(eit->second);
717 Vertex_handle vB = c->vertex(eit->third);
718
720 << "Dual face uses circumcenter defined by a "
721 << "Delaunay tetrahedron with no internal "
722 << "or boundary points. Defining Delaunay edge ends: "
723 << topoint(vA->point()) << " "
724 << topoint(vB->point()) << nl
725 << exit(FatalError);
726 }
727
728 if (cc->filterCount() > maxFC)
729 {
730 maxFC = cc->filterCount();
731 }
732
733 cc++;
734
735 } while (cc != ccStart);
736
737 return maxFC;
738}
739
740
741bool Foam::conformalVoronoiMesh::ownerAndNeighbour
742(
743 Vertex_handle vA,
744 Vertex_handle vB,
745 label& owner,
746 label& neighbour
747) const
748{
749 bool reverse = false;
750
751 owner = -1;
752 neighbour = -1;
753
754 label dualCellIndexA = vA->index();
755
756 if (!vA->internalOrBoundaryPoint() || vA->referred())
757 {
758 if (!vA->constrained())
759 {
760 dualCellIndexA = -1;
761 }
762 }
763
764 label dualCellIndexB = vB->index();
765
766 if (!vB->internalOrBoundaryPoint() || vB->referred())
767 {
768 if (!vB->constrained())
769 {
770 dualCellIndexB = -1;
771 }
772 }
773
774 if (dualCellIndexA == -1 && dualCellIndexB == -1)
775 {
777 << "Attempting to create a face joining "
778 << "two unindexed dual cells "
779 << exit(FatalError);
780 }
781 else if (dualCellIndexA == -1 || dualCellIndexB == -1)
782 {
783 // boundary face, find which is the owner
784
785 if (dualCellIndexA == -1)
786 {
787 owner = dualCellIndexB;
788
789 reverse = true;
790 }
791 else
792 {
793 owner = dualCellIndexA;
794 }
795 }
796 else
797 {
798 // internal face, find the lower cell to be the owner
799
800 if (dualCellIndexB > dualCellIndexA)
801 {
802 owner = dualCellIndexA;
803 neighbour = dualCellIndexB;
804 }
805 else
806 {
807 owner = dualCellIndexB;
808 neighbour = dualCellIndexA;
809
810 // reverse face order to correctly orientate normal
811 reverse = true;
812 }
813 }
814
815 return reverse;
816}
817
818
819// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
820
822(
823 const Time& runTime,
824 const dictionary& foamyHexMeshDict,
825 const fileName& decompDictFile
826)
827:
828 DistributedDelaunayMesh<Delaunay>(runTime),
829 runTime_(runTime),
830 rndGen_(64293*Pstream::myProcNo()),
831 foamyHexMeshControls_(foamyHexMeshDict),
832 allGeometry_
833 (
834 IOobject
835 (
836 "cvSearchableSurfaces",
837 runTime_.constant(),
838 "triSurface",
839 runTime_,
840 IOobject::MUST_READ,
841 IOobject::NO_WRITE
842 ),
843 foamyHexMeshDict.subDict("geometry"),
844 foamyHexMeshDict.getOrDefault("singleRegionName", true)
845 ),
846 geometryToConformTo_
847 (
848 runTime_,
849 rndGen_,
850 allGeometry_,
851 foamyHexMeshDict.subDict("surfaceConformation")
852 ),
853 decomposition_
854 (
855 Pstream::parRun()
856 ? new backgroundMeshDecomposition
857 (
858 runTime_,
859 rndGen_,
860 geometryToConformTo_,
861 foamyHexMeshControls().foamyHexMeshDict().subDict
862 (
863 "backgroundMeshDecomposition"
864 ),
865 decompDictFile
866 )
867 : nullptr
868 ),
869 cellShapeControl_
870 (
871 runTime_,
872 foamyHexMeshControls_,
873 allGeometry_,
874 geometryToConformTo_
875 ),
876 limitBounds_(),
877 ptPairs_(*this),
878 ftPtConformer_(*this),
879 edgeLocationTreePtr_(),
880 surfacePtLocationTreePtr_(),
881 surfaceConformationVertices_(),
882 initialPointsMethod_
883 (
884 initialPointsMethod::New
885 (
886 foamyHexMeshDict.subDict("initialPoints"),
887 runTime_,
888 rndGen_,
889 geometryToConformTo_,
890 cellShapeControl_,
891 decomposition_
892 )
893 ),
894 relaxationModel_
895 (
896 relaxationModel::New
897 (
898 foamyHexMeshDict.subDict("motionControl"),
899 runTime_
900 )
901 ),
902 faceAreaWeightModel_
903 (
904 faceAreaWeightModel::New
905 (
906 foamyHexMeshDict.subDict("motionControl")
907 )
908 )
909{}
910
911
912// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
913
915{}
916
917
918// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
919
921{
922 if (foamyHexMeshControls().objOutput())
923 {
924 geometryToConformTo_.writeFeatureObj("foamyHexMesh");
925 }
926
927 buildCellSizeAndAlignmentMesh();
928
929 insertInitialPoints();
930
931 insertFeaturePoints(true);
932
933 setVertexSizeAndAlignment();
934
935 cellSizeMeshOverlapsBackground();
936
937 // Improve the guess that the backgroundMeshDecomposition makes with the
938 // initial positions. Use before building the surface conformation to
939 // better balance the surface conformation load.
940 distributeBackground(*this);
941
942 buildSurfaceConformation();
943
944 // The introduction of the surface conformation may have distorted the
945 // balance of vertices, distribute if necessary.
946 distributeBackground(*this);
947
948 if (Pstream::parRun())
949 {
950 sync(decomposition_().procBounds());
951 }
952
953 // Do not store the surface conformation until after it has been
954 // (potentially) redistributed.
955 storeSurfaceConformation();
956
957 // Report any Delaunay vertices that do not think that they are in the
958 // domain the processor they are on.
959 // reportProcessorOccupancy();
960
961 cellSizeMeshOverlapsBackground();
962
963 if (foamyHexMeshControls().printVertexInfo())
964 {
965 printVertexInfo(Info);
966 }
967
968 if (foamyHexMeshControls().objOutput())
969 {
971 (
972 time().path()/"internalPoints_" + time().timeName() + ".obj",
973 *this,
976 );
977 }
978}
979
980
982{
983 if (Pstream::parRun())
984 {
985 decomposition_.reset
986 (
987 new backgroundMeshDecomposition
988 (
989 runTime_,
990 rndGen_,
991 geometryToConformTo_,
992 foamyHexMeshControls().foamyHexMeshDict().subDict
993 (
994 "backgroundMeshDecomposition"
995 )
996 )
997 );
998 }
999
1000 insertInitialPoints();
1001
1002 insertFeaturePoints();
1003
1004 // Improve the guess that the backgroundMeshDecomposition makes with the
1005 // initial positions. Use before building the surface conformation to
1006 // better balance the surface conformation load.
1007 distributeBackground(*this);
1008
1009 buildSurfaceConformation();
1010
1011 // The introduction of the surface conformation may have distorted the
1012 // balance of vertices, distribute if necessary.
1013 distributeBackground(*this);
1014
1015 if (Pstream::parRun())
1016 {
1017 sync(decomposition_().procBounds());
1018 }
1019
1020 cellSizeMeshOverlapsBackground();
1021
1022 if (foamyHexMeshControls().printVertexInfo())
1023 {
1024 printVertexInfo(Info);
1025 }
1026}
1027
1028
1030{
1031 timeCheck("Start of move");
1032
1033 scalar relaxation = relaxationModel_->relaxation();
1034
1035 Info<< nl << "Relaxation = " << relaxation << endl;
1036
1037 pointField dualVertices(number_of_finite_cells());
1038
1039 this->resetCellCount();
1040
1041 // Find the dual point of each tetrahedron and assign it an index.
1042 for
1043 (
1044 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1045 cit != finite_cells_end();
1046 ++cit
1047 )
1048 {
1049 cit->cellIndex() = Cb::ctUnassigned;
1050
1051 if (cit->anyInternalOrBoundaryDualVertex())
1052 {
1053 cit->cellIndex() = getNewCellIndex();
1054
1055 dualVertices[cit->cellIndex()] = cit->dual();
1056 }
1057
1058 if (cit->hasFarPoint())
1059 {
1060 cit->cellIndex() = Cb::ctFar;
1061 }
1062 }
1063
1064 dualVertices.setSize(cellCount());
1065
1066 setVertexSizeAndAlignment();
1067
1068 timeCheck("Determined sizes and alignments");
1069
1070 Info<< nl << "Determining vertex displacements" << endl;
1071
1072 vectorField cartesianDirections(3);
1073
1074 cartesianDirections[0] = vector(1, 0, 0);
1075 cartesianDirections[1] = vector(0, 1, 0);
1076 cartesianDirections[2] = vector(0, 0, 1);
1077
1078 vectorField displacementAccumulator
1079 (
1080 number_of_vertices(),
1081 Zero
1082 );
1083
1084 bitSet pointToBeRetained(number_of_vertices(), true);
1085
1086 DynamicList<Point> pointsToInsert(number_of_vertices());
1087
1088 for
1089 (
1090 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1091 eit != finite_edges_end();
1092 ++eit
1093 )
1094 {
1095 Cell_handle c = eit->first;
1096 Vertex_handle vA = c->vertex(eit->second);
1097 Vertex_handle vB = c->vertex(eit->third);
1098
1099 if
1100 (
1101 (
1102 vA->internalPoint() && !vA->referred()
1103 && vB->internalOrBoundaryPoint()
1104 )
1105 || (
1106 vB->internalPoint() && !vB->referred()
1107 && vA->internalOrBoundaryPoint()
1108 )
1109 )
1110 {
1111 pointFromPoint dVA = topoint(vA->point());
1112 pointFromPoint dVB = topoint(vB->point());
1113
1114 Field<vector> alignmentDirsA
1115 (
1116 vA->alignment().T() & cartesianDirections
1117 );
1118 Field<vector> alignmentDirsB
1119 (
1120 vB->alignment().T() & cartesianDirections
1121 );
1122
1123 Field<vector> alignmentDirs(alignmentDirsA);
1124
1125 forAll(alignmentDirsA, aA)
1126 {
1127 const vector& a = alignmentDirsA[aA];
1128
1129 scalar maxDotProduct = 0.0;
1130
1131 forAll(alignmentDirsB, aB)
1132 {
1133 const vector& b = alignmentDirsB[aB];
1134
1135 const scalar dotProduct = a & b;
1136
1137 if (mag(dotProduct) > maxDotProduct)
1138 {
1139 maxDotProduct = mag(dotProduct);
1140
1141 alignmentDirs[aA] = a + sign(dotProduct)*b;
1142
1143 alignmentDirs[aA].normalise();
1144 }
1145 }
1146 }
1147
1148 vector rAB = dVA - dVB;
1149
1150 scalar rABMag = mag(rAB);
1151
1152 if (rABMag < SMALL)
1153 {
1154 // Removal of close points
1155
1156 if
1157 (
1158 vA->internalPoint() && !vA->referred() && !vA->fixed()
1159 && vB->internalPoint() && !vB->referred() && !vB->fixed()
1160 )
1161 {
1162 // Only insert a point at the midpoint of
1163 // the short edge if neither attached
1164 // point has already been identified to be
1165 // removed.
1166
1167 if
1168 (
1169 pointToBeRetained.test(vA->index())
1170 && pointToBeRetained.test(vB->index())
1171 )
1172 {
1173 const Foam::point pt(0.5*(dVA + dVB));
1174
1175 if (internalPointIsInside(pt))
1176 {
1177 pointsToInsert.append(toPoint(pt));
1178 }
1179 }
1180 }
1181
1182 if (vA->internalPoint() && !vA->referred() && !vA->fixed())
1183 {
1184 pointToBeRetained.unset(vA->index());
1185 }
1186
1187 if (vB->internalPoint() && !vB->referred() && !vB->fixed())
1188 {
1189 pointToBeRetained.unset(vB->index());
1190 }
1191
1192 // Do not consider this Delaunay edge any further
1193
1194 continue;
1195 }
1196
1197 forAll(alignmentDirs, aD)
1198 {
1199 vector& alignmentDir = alignmentDirs[aD];
1200
1201 scalar dotProd = rAB & alignmentDir;
1202
1203 if (dotProd < 0)
1204 {
1205 // swap the direction of the alignment so that has the
1206 // same sense as rAB
1207 alignmentDir *= -1;
1208 dotProd *= -1;
1209 }
1210
1211 const scalar alignmentDotProd = dotProd/rABMag;
1212
1213 if
1214 (
1215 alignmentDotProd
1216 > foamyHexMeshControls().cosAlignmentAcceptanceAngle()
1217 )
1218 {
1219 scalar targetCellSize =
1221
1222 scalar targetFaceArea = sqr(targetCellSize);
1223
1224 const vector originalAlignmentDir = alignmentDir;
1225
1226 // Update cell size and face area
1227 cellShapeControls().aspectRatio().updateCellSizeAndFaceArea
1228 (
1229 alignmentDir,
1230 targetFaceArea,
1231 targetCellSize
1232 );
1233
1234 // Vector to move end points around middle of vector
1235 // to align edge (i.e. dual face normal) with alignment
1236 // directions.
1237 vector delta = alignmentDir - 0.5*rAB;
1238
1239 face dualFace = buildDualFace(eit);
1240
1241// Pout<< dualFace << endl;
1242// Pout<< " " << vA->info() << endl;
1243// Pout<< " " << vB->info() << endl;
1244
1245 const scalar faceArea = dualFace.mag(dualVertices);
1246
1247 // Update delta vector
1248 cellShapeControls().aspectRatio().updateDeltaVector
1249 (
1250 originalAlignmentDir,
1251 targetCellSize,
1252 rABMag,
1253 delta
1254 );
1255
1256// if (targetFaceArea == 0)
1257// {
1258// Pout<< vA->info() << vB->info();
1259//
1260// Cell_handle ch = locate(vA->point());
1261// if (is_infinite(ch))
1262// {
1263// Pout<< "vA " << vA->targetCellSize() << endl;
1264// }
1265//
1266// ch = locate(vB->point());
1267// if (is_infinite(ch))
1268// {
1269// Pout<< "vB " << vB->targetCellSize() << endl;
1270// }
1271// }
1272
1273 delta *= faceAreaWeightModel_->faceAreaWeight
1274 (
1275 faceArea/targetFaceArea
1276 );
1277
1278 if
1279 (
1280 (
1281 (vA->internalPoint() && vB->internalPoint())
1282 && (!vA->referred() || !vB->referred())
1283// ||
1284// (
1285// vA->referredInternalPoint()
1286// && vB->referredInternalPoint()
1287// )
1288 )
1289 && rABMag
1290 > foamyHexMeshControls().insertionDistCoeff()
1291 *targetCellSize
1292 && faceArea
1293 > foamyHexMeshControls().faceAreaRatioCoeff()
1294 *targetFaceArea
1295 && alignmentDotProd
1296 > foamyHexMeshControls().cosInsertionAcceptanceAngle()
1297 )
1298 {
1299 // Point insertion
1300 if
1301 (
1302 !geometryToConformTo_.findSurfaceAnyIntersection
1303 (
1304 dVA,
1305 dVB
1306 )
1307 )
1308 {
1309 const Foam::point newPt(0.5*(dVA + dVB));
1310
1311 // Prevent insertions spanning surfaces
1312 if (internalPointIsInside(newPt))
1313 {
1314 if (Pstream::parRun())
1315 {
1316 if
1317 (
1318 decomposition().positionOnThisProcessor
1319 (
1320 newPt
1321 )
1322 )
1323 {
1324 pointsToInsert.append(toPoint(newPt));
1325 }
1326 }
1327 else
1328 {
1329 pointsToInsert.append(toPoint(newPt));
1330 }
1331 }
1332 }
1333 }
1334 else if
1335 (
1336 (
1337 (vA->internalPoint() && !vA->referred())
1338 || (vB->internalPoint() && !vB->referred())
1339 )
1340 && rABMag
1341 < foamyHexMeshControls().removalDistCoeff()
1342 *targetCellSize
1343 )
1344 {
1345 // Point removal
1346 if
1347 (
1348 (
1349 vA->internalPoint()
1350 && !vA->referred()
1351 && !vA->fixed()
1352 )
1353 &&
1354 (
1355 vB->internalPoint()
1356 && !vB->referred()
1357 && !vB->fixed()
1358 )
1359 )
1360 {
1361 // Only insert a point at the midpoint of
1362 // the short edge if neither attached
1363 // point has already been identified to be
1364 // removed.
1365 if
1366 (
1367 pointToBeRetained.test(vA->index())
1368 && pointToBeRetained.test(vB->index())
1369 )
1370 {
1371 const Foam::point pt(0.5*(dVA + dVB));
1372
1373 if (internalPointIsInside(pt))
1374 {
1375 pointsToInsert.append(toPoint(pt));
1376 }
1377 }
1378 }
1379
1380 if
1381 (
1382 vA->internalPoint()
1383 && !vA->referred()
1384 && !vA->fixed()
1385 )
1386 {
1387 pointToBeRetained.unset(vA->index());
1388 }
1389
1390 if
1391 (
1392 vB->internalPoint()
1393 && !vB->referred()
1394 && !vB->fixed()
1395 )
1396 {
1397 pointToBeRetained.unset(vB->index());
1398 }
1399 }
1400 else
1401 {
1402 if
1403 (
1404 vA->internalPoint()
1405 && !vA->referred()
1406 && !vA->fixed()
1407 )
1408 {
1409 if (vB->fixed())
1410 {
1411 displacementAccumulator[vA->index()] += 2*delta;
1412 }
1413 else
1414 {
1415 displacementAccumulator[vA->index()] += delta;
1416 }
1417 }
1418
1419 if
1420 (
1421 vB->internalPoint()
1422 && !vB->referred()
1423 && !vB->fixed()
1424 )
1425 {
1426 if (vA->fixed())
1427 {
1428 displacementAccumulator[vB->index()] -= 2*delta;
1429 }
1430 else
1431 {
1432 displacementAccumulator[vB->index()] -= delta;
1433 }
1434 }
1435 }
1436 }
1437 }
1438 }
1439 }
1440
1441 Info<< "Limit displacements" << endl;
1442
1443 // Limit displacements that pierce, or get too close to the surface
1444 for
1445 (
1446 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1447 vit != finite_vertices_end();
1448 ++vit
1449 )
1450 {
1451 if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1452 {
1453 if (pointToBeRetained.test(vit->index()))
1454 {
1455 limitDisplacement
1456 (
1457 vit,
1458 displacementAccumulator[vit->index()]
1459 );
1460 }
1461 }
1462 }
1463
1464 vector totalDisp = gSum(displacementAccumulator);
1465 scalar totalDist = gSum(mag(displacementAccumulator));
1466
1467 displacementAccumulator *= relaxation;
1468
1469 Info<< "Sum displacements" << endl;
1470
1471 label nPointsToRetain = 0;
1472 label nPointsToRemove = 0;
1473
1474 for
1475 (
1476 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1477 vit != finite_vertices_end();
1478 ++vit
1479 )
1480 {
1481 if (vit->internalPoint() && !vit->referred() && !vit->fixed())
1482 {
1483 if (pointToBeRetained.test(vit->index()))
1484 {
1485 // Convert vit->point() to FOAM vector (double) to do addition,
1486 // avoids memory increase because a record of the constructions
1487 // would be kept otherwise.
1488 // See cgal-discuss@lists-sop.inria.fr:
1489 // "Memory issue with openSUSE 11.3, exact kernel, adding
1490 // points/vectors"
1491 // 14/1/2011.
1492 // Only necessary if using an exact constructions kernel
1493 // (extended precision)
1494 Foam::point pt
1495 (
1496 topoint(vit->point())
1497 + displacementAccumulator[vit->index()]
1498 );
1499
1500 if (internalPointIsInside(pt))
1501 {
1502 pointsToInsert.append(toPoint(pt));
1503 nPointsToRemove++;
1504 }
1505
1506 nPointsToRetain++;
1507 }
1508 }
1509 }
1510
1511 pointsToInsert.shrink();
1512
1514 (
1515 nPointsToRetain - nPointsToRemove,
1516 sumOp<label>()
1517 )
1518 << " internal points are outside the domain. "
1519 << "They will not be inserted." << endl;
1520
1521 // Save displacements to file.
1522 if (foamyHexMeshControls().objOutput() && time().writeTime())
1523 {
1524 Info<< "Writing point displacement vectors to file." << endl;
1525 OFstream str
1526 (
1527 time().path()/"displacements_" + runTime_.timeName() + ".obj"
1528 );
1529
1530 for
1531 (
1532 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1533 vit != finite_vertices_end();
1534 ++vit
1535 )
1536 {
1537 if (vit->internalPoint() && !vit->referred())
1538 {
1539 if (pointToBeRetained.test(vit->index()))
1540 {
1541 meshTools::writeOBJ(str, topoint(vit->point()));
1542
1543 str << "vn "
1544 << displacementAccumulator[vit->index()][0] << " "
1545 << displacementAccumulator[vit->index()][1] << " "
1546 << displacementAccumulator[vit->index()][2] << " "
1547 << endl;
1548 }
1549 }
1550 }
1551 }
1552
1553 // Remove the entire tessellation
1555
1556 timeCheck("Displacement calculated");
1557
1558 Info<< nl<< "Inserting displaced tessellation" << endl;
1559
1560 insertFeaturePoints(true);
1561
1562 insertInternalPoints(pointsToInsert, true);
1563
1564 // Fix points that have not been significantly displaced
1565// for
1566// (
1567// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1568// vit != finite_vertices_end();
1569// ++vit
1570// )
1571// {
1572// if (vit->internalPoint())
1573// {
1574// if
1575// (
1576// mag(displacementAccumulator[vit->index()])
1577// < 0.1*targetCellSize(topoint(vit->point()))
1578// )
1579// {
1580// vit->setVertexFixed();
1581// }
1582// }
1583// }
1584
1585 timeCheck("Internal points inserted");
1586
1587 {
1588 // Check that no index is shared between any of the local points
1589 labelHashSet usedIndices;
1590 for
1591 (
1592 Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
1593 vit != finite_vertices_end();
1594 ++vit
1595 )
1596 {
1597 if (!vit->referred() && !usedIndices.insert(vit->index()))
1598 {
1600 << "Index already used! Could not insert: " << nl
1601 << vit->info()
1602 << abort(FatalError);
1603 }
1604 }
1605 }
1606
1607 conformToSurface();
1608
1609 if (foamyHexMeshControls().objOutput())
1610 {
1612 (
1613 time().path()/"internalPoints_" + time().timeName() + ".obj",
1614 *this,
1616 );
1617
1618 if (reconformToSurface())
1619 {
1621 (
1622 time().path()/"boundaryPoints_" + time().timeName() + ".obj",
1623 *this
1624 );
1625
1627 (
1628 time().path()/"internalBoundaryPoints_" + time().timeName()
1629 + ".obj",
1630 *this,
1633 );
1634
1636 (
1637 time().path()/"externalBoundaryPoints_" + time().timeName()
1638 + ".obj",
1639 *this,
1642 );
1643
1644 OBJstream multipleIntersections
1645 (
1646 "multipleIntersections_"
1647 + time().timeName()
1648 + ".obj"
1649 );
1650
1651 for
1652 (
1653 Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1654 eit != finite_edges_end();
1655 ++eit
1656 )
1657 {
1658 Cell_handle c = eit->first;
1659 Vertex_handle vA = c->vertex(eit->second);
1660 Vertex_handle vB = c->vertex(eit->third);
1661
1662 Foam::point ptA(topoint(vA->point()));
1663 Foam::point ptB(topoint(vB->point()));
1664
1665 List<pointIndexHit> surfHits;
1666 labelList hitSurfaces;
1667
1668 geometryToConformTo().findSurfaceAllIntersections
1669 (
1670 ptA,
1671 ptB,
1672 surfHits,
1673 hitSurfaces
1674 );
1675
1676 if
1677 (
1678 surfHits.size() >= 2
1679 || (
1680 surfHits.size() == 0
1681 && (
1682 (vA->externalBoundaryPoint()
1683 && vB->internalBoundaryPoint())
1684 || (vB->externalBoundaryPoint()
1685 && vA->internalBoundaryPoint())
1686 )
1687 )
1688 )
1689 {
1690 multipleIntersections.write(linePointRef(ptA, ptB));
1691 }
1692 }
1693 }
1694 }
1695
1696 timeCheck("After conformToSurface");
1697
1698 if (foamyHexMeshControls().printVertexInfo())
1699 {
1700 printVertexInfo(Info);
1701 }
1702
1703 if (time().writeTime())
1704 {
1705 writeMesh(time().timeName());
1706 }
1707
1708 Info<< nl
1709 << "Total displacement = " << totalDisp << nl
1710 << "Total distance = " << totalDist << nl
1711 << endl;
1712}
1713
1714
1715void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
1716{
1717 typedef CGAL::Exact_predicates_exact_constructions_kernel Kexact;
1718 typedef CGAL::Point_3<Kexact> PointExact;
1719
1720 if (!is_valid())
1721 {
1722 Pout<< "Triangulation is invalid!" << endl;
1723 }
1724
1725 OFstream str("badCells.obj");
1726
1727 label badCells = 0;
1728
1729 for
1730 (
1731 Delaunay::Finite_cells_iterator cit = finite_cells_begin();
1732 cit != finite_cells_end();
1733 ++cit
1734 )
1735 {
1736 const scalar quality = foamyHexMeshChecks::coplanarTet(cit, 1e-16);
1737
1738 if (quality == 0)
1739 {
1740 Pout<< "COPLANAR: " << cit->info() << nl
1741 << " quality = " << quality << nl
1742 << " dual = " << topoint(cit->dual()) << endl;
1743
1744 DelaunayMeshTools::drawDelaunayCell(str, cit, badCells++);
1745
1746 FixedList<PointExact, 4> cellVerticesExact(PointExact(0,0,0));
1747 forAll(cellVerticesExact, vI)
1748 {
1749 cellVerticesExact[vI] = PointExact
1750 (
1751 cit->vertex(vI)->point().x(),
1752 cit->vertex(vI)->point().y(),
1753 cit->vertex(vI)->point().z()
1754 );
1755 }
1756
1757 PointExact synchronisedDual = CGAL::circumcenter<Kexact>
1758 (
1759 cellVerticesExact[0],
1760 cellVerticesExact[1],
1761 cellVerticesExact[2],
1762 cellVerticesExact[3]
1763 );
1764
1765 Foam::point exactPt
1766 (
1767 CGAL::to_double(synchronisedDual.x()),
1768 CGAL::to_double(synchronisedDual.y()),
1769 CGAL::to_double(synchronisedDual.z())
1770 );
1771
1772 Info<< "inexact = " << cit->dual() << nl
1773 << "exact = " << exactPt << endl;
1774 }
1775 }
1776
1777 Pout<< "There are " << badCells << " bad cells out of "
1778 << number_of_finite_cells() << endl;
1779
1780
1781 label nNonGabriel = 0;
1782 for
1783 (
1784 Delaunay::Finite_facets_iterator fit = finite_facets_begin();
1785 fit != finite_facets_end();
1786 ++fit
1787 )
1788 {
1789 if (!is_Gabriel(*fit))
1790 {
1791 nNonGabriel++;//Pout<< "Non-gabriel face" << endl;
1792 }
1793 }
1794
1795 Pout<< "There are " << nNonGabriel << " non-Gabriel faces out of "
1796 << number_of_finite_facets() << endl;
1797}
1798
1799
1800// ************************************************************************* //
CGAL::Delaunay_triangulation_3< K, Tds, CompactLocator > Delaunay
Istream and Ostream manipulators taking arguments.
scalar delta
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep track of the Delaunay vertices i...
Definition: indexedVertex.H:88
void insertSurfacePointPairs()
Insert all surface point-pairs from.
Map< label > insertPoints(const List< Vb > &vertices, const bool reIndex)
Insert the list of vertices (calls rangeInsertWithInfo)
void reset()
Clear the entire triangulation.
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: Enum.H:61
void normalise()
Normalise the field inplace. See notes in Field.
InfoProxy< IOstream > info() const
Return info proxy.
Definition: IOstream.H:413
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
iterator begin() noexcept
Return an iterator to begin traversing the UList.
Definition: UListI.H:329
iterator end() noexcept
Return an iterator to end traversing the UList.
Definition: UListI.H:350
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
cellShapeControlMesh & shapeControlMesh()
static const Enum< dualMeshPointType > dualMeshPointTypeNames_
~conformalVoronoiMesh()
Destructor.
void move()
Move the vertices according to the controller, re-conforming to the.
sideVolumeType
Normals point to the outside.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
@ BOTH
limit by both minimum and maximum values
Definition: limitFields.H:189
constant condensation/saturation model.
int myProcNo() const noexcept
Return processor number.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
engineTime & runTime
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
word timeName
Definition: getTimeIndex.H:3
#define WarningInFunction
Report a warning using Foam::Warning.
Foam::scalar averageCellSize(const VertexType &vA, const VertexType &vB)
Return the target cell size from that stored on a pair of Delaunay vertices,.
bool uninitialised(const VertexType &v)
void drawDelaunayCell(Ostream &os, const CellHandle &c, label offset=0)
Draws a tet cell to an output stream. The offset is supplied as the tet.
void writeOBJ(const fileName &fName, const List< Foam::point > &points)
Write list of points to file.
void writeBoundaryPoints(const fileName &fName, const Triangulation &t)
Write the boundary Delaunay points to an OBJ file.
const dimensionedScalar c
Speed of light in a vacuum.
scalar coplanarTet(Cell &c, const scalar tol=1e-12)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
Type gSum(const FieldField< Field, Type > &f)
dimensionedScalar sign(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
dimensionedSymmTensor sqr(const dimensionedVector &dv)
PointFrompoint toPoint(const Foam::point &p)
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
pointFromPoint topoint(const Point &P)
void reverse(UList< T > &list, const label n)
Reverse the first n elements of the list.
Definition: UListI.H:449
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
pointField vertices(const blockVertexList &bvl)
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
error FatalError
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce (copy) and return value.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Vector< scalar > vector
Definition: vector.H:61
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
srcOptions insert("case", fileName(rootDirSource/caseDirSource))
volScalarField & b
Definition: createFields.H:27
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333