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