FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/dictionary.H>
29 #include <OpenFOAM/SubField.H>
31 #include <OpenFOAM/matchPoints.H>
32 #include <OpenFOAM/OFstream.H>
33 #include <OpenFOAM/polyMesh.H>
34 #include <OpenFOAM/Time.H>
35 #include <OpenFOAM/transformList.H>
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  defineTypeNameAndDebug(processorPolyPatch, 0);
42  addToRunTimeSelectionTable(polyPatch, processorPolyPatch, dictionary);
43 }
44 
45 
46 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
47 
49 (
50  const word& name,
51  const label size,
52  const label start,
53  const label index,
54  const polyBoundaryMesh& bm,
55  const int myProcNo,
56  const int neighbProcNo
57 )
58 :
59  coupledPolyPatch(name, size, start, index, bm),
60  myProcNo_(myProcNo),
61  neighbProcNo_(neighbProcNo),
62  neighbFaceCentres_(),
63  neighbFaceAreas_(),
64  neighbFaceCellCentres_(),
65  neighbPointsPtr_(NULL),
66  neighbEdgesPtr_(NULL)
67 {}
68 
69 
71 (
72  const word& name,
73  const dictionary& dict,
74  const label index,
75  const polyBoundaryMesh& bm
76 )
77 :
78  coupledPolyPatch(name, dict, index, bm),
79  myProcNo_(readLabel(dict.lookup("myProcNo"))),
80  neighbProcNo_(readLabel(dict.lookup("neighbProcNo"))),
81  neighbFaceCentres_(),
82  neighbFaceAreas_(),
83  neighbFaceCellCentres_(),
84  neighbPointsPtr_(NULL),
85  neighbEdgesPtr_(NULL)
86 {}
87 
88 
90 (
91  const processorPolyPatch& pp,
92  const polyBoundaryMesh& bm
93 )
94 :
95  coupledPolyPatch(pp, bm),
96  myProcNo_(pp.myProcNo_),
97  neighbProcNo_(pp.neighbProcNo_),
98  neighbFaceCentres_(),
99  neighbFaceAreas_(),
100  neighbFaceCellCentres_(),
101  neighbPointsPtr_(NULL),
102  neighbEdgesPtr_(NULL)
103 {}
104 
105 
107 (
108  const processorPolyPatch& pp,
109  const polyBoundaryMesh& bm,
110  const label index,
111  const label newSize,
112  const label newStart
113 )
114 :
115  coupledPolyPatch(pp, bm, index, newSize, newStart),
116  myProcNo_(pp.myProcNo_),
117  neighbProcNo_(pp.neighbProcNo_),
118  neighbFaceCentres_(),
119  neighbFaceAreas_(),
120  neighbFaceCellCentres_(),
121  neighbPointsPtr_(NULL),
122  neighbEdgesPtr_(NULL)
123 {}
124 
125 
126 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
127 
129 {
130  deleteDemandDrivenData(neighbPointsPtr_);
131  deleteDemandDrivenData(neighbEdgesPtr_);
132 }
133 
134 
135 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
136 
138 {
139  if (Pstream::parRun())
140  {
141  OPstream toNeighbProc
142  (
144  neighbProcNo(),
145  3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
146  );
147 
148  toNeighbProc
149  << faceCentres()
150  << faceAreas()
151  << faceCellCentres();
152  }
153 }
154 
155 
157 {
158  if (Pstream::parRun())
159  {
160  {
161  IPstream fromNeighbProc
162  (
164  neighbProcNo(),
165  3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
166  );
167  fromNeighbProc
168  >> neighbFaceCentres_
169  >> neighbFaceAreas_
170  >> neighbFaceCellCentres_;
171  }
172 
173  // My normals
174  vectorField faceNormals(size());
175 
176  // Neighbour normals
177  vectorField nbrFaceNormals(neighbFaceAreas_.size());
178 
179  // Calculate normals from areas and check
180  forAll(faceNormals, facei)
181  {
182  scalar magSf = mag(faceAreas()[facei]);
183  scalar nbrMagSf = mag(neighbFaceAreas_[facei]);
184  scalar avSf = (magSf + nbrMagSf)/2.0;
185 
186  if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
187  {
188  // Undetermined normal. Use dummy normal to force separation
189  // check. (note use of sqrt(VSMALL) since that is how mag
190  // scales)
191  faceNormals[facei] = point(1, 0, 0);
192  nbrFaceNormals[facei] = faceNormals[facei];
193  }
194  else if (mag(magSf - nbrMagSf)/avSf > coupledPolyPatch::matchTol)
195  {
197  (
198  "processorPolyPatch::calcGeometry()"
199  ) << "face " << facei << " area does not match neighbour by "
200  << 100*mag(magSf - nbrMagSf)/avSf
201  << "% -- possible face ordering problem." << endl
202  << "patch:" << name()
203  << " my area:" << magSf
204  << " neighbour area:" << nbrMagSf
205  << " matching tolerance:" << coupledPolyPatch::matchTol
206  << endl
207  << "Mesh face:" << start()+facei
208  << " vertices:"
209  << UIndirectList<point>(points(), operator[](facei))()
210  << endl
211  << "Rerun with processor debug flag set for"
212  << " more information." << exit(FatalError);
213  }
214  else
215  {
216  faceNormals[facei] = faceAreas()[facei]/magSf;
217  nbrFaceNormals[facei] = neighbFaceAreas_[facei]/nbrMagSf;
218  }
219  }
220 
221  calcTransformTensors
222  (
223  faceCentres(),
224  neighbFaceCentres_,
225  faceNormals,
226  nbrFaceNormals,
227  calcFaceTol(*this, points(), faceCentres())
228  );
229  }
230 }
231 
232 
234 {
237 }
238 
239 
241 {
243 }
244 
245 
247 {
249 
250  deleteDemandDrivenData(neighbPointsPtr_);
251  deleteDemandDrivenData(neighbEdgesPtr_);
252 
253  if (Pstream::parRun())
254  {
255  // Express all points as patch face and index in face.
256  labelList pointFace(nPoints());
257  labelList pointIndex(nPoints());
258 
259  for (label patchPointI = 0; patchPointI < nPoints(); patchPointI++)
260  {
261  label faceI = pointFaces()[patchPointI][0];
262 
263  pointFace[patchPointI] = faceI;
264 
265  const face& f = localFaces()[faceI];
266 
267  pointIndex[patchPointI] = findIndex(f, patchPointI);
268  }
269 
270  // Express all edges as patch face and index in face.
271  labelList edgeFace(nEdges());
272  labelList edgeIndex(nEdges());
273 
274  for (label patchEdgeI = 0; patchEdgeI < nEdges(); patchEdgeI++)
275  {
276  label faceI = edgeFaces()[patchEdgeI][0];
277 
278  edgeFace[patchEdgeI] = faceI;
279 
280  const labelList& fEdges = faceEdges()[faceI];
281 
282  edgeIndex[patchEdgeI] = findIndex(fEdges, patchEdgeI);
283  }
284 
285  OPstream toNeighbProc
286  (
288  neighbProcNo(),
289  8*sizeof(label) // four headers of labelList
290  + 2*nPoints()*sizeof(label) // two point-based labellists
291  + 2*nEdges()*sizeof(label) // two edge-based labelLists
292  );
293 
294  toNeighbProc
295  << pointFace
296  << pointIndex
297  << edgeFace
298  << edgeIndex;
299  }
300 }
301 
302 
304 {
305  // For completeness
307 
308  if (Pstream::parRun())
309  {
310  labelList nbrPointFace;
311  labelList nbrPointIndex;
312  labelList nbrEdgeFace;
313  labelList nbrEdgeIndex;
314 
315  {
316  // Note cannot predict exact size since opposite nPoints might
317  // be different from one over here.
318  IPstream fromNeighbProc(Pstream::blocking, neighbProcNo());
319 
320  fromNeighbProc
321  >> nbrPointFace
322  >> nbrPointIndex
323  >> nbrEdgeFace
324  >> nbrEdgeIndex;
325  }
326 
327  // Convert neighbour faces and indices into face back into
328  // my edges and points.
329 
330  // Convert points.
331  // ~~~~~~~~~~~~~~~
332 
333  neighbPointsPtr_ = new labelList(nPoints(), -1);
334  labelList& neighbPoints = *neighbPointsPtr_;
335 
336  forAll(nbrPointFace, nbrPointI)
337  {
338  // Find face and index in face on this side.
339  const face& f = localFaces()[nbrPointFace[nbrPointI]];
340  label index = (f.size() - nbrPointIndex[nbrPointI]) % f.size();
341  label patchPointI = f[index];
342 
343  if (neighbPoints[patchPointI] == -1)
344  {
345  // First reference of point
346  neighbPoints[patchPointI] = nbrPointI;
347  }
348  else if (neighbPoints[patchPointI] >= 0)
349  {
350  // Point already visited. Mark as duplicate.
351  neighbPoints[patchPointI] = -2;
352  }
353  }
354 
355  // Reset all duplicate entries to -1.
356  forAll(neighbPoints, patchPointI)
357  {
358  if (neighbPoints[patchPointI] == -2)
359  {
360  neighbPoints[patchPointI] = -1;
361  }
362  }
363 
364  // Convert edges.
365  // ~~~~~~~~~~~~~~
366 
367  neighbEdgesPtr_ = new labelList(nEdges(), -1);
368  labelList& neighbEdges = *neighbEdgesPtr_;
369 
370  forAll(nbrEdgeFace, nbrEdgeI)
371  {
372  // Find face and index in face on this side.
373  const labelList& f = faceEdges()[nbrEdgeFace[nbrEdgeI]];
374  label index = (f.size() - nbrEdgeIndex[nbrEdgeI] - 1) % f.size();
375  label patchEdgeI = f[index];
376 
377  if (neighbEdges[patchEdgeI] == -1)
378  {
379  // First reference of edge
380  neighbEdges[patchEdgeI] = nbrEdgeI;
381  }
382  else if (neighbEdges[patchEdgeI] >= 0)
383  {
384  // Edge already visited. Mark as duplicate.
385  neighbEdges[patchEdgeI] = -2;
386  }
387  }
388 
389  // Reset all duplicate entries to -1.
390  forAll(neighbEdges, patchEdgeI)
391  {
392  if (neighbEdges[patchEdgeI] == -2)
393  {
394  neighbEdges[patchEdgeI] = -1;
395  }
396  }
397 
398  // Remove any addressing used for shared points/edges calculation
400  }
401 }
402 
403 
405 {
406  if (!neighbPointsPtr_)
407  {
408  FatalErrorIn("processorPolyPatch::neighbPoints() const")
409  << "No extended addressing calculated for patch " << name()
410  << abort(FatalError);
411  }
412  return *neighbPointsPtr_;
413 }
414 
415 
417 {
418  if (!neighbEdgesPtr_)
419  {
420  FatalErrorIn("processorPolyPatch::neighbEdges() const")
421  << "No extended addressing calculated for patch " << name()
422  << abort(FatalError);
423  }
424  return *neighbEdgesPtr_;
425 }
426 
427 
429 {
430  if (!Pstream::parRun())
431  {
432  return;
433  }
434 
435  if (debug)
436  {
437  fileName nm
438  (
439  boundaryMesh().mesh().time().path()
440  /name()+"_faces.obj"
441  );
442  Pout<< "processorPolyPatch::order : Writing my " << pp.size()
443  << " faces to OBJ file " << nm << endl;
444  writeOBJ(nm, pp, pp.points());
445 
446  // Calculate my face centres
447  pointField ctrs(calcFaceCentres(pp, pp.points()));
448 
449  OFstream localStr
450  (
451  boundaryMesh().mesh().time().path()
452  /name() + "_localFaceCentres.obj"
453  );
454  Pout<< "processorPolyPatch::order : "
455  << "Dumping " << ctrs.size()
456  << " local faceCentres to " << localStr.name() << endl;
457 
458  forAll(ctrs, faceI)
459  {
460  writeOBJ(localStr, ctrs[faceI]);
461  }
462  }
463 
464  const bool isMaster = Pstream::myProcNo() < neighbProcNo();
465 
466  if (isMaster)
467  {
468  pointField ctrs(calcFaceCentres(pp, pp.points()));
469 
470  pointField anchors(getAnchorPoints(pp, pp.points()));
471 
472  // Now send all info over to the neighbour
473  OPstream toNeighbour(Pstream::blocking, neighbProcNo());
474  toNeighbour << ctrs << anchors;
475  }
476 }
477 
478 
479 // Return new ordering. Ordering is -faceMap: for every face index
480 // the new face -rotation:for every new face the clockwise shift
481 // of the original face. Return false if nothing changes (faceMap
482 // is identity, rotation is 0)
484 (
485  const primitivePatch& pp,
486  labelList& faceMap,
487  labelList& rotation
488 ) const
489 {
490  if (!Pstream::parRun())
491  {
492  return false;
493  }
494 
495  faceMap.setSize(pp.size());
496  faceMap = -1;
497 
498  rotation.setSize(pp.size());
499  rotation = 0;
500 
501  const bool isMaster = Pstream::myProcNo() < neighbProcNo();
502 
503  if (isMaster)
504  {
505  // Do nothing (i.e. identical mapping, zero rotation).
506  // See comment at top.
507  forAll(faceMap, patchFaceI)
508  {
509  faceMap[patchFaceI] = patchFaceI;
510  }
511 
512  return false;
513  }
514  else
515  {
516  vectorField masterCtrs;
517  vectorField masterAnchors;
518 
519  // Receive data from neighbour
520  {
521  IPstream fromNeighbour(Pstream::blocking, neighbProcNo());
522  fromNeighbour >> masterCtrs >> masterAnchors;
523  }
524 
525  // Calculate my face centres
526  pointField ctrs(calcFaceCentres(pp, pp.points()));
527 
528  // Calculate typical distance from face centre
529  scalarField tols(calcFaceTol(pp, pp.points(), ctrs));
530 
531  if (debug || masterCtrs.size() != pp.size())
532  {
533  {
534  OFstream nbrStr
535  (
536  boundaryMesh().mesh().time().path()
537  /name() + "_nbrFaceCentres.obj"
538  );
539  Pout<< "processorPolyPatch::order : "
540  << "Dumping neighbour faceCentres to " << nbrStr.name()
541  << endl;
542  forAll(masterCtrs, faceI)
543  {
544  writeOBJ(nbrStr, masterCtrs[faceI]);
545  }
546  }
547 
548  if (masterCtrs.size() != pp.size())
549  {
551  (
552  "processorPolyPatch::order(const primitivePatch&"
553  ", labelList&, labelList&) const"
554  ) << "in patch:" << name() << " : "
555  << "Local size of patch is " << pp.size() << " (faces)."
556  << endl
557  << "Received from neighbour " << masterCtrs.size()
558  << " faceCentres!"
559  << abort(FatalError);
560  }
561  }
562 
563  // Geometric match of face centre vectors
564  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
565 
566  // 1. Try existing ordering and transformation
567  bool matchedAll = false;
568 
569  if
570  (
571  separated()
572  && (separation().size() == 1 || separation().size() == pp.size())
573  )
574  {
575  vectorField transformedCtrs;
576 
577  const vectorField& v = separation();
578 
579  if (v.size() == 1)
580  {
581  transformedCtrs = masterCtrs-v[0];
582  }
583  else
584  {
585  transformedCtrs = masterCtrs-v;
586  }
587  matchedAll = matchPoints
588  (
589  ctrs,
590  transformedCtrs,
591  tols,
592  true,
593  faceMap
594  );
595 
596  if (matchedAll)
597  {
598  // Use transformed centers from now on
599  masterCtrs = transformedCtrs;
600 
601  // Transform anchors
602  if (v.size() == 1)
603  {
604  masterAnchors -= v[0];
605  }
606  else
607  {
608  masterAnchors -= v;
609  }
610  }
611  }
612  else if
613  (
614  !parallel()
615  && (forwardT().size() == 1 || forwardT().size() == pp.size())
616  )
617  {
618  vectorField transformedCtrs = masterCtrs;
619  transformList(forwardT(), transformedCtrs);
620  matchedAll = matchPoints
621  (
622  ctrs,
623  transformedCtrs,
624  tols,
625  true,
626  faceMap
627  );
628 
629  if (matchedAll)
630  {
631  // Use transformed centers from now on
632  masterCtrs = transformedCtrs;
633 
634  // Transform anchors
635  transformList(forwardT(), masterAnchors);
636  }
637  }
638 
639 
640  // 2. Try zero separation automatic matching
641  if (!matchedAll)
642  {
643  matchedAll = matchPoints(ctrs, masterCtrs, tols, true, faceMap);
644  }
645 
646  if (!matchedAll || debug)
647  {
648  // Dump faces
649  fileName str
650  (
651  boundaryMesh().mesh().time().path()
652  /name()/name()+"_faces.obj"
653  );
654  Pout<< "processorPolyPatch::order :"
655  << " Writing faces to OBJ file " << str.name() << endl;
656  writeOBJ(str, pp, pp.points());
657 
658  OFstream ccStr
659  (
660  boundaryMesh().mesh().time().path()
661  /name() + "_faceCentresConnections.obj"
662  );
663 
664  Pout<< "processorPolyPatch::order :"
665  << " Dumping newly found match as lines between"
666  << " corresponding face centres to OBJ file " << ccStr.name()
667  << endl;
668 
669  label vertI = 0;
670 
671  forAll(ctrs, faceI)
672  {
673  label masterFaceI = faceMap[faceI];
674 
675  if (masterFaceI != -1)
676  {
677  const point& c0 = masterCtrs[masterFaceI];
678  const point& c1 = ctrs[faceI];
679  writeOBJ(ccStr, c0, c1, vertI);
680  }
681  }
682  }
683 
684  if (!matchedAll)
685  {
687  (
688  "processorPolyPatch::order(const primitivePatch&"
689  ", labelList&, labelList&) const"
690  ) << "in patch:" << name() << " : "
691  << "Cannot match vectors to faces on both sides of patch"
692  << endl
693  << " masterCtrs[0]:" << masterCtrs[0] << endl
694  << " ctrs[0]:" << ctrs[0] << endl
695  << " Please check your topology changes or maybe you have"
696  << " multiple separated (from cyclics) processor patches"
697  << endl
698  << " Continuing with incorrect face ordering from now on!"
699  << endl;
700 
701  return false;
702  }
703 
704  // Set rotation.
705  forAll(faceMap, oldFaceI)
706  {
707  // The face f will be at newFaceI (after morphing) and we want its
708  // anchorPoint (= f[0]) to align with the anchorpoint for the
709  // corresponding face on the other side.
710 
711  label newFaceI = faceMap[oldFaceI];
712 
713  const point& wantedAnchor = masterAnchors[newFaceI];
714 
715  rotation[newFaceI] = getRotation
716  (
717  pp.points(),
718  pp[oldFaceI],
719  wantedAnchor,
720  tols[oldFaceI]
721  );
722 
723  if (rotation[newFaceI] == -1)
724  {
726  (
727  "processorPolyPatch::order(const primitivePatch&"
728  ", labelList&, labelList&) const"
729  ) << "in patch " << name()
730  << " : "
731  << "Cannot find point on face " << pp[oldFaceI]
732  << " with vertices "
733  << UIndirectList<point>(pp.points(), pp[oldFaceI])()
734  << " that matches point " << wantedAnchor
735  << " when matching the halves of processor patch " << name()
736  << "Continuing with incorrect face ordering from now on!"
737  << endl;
738 
739  return false;
740  }
741  }
742 
743  forAll(faceMap, faceI)
744  {
745  if (faceMap[faceI] != faceI || rotation[faceI] != 0)
746  {
747  return true;
748  }
749  }
750 
751  return false;
752  }
753 }
754 
755 
757 {
758  polyPatch::write(os);
759  os.writeKeyword("myProcNo") << myProcNo_
760  << token::END_STATEMENT << nl;
761  os.writeKeyword("neighbProcNo") << neighbProcNo_
762  << token::END_STATEMENT << nl;
763 }
764 
765 
766 // ************************ vim: set sw=4 sts=4 et: ************************ //