oldCyclicPolyPatch.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 "oldCyclicPolyPatch.H"
28 #include "polyBoundaryMesh.H"
29 #include "polyMesh.H"
30 #include "demandDrivenData.H"
31 #include "OFstream.H"
32 #include "patchZones.H"
33 #include "matchPoints.H"
34 #include "Time.H"
35 #include "transformList.H"
36 #include "cyclicPolyPatch.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(oldCyclicPolyPatch, 0);
43 
44  addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, word);
45  addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, dictionary);
46 }
47 
48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
49 
50 Foam::pointField Foam::oldCyclicPolyPatch::calcFaceCentres
51 (
52  const UList<face>& faces,
53  const pointField& points
54 )
55 {
56  pointField ctrs(faces.size());
57 
58  forAll(faces, faceI)
59  {
60  ctrs[faceI] = faces[faceI].centre(points);
61  }
62 
63  return ctrs;
64 }
65 
66 
67 Foam::pointField Foam::oldCyclicPolyPatch::getAnchorPoints
68 (
69  const UList<face>& faces,
70  const pointField& points
71 )
72 {
73  pointField anchors(faces.size());
74 
75  forAll(faces, faceI)
76  {
77  anchors[faceI] = points[faces[faceI][0]];
78  }
79 
80  return anchors;
81 }
82 
83 
84 Foam::label Foam::oldCyclicPolyPatch::findMaxArea
85 (
86  const pointField& points,
87  const faceList& faces
88 )
89 {
90  label maxI = -1;
91  scalar maxAreaSqr = -GREAT;
92 
93  forAll(faces, faceI)
94  {
95  scalar areaSqr = magSqr(faces[faceI].normal(points));
96 
97  if (areaSqr > maxAreaSqr)
98  {
99  maxAreaSqr = areaSqr;
100  maxI = faceI;
101  }
102  }
103  return maxI;
104 }
105 
106 
107 // Get geometric zones of patch by looking at normals.
108 // Method 1: any edge with sharpish angle is edge between two halves.
109 // (this will handle e.g. wedge geometries).
110 // Also two fully disconnected regions will be handled this way.
111 // Method 2: sort faces into two halves based on face normal.
112 bool Foam::oldCyclicPolyPatch::getGeometricHalves
113 (
114  const primitivePatch& pp,
115  labelList& half0ToPatch,
116  labelList& half1ToPatch
117 ) const
118 {
119  // Calculate normals
120  const vectorField& faceNormals = pp.faceNormals();
121 
122  // Find edges with sharp angles.
123  boolList regionEdge(pp.nEdges(), false);
124 
125  const labelListList& edgeFaces = pp.edgeFaces();
126 
127  label nRegionEdges = 0;
128 
129  forAll(edgeFaces, edgeI)
130  {
131  const labelList& eFaces = edgeFaces[edgeI];
132 
133  // Check manifold edges for sharp angle.
134  // (Non-manifold already handled by patchZones)
135  if (eFaces.size() == 2)
136  {
137  if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
138  {
139  regionEdge[edgeI] = true;
140 
141  nRegionEdges++;
142  }
143  }
144  }
145 
146 
147  // For every face determine zone it is connected to (without crossing
148  // any regionEdge)
149  patchZones ppZones(pp, regionEdge);
150 
151  if (debug)
152  {
153  Pout<< "oldCyclicPolyPatch::getGeometricHalves : "
154  << "Found " << nRegionEdges << " edges on patch " << name()
155  << " where the cos of the angle between two connected faces"
156  << " was less than " << featureCos_ << nl
157  << "Patch divided by these and by single sides edges into "
158  << ppZones.nZones() << " parts." << endl;
159 
160 
161  // Dumping zones to obj files.
162 
163  labelList nZoneFaces(ppZones.nZones());
164 
165  for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
166  {
167  OFstream stream
168  (
169  boundaryMesh().mesh().time().path()
170  /name()+"_zone_"+Foam::name(zoneI)+".obj"
171  );
172  Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing zone "
173  << zoneI << " face centres to OBJ file " << stream.name()
174  << endl;
175 
176  labelList zoneFaces(findIndices(ppZones, zoneI));
177 
178  forAll(zoneFaces, i)
179  {
180  writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
181  }
182 
183  nZoneFaces[zoneI] = zoneFaces.size();
184  }
185  }
186 
187 
188  if (ppZones.nZones() == 2)
189  {
190  half0ToPatch = findIndices(ppZones, 0);
191  half1ToPatch = findIndices(ppZones, 1);
192  }
193  else
194  {
195  if (debug)
196  {
197  Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
198  << " falling back to face-normal comparison" << endl;
199  }
200  label n0Faces = 0;
201  half0ToPatch.setSize(pp.size());
202 
203  label n1Faces = 0;
204  half1ToPatch.setSize(pp.size());
205 
206  // Compare to face 0 normal.
207  forAll(faceNormals, faceI)
208  {
209  if ((faceNormals[faceI] & faceNormals[0]) > 0)
210  {
211  half0ToPatch[n0Faces++] = faceI;
212  }
213  else
214  {
215  half1ToPatch[n1Faces++] = faceI;
216  }
217  }
218  half0ToPatch.setSize(n0Faces);
219  half1ToPatch.setSize(n1Faces);
220 
221  if (debug)
222  {
223  Pout<< "oldCyclicPolyPatch::getGeometricHalves :"
224  << " Number of faces per zone:("
225  << n0Faces << ' ' << n1Faces << ')' << endl;
226  }
227  }
228 
229  if (half0ToPatch.size() != half1ToPatch.size())
230  {
231  fileName casePath(boundaryMesh().mesh().time().path());
232 
233  // Dump halves
234  {
235  fileName nm0(casePath/name()+"_half0_faces.obj");
236  Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
237  << " faces to OBJ file " << nm0 << endl;
238  writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
239 
240  fileName nm1(casePath/name()+"_half1_faces.obj");
241  Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
242  << " faces to OBJ file " << nm1 << endl;
243  writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
244  }
245 
246  // Dump face centres
247  {
248  OFstream str0(casePath/name()+"_half0.obj");
249  Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0"
250  << " face centres to OBJ file " << str0.name() << endl;
251 
252  forAll(half0ToPatch, i)
253  {
254  writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
255  }
256 
257  OFstream str1(casePath/name()+"_half1.obj");
258  Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1"
259  << " face centres to OBJ file " << str1.name() << endl;
260  forAll(half1ToPatch, i)
261  {
262  writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
263  }
264  }
265 
267  (
268  "oldCyclicPolyPatch::getGeometricHalves"
269  "(const primitivePatch&, labelList&, labelList&) const"
270  ) << "Patch " << name() << " gets decomposed in two zones of"
271  << "inequal size: " << half0ToPatch.size()
272  << " and " << half1ToPatch.size() << endl
273  << "This means that the patch is either not two separate regions"
274  << " or one region where the angle between the different regions"
275  << " is not sufficiently sharp." << endl
276  << "Please adapt the featureCos setting." << endl
277  << "Continuing with incorrect face ordering from now on!" << endl;
278 
279  return false;
280  }
281  else
282  {
283  return true;
284  }
285 }
286 
287 
288 // Given a split of faces into left and right half calculate the centres
289 // and anchor points. Transform the left points so they align with the
290 // right ones.
291 void Foam::oldCyclicPolyPatch::getCentresAndAnchors
292 (
293  const primitivePatch& pp,
294  const faceList& half0Faces,
295  const faceList& half1Faces,
296 
297  pointField& ppPoints,
298  pointField& half0Ctrs,
299  pointField& half1Ctrs,
300  pointField& anchors0,
301  scalarField& tols
302 ) const
303 {
304  // Get geometric data on both halves.
305  half0Ctrs = calcFaceCentres(half0Faces, pp.points());
306  anchors0 = getAnchorPoints(half0Faces, pp.points());
307  half1Ctrs = calcFaceCentres(half1Faces, pp.points());
308 
309  switch (transform())
310  {
311  case ROTATIONAL:
312  {
313  label face0 = getConsistentRotationFace(half0Ctrs);
314  label face1 = getConsistentRotationFace(half1Ctrs);
315 
316  vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
317  vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
318  n0 /= mag(n0) + VSMALL;
319  n1 /= mag(n1) + VSMALL;
320 
321  if (debug)
322  {
323  Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
324  << " Specified rotation :"
325  << " n0:" << n0 << " n1:" << n1 << endl;
326  }
327 
328  // Rotation (around origin)
329  const tensor reverseT(rotationTensor(n0, -n1));
330 
331  // Rotation
332  forAll(half0Ctrs, faceI)
333  {
334  half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]);
335  anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]);
336  }
337 
338  ppPoints = Foam::transform(reverseT, pp.points());
339 
340  break;
341  }
342  //- Problem: usually specified translation is not accurate enough
343  //- To get proper match so keep automatic determination over here.
344  //case TRANSLATIONAL:
345  //{
346  // // Transform 0 points.
347  //
348  // if (debug)
349  // {
350  // Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
351  // << "Specified translation : " << separationVector_
352  // << endl;
353  // }
354  //
355  // half0Ctrs += separationVector_;
356  // anchors0 += separationVector_;
357  // break;
358  //}
359  default:
360  {
361  // Assumes that cyclic is planar. This is also the initial
362  // condition for patches without faces.
363 
364  // Determine the face with max area on both halves. These
365  // two faces are used to determine the transformation tensors
366  label max0I = findMaxArea(pp.points(), half0Faces);
367  vector n0 = half0Faces[max0I].normal(pp.points());
368  n0 /= mag(n0) + VSMALL;
369 
370  label max1I = findMaxArea(pp.points(), half1Faces);
371  vector n1 = half1Faces[max1I].normal(pp.points());
372  n1 /= mag(n1) + VSMALL;
373 
374  if (mag(n0 & n1) < 1-matchTolerance())
375  {
376  if (debug)
377  {
378  Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
379  << " Detected rotation :"
380  << " n0:" << n0 << " n1:" << n1 << endl;
381  }
382 
383  // Rotation (around origin)
384  const tensor reverseT(rotationTensor(n0, -n1));
385 
386  // Rotation
387  forAll(half0Ctrs, faceI)
388  {
389  half0Ctrs[faceI] = Foam::transform
390  (
391  reverseT,
392  half0Ctrs[faceI]
393  );
394  anchors0[faceI] = Foam::transform
395  (
396  reverseT,
397  anchors0[faceI]
398  );
399  }
400  ppPoints = Foam::transform(reverseT, pp.points());
401  }
402  else
403  {
404  // Parallel translation. Get average of all used points.
405 
406  primitiveFacePatch half0(half0Faces, pp.points());
407  const pointField& half0Pts = half0.localPoints();
408  const point ctr0(sum(half0Pts)/half0Pts.size());
409 
410  primitiveFacePatch half1(half1Faces, pp.points());
411  const pointField& half1Pts = half1.localPoints();
412  const point ctr1(sum(half1Pts)/half1Pts.size());
413 
414  if (debug)
415  {
416  Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :"
417  << " Detected translation :"
418  << " n0:" << n0 << " n1:" << n1
419  << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
420  }
421 
422  half0Ctrs += ctr1 - ctr0;
423  anchors0 += ctr1 - ctr0;
424  ppPoints = pp.points() + ctr1 - ctr0;
425  }
426  break;
427  }
428  }
429 
430 
431  // Calculate typical distance per face
432  tols = matchTolerance()*calcFaceTol(half1Faces, pp.points(), half1Ctrs);
433 }
434 
435 
436 // Calculates faceMap and rotation. Returns true if all ok.
437 bool Foam::oldCyclicPolyPatch::matchAnchors
438 (
439  const bool report,
440  const primitivePatch& pp,
441  const labelList& half0ToPatch,
442  const pointField& anchors0,
443 
444  const labelList& half1ToPatch,
445  const faceList& half1Faces,
446  const labelList& from1To0,
447 
448  const scalarField& tols,
449 
450  labelList& faceMap,
451  labelList& rotation
452 ) const
453 {
454  // Set faceMap such that half0 faces get first and corresponding half1
455  // faces last.
456 
457  forAll(half0ToPatch, half0FaceI)
458  {
459  // Label in original patch
460  label patchFaceI = half0ToPatch[half0FaceI];
461 
462  faceMap[patchFaceI] = half0FaceI;
463 
464  // No rotation
465  rotation[patchFaceI] = 0;
466  }
467 
468  bool fullMatch = true;
469 
470  forAll(from1To0, half1FaceI)
471  {
472  label patchFaceI = half1ToPatch[half1FaceI];
473 
474  // This face has to match the corresponding one on half0.
475  label half0FaceI = from1To0[half1FaceI];
476 
477  label newFaceI = half0FaceI + pp.size()/2;
478 
479  faceMap[patchFaceI] = newFaceI;
480 
481  // Rotate patchFaceI such that its f[0] aligns with that of
482  // the corresponding face
483  // (which after shuffling will be at position half0FaceI)
484 
485  const point& wantedAnchor = anchors0[half0FaceI];
486 
487  rotation[newFaceI] = getRotation
488  (
489  pp.points(),
490  half1Faces[half1FaceI],
491  wantedAnchor,
492  tols[half1FaceI]
493  );
494 
495  if (rotation[newFaceI] == -1)
496  {
497  fullMatch = false;
498 
499  if (report)
500  {
501  const face& f = half1Faces[half1FaceI];
503  (
504  "oldCyclicPolyPatch::matchAnchors(..)"
505  ) << "Patch:" << name() << " : "
506  << "Cannot find point on face " << f
507  << " with vertices:"
508  << UIndirectList<point>(pp.points(), f)()
509  << " that matches point " << wantedAnchor
510  << " when matching the halves of cyclic patch " << name()
511  << endl
512  << "Continuing with incorrect face ordering from now on!"
513  << endl;
514  }
515  }
516  }
517  return fullMatch;
518 }
519 
520 
521 Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace
522 (
523  const pointField& faceCentres
524 ) const
525 {
526  const scalarField magRadSqr
527  (
528  magSqr((faceCentres - rotationCentre_) ^ rotationAxis_)
529  );
530  scalarField axisLen
531  (
532  (faceCentres - rotationCentre_) & rotationAxis_
533  );
534  axisLen = axisLen - min(axisLen);
535  const scalarField magLenSqr
536  (
537  magRadSqr + axisLen*axisLen
538  );
539 
540  label rotFace = -1;
541  scalar maxMagLenSqr = -GREAT;
542  scalar maxMagRadSqr = -GREAT;
543  forAll(faceCentres, i)
544  {
545  if (magLenSqr[i] >= maxMagLenSqr)
546  {
547  if (magRadSqr[i] > maxMagRadSqr)
548  {
549  rotFace = i;
550  maxMagLenSqr = magLenSqr[i];
551  maxMagRadSqr = magRadSqr[i];
552  }
553  }
554  }
555 
556  if (debug)
557  {
558  Info<< "getConsistentRotationFace(const pointField&)" << nl
559  << " rotFace = " << rotFace << nl
560  << " point = " << faceCentres[rotFace] << endl;
561  }
562 
563  return rotFace;
564 }
565 
566 
567 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
568 
570 (
571  const word& name,
572  const label size,
573  const label start,
574  const label index,
575  const polyBoundaryMesh& bm,
576  const word& patchType,
577  const transformType transform
578 )
579 :
580  coupledPolyPatch(name, size, start, index, bm, patchType, transform),
581  featureCos_(0.9),
582  rotationAxis_(vector::zero),
583  rotationCentre_(point::zero),
584  separationVector_(vector::zero)
585 {}
586 
587 
589 (
590  const word& name,
591  const dictionary& dict,
592  const label index,
593  const polyBoundaryMesh& bm,
594  const word& patchType
595 )
596 :
597  coupledPolyPatch(name, dict, index, bm, patchType),
598  featureCos_(0.9),
599  rotationAxis_(vector::zero),
600  rotationCentre_(point::zero),
601  separationVector_(vector::zero)
602 {
603  if (dict.found("neighbourPatch"))
604  {
606  (
607  "oldCyclicPolyPatch::oldCyclicPolyPatch\n"
608  "(\n"
609  " const word& name,\n"
610  " const dictionary& dict,\n"
611  " const label index,\n"
612  " const polyBoundaryMesh& bm\n"
613  ")",
614  dict
615  ) << "Found \"neighbourPatch\" entry when reading cyclic patch "
616  << name << endl
617  << "Is this mesh already with split cyclics?" << endl
618  << "If so run a newer version that supports it"
619  << ", if not comment out the \"neighbourPatch\" entry and re-run"
620  << exit(FatalIOError);
621  }
622 
623  dict.readIfPresent("featureCos", featureCos_);
624 
625  switch (transform())
626  {
627  case ROTATIONAL:
628  {
629  dict.lookup("rotationAxis") >> rotationAxis_;
630  dict.lookup("rotationCentre") >> rotationCentre_;
631  break;
632  }
633  case TRANSLATIONAL:
634  {
635  dict.lookup("separationVector") >> separationVector_;
636  break;
637  }
638  default:
639  {
640  // no additional info required
641  }
642  }
643 }
644 
645 
647 (
648  const oldCyclicPolyPatch& pp,
649  const polyBoundaryMesh& bm
650 )
651 :
652  coupledPolyPatch(pp, bm),
653  featureCos_(pp.featureCos_),
654  rotationAxis_(pp.rotationAxis_),
655  rotationCentre_(pp.rotationCentre_),
656  separationVector_(pp.separationVector_)
657 {}
658 
659 
661 (
662  const oldCyclicPolyPatch& pp,
663  const polyBoundaryMesh& bm,
664  const label index,
665  const label newSize,
666  const label newStart
667 )
668 :
669  coupledPolyPatch(pp, bm, index, newSize, newStart),
670  featureCos_(pp.featureCos_),
671  rotationAxis_(pp.rotationAxis_),
672  rotationCentre_(pp.rotationCentre_),
673  separationVector_(pp.separationVector_)
674 {}
675 
676 
677 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
678 
680 {}
681 
682 
683 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
684 
686 {
688 }
689 
690 
692 (
693  const primitivePatch& referPatch,
694  const pointField& thisCtrs,
695  const vectorField& thisAreas,
696  const pointField& thisCc,
697  const pointField& nbrCtrs,
698  const vectorField& nbrAreas,
699  const pointField& nbrCc
700 )
701 {}
702 
703 
705 {}
706 
707 
709 (
710  PstreamBuffers& pBufs,
711  const pointField& p
712 )
713 {
714  polyPatch::initMovePoints(pBufs, p);
715 }
716 
717 
719 (
720  PstreamBuffers& pBufs,
721  const pointField& p
722 )
723 {
724  polyPatch::movePoints(pBufs, p);
725 }
726 
727 
729 {
731 }
732 
733 
735 {
736  polyPatch::updateMesh(pBufs);
737 }
738 
739 
741 (
743  const primitivePatch& pp
744 ) const
745 {}
746 
747 
748 // Return new ordering. Ordering is -faceMap: for every face index
749 // the new face -rotation:for every new face the clockwise shift
750 // of the original face. Return false if nothing changes (faceMap
751 // is identity, rotation is 0)
753 (
755  const primitivePatch& pp,
756  labelList& faceMap,
757  labelList& rotation
758 ) const
759 {
760  faceMap.setSize(pp.size());
761  faceMap = -1;
762 
763  rotation.setSize(pp.size());
764  rotation = 0;
765 
766  if (pp.empty())
767  {
768  // No faces, nothing to change.
769  return false;
770  }
771 
772  if (pp.size()&1)
773  {
774  FatalErrorIn("oldCyclicPolyPatch::order(..)")
775  << "Size of cyclic " << name() << " should be a multiple of 2"
776  << ". It is " << pp.size() << abort(FatalError);
777  }
778 
779  label halfSize = pp.size()/2;
780 
781  // Supplied primitivePatch already with new points.
782  // Cyclics are limited to one transformation tensor
783  // currently anyway (i.e. straight plane) so should not be too big a
784  // problem.
785 
786 
787  // Indices of faces on half0
788  labelList half0ToPatch;
789  // Indices of faces on half1
790  labelList half1ToPatch;
791 
792 
793  // 1. Test if already correctly oriented by starting from trivial ordering.
794  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
795 
796  half0ToPatch = identity(halfSize);
797  half1ToPatch = half0ToPatch + halfSize;
798 
799  // Get faces
800  faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
801  faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
802 
803  // Get geometric quantities
804  pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
805  scalarField tols;
806  getCentresAndAnchors
807  (
808  pp,
809  half0Faces,
810  half1Faces,
811 
812  ppPoints,
813  half0Ctrs,
814  half1Ctrs,
815  anchors0,
816  tols
817  );
818 
819  // Geometric match of face centre vectors
820  labelList from1To0;
821  bool matchedAll = matchPoints
822  (
823  half1Ctrs,
824  half0Ctrs,
825  tols,
826  false,
827  from1To0
828  );
829 
830  if (debug)
831  {
832  Pout<< "oldCyclicPolyPatch::order : test if already ordered:"
833  << matchedAll << endl;
834 
835  // Dump halves
836  fileName nm0("match1_"+name()+"_half0_faces.obj");
837  Pout<< "oldCyclicPolyPatch::order : Writing half0"
838  << " faces to OBJ file " << nm0 << endl;
839  writeOBJ(nm0, half0Faces, ppPoints);
840 
841  fileName nm1("match1_"+name()+"_half1_faces.obj");
842  Pout<< "oldCyclicPolyPatch::order : Writing half1"
843  << " faces to OBJ file " << nm1 << endl;
844  writeOBJ(nm1, half1Faces, pp.points());
845 
846  OFstream ccStr
847  (
848  boundaryMesh().mesh().time().path()
849  /"match1_"+ name() + "_faceCentres.obj"
850  );
851  Pout<< "oldCyclicPolyPatch::order : "
852  << "Dumping currently found cyclic match as lines between"
853  << " corresponding face centres to file " << ccStr.name()
854  << endl;
855 
856  // Recalculate untransformed face centres
857  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
858  label vertI = 0;
859 
860  forAll(half1Ctrs, i)
861  {
862  //if (from1To0[i] != -1)
863  {
864  // Write edge between c1 and c0
865  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
866  //const point& c0 = half0Ctrs[from1To0[i]];
867  const point& c0 = half0Ctrs[i];
868  const point& c1 = half1Ctrs[i];
869  writeOBJ(ccStr, c0, c1, vertI);
870  }
871  }
872  }
873 
874 
875  // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
876  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
877 
878  if (!matchedAll)
879  {
880  label faceI = 0;
881  for (label i = 0; i < halfSize; i++)
882  {
883  half0ToPatch[i] = faceI++;
884  half1ToPatch[i] = faceI++;
885  }
886 
887  // And redo all matching
888  half0Faces = UIndirectList<face>(pp, half0ToPatch);
889  half1Faces = UIndirectList<face>(pp, half1ToPatch);
890 
891  getCentresAndAnchors
892  (
893  pp,
894  half0Faces,
895  half1Faces,
896 
897  ppPoints,
898  half0Ctrs,
899  half1Ctrs,
900  anchors0,
901  tols
902  );
903 
904  // Geometric match of face centre vectors
905  matchedAll = matchPoints
906  (
907  half1Ctrs,
908  half0Ctrs,
909  tols,
910  false,
911  from1To0
912  );
913 
914  if (debug)
915  {
916  Pout<< "oldCyclicPolyPatch::order : test if pairwise ordered:"
917  << matchedAll << endl;
918  // Dump halves
919  fileName nm0("match2_"+name()+"_half0_faces.obj");
920  Pout<< "oldCyclicPolyPatch::order : Writing half0"
921  << " faces to OBJ file " << nm0 << endl;
922  writeOBJ(nm0, half0Faces, ppPoints);
923 
924  fileName nm1("match2_"+name()+"_half1_faces.obj");
925  Pout<< "oldCyclicPolyPatch::order : Writing half1"
926  << " faces to OBJ file " << nm1 << endl;
927  writeOBJ(nm1, half1Faces, pp.points());
928 
929  OFstream ccStr
930  (
931  boundaryMesh().mesh().time().path()
932  /"match2_"+name()+"_faceCentres.obj"
933  );
934  Pout<< "oldCyclicPolyPatch::order : "
935  << "Dumping currently found cyclic match as lines between"
936  << " corresponding face centres to file " << ccStr.name()
937  << endl;
938 
939  // Recalculate untransformed face centres
940  label vertI = 0;
941 
942  forAll(half1Ctrs, i)
943  {
944  if (from1To0[i] != -1)
945  {
946  // Write edge between c1 and c0
947  const point& c0 = half0Ctrs[from1To0[i]];
948  const point& c1 = half1Ctrs[i];
949  writeOBJ(ccStr, c0, c1, vertI);
950  }
951  }
952  }
953  }
954 
955 
956  // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
957  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
958 
959  if (!matchedAll)
960  {
961  label baffleI = 0;
962 
963  forAll(pp, faceI)
964  {
965  const face& f = pp.localFaces()[faceI];
966  const labelList& pFaces = pp.pointFaces()[f[0]];
967 
968  label matchedFaceI = -1;
969 
970  forAll(pFaces, i)
971  {
972  label otherFaceI = pFaces[i];
973 
974  if (otherFaceI > faceI)
975  {
976  const face& otherF = pp.localFaces()[otherFaceI];
977 
978  // Note: might pick up two similar oriented faces
979  // (but that is illegal anyway)
980  if (f == otherF)
981  {
982  matchedFaceI = otherFaceI;
983  break;
984  }
985  }
986  }
987 
988  if (matchedFaceI != -1)
989  {
990  half0ToPatch[baffleI] = faceI;
991  half1ToPatch[baffleI] = matchedFaceI;
992  baffleI++;
993  }
994  }
995 
996  if (baffleI == halfSize)
997  {
998  // And redo all matching
999  half0Faces = UIndirectList<face>(pp, half0ToPatch);
1000  half1Faces = UIndirectList<face>(pp, half1ToPatch);
1001 
1002  getCentresAndAnchors
1003  (
1004  pp,
1005  half0Faces,
1006  half1Faces,
1007 
1008  ppPoints,
1009  half0Ctrs,
1010  half1Ctrs,
1011  anchors0,
1012  tols
1013  );
1014 
1015  // Geometric match of face centre vectors
1016  matchedAll = matchPoints
1017  (
1018  half1Ctrs,
1019  half0Ctrs,
1020  tols,
1021  false,
1022  from1To0
1023  );
1024 
1025  if (debug)
1026  {
1027  Pout<< "oldCyclicPolyPatch::order : test if baffles:"
1028  << matchedAll << endl;
1029  // Dump halves
1030  fileName nm0("match3_"+name()+"_half0_faces.obj");
1031  Pout<< "oldCyclicPolyPatch::order : Writing half0"
1032  << " faces to OBJ file " << nm0 << endl;
1033  writeOBJ(nm0, half0Faces, ppPoints);
1034 
1035  fileName nm1("match3_"+name()+"_half1_faces.obj");
1036  Pout<< "oldCyclicPolyPatch::order : Writing half1"
1037  << " faces to OBJ file " << nm1 << endl;
1038  writeOBJ(nm1, half1Faces, pp.points());
1039 
1040  OFstream ccStr
1041  (
1042  boundaryMesh().mesh().time().path()
1043  /"match3_"+ name() + "_faceCentres.obj"
1044  );
1045  Pout<< "oldCyclicPolyPatch::order : "
1046  << "Dumping currently found cyclic match as lines between"
1047  << " corresponding face centres to file " << ccStr.name()
1048  << endl;
1049 
1050  // Recalculate untransformed face centres
1051  label vertI = 0;
1052 
1053  forAll(half1Ctrs, i)
1054  {
1055  if (from1To0[i] != -1)
1056  {
1057  // Write edge between c1 and c0
1058  const point& c0 = half0Ctrs[from1To0[i]];
1059  const point& c1 = half1Ctrs[i];
1060  writeOBJ(ccStr, c0, c1, vertI);
1061  }
1062  }
1063  }
1064  }
1065  }
1066 
1067 
1068  // 4. Automatic geometric ordering
1069  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1070 
1071  if (!matchedAll)
1072  {
1073  // Split faces according to feature angle or topology
1074  bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1075 
1076  if (!okSplit)
1077  {
1078  // Did not split into two equal parts.
1079  return false;
1080  }
1081 
1082  // And redo all matching
1083  half0Faces = UIndirectList<face>(pp, half0ToPatch);
1084  half1Faces = UIndirectList<face>(pp, half1ToPatch);
1085 
1086  getCentresAndAnchors
1087  (
1088  pp,
1089  half0Faces,
1090  half1Faces,
1091 
1092  ppPoints,
1093  half0Ctrs,
1094  half1Ctrs,
1095  anchors0,
1096  tols
1097  );
1098 
1099  // Geometric match of face centre vectors
1100  matchedAll = matchPoints
1101  (
1102  half1Ctrs,
1103  half0Ctrs,
1104  tols,
1105  false,
1106  from1To0
1107  );
1108 
1109  if (debug)
1110  {
1111  Pout<< "oldCyclicPolyPatch::order : automatic ordering result:"
1112  << matchedAll << endl;
1113  // Dump halves
1114  fileName nm0("match4_"+name()+"_half0_faces.obj");
1115  Pout<< "oldCyclicPolyPatch::order : Writing half0"
1116  << " faces to OBJ file " << nm0 << endl;
1117  writeOBJ(nm0, half0Faces, ppPoints);
1118 
1119  fileName nm1("match4_"+name()+"_half1_faces.obj");
1120  Pout<< "oldCyclicPolyPatch::order : Writing half1"
1121  << " faces to OBJ file " << nm1 << endl;
1122  writeOBJ(nm1, half1Faces, pp.points());
1123 
1124  OFstream ccStr
1125  (
1126  boundaryMesh().mesh().time().path()
1127  /"match4_"+ name() + "_faceCentres.obj"
1128  );
1129  Pout<< "oldCyclicPolyPatch::order : "
1130  << "Dumping currently found cyclic match as lines between"
1131  << " corresponding face centres to file " << ccStr.name()
1132  << endl;
1133 
1134  // Recalculate untransformed face centres
1135  label vertI = 0;
1136 
1137  forAll(half1Ctrs, i)
1138  {
1139  if (from1To0[i] != -1)
1140  {
1141  // Write edge between c1 and c0
1142  const point& c0 = half0Ctrs[from1To0[i]];
1143  const point& c1 = half1Ctrs[i];
1144  writeOBJ(ccStr, c0, c1, vertI);
1145  }
1146  }
1147  }
1148  }
1149 
1150 
1151  if (!matchedAll || debug)
1152  {
1153  // Dump halves
1154  fileName nm0(name()+"_half0_faces.obj");
1155  Pout<< "oldCyclicPolyPatch::order : Writing half0"
1156  << " faces to OBJ file " << nm0 << endl;
1157  writeOBJ(nm0, half0Faces, pp.points());
1158 
1159  fileName nm1(name()+"_half1_faces.obj");
1160  Pout<< "oldCyclicPolyPatch::order : Writing half1"
1161  << " faces to OBJ file " << nm1 << endl;
1162  writeOBJ(nm1, half1Faces, pp.points());
1163 
1164  OFstream ccStr
1165  (
1166  boundaryMesh().mesh().time().path()
1167  /name() + "_faceCentres.obj"
1168  );
1169  Pout<< "oldCyclicPolyPatch::order : "
1170  << "Dumping currently found cyclic match as lines between"
1171  << " corresponding face centres to file " << ccStr.name()
1172  << endl;
1173 
1174  // Recalculate untransformed face centres
1175  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1176  label vertI = 0;
1177 
1178  forAll(half1Ctrs, i)
1179  {
1180  if (from1To0[i] != -1)
1181  {
1182  // Write edge between c1 and c0
1183  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1184  const point& c0 = half0Ctrs[from1To0[i]];
1185  const point& c1 = half1Ctrs[i];
1186  writeOBJ(ccStr, c0, c1, vertI);
1187  }
1188  }
1189  }
1190 
1191 
1192  if (!matchedAll)
1193  {
1195  (
1196  "oldCyclicPolyPatch::order"
1197  "(const primitivePatch&, labelList&, labelList&) const"
1198  ) << "Patch:" << name() << " : "
1199  << "Cannot match vectors to faces on both sides of patch" << endl
1200  << " Perhaps your faces do not match?"
1201  << " The obj files written contain the current match." << endl
1202  << " Continuing with incorrect face ordering from now on!"
1203  << endl;
1204 
1205  return false;
1206  }
1207 
1208 
1209  // Set faceMap such that half0 faces get first and corresponding half1
1210  // faces last.
1211  matchAnchors
1212  (
1213  true, // report if anchor matching error
1214  pp,
1215  half0ToPatch,
1216  anchors0,
1217  half1ToPatch,
1218  half1Faces,
1219  from1To0,
1220  tols,
1221  faceMap,
1222  rotation
1223  );
1224 
1225  // Return false if no change neccesary, true otherwise.
1226 
1227  forAll(faceMap, faceI)
1228  {
1229  if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1230  {
1231  return true;
1232  }
1233  }
1234 
1235  return false;
1236 }
1237 
1238 
1240 {
1241  // Replacement of polyPatch::write to write 'cyclic' instead of type():
1242  os.writeKeyword("type") << cyclicPolyPatch::typeName
1243  << token::END_STATEMENT << nl;
1245  os.writeKeyword("nFaces") << size() << token::END_STATEMENT << nl;
1246  os.writeKeyword("startFace") << start() << token::END_STATEMENT << nl;
1247 
1248 
1249  os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
1250  switch (transform())
1251  {
1252  case ROTATIONAL:
1253  {
1254  os.writeKeyword("rotationAxis") << rotationAxis_
1255  << token::END_STATEMENT << nl;
1256  os.writeKeyword("rotationCentre") << rotationCentre_
1257  << token::END_STATEMENT << nl;
1258  break;
1259  }
1260  case TRANSLATIONAL:
1261  {
1262  os.writeKeyword("separationVector") << separationVector_
1263  << token::END_STATEMENT << nl;
1264  break;
1265  }
1266  default:
1267  {
1268  // no additional info to write
1269  }
1270  }
1271 
1272  WarningIn("oldCyclicPolyPatch::write(Ostream& os) const")
1273  << "Please run foamUpgradeCyclics to convert these old-style"
1274  << " cyclics into two separate cyclics patches."
1275  << endl;
1276 }
1277 
1278 
1279 // ************************************************************************* //
#define SeriousErrorIn(functionName)
Report an error message using Foam::SeriousError.
Output to file stream.
Definition: OFstream.H:81
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:306
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Buffers for inter-processor communications streams (UOPstream, UIPstream).
virtual bool order(PstreamBuffers &, const primitivePatch &, labelList &faceMap, labelList &rotation) const
Return new ordering for primitivePatch.
virtual void initOrder(PstreamBuffers &, const primitivePatch &) const
Initialize ordering for primitivePatch. Does not.
vector point
Point is a vector.
Definition: point.H:41
dimensioned< scalar > mag(const dimensioned< Type > &)
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
labelList f(nPoints)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const labelListList & pointFaces() const
Return point-face addressing.
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
virtual void initMovePoints(PstreamBuffers &, const pointField &)
Initialise the patches for moving points.
Definition: polyPatch.H:108
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
Definition: polyPatch.C:62
messageStream Info
List< face > faceList
Definition: faceListFwd.H:43
PrimitivePatch< face, List, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
dynamicFvMesh & mesh
virtual void calcGeometry(PstreamBuffers &)
Calculate the patch geometry.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
void write(Ostream &) const
Write patchIdentifier as a dictionary.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static const char nl
Definition: Ostream.H:260
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
IOerror FatalIOError
virtual void updateMesh(PstreamBuffers &)
Update of the patch topology.
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
Definition: polyPatch.H:100
A List with indirect addressing.
Definition: fvMatrix.H:106
Foam::polyBoundaryMesh.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
Field< vector > vectorField
Specialisation of Field<T> for vector.
#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].
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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.
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
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
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
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
virtual void initGeometry(PstreamBuffers &)
Initialise the calculation of the patch geometry.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
PrimitivePatch< face, SubList, const pointField & > primitivePatch
Addressing for a faceList slice.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
tensor rotationTensor(const vector &n1, const vector &n2)
Definition: transform.H:46
virtual void initUpdateMesh(PstreamBuffers &)
Initialise the update of the patch topology.
Definition: polyPatch.H:115
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
oldCyclicPolyPatch(const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm, const word &patchType, const transformType transform=UNKNOWN)
Construct from components.
#define FatalIOErrorIn(functionName, ios)
Report an error message using Foam::FatalIOError.
Definition: error.H:325
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurences of given element. Linear search.
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
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
virtual void movePoints(PstreamBuffers &, const pointField &p)
Correct patches after moving points.
Definition: polyPatch.C:57
&#39;old&#39; style cyclic polyPatch with all faces in single patch. Does ordering but cannot be used to run...
virtual void movePoints(PstreamBuffers &, const pointField &)
Correct patches after moving points.