checkGeometry.C
Go to the documentation of this file.
1 #include "PatchTools.H"
2 #include "checkGeometry.H"
3 #include "polyMesh.H"
4 #include "cellSet.H"
5 #include "faceSet.H"
6 #include "pointSet.H"
7 #include "edgeHashes.H"
8 #include "wedgePolyPatch.H"
9 #include "unitConversion.H"
10 #include "polyMeshTetDecomposition.H"
11 #include "checkTools.H"
12 #include "functionObject.H"
13 
14 #include "vtkSurfaceWriter.H"
15 #include "writer.H"
16 
17 #include "cyclicACMIPolyPatch.H"
18 #include "Time.H"
19 
20 // Find wedge with opposite orientation. Note: does not actually check that
21 // it is opposite, only that it has opposite normal and same axis
22 Foam::label Foam::findOppositeWedge
23 (
24  const polyMesh& mesh,
25  const wedgePolyPatch& wpp
26 )
27 {
28  const polyBoundaryMesh& patches = mesh.boundaryMesh();
29 
30  scalar wppCosAngle = wpp.cosAngle();
31 
32  forAll(patches, patchi)
33  {
34  if
35  (
36  patchi != wpp.index()
37  && patches[patchi].size()
38  && isA<wedgePolyPatch>(patches[patchi])
39  )
40  {
41  const wedgePolyPatch& pp =
42  refCast<const wedgePolyPatch>(patches[patchi]);
43 
44  // Calculate (cos of) angle to wpp (not pp!) centre normal
45  scalar ppCosAngle = wpp.centreNormal() & pp.n();
46 
47  if
48  (
49  pp.size() == wpp.size()
50  && mag(pp.axis() & wpp.axis()) >= (1-1e-3)
51  && mag(ppCosAngle - wppCosAngle) >= 1e-3
52  )
53  {
54  return patchi;
55  }
56  }
57  }
58  return -1;
59 }
60 
61 
62 bool Foam::checkWedges
63 (
64  const polyMesh& mesh,
65  const bool report,
66  const Vector<label>& directions,
67  labelHashSet* setPtr
68 )
69 {
70  // To mark edges without calculating edge addressing
71  EdgeMap<label> edgesInError;
72 
73  const pointField& p = mesh.points();
74  const faceList& fcs = mesh.faces();
75 
76 
77  const polyBoundaryMesh& patches = mesh.boundaryMesh();
78  forAll(patches, patchi)
79  {
80  if (patches[patchi].size() && isA<wedgePolyPatch>(patches[patchi]))
81  {
82  const wedgePolyPatch& pp =
83  refCast<const wedgePolyPatch>(patches[patchi]);
84 
85  scalar wedgeAngle = acos(pp.cosAngle());
86 
87  if (report)
88  {
89  Info<< " Wedge " << pp.name() << " with angle "
90  << radToDeg(wedgeAngle) << " degrees"
91  << endl;
92  }
93 
94  // Find opposite
95  label oppositePatchi = findOppositeWedge(mesh, pp);
96 
97  if (oppositePatchi == -1)
98  {
99  if (report)
100  {
101  Info<< " ***Cannot find opposite wedge for wedge "
102  << pp.name() << endl;
103  }
104  return true;
105  }
106 
107  const wedgePolyPatch& opp =
108  refCast<const wedgePolyPatch>(patches[oppositePatchi]);
109 
110 
111  if (mag(opp.axis() & pp.axis()) < (1-1e-3))
112  {
113  if (report)
114  {
115  Info<< " ***Wedges do not have the same axis."
116  << " Encountered " << pp.axis()
117  << " on patch " << pp.name()
118  << " which differs from " << opp.axis()
119  << " on opposite wedge patch" << opp.axis()
120  << endl;
121  }
122  return true;
123  }
124 
125 
126 
127  // Mark edges on wedgePatches
128  forAll(pp, i)
129  {
130  const face& f = pp[i];
131  forAll(f, fp)
132  {
133  label p0 = f[fp];
134  label p1 = f.nextLabel(fp);
135  edgesInError.insert(edge(p0, p1), -1); // non-error value
136  }
137  }
138 
139 
140  // Check that wedge patch is flat
141  const point& p0 = p[pp.meshPoints()[0]];
142  forAll(pp.meshPoints(), i)
143  {
144  const point& pt = p[pp.meshPoints()[i]];
145  scalar d = mag((pt - p0) & pp.n());
146 
147  if (d > ROOTSMALL)
148  {
149  if (report)
150  {
151  Info<< " ***Wedge patch " << pp.name() << " not planar."
152  << " Point " << pt << " is not in patch plane by "
153  << d << " metre."
154  << endl;
155  }
156  return true;
157  }
158  }
159  }
160  }
161 
162 
163 
164  // Check all non-wedge faces
165  label nEdgesInError = 0;
166 
167  forAll(fcs, facei)
168  {
169  const face& f = fcs[facei];
170 
171  forAll(f, fp)
172  {
173  label p0 = f[fp];
174  label p1 = f.nextLabel(fp);
175  if (p0 < p1)
176  {
177  vector d(p[p1]-p[p0]);
178  scalar magD = mag(d);
179 
180  if (magD > ROOTVSMALL)
181  {
182  d /= magD;
183 
184  // Check how many empty directions are used by the edge.
185  label nEmptyDirs = 0;
186  label nNonEmptyDirs = 0;
187  for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
188  {
189  if (mag(d[cmpt]) > 1e-6)
190  {
191  if (directions[cmpt] == 0)
192  {
193  nEmptyDirs++;
194  }
195  else
196  {
197  nNonEmptyDirs++;
198  }
199  }
200  }
201 
202  if (nEmptyDirs == 0)
203  {
204  // Purely in ok directions.
205  }
206  else if (nEmptyDirs == 1)
207  {
208  // Ok if purely in empty directions.
209  if (nNonEmptyDirs > 0)
210  {
211  if (edgesInError.insert(edge(p0, p1), facei))
212  {
213  nEdgesInError++;
214  }
215  }
216  }
217  else if (nEmptyDirs > 1)
218  {
219  // Always an error
220  if (edgesInError.insert(edge(p0, p1), facei))
221  {
222  nEdgesInError++;
223  }
224  }
225  }
226  }
227  }
228  }
229 
230  label nErrorEdges = returnReduce(nEdgesInError, sumOp<label>());
231 
232  if (nErrorEdges > 0)
233  {
234  if (report)
235  {
236  Info<< " ***Number of edges not aligned with or perpendicular to "
237  << "non-empty directions: " << nErrorEdges << endl;
238  }
239 
240  if (setPtr)
241  {
242  setPtr->resize(2*nEdgesInError);
243  forAllConstIters(edgesInError, iter)
244  {
245  if (iter() >= 0)
246  {
247  setPtr->insert(iter.key()[0]);
248  setPtr->insert(iter.key()[1]);
249  }
250  }
251  }
252 
253  return true;
254  }
255 
256  if (report)
257  {
258  Info<< " All edges aligned with or perpendicular to "
259  << "non-empty directions." << endl;
260  }
261 
262  return false;
263 }
264 
265 
266 namespace Foam
267 {
268  //- Default transformation behaviour for position
269  class transformPositionList
270  {
271  public:
272 
273  //- Transform patch-based field
274  void operator()
275  (
276  const coupledPolyPatch& cpp,
277  List<pointField>& pts
278  ) const
279  {
280  // Each element of pts is all the points in the face. Convert into
281  // lists of size cpp to transform.
282 
283  List<pointField> newPts(pts.size());
284  forAll(pts, facei)
285  {
286  newPts[facei].setSize(pts[facei].size());
287  }
288 
289  label index = 0;
290  while (true)
291  {
292  label n = 0;
293 
294  // Extract for every face the i'th position
295  pointField ptsAtIndex(pts.size(), Zero);
296  forAll(cpp, facei)
297  {
298  const pointField& facePts = pts[facei];
299  if (facePts.size() > index)
300  {
301  ptsAtIndex[facei] = facePts[index];
302  n++;
303  }
304  }
305 
306  if (n == 0)
307  {
308  break;
309  }
310 
311  // Now ptsAtIndex will have for every face either zero or
312  // the position of the i'th vertex. Transform.
313  cpp.transformPosition(ptsAtIndex);
314 
315  // Extract back from ptsAtIndex into newPts
316  forAll(cpp, facei)
317  {
318  pointField& facePts = newPts[facei];
319  if (facePts.size() > index)
320  {
321  facePts[index] = ptsAtIndex[facei];
322  }
323  }
324 
325  index++;
326  }
327 
328  pts.transfer(newPts);
329  }
330  };
331 }
332 
333 
334 bool Foam::checkCoupledPoints
335 (
336  const polyMesh& mesh,
337  const bool report,
338  labelHashSet* setPtr
339 )
340 {
341  const pointField& p = mesh.points();
342  const faceList& fcs = mesh.faces();
343  const polyBoundaryMesh& patches = mesh.boundaryMesh();
344 
345  // Zero'th point on coupled faces
346  //pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
347  List<pointField> nbrPoints(fcs.size() - mesh.nInternalFaces());
348 
349  // Exchange zero point
350  forAll(patches, patchi)
351  {
352  if (patches[patchi].coupled())
353  {
354  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
355  (
356  patches[patchi]
357  );
358 
359  forAll(cpp, i)
360  {
361  label bFacei = cpp.start() + i - mesh.nInternalFaces();
362  const face& f = cpp[i];
363  nbrPoints[bFacei].setSize(f.size());
364  forAll(f, fp)
365  {
366  const point& p0 = p[f[fp]];
367  nbrPoints[bFacei][fp] = p0;
368  }
369  }
370  }
371  }
372  syncTools::syncBoundaryFaceList
373  (
374  mesh,
375  nbrPoints,
376  eqOp<pointField>(),
377  transformPositionList()
378  );
379 
380  // Compare to local ones. Use same tolerance as for matching
381  label nErrorFaces = 0;
382  scalar avgMismatch = 0;
383  label nCoupledPoints = 0;
384 
385  forAll(patches, patchi)
386  {
387  if (patches[patchi].coupled())
388  {
389  const coupledPolyPatch& cpp =
390  refCast<const coupledPolyPatch>(patches[patchi]);
391 
392  if (cpp.owner())
393  {
394  scalarField smallDist
395  (
396  cpp.calcFaceTol
397  (
398  //cpp.matchTolerance(),
399  cpp,
400  cpp.points(),
401  cpp.faceCentres()
402  )
403  );
404 
405  forAll(cpp, i)
406  {
407  label bFacei = cpp.start() + i - mesh.nInternalFaces();
408  const face& f = cpp[i];
409 
410  if (f.size() != nbrPoints[bFacei].size())
411  {
412  FatalErrorInFunction
413  << "Local face size : " << f.size()
414  << " does not equal neighbour face size : "
415  << nbrPoints[bFacei].size()
416  << abort(FatalError);
417  }
418 
419  label fp = 0;
420  forAll(f, j)
421  {
422  const point& p0 = p[f[fp]];
423  scalar d = mag(p0 - nbrPoints[bFacei][j]);
424 
425  if (d > smallDist[i])
426  {
427  if (setPtr)
428  {
429  setPtr->insert(cpp.start()+i);
430  }
431  nErrorFaces++;
432 
433  break;
434  }
435 
436  avgMismatch += d;
437  nCoupledPoints++;
438 
439  fp = f.rcIndex(fp);
440  }
441  }
442  }
443  }
444  }
445 
446  reduce(nErrorFaces, sumOp<label>());
447  reduce(avgMismatch, maxOp<scalar>());
448  reduce(nCoupledPoints, sumOp<label>());
449 
450  if (nCoupledPoints > 0)
451  {
452  avgMismatch /= nCoupledPoints;
453  }
454 
455  if (nErrorFaces > 0)
456  {
457  if (report)
458  {
459  Info<< " **Error in coupled point location: "
460  << nErrorFaces
461  << " faces have their 0th or consecutive vertex not opposite"
462  << " their coupled equivalent. Average mismatch "
463  << avgMismatch << "."
464  << endl;
465  }
466 
467  return true;
468  }
469 
470  if (report)
471  {
472  Info<< " Coupled point location match (average "
473  << avgMismatch << ") OK." << endl;
474  }
475 
476  return false;
477 }
478 
479 
480 Foam::label Foam::checkGeometry
481 (
482  const polyMesh& mesh,
483  const bool allGeometry,
484  autoPtr<surfaceWriter>& surfWriter,
485  const autoPtr<writer<scalar>>& setWriter
486 )
487 {
488  label noFailedChecks = 0;
489 
490  Info<< "\nChecking geometry..." << endl;
491 
492  // Get a small relative length from the bounding box
493  const boundBox& globalBb = mesh.bounds();
494 
495  Info<< " Overall domain bounding box "
496  << globalBb.min() << " " << globalBb.max() << endl;
497 
498 
499  // Min length
500  scalar minDistSqr = magSqr(1e-6 * globalBb.span());
501 
502  // Geometric directions
503  const Vector<label> validDirs = (mesh.geometricD() + Vector<label>::one)/2;
504  Info<< " Mesh has " << mesh.nGeometricD()
505  << " geometric (non-empty/wedge) directions " << validDirs << endl;
506 
507  // Solution directions
508  const Vector<label> solDirs = (mesh.solutionD() + Vector<label>::one)/2;
509  Info<< " Mesh has " << mesh.nSolutionD()
510  << " solution (non-empty) directions " << solDirs << endl;
511 
512  if (mesh.nGeometricD() < 3)
513  {
514  pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
515 
516  if
517  (
518  (
519  validDirs != solDirs
520  && checkWedges(mesh, true, validDirs, &nonAlignedPoints)
521  )
522  || (
523  validDirs == solDirs
524  && mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints)
525  )
526  )
527  {
528  noFailedChecks++;
529  label nNonAligned = returnReduce
530  (
531  nonAlignedPoints.size(),
532  sumOp<label>()
533  );
534 
535  if (nNonAligned > 0)
536  {
537  Info<< " <<Writing " << nNonAligned
538  << " points on non-aligned edges to set "
539  << nonAlignedPoints.name() << endl;
540  nonAlignedPoints.instance() = mesh.pointsInstance();
541  nonAlignedPoints.write();
542  if (setWriter)
543  {
544  mergeAndWrite(*setWriter, nonAlignedPoints);
545  }
546  }
547  }
548  }
549 
550  if (mesh.checkClosedBoundary(true)) noFailedChecks++;
551 
552  {
553  cellSet cells(mesh, "nonClosedCells", mesh.nCells()/100+1);
554  cellSet aspectCells(mesh, "highAspectRatioCells", mesh.nCells()/100+1);
555  if
556  (
557  mesh.checkClosedCells
558  (
559  true,
560  &cells,
561  &aspectCells,
562  mesh.geometricD()
563  )
564  )
565  {
566  noFailedChecks++;
567 
568  label nNonClosed = returnReduce(cells.size(), sumOp<label>());
569 
570  if (nNonClosed > 0)
571  {
572  Info<< " <<Writing " << nNonClosed
573  << " non closed cells to set " << cells.name() << endl;
574  cells.instance() = mesh.pointsInstance();
575  cells.write();
576  if (surfWriter)
577  {
578  mergeAndWrite(*surfWriter, cells);
579  }
580  }
581  }
582 
583  label nHighAspect = returnReduce(aspectCells.size(), sumOp<label>());
584 
585  if (nHighAspect > 0)
586  {
587  Info<< " <<Writing " << nHighAspect
588  << " cells with high aspect ratio to set "
589  << aspectCells.name() << endl;
590  aspectCells.instance() = mesh.pointsInstance();
591  aspectCells.write();
592  if (surfWriter)
593  {
594  mergeAndWrite(*surfWriter, aspectCells);
595  }
596  }
597  }
598 
599  {
600  faceSet faces(mesh, "zeroAreaFaces", mesh.nFaces()/100+1);
601  if (mesh.checkFaceAreas(true, &faces))
602  {
603  noFailedChecks++;
604 
605  label nFaces = returnReduce(faces.size(), sumOp<label>());
606 
607  if (nFaces > 0)
608  {
609  Info<< " <<Writing " << nFaces
610  << " zero area faces to set " << faces.name() << endl;
611  faces.instance() = mesh.pointsInstance();
612  faces.write();
613  if (surfWriter)
614  {
615  mergeAndWrite(*surfWriter, faces);
616  }
617  }
618  }
619  }
620 
621  {
622  cellSet cells(mesh, "zeroVolumeCells", mesh.nCells()/100+1);
623  if (mesh.checkCellVolumes(true, &cells))
624  {
625  noFailedChecks++;
626 
627  label nCells = returnReduce(cells.size(), sumOp<label>());
628 
629  if (nCells > 0)
630  {
631  Info<< " <<Writing " << nCells
632  << " zero volume cells to set " << cells.name() << endl;
633  cells.instance() = mesh.pointsInstance();
634  cells.write();
635  if (surfWriter)
636  {
637  mergeAndWrite(*surfWriter, cells);
638  }
639  }
640  }
641  }
642 
643  {
644  faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
645  if (mesh.checkFaceOrthogonality(true, &faces))
646  {
647  noFailedChecks++;
648  }
649 
650  label nFaces = returnReduce(faces.size(), sumOp<label>());
651 
652  if (nFaces > 0)
653  {
654  Info<< " <<Writing " << nFaces
655  << " non-orthogonal faces to set " << faces.name() << endl;
656  faces.instance() = mesh.pointsInstance();
657  faces.write();
658  if (surfWriter)
659  {
660  mergeAndWrite(*surfWriter, faces);
661  }
662  }
663  }
664 
665  {
666  faceSet faces(mesh, "wrongOrientedFaces", mesh.nFaces()/100 + 1);
667  if (mesh.checkFacePyramids(true, -SMALL, &faces))
668  {
669  noFailedChecks++;
670 
671  label nFaces = returnReduce(faces.size(), sumOp<label>());
672 
673  if (nFaces > 0)
674  {
675  Info<< " <<Writing " << nFaces
676  << " faces with incorrect orientation to set "
677  << faces.name() << endl;
678  faces.instance() = mesh.pointsInstance();
679  faces.write();
680  if (surfWriter)
681  {
682  mergeAndWrite(*surfWriter, faces);
683  }
684  }
685  }
686  }
687 
688  {
689  faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
690  if (mesh.checkFaceSkewness(true, &faces))
691  {
692  noFailedChecks++;
693 
694  label nFaces = returnReduce(faces.size(), sumOp<label>());
695 
696  if (nFaces > 0)
697  {
698  Info<< " <<Writing " << nFaces
699  << " skew faces to set " << faces.name() << endl;
700  faces.instance() = mesh.pointsInstance();
701  faces.write();
702  if (surfWriter)
703  {
704  mergeAndWrite(*surfWriter, faces);
705  }
706  }
707  }
708  }
709 
710  {
711  faceSet faces(mesh, "coupledFaces", mesh.nFaces()/100 + 1);
712  if (checkCoupledPoints(mesh, true, &faces))
713  {
714  noFailedChecks++;
715 
716  label nFaces = returnReduce(faces.size(), sumOp<label>());
717 
718  if (nFaces > 0)
719  {
720  Info<< " <<Writing " << nFaces
721  << " faces with incorrectly matched 0th (or consecutive)"
722  << " vertex to set "
723  << faces.name() << endl;
724  faces.instance() = mesh.pointsInstance();
725  faces.write();
726  if (surfWriter)
727  {
728  mergeAndWrite(*surfWriter, faces);
729  }
730  }
731  }
732  }
733 
734  if (allGeometry)
735  {
736  faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
737  if
738  (
739  polyMeshTetDecomposition::checkFaceTets
740  (
741  mesh,
742  polyMeshTetDecomposition::minTetQuality,
743  true,
744  &faces
745  )
746  )
747  {
748  noFailedChecks++;
749 
750  label nFaces = returnReduce(faces.size(), sumOp<label>());
751 
752  if (nFaces > 0)
753  {
754  Info<< " <<Writing " << nFaces
755  << " faces with low quality or negative volume "
756  << "decomposition tets to set " << faces.name() << endl;
757  faces.instance() = mesh.pointsInstance();
758  faces.write();
759  if (surfWriter)
760  {
761  mergeAndWrite(*surfWriter, faces);
762  }
763  }
764  }
765  }
766 
767  if (allGeometry)
768  {
769  // Note use of nPoints since don't want edge construction.
770  pointSet points(mesh, "shortEdges", mesh.nPoints()/1000 + 1);
771  if (mesh.checkEdgeLength(true, minDistSqr, &points))
772  {
773  //noFailedChecks++;
774 
775  label nPoints = returnReduce(points.size(), sumOp<label>());
776 
777  if (nPoints > 0)
778  {
779  Info<< " <<Writing " << nPoints
780  << " points on short edges to set " << points.name()
781  << endl;
782  points.instance() = mesh.pointsInstance();
783  points.write();
784  if (setWriter)
785  {
786  mergeAndWrite(*setWriter, points);
787  }
788  }
789  }
790 
791  label nEdgeClose = returnReduce(points.size(), sumOp<label>());
792 
793  if (mesh.checkPointNearness(false, minDistSqr, &points))
794  {
795  //noFailedChecks++;
796 
797  label nPoints = returnReduce(points.size(), sumOp<label>());
798 
799  if (nPoints > nEdgeClose)
800  {
801  pointSet nearPoints(mesh, "nearPoints", points);
802  Info<< " <<Writing " << nPoints
803  << " near (closer than " << Foam::sqrt(minDistSqr)
804  << " apart) points to set " << nearPoints.name() << endl;
805  nearPoints.instance() = mesh.pointsInstance();
806  nearPoints.write();
807  if (setWriter)
808  {
809  mergeAndWrite(*setWriter, nearPoints);
810  }
811  }
812  }
813  }
814 
815  if (allGeometry)
816  {
817  faceSet faces(mesh, "concaveFaces", mesh.nFaces()/100 + 1);
818  if (mesh.checkFaceAngles(true, 10, &faces))
819  {
820  //noFailedChecks++;
821 
822  label nFaces = returnReduce(faces.size(), sumOp<label>());
823 
824  if (nFaces > 0)
825  {
826  Info<< " <<Writing " << nFaces
827  << " faces with concave angles to set " << faces.name()
828  << endl;
829  faces.instance() = mesh.pointsInstance();
830  faces.write();
831  if (surfWriter)
832  {
833  mergeAndWrite(*surfWriter, faces);
834  }
835  }
836  }
837  }
838 
839  if (allGeometry)
840  {
841  faceSet faces(mesh, "warpedFaces", mesh.nFaces()/100 + 1);
842  if (mesh.checkFaceFlatness(true, 0.8, &faces))
843  {
844  //noFailedChecks++;
845 
846  label nFaces = returnReduce(faces.size(), sumOp<label>());
847 
848  if (nFaces > 0)
849  {
850  Info<< " <<Writing " << nFaces
851  << " warped faces to set " << faces.name() << endl;
852  faces.instance() = mesh.pointsInstance();
853  faces.write();
854  if (surfWriter)
855  {
856  mergeAndWrite(*surfWriter, faces);
857  }
858  }
859  }
860  }
861 
862  if (allGeometry)
863  {
864  cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
865  if (mesh.checkCellDeterminant(true, &cells))
866  {
867  noFailedChecks++;
868 
869  label nCells = returnReduce(cells.size(), sumOp<label>());
870 
871  Info<< " <<Writing " << nCells
872  << " under-determined cells to set " << cells.name() << endl;
873  cells.instance() = mesh.pointsInstance();
874  cells.write();
875  if (surfWriter)
876  {
877  mergeAndWrite(*surfWriter, cells);
878  }
879  }
880  }
881 
882  if (allGeometry)
883  {
884  cellSet cells(mesh, "concaveCells", mesh.nCells()/100);
885  if (mesh.checkConcaveCells(true, &cells))
886  {
887  noFailedChecks++;
888 
889  label nCells = returnReduce(cells.size(), sumOp<label>());
890 
891  Info<< " <<Writing " << nCells
892  << " concave cells to set " << cells.name() << endl;
893  cells.instance() = mesh.pointsInstance();
894  cells.write();
895  if (surfWriter)
896  {
897  mergeAndWrite(*surfWriter, cells);
898  }
899  }
900  }
901 
902  if (allGeometry)
903  {
904  faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
905  if (mesh.checkFaceWeight(true, 0.05, &faces))
906  {
907  noFailedChecks++;
908 
909  label nFaces = returnReduce(faces.size(), sumOp<label>());
910 
911  Info<< " <<Writing " << nFaces
912  << " faces with low interpolation weights to set "
913  << faces.name() << endl;
914  faces.instance() = mesh.pointsInstance();
915  faces.write();
916  if (surfWriter)
917  {
918  mergeAndWrite(*surfWriter, faces);
919  }
920  }
921  }
922 
923  if (allGeometry)
924  {
925  faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
926  if (mesh.checkVolRatio(true, 0.01, &faces))
927  {
928  noFailedChecks++;
929 
930  label nFaces = returnReduce(faces.size(), sumOp<label>());
931 
932  Info<< " <<Writing " << nFaces
933  << " faces with low volume ratio cells to set "
934  << faces.name() << endl;
935  faces.instance() = mesh.pointsInstance();
936  faces.write();
937  if (surfWriter)
938  {
939  mergeAndWrite(*surfWriter, faces);
940  }
941  }
942  }
943 
944  if (allGeometry)
945  {
946  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
947 
948  const word tmName(mesh.time().timeName());
949  const word procAndTime(Foam::name(Pstream::myProcNo()) + "_" + tmName);
950 
951  autoPtr<surfaceWriter> patchWriter;
952  if (!surfWriter)
953  {
954  patchWriter.reset(new surfaceWriters::vtkWriter());
955  }
956 
957  surfaceWriter& wr = (surfWriter ? *surfWriter : *patchWriter);
958 
959  // Currently only do AMI checks
960 
961  const fileName outputDir
962  (
963  mesh.time().globalPath()/functionObject::outputPrefix/"checkMesh"
964  );
965 
966  forAll(pbm, patchi)
967  {
968  if (isA<cyclicAMIPolyPatch>(pbm[patchi]))
969  {
970  const cyclicAMIPolyPatch& cpp =
971  refCast<const cyclicAMIPolyPatch>(pbm[patchi]);
972 
973  if (cpp.owner())
974  {
975  Info<< "Calculating AMI weights between owner patch: "
976  << cpp.name() << " and neighbour patch: "
977  << cpp.neighbPatch().name() << endl;
978 
979  const AMIPatchToPatchInterpolation& ami =
980  cpp.AMI();
981 
982  {
983  // Collect geometry
984  labelList pointToGlobal;
985  labelList uniqueMeshPointLabels;
986  autoPtr<globalIndex> globalPoints;
987  autoPtr<globalIndex> globalFaces;
988  faceList mergedFaces;
989  pointField mergedPoints;
990  Foam::PatchTools::gatherAndMerge
991  (
992  mesh,
993  cpp.localFaces(),
994  cpp.meshPoints(),
995  cpp.meshPointMap(),
996 
997  pointToGlobal,
998  uniqueMeshPointLabels,
999  globalPoints,
1000  globalFaces,
1001 
1002  mergedFaces,
1003  mergedPoints
1004  );
1005  // Collect field
1006  scalarField mergedWeights;
1007  globalFaces().gather
1008  (
1009  ami.srcWeightsSum(),
1010  mergedWeights
1011  );
1012 
1013  if (Pstream::master())
1014  {
1015  const word fName
1016  (
1017  "patch" + Foam::name(cpp.index())
1018  + "-src_" + tmName
1019  );
1020 
1021  wr.open
1022  (
1023  mergedPoints,
1024  mergedFaces,
1025  (outputDir / fName),
1026  false // serial - already merged
1027  );
1028 
1029  wr.write("weightsSum", mergedWeights);
1030  wr.clear();
1031  }
1032 
1033  if (isA<cyclicACMIPolyPatch>(pbm[patchi]))
1034  {
1035  const cyclicACMIPolyPatch& pp =
1036  refCast<const cyclicACMIPolyPatch>(pbm[patchi]);
1037 
1038  scalarField mergedMask;
1039  globalFaces().gather
1040  (
1041  pp.mask(),
1042  mergedMask
1043  );
1044 
1045  if (Pstream::master())
1046  {
1047  const word fName
1048  (
1049  "patch" + Foam::name(cpp.index())
1050  + "-src_" + tmName
1051  );
1052 
1053  wr.open
1054  (
1055  mergedPoints,
1056  mergedFaces,
1057  (outputDir / fName),
1058  false // serial - already merged
1059  );
1060 
1061  wr.write("mask", mergedMask);
1062  wr.clear();
1063  }
1064  }
1065  }
1066  {
1067  // Collect geometry
1068  labelList pointToGlobal;
1069  labelList uniqueMeshPointLabels;
1070  autoPtr<globalIndex> globalPoints;
1071  autoPtr<globalIndex> globalFaces;
1072  faceList mergedFaces;
1073  pointField mergedPoints;
1074  Foam::PatchTools::gatherAndMerge
1075  (
1076  mesh,
1077  cpp.neighbPatch().localFaces(),
1078  cpp.neighbPatch().meshPoints(),
1079  cpp.neighbPatch().meshPointMap(),
1080 
1081  pointToGlobal,
1082  uniqueMeshPointLabels,
1083  globalPoints,
1084  globalFaces,
1085 
1086  mergedFaces,
1087  mergedPoints
1088  );
1089  // Collect field
1090  scalarField mergedWeights;
1091  globalFaces().gather
1092  (
1093  ami.tgtWeightsSum(),
1094  mergedWeights
1095  );
1096 
1097  if (Pstream::master())
1098  {
1099  const word fName
1100  (
1101  "patch" + Foam::name(cpp.index())
1102  + "-tgt_" + tmName
1103  );
1104 
1105  wr.open
1106  (
1107  mergedPoints,
1108  mergedFaces,
1109  (outputDir / fName),
1110  false // serial - already merged
1111  );
1112 
1113  wr.write("weightsSum", mergedWeights);
1114  wr.clear();
1115  }
1116 
1117  if (isA<cyclicACMIPolyPatch>(pbm[patchi]))
1118  {
1119  const cyclicACMIPolyPatch& pp =
1120  refCast<const cyclicACMIPolyPatch>(pbm[patchi]);
1121  scalarField mergedMask;
1122  globalFaces().gather
1123  (
1124  pp.neighbPatch().mask(),
1125  mergedMask
1126  );
1127 
1128  if (Pstream::master())
1129  {
1130  const word fName
1131  (
1132  "patch" + Foam::name(cpp.index())
1133  + "-tgt_" + tmName
1134  );
1135 
1136  wr.open
1137  (
1138  mergedPoints,
1139  mergedFaces,
1140  (outputDir / fName),
1141  false // serial - already merged
1142  );
1143 
1144  wr.write("mask", mergedMask);
1145  wr.clear();
1146  }
1147  }
1148  }
1149  }
1150  }
1151  }
1152  }
1153 
1154  return noFailedChecks;
1155 }
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:63
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::VectorSpace< Vector< label >, label, 3 >::one
static const Form one
Definition: VectorSpace.H:112
Foam::checkWedges
bool checkWedges(const polyMesh &, const bool report, const Vector< label > &, labelHashSet *)
cyclicACMIPolyPatch.H
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:88
Foam::Zero
static constexpr const zero Zero
Definition: zero.H:131
PatchTools.H
checkGeometry.H
wedgePolyPatch.H
Foam::PatchTools::gatherAndMerge
static void gatherAndMerge(const scalar mergeDist, const PrimitivePatch< FaceList, PointField > &p, Field< typename PrimitivePatch< FaceList, PointField >::point_type > &mergedPoints, List< typename PrimitivePatch< FaceList, PointField >::face_type > &mergedFaces, labelList &pointMergeMap)
Definition: PatchToolsGatherAndMerge.C:31
Foam::findOppositeWedge
label findOppositeWedge(const polyMesh &, const wedgePolyPatch &)
Foam::dimensioned::name
const word & name() const
Definition: dimensionedType.C:399
unitConversion.H
Unit conversion functions.
Foam::endl
Ostream & endl(Ostream &os)
Definition: Ostream.H:381
polyMesh.H
Foam::checkCoupledPoints
bool checkCoupledPoints(const polyMesh &, const bool report, labelHashSet *)
polyMeshTetDecomposition.H
forAll
#define forAll(list, i)
Definition: stdFoam.H:349
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::Info
messageStream Info
faceSet.H
Foam::radToDeg
constexpr scalar radToDeg(const scalar rad) noexcept
Definition: unitConversion.H:52
Foam::checkGeometry
label checkGeometry(const polyMesh &mesh, const bool allGeometry, autoPtr< surfaceWriter > &surfWriter, const autoPtr< writer< scalar >> &setWriter)
Foam::FatalError
error FatalError
reduce
reduce(hasMovingMesh, orOp< bool >())
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam
Definition: atmBoundaryLayer.C:26
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:139
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:47
Time.H
FatalErrorInFunction
#define FatalErrorInFunction
Definition: error.H:465
edgeHashes.H
forAllConstIters
forAllConstIters(mixture.phases(), phase)
Definition: pEqn.H:28
f
labelList f(nPoints)
Foam::faceList
List< face > faceList
A List of faces.
Definition: faceListFwd.H:41
vtkSurfaceWriter.H
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:137
Foam::acos
dimensionedScalar acos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:261
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::AMIPatchToPatchInterpolation
AMIInterpolation AMIPatchToPatchInterpolation
Definition: AMIPatchToPatchInterpolation.H:31
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Definition: createFields.H:11
Foam::direction
uint8_t direction
Definition: direction.H:46
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:59
Foam::mergeAndWrite
void mergeAndWrite(const polyMesh &mesh, surfaceWriter &writer, const word &name, const indirectPrimitivePatch &setPatch, const fileName &outputDir)
Foam::name
word name(const expressions::valueTypeCode typeCode)
Definition: exprTraits.C:52
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
cellSet.H
coupled
bool coupled(solutionDict.getOrDefault("coupledEnergyField", false))
p0
const volScalarField & p0
Definition: EEqn.H:36
Foam::point
vector point
Point is a vector.
Definition: point.H:37
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys and label hasher.
Definition: HashSet.H:85
writer.H
functionObject.H
checkTools.H
pointSet.H

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17