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 | Copyright (C) 2011-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "polyMesh.H"
27 #include "Time.H"
28 #include "cellIOList.H"
29 #include "wedgePolyPatch.H"
30 #include "emptyPolyPatch.H"
31 #include "globalMeshData.H"
32 #include "processorPolyPatch.H"
34 #include "indexedOctree.H"
35 #include "treeDataCell.H"
36 #include "MeshObject.H"
37 #include "pointMesh.H"
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43  defineTypeNameAndDebug(polyMesh, 0);
44 
45  word polyMesh::defaultRegion = "region0";
46  word polyMesh::meshSubDir = "polyMesh";
47 }
48 
49 
50 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51 
52 void Foam::polyMesh::calcDirections() const
53 {
54  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
55  {
56  solutionD_[cmpt] = 1;
57  }
58 
59  // Knock out empty and wedge directions. Note:they will be present on all
60  // domains.
61 
62  label nEmptyPatches = 0;
63  label nWedgePatches = 0;
64 
65  vector emptyDirVec = Zero;
66  vector wedgeDirVec = Zero;
67 
68  forAll(boundaryMesh(), patchi)
69  {
70  if (boundaryMesh()[patchi].size())
71  {
72  if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
73  {
74  nEmptyPatches++;
75  emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
76  }
77  else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
78  {
79  const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
80  (
81  boundaryMesh()[patchi]
82  );
83 
84  nWedgePatches++;
85  wedgeDirVec += cmptMag(wpp.centreNormal());
86  }
87  }
88  }
89 
90  reduce(nEmptyPatches, maxOp<label>());
91  reduce(nWedgePatches, maxOp<label>());
92 
93  if (nEmptyPatches)
94  {
95  reduce(emptyDirVec, sumOp<vector>());
96 
97  emptyDirVec /= mag(emptyDirVec);
98 
99  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
100  {
101  if (emptyDirVec[cmpt] > 1e-6)
102  {
103  solutionD_[cmpt] = -1;
104  }
105  else
106  {
107  solutionD_[cmpt] = 1;
108  }
109  }
110  }
111 
112 
113  // Knock out wedge directions
114 
115  geometricD_ = solutionD_;
116 
117  if (nWedgePatches)
118  {
119  reduce(wedgeDirVec, sumOp<vector>());
120 
121  wedgeDirVec /= mag(wedgeDirVec);
122 
123  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
124  {
125  if (wedgeDirVec[cmpt] > 1e-6)
126  {
127  geometricD_[cmpt] = -1;
128  }
129  else
130  {
131  geometricD_[cmpt] = 1;
132  }
133  }
134  }
135 }
136 
137 
138 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
139 
140 Foam::polyMesh::polyMesh(const IOobject& io)
141 :
142  objectRegistry(io),
143  primitiveMesh(),
144  points_
145  (
146  IOobject
147  (
148  "points",
149  time().findInstance(meshDir(), "points"),
150  meshSubDir,
151  *this,
152  IOobject::MUST_READ,
153  IOobject::NO_WRITE
154  )
155  ),
156  faces_
157  (
158  IOobject
159  (
160  "faces",
161  time().findInstance(meshDir(), "faces"),
162  meshSubDir,
163  *this,
164  IOobject::MUST_READ,
165  IOobject::NO_WRITE
166  )
167  ),
168  owner_
169  (
170  IOobject
171  (
172  "owner",
173  faces_.instance(),
174  meshSubDir,
175  *this,
176  IOobject::READ_IF_PRESENT,
177  IOobject::NO_WRITE
178  )
179  ),
180  neighbour_
181  (
182  IOobject
183  (
184  "neighbour",
185  faces_.instance(),
186  meshSubDir,
187  *this,
188  IOobject::READ_IF_PRESENT,
189  IOobject::NO_WRITE
190  )
191  ),
192  clearedPrimitives_(false),
193  boundary_
194  (
195  IOobject
196  (
197  "boundary",
198  time().findInstance(meshDir(), "boundary"),
199  meshSubDir,
200  *this,
201  IOobject::MUST_READ,
202  IOobject::NO_WRITE
203  ),
204  *this
205  ),
206  bounds_(points_),
207  comm_(UPstream::worldComm),
208  geometricD_(Zero),
209  solutionD_(Zero),
210  tetBasePtIsPtr_(NULL),
211  cellTreePtr_(NULL),
212  pointZones_
213  (
214  IOobject
215  (
216  "pointZones",
217  time().findInstance
218  (
219  meshDir(),
220  "pointZones",
221  IOobject::READ_IF_PRESENT
222  ),
223  meshSubDir,
224  *this,
225  IOobject::READ_IF_PRESENT,
226  IOobject::NO_WRITE
227  ),
228  *this
229  ),
230  faceZones_
231  (
232  IOobject
233  (
234  "faceZones",
235  time().findInstance
236  (
237  meshDir(),
238  "faceZones",
239  IOobject::READ_IF_PRESENT
240  ),
241  meshSubDir,
242  *this,
243  IOobject::READ_IF_PRESENT,
244  IOobject::NO_WRITE
245  ),
246  *this
247  ),
248  cellZones_
249  (
250  IOobject
251  (
252  "cellZones",
253  time().findInstance
254  (
255  meshDir(),
256  "cellZones",
257  IOobject::READ_IF_PRESENT
258  ),
259  meshSubDir,
260  *this,
261  IOobject::READ_IF_PRESENT,
262  IOobject::NO_WRITE
263  ),
264  *this
265  ),
266  globalMeshDataPtr_(NULL),
267  moving_(false),
268  topoChanging_(false),
269  curMotionTimeIndex_(time().timeIndex()),
270  oldPointsPtr_(NULL)
271 {
272  if (exists(owner_.objectPath()))
273  {
274  initMesh();
275  }
276  else
277  {
278  cellCompactIOList cLst
279  (
280  IOobject
281  (
282  "cells",
283  time().findInstance(meshDir(), "cells"),
284  meshSubDir,
285  *this,
288  )
289  );
290 
291  // Set the primitive mesh
292  initMesh(cLst);
293 
294  owner_.write();
295  neighbour_.write();
296  }
297 
298  // Calculate topology for the patches (processor-processor comms etc.)
299  boundary_.updateMesh();
300 
301  // Calculate the geometry for the patches (transformation tensors etc.)
302  boundary_.calcGeometry();
303 
304  // Warn if global empty mesh
305  if (returnReduce(nPoints(), sumOp<label>()) == 0)
306  {
308  << "no points in mesh" << endl;
309  }
310  if (returnReduce(nCells(), sumOp<label>()) == 0)
311  {
313  << "no cells in mesh" << endl;
314  }
315 
316  // Initialise demand-driven data
317  calcDirections();
318 }
319 
320 
321 Foam::polyMesh::polyMesh
322 (
323  const IOobject& io,
324  const Xfer<pointField>& points,
325  const Xfer<faceList>& faces,
326  const Xfer<labelList>& owner,
327  const Xfer<labelList>& neighbour,
328  const bool syncPar
329 )
330 :
331  objectRegistry(io),
332  primitiveMesh(),
333  points_
334  (
335  IOobject
336  (
337  "points",
338  instance(),
339  meshSubDir,
340  *this,
341  io.readOpt(),
343  ),
344  points
345  ),
346  faces_
347  (
348  IOobject
349  (
350  "faces",
351  instance(),
352  meshSubDir,
353  *this,
354  io.readOpt(),
356  ),
357  faces
358  ),
359  owner_
360  (
361  IOobject
362  (
363  "owner",
364  instance(),
365  meshSubDir,
366  *this,
367  io.readOpt(),
369  ),
370  owner
371  ),
372  neighbour_
373  (
374  IOobject
375  (
376  "neighbour",
377  instance(),
378  meshSubDir,
379  *this,
380  io.readOpt(),
382  ),
383  neighbour
384  ),
385  clearedPrimitives_(false),
386  boundary_
387  (
388  IOobject
389  (
390  "boundary",
391  instance(),
392  meshSubDir,
393  *this,
394  io.readOpt(),
396  ),
397  *this,
398  polyPatchList()
399  ),
400  bounds_(points_, syncPar),
401  comm_(UPstream::worldComm),
402  geometricD_(Zero),
403  solutionD_(Zero),
404  tetBasePtIsPtr_(NULL),
405  cellTreePtr_(NULL),
406  pointZones_
407  (
408  IOobject
409  (
410  "pointZones",
411  instance(),
412  meshSubDir,
413  *this,
414  io.readOpt(),
416  ),
417  *this,
419  ),
420  faceZones_
421  (
422  IOobject
423  (
424  "faceZones",
425  instance(),
426  meshSubDir,
427  *this,
428  io.readOpt(),
430  ),
431  *this,
433  ),
434  cellZones_
435  (
436  IOobject
437  (
438  "cellZones",
439  instance(),
440  meshSubDir,
441  *this,
442  io.readOpt(),
444  ),
445  *this,
447  ),
448  globalMeshDataPtr_(NULL),
449  moving_(false),
450  topoChanging_(false),
451  curMotionTimeIndex_(time().timeIndex()),
452  oldPointsPtr_(NULL)
453 {
454  // Check if the faces and cells are valid
455  forAll(faces_, facei)
456  {
457  const face& curFace = faces_[facei];
458 
459  if (min(curFace) < 0 || max(curFace) > points_.size())
460  {
462  << "Face " << facei << "contains vertex labels out of range: "
463  << curFace << " Max point index = " << points_.size()
464  << abort(FatalError);
465  }
466  }
467 
468  // Set the primitive mesh
469  initMesh();
470 }
471 
472 
473 Foam::polyMesh::polyMesh
474 (
475  const IOobject& io,
476  const Xfer<pointField>& points,
477  const Xfer<faceList>& faces,
478  const Xfer<cellList>& cells,
479  const bool syncPar
480 )
481 :
482  objectRegistry(io),
483  primitiveMesh(),
484  points_
485  (
486  IOobject
487  (
488  "points",
489  instance(),
490  meshSubDir,
491  *this,
494  ),
495  points
496  ),
497  faces_
498  (
499  IOobject
500  (
501  "faces",
502  instance(),
503  meshSubDir,
504  *this,
507  ),
508  faces
509  ),
510  owner_
511  (
512  IOobject
513  (
514  "owner",
515  instance(),
516  meshSubDir,
517  *this,
520  ),
521  0
522  ),
523  neighbour_
524  (
525  IOobject
526  (
527  "neighbour",
528  instance(),
529  meshSubDir,
530  *this,
533  ),
534  0
535  ),
536  clearedPrimitives_(false),
537  boundary_
538  (
539  IOobject
540  (
541  "boundary",
542  instance(),
543  meshSubDir,
544  *this,
547  ),
548  *this,
549  0
550  ),
551  bounds_(points_, syncPar),
552  comm_(UPstream::worldComm),
553  geometricD_(Zero),
554  solutionD_(Zero),
555  tetBasePtIsPtr_(NULL),
556  cellTreePtr_(NULL),
557  pointZones_
558  (
559  IOobject
560  (
561  "pointZones",
562  instance(),
563  meshSubDir,
564  *this,
567  ),
568  *this,
569  0
570  ),
571  faceZones_
572  (
573  IOobject
574  (
575  "faceZones",
576  instance(),
577  meshSubDir,
578  *this,
581  ),
582  *this,
583  0
584  ),
585  cellZones_
586  (
587  IOobject
588  (
589  "cellZones",
590  instance(),
591  meshSubDir,
592  *this,
595  ),
596  *this,
597  0
598  ),
599  globalMeshDataPtr_(NULL),
600  moving_(false),
601  topoChanging_(false),
602  curMotionTimeIndex_(time().timeIndex()),
603  oldPointsPtr_(NULL)
604 {
605  // Check if faces are valid
606  forAll(faces_, facei)
607  {
608  const face& curFace = faces_[facei];
609 
610  if (min(curFace) < 0 || max(curFace) > points_.size())
611  {
613  << "Face " << facei << "contains vertex labels out of range: "
614  << curFace << " Max point index = " << points_.size()
615  << abort(FatalError);
616  }
617  }
618 
619  // transfer in cell list
620  cellList cLst(cells);
621 
622  // Check if cells are valid
623  forAll(cLst, celli)
624  {
625  const cell& curCell = cLst[celli];
626 
627  if (min(curCell) < 0 || max(curCell) > faces_.size())
628  {
630  << "Cell " << celli << "contains face labels out of range: "
631  << curCell << " Max face index = " << faces_.size()
632  << abort(FatalError);
633  }
634  }
635 
636  // Set the primitive mesh
637  initMesh(cLst);
638 }
639 
640 
642 (
643  const Xfer<pointField>& points,
644  const Xfer<faceList>& faces,
645  const Xfer<labelList>& owner,
646  const Xfer<labelList>& neighbour,
647  const labelList& patchSizes,
648  const labelList& patchStarts,
649  const bool validBoundary
650 )
651 {
652  // Clear addressing. Keep geometric props and updateable props for mapping.
653  clearAddressing(true);
654 
655  // Take over new primitive data.
656  // Optimized to avoid overwriting data at all
657  if (notNull(points))
658  {
659  points_.transfer(points());
660  bounds_ = boundBox(points_, validBoundary);
661  }
662 
663  if (notNull(faces))
664  {
665  faces_.transfer(faces());
666  }
667 
668  if (notNull(owner))
669  {
670  owner_.transfer(owner());
671  }
672 
673  if (notNull(neighbour))
674  {
675  neighbour_.transfer(neighbour());
676  }
677 
678 
679  // Reset patch sizes and starts
680  forAll(boundary_, patchi)
681  {
682  boundary_[patchi] = polyPatch
683  (
684  boundary_[patchi],
685  boundary_,
686  patchi,
687  patchSizes[patchi],
688  patchStarts[patchi]
689  );
690  }
691 
692 
693  // Flags the mesh files as being changed
695 
696  // Check if the faces and cells are valid
697  forAll(faces_, facei)
698  {
699  const face& curFace = faces_[facei];
700 
701  if (min(curFace) < 0 || max(curFace) > points_.size())
702  {
704  << "Face " << facei << " contains vertex labels out of range: "
705  << curFace << " Max point index = " << points_.size()
706  << abort(FatalError);
707  }
708  }
709 
710 
711  // Set the primitive mesh from the owner_, neighbour_.
712  // Works out from patch end where the active faces stop.
713  initMesh();
714 
715 
716  if (validBoundary)
717  {
718  // Note that we assume that all the patches stay the same and are
719  // correct etc. so we can already use the patches to do
720  // processor-processor comms.
721 
722  // Calculate topology for the patches (processor-processor comms etc.)
723  boundary_.updateMesh();
724 
725  // Calculate the geometry for the patches (transformation tensors etc.)
726  boundary_.calcGeometry();
727 
728  // Warn if global empty mesh
729  if
730  (
731  (returnReduce(nPoints(), sumOp<label>()) == 0)
732  || (returnReduce(nCells(), sumOp<label>()) == 0)
733  )
734  {
736  << "no points or no cells in mesh" << endl;
737  }
738  }
739 }
740 
741 
742 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
743 
745 {
746  clearOut();
747  resetMotion();
748 }
749 
750 
751 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
752 
754 {
756  {
757  return parent().dbDir();
758  }
759  else
760  {
761  return objectRegistry::dbDir();
762  }
763 }
764 
765 
767 {
768  return dbDir()/meshSubDir;
769 }
770 
771 
773 {
774  return points_.instance();
775 }
776 
777 
779 {
780  return faces_.instance();
781 }
782 
783 
785 {
786  if (geometricD_.x() == 0)
787  {
788  calcDirections();
789  }
790 
791  return geometricD_;
792 }
793 
794 
796 {
797  return cmptSum(geometricD() + Vector<label>::one)/2;
798 }
799 
800 
802 {
803  if (solutionD_.x() == 0)
804  {
805  calcDirections();
806  }
807 
808  return solutionD_;
809 }
810 
811 
813 {
814  return cmptSum(solutionD() + Vector<label>::one)/2;
815 }
816 
817 
819 {
820  if (tetBasePtIsPtr_.empty())
821  {
822  if (debug)
823  {
825  << "Forcing storage of base points."
826  << endl;
827  }
828 
829  tetBasePtIsPtr_.reset
830  (
831  new labelList
832  (
834  )
835  );
836  }
837 
838  return tetBasePtIsPtr_();
839 }
840 
841 
844 {
845  if (cellTreePtr_.empty())
846  {
847  treeBoundBox overallBb(points());
848 
849  Random rndGen(261782);
850 
851  overallBb = overallBb.extend(rndGen, 1e-4);
852  overallBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
853  overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
854 
855  cellTreePtr_.reset
856  (
858  (
860  (
861  false, // not cache bb
862  *this,
863  CELL_TETS // use tet-decomposition for any inside test
864  ),
865  overallBb,
866  8, // maxLevel
867  10, // leafsize
868  5.0 // duplicity
869  )
870  );
871  }
872 
873  return cellTreePtr_();
874 }
875 
876 
878 (
879  const List<polyPatch*>& p,
880  const bool validBoundary
881 )
882 {
883  if (boundaryMesh().size())
884  {
886  << "boundary already exists"
887  << abort(FatalError);
888  }
889 
890  // Reset valid directions
891  geometricD_ = Zero;
892  solutionD_ = Zero;
893 
894  boundary_.setSize(p.size());
895 
896  // Copy the patch pointers
897  forAll(p, pI)
898  {
899  boundary_.set(pI, p[pI]);
900  }
901 
902  // parallelData depends on the processorPatch ordering so force
903  // recalculation. Problem: should really be done in removeBoundary but
904  // there is some info in parallelData which might be interesting inbetween
905  // removeBoundary and addPatches.
906  globalMeshDataPtr_.clear();
907 
908  if (validBoundary)
909  {
910  // Calculate topology for the patches (processor-processor comms etc.)
911  boundary_.updateMesh();
912 
913  // Calculate the geometry for the patches (transformation tensors etc.)
914  boundary_.calcGeometry();
915 
916  boundary_.checkDefinition();
917  }
918 }
919 
920 
922 (
923  const List<pointZone*>& pz,
924  const List<faceZone*>& fz,
925  const List<cellZone*>& cz
926 )
927 {
928  if (pointZones().size() || faceZones().size() || cellZones().size())
929  {
931  << "point, face or cell zone already exists"
932  << abort(FatalError);
933  }
934 
935  // Point zones
936  if (pz.size())
937  {
938  pointZones_.setSize(pz.size());
939 
940  // Copy the zone pointers
941  forAll(pz, pI)
942  {
943  pointZones_.set(pI, pz[pI]);
944  }
945 
946  pointZones_.writeOpt() = IOobject::AUTO_WRITE;
947  }
948 
949  // Face zones
950  if (fz.size())
951  {
952  faceZones_.setSize(fz.size());
953 
954  // Copy the zone pointers
955  forAll(fz, fI)
956  {
957  faceZones_.set(fI, fz[fI]);
958  }
959 
960  faceZones_.writeOpt() = IOobject::AUTO_WRITE;
961  }
962 
963  // Cell zones
964  if (cz.size())
965  {
966  cellZones_.setSize(cz.size());
967 
968  // Copy the zone pointers
969  forAll(cz, cI)
970  {
971  cellZones_.set(cI, cz[cI]);
972  }
973 
974  cellZones_.writeOpt() = IOobject::AUTO_WRITE;
975  }
976 }
977 
978 
980 {
981  if (clearedPrimitives_)
982  {
984  << "points deallocated"
985  << abort(FatalError);
986  }
987 
988  return points_;
989 }
990 
991 
993 {
994  return io.upToDate(points_);
995 }
996 
997 
999 {
1000  io.eventNo() = points_.eventNo()+1;
1001 }
1002 
1003 
1005 {
1006  if (clearedPrimitives_)
1007  {
1009  << "faces deallocated"
1010  << abort(FatalError);
1011  }
1012 
1013  return faces_;
1014 }
1015 
1016 
1018 {
1019  return owner_;
1020 }
1021 
1022 
1024 {
1025  return neighbour_;
1026 }
1027 
1028 
1030 {
1031  if (oldPointsPtr_.empty())
1032  {
1033  if (debug)
1034  {
1036  << endl;
1037  }
1038 
1039  oldPointsPtr_.reset(new pointField(points_));
1040  curMotionTimeIndex_ = time().timeIndex();
1041  }
1042 
1043  return oldPointsPtr_();
1044 }
1045 
1046 
1049  const pointField& newPoints
1050 )
1051 {
1052  if (debug)
1053  {
1055  << "Moving points for time " << time().value()
1056  << " index " << time().timeIndex() << endl;
1057  }
1058 
1059  moving(true);
1060 
1061  // Pick up old points
1062  if (curMotionTimeIndex_ != time().timeIndex())
1063  {
1064  // Mesh motion in the new time step
1065  oldPointsPtr_.clear();
1066  oldPointsPtr_.reset(new pointField(points_));
1067  curMotionTimeIndex_ = time().timeIndex();
1068  }
1069 
1070  points_ = newPoints;
1071 
1072  bool moveError = false;
1073  if (debug)
1074  {
1075  // Check mesh motion
1076  if (checkMeshMotion(points_, true))
1077  {
1078  moveError = true;
1079 
1081  << "Moving the mesh with given points will "
1082  << "invalidate the mesh." << nl
1083  << "Mesh motion should not be executed." << endl;
1084  }
1085  }
1086 
1087  points_.writeOpt() = IOobject::AUTO_WRITE;
1088  points_.instance() = time().timeName();
1089  points_.eventNo() = getEvent();
1090 
1092  (
1093  points_,
1094  oldPoints()
1095  );
1096 
1097  // Adjust parallel shared points
1098  if (globalMeshDataPtr_.valid())
1099  {
1100  globalMeshDataPtr_().movePoints(points_);
1101  }
1102 
1103  // Force recalculation of all geometric data with new points
1104 
1105  bounds_ = boundBox(points_);
1106  boundary_.movePoints(points_);
1107 
1108  pointZones_.movePoints(points_);
1109  faceZones_.movePoints(points_);
1110  cellZones_.movePoints(points_);
1111 
1112  // Reset valid directions (could change with rotation)
1113  geometricD_ = Zero;
1114  solutionD_ = Zero;
1115 
1116  meshObject::movePoints<polyMesh>(*this);
1117  meshObject::movePoints<pointMesh>(*this);
1118 
1119  const_cast<Time&>(time()).functionObjects().movePoints(*this);
1120 
1121 
1122  if (debug && moveError)
1123  {
1124  // Write mesh to ease debugging. Note we want to avoid calling
1125  // e.g. fvMesh::write since meshPhi not yet complete.
1126  polyMesh::write();
1127  }
1128 
1129  return sweptVols;
1130 }
1131 
1132 
1134 {
1135  curMotionTimeIndex_ = 0;
1136  oldPointsPtr_.clear();
1137 }
1138 
1139 
1141 {
1142  if (globalMeshDataPtr_.empty())
1143  {
1144  if (debug)
1145  {
1146  Pout<< "polyMesh::globalData() const : "
1147  << "Constructing parallelData from processor topology"
1148  << endl;
1149  }
1150  // Construct globalMeshData using processorPatch information only.
1151  globalMeshDataPtr_.reset(new globalMeshData(*this));
1152  }
1153 
1154  return globalMeshDataPtr_();
1155 }
1156 
1157 
1159 {
1160  return comm_;
1161 }
1162 
1163 
1165 {
1166  return comm_;
1167 }
1168 
1169 
1170 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
1171 {
1172  fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir();
1173 
1174  rm(meshFilesPath/"points");
1175  rm(meshFilesPath/"faces");
1176  rm(meshFilesPath/"owner");
1177  rm(meshFilesPath/"neighbour");
1178  rm(meshFilesPath/"cells");
1179  rm(meshFilesPath/"boundary");
1180  rm(meshFilesPath/"pointZones");
1181  rm(meshFilesPath/"faceZones");
1182  rm(meshFilesPath/"cellZones");
1183  rm(meshFilesPath/"meshModifiers");
1184  rm(meshFilesPath/"parallelData");
1185 
1186  // remove subdirectories
1187  if (isDir(meshFilesPath/"sets"))
1188  {
1189  rmDir(meshFilesPath/"sets");
1190  }
1191 }
1192 
1193 
1195 {
1196  removeFiles(instance());
1197 }
1198 
1199 
1202  const point& p,
1203  label& celli,
1204  label& tetFacei,
1205  label& tetPti
1206 ) const
1207 {
1208  celli = -1;
1209  tetFacei = -1;
1210  tetPti = -1;
1211 
1212  const indexedOctree<treeDataCell>& tree = cellTree();
1213 
1214  // Find nearest cell to the point
1215  pointIndexHit info = tree.findNearest(p, sqr(GREAT));
1216 
1217  if (info.hit())
1218  {
1219  label nearestCelli = tree.shapes().cellLabels()[info.index()];
1220 
1221  // Check the nearest cell to see if the point is inside.
1222  findTetFacePt(nearestCelli, p, tetFacei, tetPti);
1223 
1224  if (tetFacei != -1)
1225  {
1226  // Point was in the nearest cell
1227 
1228  celli = nearestCelli;
1229 
1230  return;
1231  }
1232  else
1233  {
1234  // Check the other possible cells that the point may be in
1235 
1236  labelList testCells = tree.findIndices(p);
1237 
1238  forAll(testCells, pCI)
1239  {
1240  label testCelli = tree.shapes().cellLabels()[testCells[pCI]];
1241 
1242  if (testCelli == nearestCelli)
1243  {
1244  // Don't retest the nearest cell
1245 
1246  continue;
1247  }
1248 
1249  // Check the test cell to see if the point is inside.
1250  findTetFacePt(testCelli, p, tetFacei, tetPti);
1251 
1252  if (tetFacei != -1)
1253  {
1254  // Point was in the test cell
1255 
1256  celli = testCelli;
1257 
1258  return;
1259  }
1260  }
1261  }
1262  }
1263  else
1264  {
1266  << "Did not find nearest cell in search tree."
1267  << abort(FatalError);
1268  }
1269 }
1270 
1271 
1274  const label celli,
1275  const point& p,
1276  label& tetFacei,
1277  label& tetPti
1278 ) const
1279 {
1280  const polyMesh& mesh = *this;
1281 
1282  tetIndices tet(polyMeshTetDecomposition::findTet(mesh, celli, p));
1283  tetFacei = tet.face();
1284  tetPti = tet.tetPt();
1285 }
1286 
1287 
1290  const point& p,
1291  label celli,
1292  const cellDecomposition decompMode
1293 ) const
1294 {
1295  switch (decompMode)
1296  {
1297  case FACE_PLANES:
1298  {
1299  return primitiveMesh::pointInCell(p, celli);
1300  }
1301  break;
1302 
1303  case FACE_CENTRE_TRIS:
1304  {
1305  // only test that point is on inside of plane defined by cell face
1306  // triangles
1307  const cell& cFaces = cells()[celli];
1308 
1309  forAll(cFaces, cFacei)
1310  {
1311  label facei = cFaces[cFacei];
1312  const face& f = faces_[facei];
1313  const point& fc = faceCentres()[facei];
1314  bool isOwn = (owner_[facei] == celli);
1315 
1316  forAll(f, fp)
1317  {
1318  label pointi;
1319  label nextPointi;
1320 
1321  if (isOwn)
1322  {
1323  pointi = f[fp];
1324  nextPointi = f.nextLabel(fp);
1325  }
1326  else
1327  {
1328  pointi = f.nextLabel(fp);
1329  nextPointi = f[fp];
1330  }
1331 
1332  triPointRef faceTri
1333  (
1334  points()[pointi],
1335  points()[nextPointi],
1336  fc
1337  );
1338 
1339  vector proj = p - faceTri.centre();
1340 
1341  if ((faceTri.normal() & proj) > 0)
1342  {
1343  return false;
1344  }
1345  }
1346  }
1347  return true;
1348  }
1349  break;
1350 
1351  case FACE_DIAG_TRIS:
1352  {
1353  // only test that point is on inside of plane defined by cell face
1354  // triangles
1355  const cell& cFaces = cells()[celli];
1356 
1357  forAll(cFaces, cFacei)
1358  {
1359  label facei = cFaces[cFacei];
1360  const face& f = faces_[facei];
1361 
1362  for (label tetPti = 1; tetPti < f.size() - 1; tetPti++)
1363  {
1364  // Get tetIndices of face triangle
1365  tetIndices faceTetIs
1366  (
1368  (
1369  *this,
1370  facei,
1371  celli,
1372  tetPti
1373  )
1374  );
1375 
1376  triPointRef faceTri = faceTetIs.faceTri(*this);
1377 
1378  vector proj = p - faceTri.centre();
1379 
1380  if ((faceTri.normal() & proj) > 0)
1381  {
1382  return false;
1383  }
1384  }
1385  }
1386 
1387  return true;
1388  }
1389  break;
1390 
1391  case CELL_TETS:
1392  {
1393  label tetFacei;
1394  label tetPti;
1395 
1396  findTetFacePt(celli, p, tetFacei, tetPti);
1397 
1398  return tetFacei != -1;
1399  }
1400  break;
1401  }
1402 
1403  return false;
1404 }
1405 
1406 
1409  const point& p,
1410  const cellDecomposition decompMode
1411 ) const
1412 {
1413  if
1414  (
1415  Pstream::parRun()
1416  && (decompMode == FACE_DIAG_TRIS || decompMode == CELL_TETS)
1417  )
1418  {
1419  // Force construction of face-diagonal decomposition before testing
1420  // for zero cells.
1421  //
1422  // If parallel running a local domain might have zero cells so never
1423  // construct the face-diagonal decomposition which uses parallel
1424  // transfers.
1425  (void)tetBasePtIs();
1426  }
1427 
1428  if (nCells() == 0)
1429  {
1430  return -1;
1431  }
1432 
1433  if (decompMode == CELL_TETS)
1434  {
1435  // Advanced search method utilizing an octree
1436  // and tet-decomposition of the cells
1437 
1438  label celli;
1439  label tetFacei;
1440  label tetPti;
1441 
1442  findCellFacePt(p, celli, tetFacei, tetPti);
1443 
1444  return celli;
1445  }
1446  else
1447  {
1448  // Approximate search avoiding the construction of an octree
1449  // and cell decomposition
1450 
1451  // Find the nearest cell centre to this location
1452  label celli = findNearestCell(p);
1453 
1454  // If point is in the nearest cell return
1455  if (pointInCell(p, celli, decompMode))
1456  {
1457  return celli;
1458  }
1459  else
1460  {
1461  // Point is not in the nearest cell so search all cells
1462 
1463  for (label celli = 0; celli < nCells(); celli++)
1464  {
1465  if (pointInCell(p, celli, decompMode))
1466  {
1467  return celli;
1468  }
1469  }
1470 
1471  return -1;
1472  }
1473  }
1474 }
1475 
1476 
1477 // ************************************************************************* //
label timeIndex() const
Return current time index.
Definition: TimeStateI.H:35
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Definition: polyPatchList.H:45
cachedRandom rndGen(label(0),-1)
const Time & time() const
Return time.
A simple container for copying or transferring objects of type <T>.
Definition: Xfer.H:85
fileName objectPath() const
Return complete path + object name.
Definition: IOobject.H:363
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: polyMesh.C:1048
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
const indexedOctree< treeDataCell > & cellTree() const
Return the cell search tree.
Definition: polyMesh.C:843
void clearAddressing()
Clear topological data.
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:57
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
uint8_t direction
Definition: direction.H:46
label comm() const
Return communicator used for parallel communication.
Definition: polyMesh.C:1158
A class for handling file names.
Definition: fileName.H:69
static labelList findFaceBasePts(const polyMesh &mesh, scalar tol=minTetQuality, bool report=false)
Find a suitable base point for each face for decomposition.
void removeFiles() const
Remove all files from mesh instance()
Definition: polyMesh.C:1194
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Cmpt & x() const
Definition: VectorI.H:75
void resetMotion() const
Reset motion.
Definition: polyMesh.C:1133
const double e
Elementary charge.
Definition: doubleFloat.H:78
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
bool checkDefinition(const bool report=false) const
Check boundary definition. Return true if in error.
error FatalError
const fileName & facesInstance() const
Return the current instance directory for faces.
Definition: polyMesh.C:778
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:801
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:109
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
bool hit() const
Is there a hit.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
void clearOut()
Clear all geometry and addressing unnecessary for CFD.
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:309
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Cmpt cmptSum(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:509
const Vector< label > & geometricD() const
Return the vector of geometric directions in mesh.
Definition: polyMesh.C:784
treeBoundBox extend(Random &, const scalar s) const
Return slightly wider bounding box.
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:306
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
primitiveMesh()
Construct null.
Definition: primitiveMesh.C:39
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
const Type & value() const
Return const reference to value.
label nSolutionD() const
Return the number of valid solved-for dimensions in the mesh.
Definition: polyMesh.C:812
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
bool isDir(const fileName &)
Does the name exist as a DIRECTORY in the file system?
Definition: POSIX.C:486
label size() const
Return number of elements in table.
label eventNo() const
Event number at last update.
Definition: regIOobjectI.H:80
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:715
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const objectRegistry & parent() const
Return the parent objectRegistry.
static label worldComm
Default communicator (all processors)
Definition: UPstream.H:274
const vectorField & faceCentres() const
label tetPt() const
Return the characterising tetPtI.
Definition: tetIndicesI.H:60
bool upToDate(const regIOobject &) const
Return true if up-to-date with respect to given object.
Definition: regIOobject.C:320
const cellList & cells() const
writeOption writeOpt() const
Definition: IOobject.H:314
virtual bool checkMeshMotion(const pointField &newPoints, const bool report=false, const bool detailedReport=false) const
Check mesh motion for correctness given motion points.
label nGeometricD() const
Return the number of valid geometric dimensions in the mesh.
Definition: polyMesh.C:795
cellDecomposition
Enumeration defining the decomposition of the cell for.
Definition: polyMesh.H:98
void movePoints(const pointField &)
Correct zone mesh after moving points.
Definition: ZoneMesh.C:505
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
fileName meshDir() const
Return the local mesh directory (dbDir()/meshSubDir)
Definition: polyMesh.C:766
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
bool moving() const
Is mesh moving.
Definition: polyMesh.H:493
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
static tetIndices findTet(const polyMesh &mesh, label cI, const point &pt)
Find the tet decomposition of the cell containing the given point.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:96
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
tmp< scalarField > movePoints(const pointField &p, const pointField &oldP)
Move points, returns volumes swept by faces in motion.
dynamicFvMesh & mesh
virtual void setUpToDatePoints(regIOobject &io) const
Set io to be up-to-date with points.
Definition: polyMesh.C:998
virtual const pointField & oldPoints() const
Return old points for mesh motion.
Definition: polyMesh.C:1029
void findCellFacePt(const point &p, label &celli, label &tetFacei, label &tetPti) const
Find the cell, tetFacei and tetPti for point p.
Definition: polyMesh.C:1201
void movePoints(const pointField &)
Correct polyBoundaryMesh after moving points.
void resetPrimitives(const Xfer< pointField > &points, const Xfer< faceList > &faces, const Xfer< labelList > &owner, const Xfer< labelList > &neighbour, const labelList &patchSizes, const labelList &patchStarts, const bool validBoundary=true)
Reset mesh primitive data. Assumes all patch info correct.
Definition: polyMesh.C:642
virtual bool upToDatePoints(const regIOobject &io) const
Return true if io is up-to-date with points.
Definition: polyMesh.C:992
virtual const fileName & dbDir() const
Override the objectRegistry dbDir for a single-region case.
Definition: polyMesh.C:753
InfoProxy< IOobject > info() const
Return info proxy.
Definition: IOobject.H:421
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
const objectRegistry & thisDb() const
Return the object registry.
Definition: polyMesh.H:484
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:91
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:1273
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
Point centre() const
Return centre (centroid)
Definition: triangleI.H:89
label face() const
Return the face.
Definition: tetIndicesI.H:36
errorManip< error > abort(error &err)
Definition: errorManip.H:131
Simple random number generator.
Definition: Random.H:49
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1140
Templated 3D Vector derived from VectorSpace adding construction from 3 components, element access using x(), y() and z() member functions and the inner-product (dot-product) and cross product operators.
Definition: Vector.H:57
bool rmDir(const fileName &)
Remove a dirctory and its contents.
Definition: POSIX.C:839
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
static tetIndices triangleTetIndices(const polyMesh &mesh, label fI, label cI, const label tetPtI)
Return the tet decomposition of the given triangle of the given face.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool pointInCell(const point &p, label celli) const
Return true if the point is in the cell.
static const char nl
Definition: Ostream.H:262
label getEvent() const
Return new event number.
void updateMesh()
Correct polyBoundaryMesh after topology update.
defineTypeNameAndDebug(combustionModel, 0)
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:818
labelList f(nPoints)
A List of objects of type <T> with automated input and output using a compact storage. Behaves like IOList except when binary output in case it writes a CompactListList.
Definition: CompactIOList.H:53
Encapsulation of data needed to search in/for cells. Used to find the cell containing a point (e...
Definition: treeDataCell.H:54
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const pointZoneMesh & pointZones() const
Return point zone mesh.
Definition: polyMesh.H:457
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
bool notNull(const T &t)
Return true if t is not a reference to the nullObject of type T.
Definition: nullObjectI.H:46
void addPatches(const List< polyPatch * > &, const bool validBoundary=true)
Add boundary patches.
Definition: polyMesh.C:878
bool pointInCell(const point &p, label celli, const cellDecomposition=CELL_TETS) const
Test if point p is in the celli.
Definition: polyMesh.C:1289
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:393
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
bool exists(const fileName &, const bool checkGzip=true)
Does the name exist (as DIRECTORY or FILE) in the file system?
Definition: POSIX.C:480
label patchi
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
vector point
Point is a vector.
Definition: point.H:41
label findNearestCell(const point &location) const
Find the cell with the nearest cell centre to location.
#define WarningInFunction
Report a warning using Foam::Warning.
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:772
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
const Type & shapes() const
Reference to shape.
label index() const
Return index.
label timeIndex
Definition: getTimeIndex.H:4
readOption readOpt() const
Definition: IOobject.H:304
virtual const fileName & dbDir() const
Local directory path of this objectRegistry relative to the time.
virtual bool write() const
Write using setting from DB.
label nPoints() const
virtual ~polyMesh()
Destructor.
Definition: polyMesh.C:744
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
Definition: polyMesh.C:1408
regIOobject is an abstract class derived from IOobject to handle automatic object registration with t...
Definition: regIOobject.H:60
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
void addZones(const List< pointZone * > &pz, const List< faceZone * > &fz, const List< cellZone * > &cz)
Add mesh zones.
Definition: polyMesh.C:922
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
const labelList & findIndices(const point &) const
Find the shape indices that occupy the result of findNode.
A class for managing temporary objects.
Definition: PtrList.H:54
bool rm(const fileName &)
Remove a file, returning true if successful otherwise false.
Definition: POSIX.C:819
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
Registry of regIOobjects.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
fileName path() const
Return path.
Definition: Time.H:269
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const fileName & instance() const
Definition: IOobject.H:337
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:365
Inter-processor communications stream.
Definition: UPstream.H:58
Namespace for OpenFOAM.
vector normal() const
Return vector normal.
Definition: triangleI.H:103
#define InfoInFunction
Report an information message using Foam::Info.