polyMesh.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) 2011-2017, 2020 OpenFOAM Foundation
9 Copyright (C) 2016-2022 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
29#include "polyMesh.H"
30#include "Time.H"
31#include "cellIOList.H"
32#include "wedgePolyPatch.H"
33#include "emptyPolyPatch.H"
34#include "globalMeshData.H"
35#include "processorPolyPatch.H"
37#include "indexedOctree.H"
38#include "treeDataCell.H"
39#include "MeshObject.H"
40#include "pointMesh.H"
41
42// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43
44namespace Foam
45{
47}
48
51
52
53// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54
55void Foam::polyMesh::calcDirections() const
56{
57 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
58 {
59 solutionD_[cmpt] = 1;
60 }
61
62 // Knock out empty and wedge directions. Note:they will be present on all
63 // domains.
64
65 label nEmptyPatches = 0;
66 label nWedgePatches = 0;
67
68 vector emptyDirVec = Zero;
69 vector wedgeDirVec = Zero;
70
71 forAll(boundaryMesh(), patchi)
72 {
73 const polyPatch& pp = boundaryMesh()[patchi];
74 if (isA<emptyPolyPatch>(pp))
75 {
76 // Force calculation of geometric properties, independent of
77 // size. This avoids parallel synchronisation problems.
78 const vectorField::subField fa(pp.faceAreas());
79
80 if (pp.size())
81 {
82 nEmptyPatches++;
83 emptyDirVec += sum(cmptMag(fa));
84 }
85 }
86 else if (isA<wedgePolyPatch>(pp))
87 {
88 const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>(pp);
89
90 // Force calculation of geometric properties, independent of
91 // size. This avoids parallel synchronisation problems.
92 (void)wpp.faceNormals();
93
94 if (pp.size())
95 {
96 nWedgePatches++;
97 wedgeDirVec += cmptMag(wpp.centreNormal());
98 }
99 }
100 }
101
102 reduce(nEmptyPatches, maxOp<label>());
103 reduce(nWedgePatches, maxOp<label>());
104
105 if (nEmptyPatches)
106 {
107 reduce(emptyDirVec, sumOp<vector>());
108
109 emptyDirVec.normalise();
110
111 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
112 {
113 if (emptyDirVec[cmpt] > 1e-6)
114 {
115 solutionD_[cmpt] = -1;
116 }
117 else
118 {
119 solutionD_[cmpt] = 1;
120 }
121 }
122 }
123
124
125 // Knock out wedge directions
126
127 geometricD_ = solutionD_;
128
129 if (nWedgePatches)
130 {
131 reduce(wedgeDirVec, sumOp<vector>());
132
133 wedgeDirVec.normalise();
134
135 for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
136 {
137 if (wedgeDirVec[cmpt] > 1e-6)
138 {
139 geometricD_[cmpt] = -1;
140 }
141 else
142 {
143 geometricD_[cmpt] = 1;
144 }
145 }
146 }
147}
148
149
150Foam::autoPtr<Foam::labelIOList> Foam::polyMesh::readTetBasePtIs() const
151{
152 IOobject io
153 (
154 "tetBasePtIs",
155 instance(),
156 meshSubDir,
157 *this,
160 );
161
162 if (io.typeHeaderOk<labelIOList>(true))
163 {
165 }
166
167 return nullptr;
168}
169
170
171// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
172
173Foam::polyMesh::polyMesh(const IOobject& io, const bool doInit)
174:
177 points_
178 (
180 (
181 "points",
182 time().findInstance(meshDir(), "points"),
183 meshSubDir,
184 *this,
185 IOobject::MUST_READ,
186 IOobject::NO_WRITE
187 )
188 ),
189 faces_
190 (
192 (
193 "faces",
194 time().findInstance(meshDir(), "faces"),
195 meshSubDir,
196 *this,
197 IOobject::MUST_READ,
198 IOobject::NO_WRITE
199 )
200 ),
201 owner_
202 (
204 (
205 "owner",
206 faces_.instance(),
207 meshSubDir,
208 *this,
209 IOobject::READ_IF_PRESENT,
210 IOobject::NO_WRITE
211 )
212 ),
213 neighbour_
214 (
216 (
217 "neighbour",
218 faces_.instance(),
219 meshSubDir,
220 *this,
221 IOobject::READ_IF_PRESENT,
222 IOobject::NO_WRITE
223 )
224 ),
225 clearedPrimitives_(false),
226 boundary_
227 (
229 (
230 "boundary",
231 time().findInstance // allow 'newer' boundary file
232 (
233 meshDir(),
234 "boundary",
235 IOobject::MUST_READ,
236 faces_.instance()
237 ),
238 meshSubDir,
239 *this,
240 IOobject::MUST_READ,
241 IOobject::NO_WRITE
242 ),
243 *this
244 ),
245 bounds_(points_),
246 comm_(UPstream::worldComm),
247 geometricD_(Zero),
248 solutionD_(Zero),
249 tetBasePtIsPtr_(readTetBasePtIs()),
250 cellTreePtr_(nullptr),
251 pointZones_
252 (
254 (
255 "pointZones",
256 //time().findInstance
257 //(
258 // meshDir(),
259 // "pointZones",
260 // IOobject::READ_IF_PRESENT
261 //),
262 faces_.instance(),
263 meshSubDir,
264 *this,
265 IOobject::READ_IF_PRESENT,
266 IOobject::NO_WRITE
267 ),
268 *this
269 ),
270 faceZones_
271 (
273 (
274 "faceZones",
275 //time().findInstance
276 //(
277 // meshDir(),
278 // "faceZones",
279 // IOobject::READ_IF_PRESENT
280 //),
281 faces_.instance(),
282 meshSubDir,
283 *this,
284 IOobject::READ_IF_PRESENT,
285 IOobject::NO_WRITE
286 ),
287 *this
288 ),
289 cellZones_
290 (
292 (
293 "cellZones",
294 //time().findInstance
295 //(
296 // meshDir(),
297 // "cellZones",
298 // IOobject::READ_IF_PRESENT
299 //),
300 faces_.instance(),
301 meshSubDir,
302 *this,
303 IOobject::READ_IF_PRESENT,
304 IOobject::NO_WRITE
305 ),
306 *this
307 ),
308 globalMeshDataPtr_(nullptr),
309 moving_(false),
310 topoChanging_(false),
311 storeOldCellCentres_(false),
312 curMotionTimeIndex_(time().timeIndex()),
313 oldPointsPtr_(nullptr),
314 oldCellCentresPtr_(nullptr)
315{
316 if (owner_.hasHeaderClass())
317 {
318 initMesh();
319 }
320 else
321 {
323 (
325 (
326 "cells",
327 time().findInstance(meshDir(), "cells"),
329 *this,
332 )
333 );
334
335 // Set the primitive mesh
336 initMesh(cLst);
337
338 owner_.write();
339 neighbour_.write();
340 }
341
342 // Warn if global empty mesh
343 if (returnReduce(boundary_.empty(), orOp<bool>()))
344 {
346 << "mesh missing boundary on one or more domains" << endl;
347
348 if (returnReduce(nPoints(), sumOp<label>()) == 0)
349 {
351 << "no points in mesh" << endl;
352 }
353 if (returnReduce(nCells(), sumOp<label>()) == 0)
354 {
356 << "no cells in mesh" << endl;
357 }
358 }
359
360 if (doInit)
361 {
362 polyMesh::init(false); // do not init lower levels
363 }
364}
365
366
367bool Foam::polyMesh::init(const bool doInit)
368{
369 if (doInit)
370 {
371 primitiveMesh::init(doInit);
372 }
373
374 // Calculate topology for the patches (processor-processor comms etc.)
375 boundary_.updateMesh();
376
377 // Calculate the geometry for the patches (transformation tensors etc.)
378 boundary_.calcGeometry();
379
380 // Initialise demand-driven data
381 calcDirections();
382
383 return false;
384}
385
386
388(
389 const IOobject& io,
391 faceList&& faces,
392 labelList&& owner,
393 labelList&& neighbour,
394 const bool syncPar
395)
396:
399 points_
400 (
402 (
403 "points",
404 instance(),
405 meshSubDir,
406 *this,
407 IOobject::NO_READ, //io.readOpt(),
408 io.writeOpt()
409 ),
410 std::move(points)
411 ),
412 faces_
413 (
415 (
416 "faces",
417 instance(),
418 meshSubDir,
419 *this,
420 IOobject::NO_READ, //io.readOpt(),
421 io.writeOpt()
422 ),
423 std::move(faces)
424 ),
425 owner_
426 (
428 (
429 "owner",
430 instance(),
431 meshSubDir,
432 *this,
433 IOobject::NO_READ, //io.readOpt(),
434 io.writeOpt()
435 ),
436 std::move(owner)
437 ),
438 neighbour_
439 (
441 (
442 "neighbour",
443 instance(),
444 meshSubDir,
445 *this,
446 IOobject::NO_READ, //io.readOpt(),
447 io.writeOpt()
448 ),
449 std::move(neighbour)
450 ),
451 clearedPrimitives_(false),
452 boundary_
453 (
455 (
456 "boundary",
457 instance(),
458 meshSubDir,
459 *this,
460 IOobject::NO_READ, // ignore since no alternative can be supplied
461 io.writeOpt()
462 ),
463 *this,
465 ),
466 bounds_(points_, syncPar),
467 comm_(UPstream::worldComm),
468 geometricD_(Zero),
469 solutionD_(Zero),
470 tetBasePtIsPtr_(nullptr),
471 cellTreePtr_(nullptr),
472 pointZones_
473 (
475 (
476 "pointZones",
477 instance(),
478 meshSubDir,
479 *this,
480 IOobject::NO_READ, // ignore since no alternative can be supplied
481 IOobject::NO_WRITE
482 ),
483 *this,
485 ),
486 faceZones_
487 (
489 (
490 "faceZones",
491 instance(),
492 meshSubDir,
493 *this,
494 IOobject::NO_READ,// ignore since no alternative can be supplied
495 IOobject::NO_WRITE
496 ),
497 *this,
499 ),
500 cellZones_
501 (
503 (
504 "cellZones",
505 instance(),
506 meshSubDir,
507 *this,
508 IOobject::NO_READ, // ignore since no alternative can be supplied
509 IOobject::NO_WRITE
510 ),
511 *this,
513 ),
514 globalMeshDataPtr_(nullptr),
515 moving_(false),
516 topoChanging_(false),
517 storeOldCellCentres_(false),
518 curMotionTimeIndex_(time().timeIndex()),
519 oldPointsPtr_(nullptr),
520 oldCellCentresPtr_(nullptr)
521{
522 // Check if the faces and cells are valid
523 forAll(faces_, facei)
524 {
525 const face& curFace = faces_[facei];
526
527 if (min(curFace) < 0 || max(curFace) > points_.size())
528 {
530 << "Face " << facei << "contains vertex labels out of range: "
531 << curFace << " Max point index = " << points_.size()
532 << abort(FatalError);
533 }
534 }
535
536 // Set the primitive mesh
537 initMesh();
538}
539
540
542(
543 const IOobject& io,
545 faceList&& faces,
546 cellList&& cells,
547 const bool syncPar
548)
549:
552 points_
553 (
555 (
556 "points",
557 instance(),
558 meshSubDir,
559 *this,
560 IOobject::NO_READ,
561 io.writeOpt()
562 ),
563 std::move(points)
564 ),
565 faces_
566 (
568 (
569 "faces",
570 instance(),
571 meshSubDir,
572 *this,
573 IOobject::NO_READ,
574 io.writeOpt()
575 ),
576 std::move(faces)
577 ),
578 owner_
579 (
581 (
582 "owner",
583 instance(),
584 meshSubDir,
585 *this,
586 IOobject::NO_READ,
587 io.writeOpt()
588 ),
589 0
590 ),
591 neighbour_
592 (
594 (
595 "neighbour",
596 instance(),
597 meshSubDir,
598 *this,
599 IOobject::NO_READ,
600 io.writeOpt()
601 ),
602 0
603 ),
604 clearedPrimitives_(false),
605 boundary_
606 (
608 (
609 "boundary",
610 instance(),
611 meshSubDir,
612 *this,
613 IOobject::NO_READ,
614 io.writeOpt()
615 ),
616 *this,
617 0
618 ),
619 bounds_(points_, syncPar),
620 comm_(UPstream::worldComm),
621 geometricD_(Zero),
622 solutionD_(Zero),
623 tetBasePtIsPtr_(nullptr),
624 cellTreePtr_(nullptr),
625 pointZones_
626 (
628 (
629 "pointZones",
630 instance(),
631 meshSubDir,
632 *this,
633 IOobject::NO_READ,
634 IOobject::NO_WRITE
635 ),
636 *this,
637 0
638 ),
639 faceZones_
640 (
642 (
643 "faceZones",
644 instance(),
645 meshSubDir,
646 *this,
647 IOobject::NO_READ,
648 IOobject::NO_WRITE
649 ),
650 *this,
651 0
652 ),
653 cellZones_
654 (
656 (
657 "cellZones",
658 instance(),
659 meshSubDir,
660 *this,
661 IOobject::NO_READ,
662 IOobject::NO_WRITE
663 ),
664 *this,
665 0
666 ),
667 globalMeshDataPtr_(nullptr),
668 moving_(false),
669 topoChanging_(false),
670 storeOldCellCentres_(false),
671 curMotionTimeIndex_(time().timeIndex()),
672 oldPointsPtr_(nullptr),
673 oldCellCentresPtr_(nullptr)
674{
675 // Check if faces are valid
676 forAll(faces_, facei)
677 {
678 const face& curFace = faces_[facei];
679
680 if (min(curFace) < 0 || max(curFace) > points_.size())
681 {
683 << "Face " << facei << "contains vertex labels out of range: "
684 << curFace << " Max point index = " << points_.size()
685 << abort(FatalError);
686 }
687 }
688
689 // Transfer in cell list
690 cellList cLst(std::move(cells));
691
692 // Check if cells are valid
693 forAll(cLst, celli)
694 {
695 const cell& curCell = cLst[celli];
696
697 if (min(curCell) < 0 || max(curCell) > faces_.size())
698 {
700 << "Cell " << celli << "contains face labels out of range: "
701 << curCell << " Max face index = " << faces_.size()
702 << abort(FatalError);
703 }
704 }
705
706 // Set the primitive mesh
707 initMesh(cLst);
708}
709
710
711Foam::polyMesh::polyMesh(const IOobject& io, const zero, const bool syncPar)
712:
713 polyMesh(io, pointField(), faceList(), labelList(), labelList(), syncPar)
714{}
715
716
718(
720 autoPtr<faceList>&& faces,
721 autoPtr<labelList>&& owner,
722 autoPtr<labelList>&& neighbour,
723 const labelUList& patchSizes,
724 const labelUList& patchStarts,
725 const bool validBoundary
726)
727{
728 // Clear addressing. Keep geometric props and updateable props for mapping.
729 clearAddressing(true);
730
731 // Take over new primitive data.
732 // Optimized to avoid overwriting data at all
733 if (points)
734 {
735 points_.transfer(*points);
736 bounds_ = boundBox(points_, validBoundary);
737 }
738
739 if (faces)
740 {
741 faces_.transfer(*faces);
742 }
743
744 if (owner)
745 {
746 owner_.transfer(*owner);
747 }
748
749 if (neighbour)
750 {
751 neighbour_.transfer(*neighbour);
752 }
753
754
755 // Reset patch sizes and starts
756 forAll(boundary_, patchi)
757 {
758 boundary_[patchi] = polyPatch
759 (
760 boundary_[patchi],
761 boundary_,
762 patchi,
763 patchSizes[patchi],
764 patchStarts[patchi]
765 );
766 }
767
768
769 // Flags the mesh files as being changed
770 setInstance(time().timeName());
771
772 // Check if the faces and cells are valid
773 forAll(faces_, facei)
774 {
775 const face& curFace = faces_[facei];
776
777 if (min(curFace) < 0 || max(curFace) > points_.size())
778 {
780 << "Face " << facei << " contains vertex labels out of range: "
781 << curFace << " Max point index = " << points_.size()
782 << abort(FatalError);
783 }
784 }
785
786
787 // Set the primitive mesh from the owner_, neighbour_.
788 // Works out from patch end where the active faces stop.
789 initMesh();
790
791
792 if (validBoundary)
793 {
794 // Note that we assume that all the patches stay the same and are
795 // correct etc. so we can already use the patches to do
796 // processor-processor comms.
797
798 // Calculate topology for the patches (processor-processor comms etc.)
799 boundary_.updateMesh();
800
801 // Calculate the geometry for the patches (transformation tensors etc.)
802 boundary_.calcGeometry();
803
804 // Warn if global empty mesh
805 if
806 (
808 || (returnReduce(nCells(), sumOp<label>()) == 0)
809 )
810 {
812 << "no points or no cells in mesh" << endl;
813 }
814 }
815}
816
817
818// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
819
821{
822 clearOut();
823 resetMotion();
824}
825
826
827// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
828
830{
831 return (region == polyMesh::defaultRegion ? word::null : region);
832}
833
834
835// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
836
838{
839 if (objectRegistry::dbDir() == defaultRegion)
840 {
841 return parent().dbDir();
842 }
843
844 return objectRegistry::dbDir();
845}
846
847
849{
850 return dbDir()/meshSubDir;
851}
852
853
855{
857}
858
859
861{
862 return points_.instance();
863}
864
865
867{
868 return faces_.instance();
869}
870
871
873{
874 if (geometricD_.x() == 0)
875 {
876 calcDirections();
877 }
878
879 return geometricD_;
880}
881
882
884{
885 return cmptSum(geometricD() + Vector<label>::one)/2;
886}
887
888
890{
891 if (solutionD_.x() == 0)
892 {
893 calcDirections();
894 }
895
896 return solutionD_;
897}
898
899
900Foam::label Foam::polyMesh::nSolutionD() const
901{
902 return cmptSum(solutionD() + Vector<label>::one)/2;
903}
904
905
907{
908 if (!tetBasePtIsPtr_)
909 {
910 if (debug)
911 {
913 << "Forcing storage of base points."
914 << endl;
915 }
916
917 tetBasePtIsPtr_.reset
918 (
919 new labelIOList
920 (
922 (
923 "tetBasePtIs",
924 instance(),
925 meshSubDir,
926 *this,
929 ),
931 )
932 );
933 }
934
935 return *tetBasePtIsPtr_;
936}
937
938
941{
942 if (!cellTreePtr_)
943 {
944 treeBoundBox overallBb(points());
945
946 Random rndGen(261782);
947
948 overallBb = overallBb.extend(rndGen, 1e-4);
949 overallBb.min() -= point::uniform(ROOTVSMALL);
950 overallBb.max() += point::uniform(ROOTVSMALL);
951
952 cellTreePtr_.reset
953 (
955 (
957 (
958 false, // not cache bb
959 *this,
960 CELL_TETS // use tet-decomposition for any inside test
961 ),
962 overallBb,
963 8, // maxLevel
964 10, // leafsize
965 5.0 // duplicity
966 )
967 );
968 }
969
970 return *cellTreePtr_;
971}
972
973
975(
976 polyPatchList& plist,
977 const bool validBoundary
978)
979{
980 if (boundaryMesh().size())
981 {
983 << "boundary already exists"
984 << abort(FatalError);
985 }
986
987 // Reset valid directions
988 geometricD_ = Zero;
989 solutionD_ = Zero;
990
991 boundary_.transfer(plist);
992
993 // parallelData depends on the processorPatch ordering so force
994 // recalculation. Problem: should really be done in removeBoundary but
995 // there is some info in parallelData which might be interesting inbetween
996 // removeBoundary and addPatches.
997 globalMeshDataPtr_.clear();
998
999 if (validBoundary)
1000 {
1001 // Calculate topology for the patches (processor-processor comms etc.)
1002 boundary_.updateMesh();
1003
1004 // Calculate the geometry for the patches (transformation tensors etc.)
1005 boundary_.calcGeometry();
1006
1007 boundary_.checkDefinition();
1008 }
1009}
1010
1011
1013(
1014 const List<pointZone*>& pz,
1015 const List<faceZone*>& fz,
1016 const List<cellZone*>& cz
1017)
1018{
1019 if (pointZones().size() || faceZones().size() || cellZones().size())
1020 {
1022 << "point, face or cell zone already exists"
1023 << abort(FatalError);
1024 }
1025
1026 // Point zones
1027 if (pz.size())
1028 {
1029 pointZones_.setSize(pz.size());
1030
1031 // Copy the zone pointers
1032 forAll(pz, pI)
1033 {
1034 pointZones_.set(pI, pz[pI]);
1035 }
1036
1037 pointZones_.writeOpt(IOobject::AUTO_WRITE);
1038 }
1039
1040 // Face zones
1041 if (fz.size())
1042 {
1043 faceZones_.setSize(fz.size());
1044
1045 // Copy the zone pointers
1046 forAll(fz, fI)
1047 {
1048 faceZones_.set(fI, fz[fI]);
1049 }
1050
1051 faceZones_.writeOpt(IOobject::AUTO_WRITE);
1052 }
1053
1054 // Cell zones
1055 if (cz.size())
1056 {
1057 cellZones_.setSize(cz.size());
1058
1059 // Copy the zone pointers
1060 forAll(cz, cI)
1061 {
1062 cellZones_.set(cI, cz[cI]);
1063 }
1064
1065 cellZones_.writeOpt(IOobject::AUTO_WRITE);
1066 }
1067}
1068
1069
1071(
1072 const List<polyPatch*>& p,
1073 const bool validBoundary
1074)
1075{
1076 // Acquire ownership of the pointers
1077 polyPatchList plist(const_cast<List<polyPatch*>&>(p));
1078
1079 addPatches(plist, validBoundary);
1080}
1081
1082
1084{
1085 if (clearedPrimitives_)
1086 {
1088 << "points deallocated"
1089 << abort(FatalError);
1090 }
1091
1092 return points_;
1093}
1094
1095
1097{
1098 return io.upToDate(points_);
1099}
1100
1101
1103{
1104 io.eventNo() = points_.eventNo()+1;
1105}
1106
1107
1109{
1110 if (clearedPrimitives_)
1111 {
1113 << "faces deallocated"
1114 << abort(FatalError);
1115 }
1116
1117 return faces_;
1118}
1119
1120
1122{
1123 return owner_;
1124}
1125
1126
1128{
1129 return neighbour_;
1130}
1131
1132
1134{
1135 if (!moving_)
1136 {
1137 return points_;
1138 }
1139
1140 if (!oldPointsPtr_)
1141 {
1142 if (debug)
1143 {
1145 }
1146
1147 oldPointsPtr_.reset(new pointField(points_));
1148 curMotionTimeIndex_ = time().timeIndex();
1149 }
1150
1151 return *oldPointsPtr_;
1152}
1153
1154
1156{
1157 storeOldCellCentres_ = true;
1158
1159 if (!moving_)
1160 {
1161 return cellCentres();
1162 }
1163
1164 if (!oldCellCentresPtr_)
1165 {
1166 oldCellCentresPtr_.reset(new pointField(cellCentres()));
1167 }
1168
1169 return *oldCellCentresPtr_;
1170}
1171
1172
1174{
1176 << "Moving points for time " << time().value()
1177 << " index " << time().timeIndex() << endl;
1178
1179 if (newPoints.size() != points_.size())
1180 {
1182 << "Size of newPoints " << newPoints.size()
1183 << " does not correspond to current mesh points size "
1184 << points_.size()
1185 << exit(FatalError);
1186 }
1187
1188
1189 moving(true);
1190
1191 // Pick up old points
1192 if (curMotionTimeIndex_ != time().timeIndex())
1193 {
1194 if (debug)
1195 {
1196 Info<< "tmp<scalarField> polyMesh::movePoints(const pointField&) : "
1197 << " Storing current points for time " << time().value()
1198 << " index " << time().timeIndex() << endl;
1199 }
1200
1201 if (storeOldCellCentres_)
1202 {
1203 oldCellCentresPtr_.clear();
1204 oldCellCentresPtr_.reset(new pointField(cellCentres()));
1205 }
1206
1207 // Mesh motion in the new time step
1208 oldPointsPtr_.clear();
1209 oldPointsPtr_.reset(new pointField(points_));
1210 curMotionTimeIndex_ = time().timeIndex();
1211 }
1212
1213 points_ = newPoints;
1214
1215 bool moveError = false;
1216 if (debug)
1217 {
1218 // Check mesh motion
1219 if (checkMeshMotion(points_, true))
1220 {
1221 moveError = true;
1222
1224 << "Moving the mesh with given points will "
1225 << "invalidate the mesh." << nl
1226 << "Mesh motion should not be executed." << endl;
1227 }
1228 }
1229
1230 points_.writeOpt(IOobject::AUTO_WRITE);
1231 points_.instance() = time().timeName();
1232 points_.eventNo() = getEvent();
1233
1234 if (tetBasePtIsPtr_)
1235 {
1236 tetBasePtIsPtr_->writeOpt(IOobject::AUTO_WRITE);
1237 tetBasePtIsPtr_->instance() = time().timeName();
1238 tetBasePtIsPtr_->eventNo() = getEvent();
1239 }
1240
1241 // Currently a no-op; earlier versions set meshPhi and call
1242 // primitiveMesh::clearGeom
1243 (void)primitiveMesh::movePoints(points_, oldPoints());
1244
1245 // Update the mesh geometry (via fvGeometryScheme)
1246 // - updateGeom is virtual -> calls fvMesh::updateGeom (or higher)
1247 // - fvMesh::updateGeom defers to surfaceInterpolation::updateGeom(),
1248 // which defers to fvGeometryScheme::movePoints()
1249 // - set the mesh flux
1250 // - clear out/recalculate stale geometry
1251 updateGeom();
1252
1253 // Adjust parallel shared points
1254 if (globalMeshDataPtr_)
1255 {
1256 globalMeshDataPtr_->movePoints(points_);
1257 }
1258
1259 // Force recalculation of all geometric data with new points
1260
1261 bounds_ = boundBox(points_);
1262 boundary_.movePoints(points_);
1263
1264 pointZones_.movePoints(points_);
1265 faceZones_.movePoints(points_);
1266 cellZones_.movePoints(points_);
1267
1268 // Reset cell tree - it gets built from mesh geometry so might have
1269 // wrong boxes. It is correct as long as none of the cells leaves
1270 // the boxes it is in which most likely is almost never the case except
1271 // for tiny displacements. An alternative is to check the displacements
1272 // to see if they are tiny - imagine a big windtunnel with a small rotating
1273 // object. In this case the processors without the rotating object wouldn't
1274 // have to clear any geometry. However your critical path still stays the
1275 // same so no time would be gained (unless the decomposition gets weighted).
1276 // Small benefit for lots of scope for problems so not done.
1277 cellTreePtr_.clear();
1278
1279 // Reset valid directions (could change with rotation)
1280 geometricD_ = Zero;
1281 solutionD_ = Zero;
1282
1283 // Note: tet-base decomposition does not get cleared. Ideally your face
1284 // decomposition should not change during mesh motion ...
1285
1286
1287 meshObject::movePoints<polyMesh>(*this);
1288 meshObject::movePoints<pointMesh>(*this);
1289
1290 const_cast<Time&>(time()).functionObjects().movePoints(*this);
1291
1292
1293 if (debug && moveError)
1294 {
1295 // Write mesh to ease debugging. Note we want to avoid calling
1296 // e.g. fvMesh::write since meshPhi not yet complete.
1298 }
1299}
1300
1301
1303{
1304 curMotionTimeIndex_ = 0;
1305 oldPointsPtr_.clear();
1306 oldCellCentresPtr_.clear();
1307}
1308
1309
1311{
1312 if (!globalMeshDataPtr_)
1313 {
1314 if (debug)
1315 {
1316 Pout<< "polyMesh::globalData() const : "
1317 << "Constructing parallelData from processor topology"
1318 << endl;
1319 }
1320 // Construct globalMeshData using processorPatch information only.
1321 globalMeshDataPtr_.reset(new globalMeshData(*this));
1322 }
1323
1324 return *globalMeshDataPtr_;
1325}
1326
1327
1329{
1330 return comm_;
1331}
1332
1333
1335{
1336 return comm_;
1337}
1338
1339
1340void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
1341{
1342 fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();
1343
1344 rm(meshFilesPath/"points");
1345 rm(meshFilesPath/"faces");
1346 rm(meshFilesPath/"owner");
1347 rm(meshFilesPath/"neighbour");
1348 rm(meshFilesPath/"cells");
1349 rm(meshFilesPath/"boundary");
1350 rm(meshFilesPath/"pointZones");
1351 rm(meshFilesPath/"faceZones");
1352 rm(meshFilesPath/"cellZones");
1353 rm(meshFilesPath/"meshModifiers");
1354 rm(meshFilesPath/"parallelData");
1355
1356 // remove subdirectories
1357 if (isDir(meshFilesPath/"sets"))
1358 {
1359 rmDir(meshFilesPath/"sets");
1360 }
1361}
1362
1363
1365{
1366 removeFiles(instance());
1367}
1368
1369
1371(
1372 const point& p,
1373 label& celli,
1374 label& tetFacei,
1375 label& tetPti
1376) const
1377{
1378 celli = -1;
1379 tetFacei = -1;
1380 tetPti = -1;
1381
1382 const indexedOctree<treeDataCell>& tree = cellTree();
1383
1384 // Find point inside cell
1385 celli = tree.findInside(p);
1386
1387 if (celli != -1)
1388 {
1389 // Check the nearest cell to see if the point is inside.
1390 findTetFacePt(celli, p, tetFacei, tetPti);
1391 }
1392}
1393
1394
1396(
1397 const label celli,
1398 const point& p,
1399 label& tetFacei,
1400 label& tetPti
1401) const
1402{
1403 const polyMesh& mesh = *this;
1404
1406 tetFacei = tet.face();
1407 tetPti = tet.tetPt();
1408}
1409
1410
1412(
1413 const point& p,
1414 label celli,
1415 const cellDecomposition decompMode
1416) const
1417{
1418 switch (decompMode)
1419 {
1420 case FACE_PLANES:
1421 {
1422 return primitiveMesh::pointInCell(p, celli);
1423 }
1424 break;
1425
1426 case FACE_CENTRE_TRIS:
1427 {
1428 // only test that point is on inside of plane defined by cell face
1429 // triangles
1430 const cell& cFaces = cells()[celli];
1431
1432 forAll(cFaces, cFacei)
1433 {
1434 label facei = cFaces[cFacei];
1435 const face& f = faces_[facei];
1436 const point& fc = faceCentres()[facei];
1437 bool isOwn = (owner_[facei] == celli);
1438
1439 forAll(f, fp)
1440 {
1441 label pointi;
1442 label nextPointi;
1443
1444 if (isOwn)
1445 {
1446 pointi = f[fp];
1447 nextPointi = f.nextLabel(fp);
1448 }
1449 else
1450 {
1451 pointi = f.nextLabel(fp);
1452 nextPointi = f[fp];
1453 }
1454
1455 triPointRef faceTri
1456 (
1457 points()[pointi],
1458 points()[nextPointi],
1459 fc
1460 );
1461
1462 vector proj = p - faceTri.centre();
1463
1464 if ((faceTri.areaNormal() & proj) > 0)
1465 {
1466 return false;
1467 }
1468 }
1469 }
1470 return true;
1471 }
1472 break;
1473
1474 case FACE_DIAG_TRIS:
1475 {
1476 // only test that point is on inside of plane defined by cell face
1477 // triangles
1478 const cell& cFaces = cells()[celli];
1479
1480 forAll(cFaces, cFacei)
1481 {
1482 label facei = cFaces[cFacei];
1483 const face& f = faces_[facei];
1484
1485 for (label tetPti = 1; tetPti < f.size() - 1; tetPti++)
1486 {
1487 // Get tetIndices of face triangle
1488 tetIndices faceTetIs(celli, facei, tetPti);
1489
1490 triPointRef faceTri = faceTetIs.faceTri(*this);
1491
1492 vector proj = p - faceTri.centre();
1493
1494 if ((faceTri.areaNormal() & proj) > 0)
1495 {
1496 return false;
1497 }
1498 }
1499 }
1500
1501 return true;
1502 }
1503 break;
1504
1505 case CELL_TETS:
1506 {
1507 label tetFacei;
1508 label tetPti;
1509
1510 findTetFacePt(celli, p, tetFacei, tetPti);
1511
1512 return tetFacei != -1;
1513 }
1514 break;
1515 }
1516
1517 return false;
1518}
1519
1520
1522(
1523 const point& p,
1524 const cellDecomposition decompMode
1525) const
1526{
1527 if
1528 (
1530 && (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
1531 )
1532 {
1533 // Force construction of face-diagonal decomposition before testing
1534 // for zero cells.
1535 //
1536 // If parallel running a local domain might have zero cells so never
1537 // construct the face-diagonal decomposition which uses parallel
1538 // transfers.
1539 (void)tetBasePtIs();
1540 }
1541
1542 if (nCells() == 0)
1543 {
1544 return -1;
1545 }
1546
1547 if (decompMode == CELL_TETS)
1548 {
1549 // Advanced search method utilizing an octree
1550 // and tet-decomposition of the cells
1551
1552 label celli;
1553 label tetFacei;
1554 label tetPti;
1555
1556 findCellFacePt(p, celli, tetFacei, tetPti);
1557
1558 return celli;
1559 }
1560 else
1561 {
1562 // Approximate search avoiding the construction of an octree
1563 // and cell decomposition
1564
1565 if (Pstream::parRun() && decompMode == FACE_DIAG_TRIS)
1566 {
1567 // Force construction of face-diagonal decomposition before testing
1568 // for zero cells. If parallel running a local domain might have
1569 // zero cells so never construct the face-diagonal decomposition
1570 // (which uses parallel transfers)
1571 (void)tetBasePtIs();
1572 }
1573
1574 // Find the nearest cell centre to this location
1575 label celli = findNearestCell(p);
1576
1577 // If point is in the nearest cell return
1578 if (pointInCell(p, celli, decompMode))
1579 {
1580 return celli;
1581 }
1582 else
1583 {
1584 // Point is not in the nearest cell so search all cells
1585
1586 for (label celli = 0; celli < nCells(); celli++)
1587 {
1588 if (pointInCell(p, celli, decompMode))
1589 {
1590 return celli;
1591 }
1592 }
1593
1594 return -1;
1595 }
1596 }
1597}
1598
1599
1600// ************************************************************************* //
A List of objects of type <T> with automated input and output using a compact storage....
Definition: CompactIOList.H:77
SubField< vector > subField
Declare type of subField.
Definition: Field.H:89
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
bool typeHeaderOk(const bool checkType=true, const bool search=true, const bool verbose=true)
Read header (uses typeFilePath to find file) and check its info.
bool hasHeaderClass() const noexcept
True if headerClassName() is non-empty (after reading)
Definition: IOobjectI.H:149
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
Random number generator.
Definition: Random.H:60
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Inter-processor communications stream.
Definition: UPstream.H:59
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
bool empty() const noexcept
True if the list is empty (ie, size() is zero)
Definition: UPtrListI.H:113
Templated 3D Vector derived from VectorSpace adding construction from 3 components,...
Definition: Vector.H:65
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:91
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface....
Definition: boundaryMesh.H:63
A subset of mesh cells.
Definition: cellZone.H:65
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:57
Particle-size distribution model wherein random samples are drawn from the doubly-truncated uniform p...
Definition: uniform.H:164
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:67
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
A class for handling file names.
Definition: fileName.H:76
static std::string path(const std::string &str)
Return directory path name (part before last /)
Definition: fileNameI.H:176
virtual bool write()
Write the output fields.
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
Non-pointer based hierarchical recursive searching.
Definition: indexedOctree.H:74
label findInside(const point &) const
Find shape containing point. Only implemented for certain.
void movePoints()
Update for new mesh geometry.
Registry of regIOobjects.
virtual const fileName & dbDir() const
Local directory path of this objectRegistry relative to the time.
const Time & time() const noexcept
Return time registry.
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
A subset of mesh points.
Definition: pointZone.H:68
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual ~polyMesh()
Destructor.
Definition: polyMesh.C:820
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
Definition: polyMesh.C:1371
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:866
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
Definition: polyMesh.C:883
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:321
void addPatches(polyPatchList &plist, const bool validBoundary=true)
Add boundary patches.
Definition: polyMesh.C:975
label comm() const noexcept
Return communicator used for parallel communication.
Definition: polyMesh.C:1328
fileName meshDir() const
Return the local mesh directory (dbDir()/meshSubDir)
Definition: polyMesh.C:848
cellDecomposition
Enumeration defining the decomposition of the cell for.
Definition: polyMesh.H:101
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1108
virtual const fileName & dbDir() const
Override the objectRegistry dbDir for a single-region case.
Definition: polyMesh.C:837
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
Definition: polyMesh.C:1522
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const labelIOList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:906
void findTetFacePt(const label celli, const point &p, label &tetFacei, label &tetPti) const
Find the tetFacei and tetPti for point p in celli.
Definition: polyMesh.C:1396
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1310
bool pointInCell(const point &p, label celli, const cellDecomposition=CELL_TETS) const
Test if point p is in the celli.
Definition: polyMesh.C:1412
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
Definition: polyMesh.C:900
virtual const pointField & oldPoints() const
Return old points (mesh motion)
Definition: polyMesh.C:1133
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:860
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
void resetMotion() const
Reset motion.
Definition: polyMesh.C:1302
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
virtual bool upToDatePoints(const regIOobject &io) const
Return true if io is up-to-date with points.
Definition: polyMesh.C:1096
void addZones(const List< pointZone * > &pz, const List< faceZone * > &fz, const List< cellZone * > &cz)
Add mesh zones.
Definition: polyMesh.C:1013
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
virtual void setUpToDatePoints(regIOobject &io) const
Set io to be up-to-date with points.
Definition: polyMesh.C:1102
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:940
void removeFiles() const
Remove all files from mesh instance()
Definition: polyMesh.C:1364
virtual const pointField & oldCellCentres() const
Return old cellCentres (mesh motion)
Definition: polyMesh.C:1155
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:324
void resetPrimitives(autoPtr< pointField > &&points, autoPtr< faceList > &&faces, autoPtr< labelList > &&owner, autoPtr< labelList > &&neighbour, const labelUList &patchSizes, const labelUList &patchStarts, const bool validBoundary=true)
Reset mesh primitive data. Assumes all patch info correct.
Definition: polyMesh.C:718
const word & regionName() const
The mesh region name or word::null if polyMesh::defaultRegion.
Definition: polyMesh.C:854
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:889
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:872
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:79
bool pointInCell(const point &p, label celli) const
Return true if the point is in the cell.
label nPoints() const noexcept
Number of mesh points.
label nCells() const noexcept
Number of mesh cells.
const cellList & cells() const
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
Definition: regIOobject.H:76
virtual bool write(const bool valid=true) const
Write using setting from DB.
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:84
label face() const noexcept
Return the face index.
Definition: tetIndices.H:139
triPointRef faceTri(const polyMesh &mesh) const
Definition: tetIndicesI.H:149
label tetPt() const noexcept
Return the characterising tet point index.
Definition: tetIndices.H:145
Standard boundBox with extra functionality for use in octree.
Definition: treeBoundBox.H:89
treeBoundBox extend(Random &rndGen, const scalar s) const
Return slightly wider bounding box.
Encapsulation of data needed to search in/for cells. Used to find the cell containing a point (e....
Definition: treeDataCell.H:57
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
vector areaNormal() const
The area normal - with magnitude equal to area of triangle.
Definition: triangleI.H:112
Point centre() const
Return centre (centroid)
Definition: triangleI.H:105
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
A class representing the concept of 0 (zero) that can be used to avoid manipulating objects known to ...
Definition: zero.H:63
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
volScalarField & p
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
label nPoints
const cellShapeList & cells
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
word timeName
Definition: getTimeIndex.H:3
#define WarningInFunction
Report a warning using Foam::Warning.
#define DebugInFunction
Report an information message using Foam::Info.
#define InfoInFunction
Report an information message using Foam::Info.
Namespace for OpenFOAM.
bool rm(const fileName &file)
Remove a file (or its gz equivalent), returning true if successful.
Definition: MSwindows.C:1012
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
bool rmDir(const fileName &directory, const bool silent=false)
Remove a directory and its contents (optionally silencing warnings)
Definition: MSwindows.C:1036
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
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
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:44
uint8_t direction
Definition: direction.H:56
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
const direction noexcept
Definition: Scalar.H:223
error FatalError
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
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.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
bool isDir(const fileName &name, const bool followLink=true)
Does the name exist as a DIRECTORY in the file system?
Definition: MSwindows.C:651
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
label timeIndex
Definition: getTimeIndex.H:30
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Random rndGen
Definition: createFields.H:23