processorPolyPatch.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-2015 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 "processorPolyPatch.H"
28 #include "dictionary.H"
29 #include "SubField.H"
30 #include "demandDrivenData.H"
31 #include "matchPoints.H"
32 #include "OFstream.H"
33 #include "polyMesh.H"
34 #include "Time.H"
35 #include "transformList.H"
36 #include "PstreamBuffers.H"
37 #include "ConstCirculator.H"
38 
39 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43  defineTypeNameAndDebug(processorPolyPatch, 0);
44  addToRunTimeSelectionTable(polyPatch, processorPolyPatch, dictionary);
45 }
46 
47 
48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
49 
51 (
52  const word& name,
53  const label size,
54  const label start,
55  const label index,
56  const polyBoundaryMesh& bm,
57  const int myProcNo,
58  const int neighbProcNo,
59  const transformType transform,
60  const word& patchType
61 )
62 :
63  coupledPolyPatch(name, size, start, index, bm, patchType, transform),
64  myProcNo_(myProcNo),
65  neighbProcNo_(neighbProcNo),
66  neighbFaceCentres_(),
67  neighbFaceAreas_(),
68  neighbFaceCellCentres_()
69 {}
70 
71 
73 (
74  const word& name,
75  const dictionary& dict,
76  const label index,
77  const polyBoundaryMesh& bm,
78  const word& patchType
79 )
80 :
81  coupledPolyPatch(name, dict, index, bm, patchType),
82  myProcNo_(readLabel(dict.lookup("myProcNo"))),
83  neighbProcNo_(readLabel(dict.lookup("neighbProcNo"))),
84  neighbFaceCentres_(),
85  neighbFaceAreas_(),
86  neighbFaceCellCentres_()
87 {}
88 
89 
91 (
92  const processorPolyPatch& pp,
93  const polyBoundaryMesh& bm
94 )
95 :
96  coupledPolyPatch(pp, bm),
97  myProcNo_(pp.myProcNo_),
98  neighbProcNo_(pp.neighbProcNo_),
99  neighbFaceCentres_(),
100  neighbFaceAreas_(),
101  neighbFaceCellCentres_()
102 {}
103 
104 
106 (
107  const processorPolyPatch& pp,
108  const polyBoundaryMesh& bm,
109  const label index,
110  const label newSize,
111  const label newStart
112 )
113 :
114  coupledPolyPatch(pp, bm, index, newSize, newStart),
115  myProcNo_(pp.myProcNo_),
116  neighbProcNo_(pp.neighbProcNo_),
117  neighbFaceCentres_(),
118  neighbFaceAreas_(),
119  neighbFaceCellCentres_()
120 {}
121 
122 
124 (
125  const processorPolyPatch& pp,
126  const polyBoundaryMesh& bm,
127  const label index,
128  const labelUList& mapAddressing,
129  const label newStart
130 )
131 :
132  coupledPolyPatch(pp, bm, index, mapAddressing, newStart),
133  myProcNo_(pp.myProcNo_),
134  neighbProcNo_(pp.neighbProcNo_),
135  neighbFaceCentres_(),
136  neighbFaceAreas_(),
137  neighbFaceCellCentres_()
138 {}
139 
140 
141 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
142 
144 {
145  neighbPointsPtr_.clear();
146  neighbEdgesPtr_.clear();
147 }
148 
149 
150 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
151 
153 {
154  if (Pstream::parRun())
155  {
156  UOPstream toNeighbProc(neighbProcNo(), pBufs);
157 
158  toNeighbProc
159  << faceCentres()
160  << faceAreas()
161  << faceCellCentres();
162  }
163 }
164 
165 
167 {
168  if (Pstream::parRun())
169  {
170  {
171  UIPstream fromNeighbProc(neighbProcNo(), pBufs);
172 
173  fromNeighbProc
174  >> neighbFaceCentres_
175  >> neighbFaceAreas_
176  >> neighbFaceCellCentres_;
177  }
178 
179  // My normals
180  vectorField faceNormals(size());
181 
182  // Neighbour normals
183  vectorField nbrFaceNormals(neighbFaceAreas_.size());
184 
185  // Face match tolerances
186  scalarField tols = calcFaceTol(*this, points(), faceCentres());
187 
188  // Calculate normals from areas and check
189  forAll(faceNormals, facei)
190  {
191  scalar magSf = mag(faceAreas()[facei]);
192  scalar nbrMagSf = mag(neighbFaceAreas_[facei]);
193  scalar avSf = (magSf + nbrMagSf)/2.0;
194 
195  // For small face area calculation the results of the area
196  // calculation have been found to only be accurate to ~1e-20
197  if (magSf < SMALL || nbrMagSf < SMALL)
198  {
199  // Undetermined normal. Use dummy normal to force separation
200  // check.
201  faceNormals[facei] = point(1, 0, 0);
202  nbrFaceNormals[facei] = -faceNormals[facei];
203  tols[facei] = GREAT;
204  }
205  else if (mag(magSf - nbrMagSf) > matchTolerance()*sqr(tols[facei]))
206  {
207  fileName nm
208  (
209  boundaryMesh().mesh().time().path()
210  /name()+"_faces.obj"
211  );
212 
213  Pout<< "processorPolyPatch::calcGeometry : Writing my "
214  << size()
215  << " faces to OBJ file " << nm << endl;
216 
217  writeOBJ(nm, *this, points());
218 
219  OFstream ccStr
220  (
221  boundaryMesh().mesh().time().path()
222  /name() + "_faceCentresConnections.obj"
223  );
224 
225  Pout<< "processorPolyPatch::calcGeometry :"
226  << " Dumping cell centre lines between"
227  << " corresponding face centres to OBJ file" << ccStr.name()
228  << endl;
229 
230  label vertI = 0;
231 
232  forAll(faceCentres(), faceI)
233  {
234  const point& c0 = neighbFaceCentres_[faceI];
235  const point& c1 = faceCentres()[faceI];
236 
237  writeOBJ(ccStr, c0, c1, vertI);
238  }
239 
241  (
242  "processorPolyPatch::calcGeometry()"
243  ) << "face " << facei << " area does not match neighbour by "
244  << 100*mag(magSf - nbrMagSf)/avSf
245  << "% -- possible face ordering problem." << endl
246  << "patch:" << name()
247  << " my area:" << magSf
248  << " neighbour area:" << nbrMagSf
249  << " matching tolerance:"
250  << matchTolerance()*sqr(tols[facei])
251  << endl
252  << "Mesh face:" << start()+facei
253  << " vertices:"
254  << UIndirectList<point>(points(), operator[](facei))()
255  << endl
256  << "If you are certain your matching is correct"
257  << " you can increase the 'matchTolerance' setting"
258  << " in the patch dictionary in the boundary file."
259  << endl
260  << "Rerun with processor debug flag set for"
261  << " more information." << exit(FatalError);
262  }
263  else
264  {
265  faceNormals[facei] = faceAreas()[facei]/magSf;
266  nbrFaceNormals[facei] = neighbFaceAreas_[facei]/nbrMagSf;
267  }
268  }
269 
270  calcTransformTensors
271  (
272  faceCentres(),
273  neighbFaceCentres_,
274  faceNormals,
275  nbrFaceNormals,
276  matchTolerance()*tols,
277  matchTolerance(),
278  transform()
279  );
280  }
281 }
282 
283 
285 (
286  PstreamBuffers& pBufs,
287  const pointField& p
288 )
289 {
290  polyPatch::movePoints(pBufs, p);
292 }
293 
294 
296 (
297  PstreamBuffers& pBufs,
298  const pointField&
299 )
300 {
302 }
303 
304 
306 {
308 
309  if (Pstream::parRun())
310  {
311  // Express all points as patch face and index in face.
312  labelList pointFace(nPoints());
313  labelList pointIndex(nPoints());
314 
315  for (label patchPointI = 0; patchPointI < nPoints(); patchPointI++)
316  {
317  label faceI = pointFaces()[patchPointI][0];
318 
319  pointFace[patchPointI] = faceI;
320 
321  const face& f = localFaces()[faceI];
322 
323  pointIndex[patchPointI] = findIndex(f, patchPointI);
324  }
325 
326  // Express all edges as patch face and index in face.
327  labelList edgeFace(nEdges());
328  labelList edgeIndex(nEdges());
329 
330  for (label patchEdgeI = 0; patchEdgeI < nEdges(); patchEdgeI++)
331  {
332  label faceI = edgeFaces()[patchEdgeI][0];
333 
334  edgeFace[patchEdgeI] = faceI;
335 
336  const labelList& fEdges = faceEdges()[faceI];
337 
338  edgeIndex[patchEdgeI] = findIndex(fEdges, patchEdgeI);
339  }
340 
341  UOPstream toNeighbProc(neighbProcNo(), pBufs);
342 
343  toNeighbProc
344  << pointFace
345  << pointIndex
346  << edgeFace
347  << edgeIndex;
348  }
349 }
350 
351 
353 {
354  // For completeness
355  polyPatch::updateMesh(pBufs);
356 
357  neighbPointsPtr_.clear();
358  neighbEdgesPtr_.clear();
359 
360  if (Pstream::parRun())
361  {
362  labelList nbrPointFace;
363  labelList nbrPointIndex;
364  labelList nbrEdgeFace;
365  labelList nbrEdgeIndex;
366 
367  {
368  // !Note: there is one situation where the opposite points and
369  // do not exactly match and that is while we are distributing
370  // meshes and multiple parts come together from different
371  // processors. This can temporarily create the situation that
372  // we have points which will be merged out later but we still
373  // need the face connectivity to be correct.
374  // So: cannot check here on same points and edges.
375 
376  UIPstream fromNeighbProc(neighbProcNo(), pBufs);
377 
378  fromNeighbProc
379  >> nbrPointFace
380  >> nbrPointIndex
381  >> nbrEdgeFace
382  >> nbrEdgeIndex;
383  }
384 
385  // Convert neighbour faces and indices into face back into
386  // my edges and points.
387 
388  // Convert points.
389  // ~~~~~~~~~~~~~~~
390 
391  neighbPointsPtr_.reset(new labelList(nPoints(), -1));
392  labelList& neighbPoints = neighbPointsPtr_();
393 
394  forAll(nbrPointFace, nbrPointI)
395  {
396  // Find face and index in face on this side.
397  const face& f = localFaces()[nbrPointFace[nbrPointI]];
398 
399  label index = (f.size() - nbrPointIndex[nbrPointI]) % f.size();
400  label patchPointI = f[index];
401 
402  if (neighbPoints[patchPointI] == -1)
403  {
404  // First reference of point
405  neighbPoints[patchPointI] = nbrPointI;
406  }
407  else if (neighbPoints[patchPointI] >= 0)
408  {
409  // Point already visited. Mark as duplicate.
410  neighbPoints[patchPointI] = -2;
411  }
412  }
413 
414  // Reset all duplicate entries to -1.
415  forAll(neighbPoints, patchPointI)
416  {
417  if (neighbPoints[patchPointI] == -2)
418  {
419  neighbPoints[patchPointI] = -1;
420  }
421  }
422 
423  // Convert edges.
424  // ~~~~~~~~~~~~~~
425 
426  neighbEdgesPtr_.reset(new labelList(nEdges(), -1));
427  labelList& neighbEdges = neighbEdgesPtr_();
428 
429  forAll(nbrEdgeFace, nbrEdgeI)
430  {
431  // Find face and index in face on this side.
432  const labelList& f = faceEdges()[nbrEdgeFace[nbrEdgeI]];
433  label index = (f.size() - nbrEdgeIndex[nbrEdgeI] - 1) % f.size();
434  label patchEdgeI = f[index];
435 
436  if (neighbEdges[patchEdgeI] == -1)
437  {
438  // First reference of edge
439  neighbEdges[patchEdgeI] = nbrEdgeI;
440  }
441  else if (neighbEdges[patchEdgeI] >= 0)
442  {
443  // Edge already visited. Mark as duplicate.
444  neighbEdges[patchEdgeI] = -2;
445  }
446  }
447 
448  // Reset all duplicate entries to -1.
449  forAll(neighbEdges, patchEdgeI)
450  {
451  if (neighbEdges[patchEdgeI] == -2)
452  {
453  neighbEdges[patchEdgeI] = -1;
454  }
455  }
456 
457  // Remove any addressing used for shared points/edges calculation
458  // since mostly not needed.
460  }
461 }
462 
463 
465 {
466  if (!neighbPointsPtr_.valid())
467  {
468  FatalErrorIn("processorPolyPatch::neighbPoints() const")
469  << "No extended addressing calculated for patch " << name()
470  << abort(FatalError);
471  }
472  return neighbPointsPtr_();
473 }
474 
475 
477 {
478  if (!neighbEdgesPtr_.valid())
479  {
480  FatalErrorIn("processorPolyPatch::neighbEdges() const")
481  << "No extended addressing calculated for patch " << name()
482  << abort(FatalError);
483  }
484  return neighbEdgesPtr_();
485 }
486 
487 
489 (
490  PstreamBuffers& pBufs,
491  const primitivePatch& pp
492 ) const
493 {
494  if
495  (
496  !Pstream::parRun()
497  || transform() == NOORDERING
498  )
499  {
500  return;
501  }
502 
503  if (debug)
504  {
505  fileName nm
506  (
507  boundaryMesh().mesh().time().path()
508  /name()+"_faces.obj"
509  );
510  Pout<< "processorPolyPatch::order : Writing my " << pp.size()
511  << " faces to OBJ file " << nm << endl;
512  writeOBJ(nm, pp, pp.points());
513 
514  // Calculate my face centres
515  const pointField& fc = pp.faceCentres();
516 
517  OFstream localStr
518  (
519  boundaryMesh().mesh().time().path()
520  /name() + "_localFaceCentres.obj"
521  );
522  Pout<< "processorPolyPatch::order : "
523  << "Dumping " << fc.size()
524  << " local faceCentres to " << localStr.name() << endl;
525 
526  forAll(fc, faceI)
527  {
528  writeOBJ(localStr, fc[faceI]);
529  }
530  }
531 
532  if (owner())
533  {
534  if (transform() == COINCIDENTFULLMATCH)
535  {
536  // Pass the patch points and faces across
537  UOPstream toNeighbour(neighbProcNo(), pBufs);
538  toNeighbour << pp.localPoints()
539  << pp.localFaces();
540  }
541  else
542  {
543  const pointField& ppPoints = pp.points();
544 
545  pointField anchors(getAnchorPoints(pp, ppPoints, transform()));
546 
547  // Get the average of the points of each face. This is needed in
548  // case the face centroid calculation is incorrect due to the face
549  // having a very high aspect ratio.
550  pointField facePointAverages(pp.size(), point::zero);
551  forAll(pp, fI)
552  {
553  const labelList& facePoints = pp[fI];
554 
555  forAll(facePoints, pI)
556  {
557  facePointAverages[fI] += ppPoints[facePoints[pI]];
558  }
559 
560  facePointAverages[fI] /= facePoints.size();
561  }
562 
563  // Now send all info over to the neighbour
564  UOPstream toNeighbour(neighbProcNo(), pBufs);
565  toNeighbour << pp.faceCentres() << pp.faceNormals()
566  << anchors << facePointAverages;
567  }
568  }
569 }
570 
571 
572 // Returns rotation.
573 // + -1 : no match
574 // + 0 : match
575 // + >0 : match if rotated clockwise by this amount
577 (
578  const face& a,
579  const pointField& aPts,
580  const face& b,
581  const pointField& bPts,
582  const bool sameOrientation,
583  const scalar absTolSqr,
584  scalar& matchDistSqr
585 )
586 {
587  if (a.size() != b.size())
588  {
589  return -1;
590  }
591 
592  enum CirculatorBase::direction circulateDirection
594 
595  if (!sameOrientation)
596  {
597  circulateDirection = CirculatorBase::ANTICLOCKWISE;
598  }
599 
600  label matchFp = -1;
601 
602  scalar closestMatchDistSqr = sqr(GREAT);
603 
604  ConstCirculator<face> aCirc(a);
605  ConstCirculator<face> bCirc(b);
606 
607  do
608  {
609  const scalar diffSqr = magSqr(aPts[aCirc()] - bPts[bCirc()]);
610 
611  if (diffSqr < absTolSqr)
612  {
613  // Found a matching point. Set the fulcrum of b to the iterator
614  ConstCirculator<face> bCirc2 = bCirc;
615  ++aCirc;
616 
617  bCirc2.setFulcrumToIterator();
618 
619  if (!sameOrientation)
620  {
621  --bCirc2;
622  }
623  else
624  {
625  ++bCirc2;
626  }
627 
628  matchDistSqr = diffSqr;
629 
630  do
631  {
632  const scalar diffSqr2 = magSqr(aPts[aCirc()] - bPts[bCirc2()]);
633 
634  if (diffSqr2 > absTolSqr)
635  {
636  // No match.
637  break;
638  }
639 
640  matchDistSqr += diffSqr2;
641  }
642  while
643  (
645  bCirc2.circulate(circulateDirection)
646  );
647 
648  if (!aCirc.circulate())
649  {
650  if (matchDistSqr < closestMatchDistSqr)
651  {
652  closestMatchDistSqr = matchDistSqr;
653 
654  if (!sameOrientation)
655  {
656  matchFp = a.size() - bCirc.nRotations();
657  }
658  else
659  {
660  matchFp = bCirc.nRotations();
661  }
662 
663  if (closestMatchDistSqr == 0)
664  {
665  break;
666  }
667  }
668  }
669 
670  // Reset aCirc
671  aCirc.setIteratorToFulcrum();
672  }
673 
674  } while (bCirc.circulate(circulateDirection));
675 
676  matchDistSqr = closestMatchDistSqr;
677 
678  return matchFp;
679 }
680 
681 
682 // Return new ordering. Ordering is -faceMap: for every face index
683 // the new face -rotation:for every new face the clockwise shift
684 // of the original face. Return false if nothing changes (faceMap
685 // is identity, rotation is 0)
687 (
688  PstreamBuffers& pBufs,
689  const primitivePatch& pp,
690  labelList& faceMap,
691  labelList& rotation
692 ) const
693 {
694  // Note: we only get the faces that originate from internal faces.
695 
696  if
697  (
698  !Pstream::parRun()
699  || transform() == NOORDERING
700  )
701  {
702  return false;
703  }
704 
705  faceMap.setSize(pp.size());
706  faceMap = -1;
707 
708  rotation.setSize(pp.size());
709  rotation = 0;
710 
711  bool change = false;
712 
713  if (owner())
714  {
715  // Do nothing (i.e. identical mapping, zero rotation).
716  // See comment at top.
717  forAll(faceMap, patchFaceI)
718  {
719  faceMap[patchFaceI] = patchFaceI;
720  }
721 
722  if (transform() != COINCIDENTFULLMATCH)
723  {
724  const pointField& ppPoints = pp.points();
725 
726  pointField anchors(getAnchorPoints(pp, ppPoints, transform()));
727 
728  // Calculate typical distance from face centre
729  scalarField tols
730  (
731  matchTolerance()*calcFaceTol(pp, pp.points(), pp.faceCentres())
732  );
733 
734  forAll(faceMap, patchFaceI)
735  {
736  const point& wantedAnchor = anchors[patchFaceI];
737 
738  rotation[patchFaceI] = getRotation
739  (
740  ppPoints,
741  pp[patchFaceI],
742  wantedAnchor,
743  tols[patchFaceI]
744  );
745 
746  if (rotation[patchFaceI] > 0)
747  {
748  change = true;
749  }
750  }
751  }
752 
753  return change;
754  }
755  else
756  {
757  if (transform() == COINCIDENTFULLMATCH)
758  {
759  vectorField masterPts;
760  faceList masterFaces;
761 
762  {
763  UIPstream fromNeighbour(neighbProcNo(), pBufs);
764  fromNeighbour >> masterPts >> masterFaces;
765  }
766 
767  const pointField& localPts = pp.localPoints();
768  const faceList& localFaces = pp.localFaces();
769 
770  label nMatches = 0;
771 
772  forAll(pp, lFaceI)
773  {
774  const face& localFace = localFaces[lFaceI];
775  label faceRotation = -1;
776 
777  const scalar absTolSqr = sqr(SMALL);
778 
779  scalar closestMatchDistSqr = sqr(GREAT);
780  scalar matchDistSqr = sqr(GREAT);
781  label closestFaceMatch = -1;
782  label closestFaceRotation = -1;
783 
784  forAll(masterFaces, mFaceI)
785  {
786  const face& masterFace = masterFaces[mFaceI];
787 
788  faceRotation = matchFace
789  (
790  localFace,
791  localPts,
792  masterFace,
793  masterPts,
794  false,
795  absTolSqr,
796  matchDistSqr
797  );
798 
799  if
800  (
801  faceRotation != -1
802  && matchDistSqr < closestMatchDistSqr
803  )
804  {
805  closestMatchDistSqr = matchDistSqr;
806  closestFaceMatch = mFaceI;
807  closestFaceRotation = faceRotation;
808  }
809 
810  if (closestMatchDistSqr == 0)
811  {
812  break;
813  }
814  }
815 
816  if (closestFaceRotation != -1 && closestMatchDistSqr == 0)
817  {
818  faceMap[lFaceI] = closestFaceMatch;
819 
820  rotation[lFaceI] = closestFaceRotation;
821 
822  if (lFaceI != closestFaceMatch || closestFaceRotation > 0)
823  {
824  change = true;
825  }
826 
827  nMatches++;
828  }
829  else
830  {
831  Pout<< "Number of matches = " << nMatches << " / "
832  << pp.size() << endl;
833 
834  pointField pts(localFace.size());
835  forAll(localFace, pI)
836  {
837  const label localPtI = localFace[pI];
838  pts[pI] = localPts[localPtI];
839  }
840 
841  FatalErrorIn("Foam::processorPolyPatch::order(...) const")
842  << "No match for face " << localFace << nl << pts
843  << abort(FatalError);
844  }
845  }
846 
847  return change;
848  }
849  else
850  {
851  vectorField masterCtrs;
852  vectorField masterNormals;
853  vectorField masterAnchors;
854  vectorField masterFacePointAverages;
855 
856  // Receive data from neighbour
857  {
858  UIPstream fromNeighbour(neighbProcNo(), pBufs);
859  fromNeighbour >> masterCtrs >> masterNormals
860  >> masterAnchors >> masterFacePointAverages;
861  }
862 
863  // Calculate typical distance from face centre
864  scalarField tols
865  (
866  matchTolerance()*calcFaceTol(pp, pp.points(), pp.faceCentres())
867  );
868 
869  if (debug || masterCtrs.size() != pp.size())
870  {
871  {
872  OFstream nbrStr
873  (
874  boundaryMesh().mesh().time().path()
875  /name() + "_nbrFaceCentres.obj"
876  );
877  Pout<< "processorPolyPatch::order : "
878  << "Dumping neighbour faceCentres to " << nbrStr.name()
879  << endl;
880  forAll(masterCtrs, faceI)
881  {
882  writeOBJ(nbrStr, masterCtrs[faceI]);
883  }
884  }
885 
886  if (masterCtrs.size() != pp.size())
887  {
889  (
890  "processorPolyPatch::order(const primitivePatch&"
891  ", labelList&, labelList&) const"
892  ) << "in patch:" << name() << " : "
893  << "Local size of patch is " << pp.size() << " (faces)."
894  << endl
895  << "Received from neighbour " << masterCtrs.size()
896  << " faceCentres!"
897  << abort(FatalError);
898  }
899  }
900 
901  // Geometric match of face centre vectors
902  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
903 
904  // 1. Try existing ordering and transformation
905  bool matchedAll = matchPoints
906  (
907  pp.faceCentres(),
908  masterCtrs,
909  pp.faceNormals(),
910  masterNormals,
911  tols,
912  false,
913  faceMap
914  );
915 
916  // Try using face point average for matching
917  if (!matchedAll)
918  {
919  const pointField& ppPoints = pp.points();
920 
921  pointField facePointAverages(pp.size(), point::zero);
922  forAll(pp, fI)
923  {
924  const labelList& facePoints = pp[fI];
925 
926  forAll(facePoints, pI)
927  {
928  facePointAverages[fI] += ppPoints[facePoints[pI]];
929  }
930 
931  facePointAverages[fI] /= facePoints.size();
932  }
933 
934  scalarField tols2
935  (
936  calcFaceTol(pp, pp.points(), facePointAverages)
937  );
938 
939  labelList faceMap2(faceMap.size(), -1);
940  matchedAll = matchPoints
941  (
942  facePointAverages,
943  masterFacePointAverages,
944  pp.faceNormals(),
945  masterNormals,
946  tols2,
947  true,
948  faceMap2
949  );
950 
951  forAll(faceMap, oldFaceI)
952  {
953  if (faceMap[oldFaceI] == -1)
954  {
955  faceMap[oldFaceI] = faceMap2[oldFaceI];
956  }
957  }
958 
959  matchedAll = true;
960 
961  forAll(faceMap, oldFaceI)
962  {
963  if (faceMap[oldFaceI] == -1)
964  {
965  matchedAll = false;
966  }
967  }
968  }
969 
970  if (!matchedAll || debug)
971  {
972  // Dump faces
973  fileName str
974  (
975  boundaryMesh().mesh().time().path()
976  /name() + "_faces.obj"
977  );
978  Pout<< "processorPolyPatch::order :"
979  << " Writing faces to OBJ file " << str.name() << endl;
980  writeOBJ(str, pp, pp.points());
981 
982  OFstream ccStr
983  (
984  boundaryMesh().mesh().time().path()
985  /name() + "_faceCentresConnections.obj"
986  );
987 
988  Pout<< "processorPolyPatch::order :"
989  << " Dumping newly found match as lines between"
990  << " corresponding face centres to OBJ file "
991  << ccStr.name()
992  << endl;
993 
994  label vertI = 0;
995 
996  forAll(pp.faceCentres(), faceI)
997  {
998  label masterFaceI = faceMap[faceI];
999 
1000  if (masterFaceI != -1)
1001  {
1002  const point& c0 = masterCtrs[masterFaceI];
1003  const point& c1 = pp.faceCentres()[faceI];
1004  writeOBJ(ccStr, c0, c1, vertI);
1005  }
1006  }
1007  }
1008 
1009  if (!matchedAll)
1010  {
1012  (
1013  "processorPolyPatch::order(const primitivePatch&"
1014  ", labelList&, labelList&) const"
1015  ) << "in patch:" << name() << " : "
1016  << "Cannot match vectors to faces on both sides of patch"
1017  << endl
1018  << " masterCtrs[0]:" << masterCtrs[0] << endl
1019  << " ctrs[0]:" << pp.faceCentres()[0] << endl
1020  << " Check your topology changes or maybe you have"
1021  << " multiple separated (from cyclics) processor patches"
1022  << endl
1023  << " Continuing with incorrect face ordering from now on"
1024  << endl;
1025 
1026  return false;
1027  }
1028 
1029  // Set rotation.
1030  forAll(faceMap, oldFaceI)
1031  {
1032  // The face f will be at newFaceI (after morphing) and we want
1033  // its anchorPoint (= f[0]) to align with the anchorpoint for
1034  // the corresponding face on the other side.
1035 
1036  label newFaceI = faceMap[oldFaceI];
1037 
1038  const point& wantedAnchor = masterAnchors[newFaceI];
1039 
1040  rotation[newFaceI] = getRotation
1041  (
1042  pp.points(),
1043  pp[oldFaceI],
1044  wantedAnchor,
1045  tols[oldFaceI]
1046  );
1047 
1048  if (rotation[newFaceI] == -1)
1049  {
1051  (
1052  "processorPolyPatch::order(const primitivePatch&"
1053  ", labelList&, labelList&) const"
1054  ) << "in patch " << name()
1055  << " : "
1056  << "Cannot find point on face " << pp[oldFaceI]
1057  << " with vertices "
1058  << UIndirectList<point>(pp.points(), pp[oldFaceI])()
1059  << " that matches point " << wantedAnchor
1060  << " when matching the halves of processor patch "
1061  << name()
1062  << "Continuing with incorrect face ordering from now on"
1063  << endl;
1064 
1065  return false;
1066  }
1067  }
1068 
1069  forAll(faceMap, faceI)
1070  {
1071  if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1072  {
1073  return true;
1074  }
1075  }
1076 
1077  return false;
1078  }
1079  }
1080 }
1081 
1082 
1084 {
1086  os.writeKeyword("myProcNo") << myProcNo_
1087  << token::END_STATEMENT << nl;
1088  os.writeKeyword("neighbProcNo") << neighbProcNo_
1089  << token::END_STATEMENT << nl;
1090 }
1091 
1092 
1093 // ************************************************************************* //
#define SeriousErrorIn(functionName)
Report an error message using Foam::SeriousError.
Output to file stream.
Definition: OFstream.H:81
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:380
const pointField & points
word name() const
Return file name (part beyond last /)
Definition: fileName.C:206
void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
void movePoints(PstreamBuffers &, const pointField &)
Correct patches after moving points.
vector point
Point is a vector.
Definition: point.H:41
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
dimensioned< scalar > mag(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Output inter-processor communications stream operating on external buffer.
Definition: UOPstream.H:54
Input inter-processor communications stream operating on external buffer.
Definition: UIPstream.H:53
labelList f(nPoints)
Walks over a container as if it were circular. The container must have the following members defined:...
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
A class for handling words, derived from string.
Definition: word.H:59
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
Addressing for all faces on surface of mesh. Can either be read from polyMesh or from triSurface...
Definition: boundaryMesh.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
Definition: polyPatch.C:62
dynamicFvMesh & mesh
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
direction
Direction type enumeration.
Namespace for OpenFOAM.
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
label readLabel(Istream &is)
Definition: label.H:64
virtual ~processorPolyPatch()
Destructor.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Field< PointType > & localPoints() const
Return pointField of points in patch.
static const char nl
Definition: Ostream.H:260
const Field< PointType > & faceCentres() const
Return face centres for patch.
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Neighbour processor patch.
A List with indirect addressing.
Definition: fvMatrix.H:106
const labelList & neighbEdges() const
Return neighbour edge labels. WIP.
Foam::polyBoundaryMesh.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
void setIteratorToFulcrum()
Set the iterator to the current position of the fulcrum.
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:118
#define forAll(list, i)
Definition: UList.H:421
Macros for easy insertion into run-time selection tables.
Template functions to aid in the implementation of demand driven data.
const dimensionedScalar c1
First radiation constant: default SI units: [W/m2].
void setFulcrumToIterator()
Set the fulcrum to the current position of the iterator.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Field< PointType > & faceNormals() const
Return face normals for patch.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:59
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
A list of faces which address into the list of points.
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
Determine correspondence between points. See below.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
error FatalError
processorPolyPatch(const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm, const int myProcNo, const int neighbProcNo, const transformType transform=UNKNOWN, const word &patchType=typeName)
Construct from components.
List< label > labelList
A List of labels.
Definition: labelList.H:56
Spatial transformation functions for primitive fields.
A class for handling file names.
Definition: fileName.H:69
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
static const Vector zero
Definition: Vector.H:80
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:33
virtual void initOrder(PstreamBuffers &, const primitivePatch &) const
Initialize ordering for primitivePatch. Does not.
bool circulate(const CirculatorBase::direction dir=NONE)
Circulate around the list in the given direction.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
virtual bool order(PstreamBuffers &, const primitivePatch &, labelList &faceMap, labelList &rotation) const
Return new ordering for primitivePatch.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Definition: polyPatch.H:115
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static label matchFace(const face &localFace, const pointField &localPts, const face &masterFace, const pointField &masterPts, const bool sameOrientation, const scalar absTolSqr, scalar &matchDistSqr)
difference_type nRotations() const
Return the distance between the iterator and the fulcrum. This is.
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
label nPoints
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
const labelList & neighbPoints() const
Return neighbour point labels. WIP.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
const Field< PointType > & points() const
Return reference to global points.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
virtual void movePoints(PstreamBuffers &, const pointField &p)
Correct patches after moving points.
Definition: polyPatch.C:57