FreeFOAM The Cross-Platform CFD Toolkit
cyclicPolyPatch.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-2011 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 "cyclicPolyPatch.H"
29 #include <OpenFOAM/polyMesh.H>
31 #include <OpenFOAM/OFstream.H>
32 #include <OpenFOAM/patchZones.H>
33 #include <OpenFOAM/matchPoints.H>
34 #include <OpenFOAM/EdgeMap.H>
35 #include <OpenFOAM/Time.H>
36 #include <OpenFOAM/transformList.H>
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(cyclicPolyPatch, 0);
43 
44  addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, word);
45  addToRunTimeSelectionTable(polyPatch, cyclicPolyPatch, dictionary);
46 }
47 
48 
49 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
50 
51 Foam::label Foam::cyclicPolyPatch::findMaxArea
52 (
53  const pointField& points,
54  const faceList& faces
55 )
56 {
57  label maxI = -1;
58  scalar maxAreaSqr = -GREAT;
59 
60  forAll(faces, faceI)
61  {
62  scalar areaSqr = magSqr(faces[faceI].normal(points));
63 
64  if (areaSqr > maxAreaSqr)
65  {
66  maxAreaSqr = areaSqr;
67  maxI = faceI;
68  }
69  }
70  return maxI;
71 }
72 
73 
74 void Foam::cyclicPolyPatch::calcTransforms()
75 {
76  if (size())
77  {
78  const pointField& points = this->points();
79 
80  // Determine geometric quantities on the two halves
81  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
82 
83  primitivePatch half0
84  (
85  SubList<face>
86  (
87  *this,
88  size()/2
89  ),
90  points
91  );
92 
93  pointField half0Ctrs(calcFaceCentres(half0, half0.points()));
94 
95  scalarField half0Tols(calcFaceTol(half0, half0.points(), half0Ctrs));
96 
97  primitivePatch half1
98  (
99  SubList<face>
100  (
101  *this,
102  size()/2,
103  size()/2
104  ),
105  points
106  );
107  pointField half1Ctrs(calcFaceCentres(half1, half1.points()));
108 
109  // Dump halves
110  if (debug)
111  {
112  fileName casePath(boundaryMesh().mesh().time().path());
113 
114  fileName nm0(casePath/name()+"_half0_faces.obj");
115  Pout<< "cyclicPolyPatch::calcTransforms : Writing half0"
116  << " faces to OBJ file " << nm0 << endl;
117  writeOBJ(nm0, half0, half0.points());
118 
119  fileName nm1(casePath/name()+"_half1_faces.obj");
120  Pout<< "cyclicPolyPatch::calcTransforms : Writing half1"
121  << " faces to OBJ file " << nm1 << endl;
122  writeOBJ(nm1, half1, half1.points());
123 
124  OFstream str(casePath/name()+"_half0_to_half1.obj");
125  label vertI = 0;
126  Pout<< "cyclicPolyPatch::calcTransforms :"
127  << " Writing coupled face centres as lines to " << str.name()
128  << endl;
129  forAll(half0Ctrs, i)
130  {
131  const point& p0 = half0Ctrs[i];
132  str << "v " << p0.x() << ' ' << p0.y() << ' ' << p0.z() << nl;
133  vertI++;
134  const point& p1 = half1Ctrs[i];
135  str << "v " << p1.x() << ' ' << p1.y() << ' ' << p1.z() << nl;
136  vertI++;
137  str << "l " << vertI-1 << ' ' << vertI << nl;
138  }
139  }
140 
141  vectorField half0Normals(half0.size());
142  vectorField half1Normals(half1.size());
143 
144  scalar maxAreaDiff = -GREAT;
145  label maxAreaFacei = -1;
146 
147  for (label facei = 0; facei < size()/2; facei++)
148  {
149  half0Normals[facei] = operator[](facei).normal(points);
150  label nbrFacei = facei+size()/2;
151  half1Normals[facei] = operator[](nbrFacei).normal(points);
152 
153  scalar magSf = mag(half0Normals[facei]);
154  scalar nbrMagSf = mag(half1Normals[facei]);
155  scalar avSf = (magSf + nbrMagSf)/2.0;
156 
157  if (magSf < ROOTVSMALL && nbrMagSf < ROOTVSMALL)
158  {
159  // Undetermined normal. Use dummy normal to force separation
160  // check. (note use of sqrt(VSMALL) since that is how mag
161  // scales)
162  half0Normals[facei] = point(1, 0, 0);
163  half1Normals[facei] = half0Normals[facei];
164  }
165  else
166  {
167  scalar areaDiff = mag(magSf - nbrMagSf)/avSf;
168 
169  if (areaDiff > maxAreaDiff)
170  {
171  maxAreaDiff = areaDiff;
172  maxAreaFacei = facei;
173  }
174 
175  if (areaDiff > coupledPolyPatch::matchTol)
176  {
178  (
179  "cyclicPolyPatch::calcTransforms()"
180  ) << "face " << facei << " area does not match neighbour "
181  << nbrFacei << " by "
182  << 100*areaDiff
183  << "% -- possible face ordering problem." << endl
184  << "patch:" << name()
185  << " my area:" << magSf
186  << " neighbour area:" << nbrMagSf
187  << " matching tolerance:"
189  << endl
190  << "Mesh face:" << start()+facei
191  << " vertices:"
192  << UIndirectList<point>(points, operator[](facei))()
193  << endl
194  << "Neighbour face:" << start()+nbrFacei
195  << " vertices:"
196  << UIndirectList<point>(points, operator[](nbrFacei))()
197  << endl
198  << "Rerun with cyclic debug flag set"
199  << " for more information." << exit(FatalError);
200  }
201  else
202  {
203  half0Normals[facei] /= magSf;
204  half1Normals[facei] /= nbrMagSf;
205  }
206  }
207  }
208 
209 
210  // Print area match
211  if (debug)
212  {
213  label nbrFacei = maxAreaFacei+size()/2;
214  Pout<< "cyclicPolyPatch::calcTransforms :"
215  << " Max area error:" << 100*maxAreaDiff << "% at face:"
216  << maxAreaFacei << " at:" << half0Ctrs[maxAreaFacei]
217  << " coupled face:" << nbrFacei
218  << " at:" << half1Ctrs[maxAreaFacei]
219  << endl;
220  }
221 
222 
223 
224  // See if transformation is prescribed
225  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
226 
227  switch (transform_)
228  {
229  case ROTATIONAL:
230  {
231  // Specified single rotation tensor.
232 
233  // Get best fitting face and its opposite number
234  label face0 = getConsistentRotationFace(half0Ctrs);
235  label face1 = face0;
236 
237  vector n0 =
238  (
239  (half0Ctrs[face0] - rotationCentre_)
240  ^ rotationAxis_
241  );
242  vector n1 =
243  (
244  (half1Ctrs[face1] - rotationCentre_)
245  ^ -rotationAxis_
246  );
247  n0 /= mag(n0) + VSMALL;
248  n1 /= mag(n1) + VSMALL;
249 
250  if (debug)
251  {
252  Pout<< "cyclicPolyPatch::calcTransforms :"
253  << " Specified rotation :"
254  << " n0:" << n0 << " from face " << half0Ctrs[face0]
255  << " and n1:" << n1 << " from face " << half1Ctrs[face1]
256  << endl;
257  }
258 
259  // Calculate transformation tensors from face0,1 only.
260  // Note: can use tight tolerance now.
261  calcTransformTensors
262  (
263  pointField(1, half0Ctrs[face0]),
264  pointField(1, half1Ctrs[face1]),
265  vectorField(1, n0),
266  vectorField(1, n1),
267  scalarField(1, half0Tols[face0]),
268  1E-4,
269  ROTATIONAL
270  );
271 
272  break;
273  }
274 
275  case TRANSLATIONAL:
276  {
277  // Calculate transformation tensors from all faces just to
278  // compare against user provided tolerance.
279  calcTransformTensors
280  (
281  half0Ctrs,
282  half1Ctrs,
283  half0Normals,
284  half1Normals,
285  half0Tols,
286  matchTol,
287  transform_
288  );
289 
290  if (debug)
291  {
292  Pout<< "cyclicPolyPatch::calcTransforms :"
293  << " Specified separation vector : "
294  << separationVector_
295  << " . Calculated average separation : "
297  << endl;
298  }
299 
300  // Override computed transform with specified.
301  const scalar avgTol = average(half0Tols);
302  if
303  (
304  coupledPolyPatch::separation().size() != 1
305  || mag(separation(0) - separationVector_) > avgTol
306  )
307  {
308  WarningIn
309  (
310  "cyclicPolyPatch::calcTransforms()"
311  ) << "Specified separationVector " << separationVector_
312  << " differs from computed separation vector "
313  << coupledPolyPatch::separation() << endl
314  << "This probably means your geometry is not consistent"
315  << " with the specified separation and might lead"
316  << " to problems." << endl
317  << "Continuing with specified separation vector "
318  << separationVector_ << endl
319  << "patch:" << name()
320  << endl;
321  }
322 
323  // Set transformation tensor.
324  const_cast<tensorField&>(coupledPolyPatch::forwardT()).clear();
325  const_cast<tensorField&>(coupledPolyPatch::reverseT()).clear();
326  const_cast<vectorField&>(coupledPolyPatch::separation()) =
328  (
329  1,
330  separationVector_
331  );
332 
333  break;
334  }
335 
336  default:
337  {
338  // Calculate transformation tensors from all faces.
339  calcTransformTensors
340  (
341  half0Ctrs,
342  half1Ctrs,
343  half0Normals,
344  half1Normals,
345  half0Tols,
346  matchTol,
347  transform_
348  );
349 
350  break;
351  }
352  }
353  }
354 }
355 
356 
357 // Get geometric zones of patch by looking at normals.
358 // Method 1: any edge with sharpish angle is edge between two halves.
359 // (this will handle e.g. wedge geometries).
360 // Also two fully disconnected regions will be handled this way.
361 // Method 2: sort faces into two halves based on face normal.
362 bool Foam::cyclicPolyPatch::getGeometricHalves
363 (
364  const primitivePatch& pp,
365  labelList& half0ToPatch,
366  labelList& half1ToPatch
367 ) const
368 {
369  // Calculate normals
370  const vectorField& faceNormals = pp.faceNormals();
371 
372  // Find edges with sharp angles.
373  boolList regionEdge(pp.nEdges(), false);
374 
375  const labelListList& edgeFaces = pp.edgeFaces();
376 
377  label nRegionEdges = 0;
378 
379  forAll(edgeFaces, edgeI)
380  {
381  const labelList& eFaces = edgeFaces[edgeI];
382 
383  // Check manifold edges for sharp angle.
384  // (Non-manifold already handled by patchZones)
385  if (eFaces.size() == 2)
386  {
387  if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_)
388  {
389  regionEdge[edgeI] = true;
390 
391  nRegionEdges++;
392  }
393  }
394  }
395 
396 
397  // For every face determine zone it is connected to (without crossing
398  // any regionEdge)
399  patchZones ppZones(pp, regionEdge);
400 
401  if (debug)
402  {
403  Pout<< "cyclicPolyPatch::getGeometricHalves : "
404  << "Found " << nRegionEdges << " edges on patch " << name()
405  << " where the cos of the angle between two connected faces"
406  << " was less than " << featureCos_ << nl
407  << "Patch divided by these and by single sides edges into "
408  << ppZones.nZones() << " parts." << endl;
409 
410 
411  // Dumping zones to obj files.
412 
413  labelList nZoneFaces(ppZones.nZones());
414 
415  for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++)
416  {
417  OFstream stream
418  (
419  boundaryMesh().mesh().time().path()
420  /name()+"_zone_"+Foam::name(zoneI)+".obj"
421  );
422  Pout<< "cyclicPolyPatch::getGeometricHalves : Writing zone "
423  << zoneI << " face centres to OBJ file " << stream.name()
424  << endl;
425 
426  labelList zoneFaces(findIndices(ppZones, zoneI));
427 
428  forAll(zoneFaces, i)
429  {
430  writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points()));
431  }
432 
433  nZoneFaces[zoneI] = zoneFaces.size();
434  }
435  }
436 
437 
438  if (ppZones.nZones() == 2)
439  {
440  half0ToPatch = findIndices(ppZones, 0);
441  half1ToPatch = findIndices(ppZones, 1);
442  }
443  else
444  {
445  if (debug)
446  {
447  Pout<< "cyclicPolyPatch::getGeometricHalves :"
448  << " falling back to face-normal comparison" << endl;
449  }
450  label n0Faces = 0;
451  half0ToPatch.setSize(pp.size());
452 
453  label n1Faces = 0;
454  half1ToPatch.setSize(pp.size());
455 
456  // Compare to face 0 normal.
457  forAll(faceNormals, faceI)
458  {
459  if ((faceNormals[faceI] & faceNormals[0]) > 0)
460  {
461  half0ToPatch[n0Faces++] = faceI;
462  }
463  else
464  {
465  half1ToPatch[n1Faces++] = faceI;
466  }
467  }
468  half0ToPatch.setSize(n0Faces);
469  half1ToPatch.setSize(n1Faces);
470 
471  if (debug)
472  {
473  Pout<< "cyclicPolyPatch::getGeometricHalves :"
474  << " Number of faces per zone:("
475  << n0Faces << ' ' << n1Faces << ')' << endl;
476  }
477  }
478 
479  if (half0ToPatch.size() != half1ToPatch.size())
480  {
481  fileName casePath(boundaryMesh().mesh().time().path());
482 
483  // Dump halves
484  {
485  fileName nm0(casePath/name()+"_half0_faces.obj");
486  Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
487  << " faces to OBJ file " << nm0 << endl;
488  writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points());
489 
490  fileName nm1(casePath/name()+"_half1_faces.obj");
491  Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
492  << " faces to OBJ file " << nm1 << endl;
493  writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points());
494  }
495 
496  // Dump face centres
497  {
498  OFstream str0(casePath/name()+"_half0.obj");
499  Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half0"
500  << " face centres to OBJ file " << str0.name() << endl;
501 
502  forAll(half0ToPatch, i)
503  {
504  writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points()));
505  }
506 
507  OFstream str1(casePath/name()+"_half1.obj");
508  Pout<< "cyclicPolyPatch::getGeometricHalves : Writing half1"
509  << " face centres to OBJ file " << str1.name() << endl;
510  forAll(half1ToPatch, i)
511  {
512  writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points()));
513  }
514  }
515 
517  (
518  "cyclicPolyPatch::getGeometricHalves"
519  "(const primitivePatch&, labelList&, labelList&) const"
520  ) << "Patch " << name() << " gets decomposed in two zones of"
521  << "inequal size: " << half0ToPatch.size()
522  << " and " << half1ToPatch.size() << endl
523  << "This means that the patch is either not two separate regions"
524  << " or one region where the angle between the different regions"
525  << " is not sufficiently sharp." << endl
526  << "Please adapt the featureCos setting." << endl
527  << "Continuing with incorrect face ordering from now on!" << endl;
528 
529  return false;
530  }
531  else
532  {
533  return true;
534  }
535 }
536 
537 
538 // Given a split of faces into left and right half calculate the centres
539 // and anchor points. Transform the left points so they align with the
540 // right ones.
541 void Foam::cyclicPolyPatch::getCentresAndAnchors
542 (
543  const primitivePatch& pp,
544  const faceList& half0Faces,
545  const faceList& half1Faces,
546 
547  pointField& ppPoints,
548  pointField& half0Ctrs,
549  pointField& half1Ctrs,
550  pointField& anchors0,
551  scalarField& tols
552 ) const
553 {
554  // Get geometric data on both halves.
555  half0Ctrs = calcFaceCentres(half0Faces, pp.points());
556  anchors0 = getAnchorPoints(half0Faces, pp.points());
557  half1Ctrs = calcFaceCentres(half1Faces, pp.points());
558 
559  switch (transform_)
560  {
561  case ROTATIONAL:
562  {
563  label face0 = getConsistentRotationFace(half0Ctrs);
564  label face1 = getConsistentRotationFace(half1Ctrs);
565 
566  vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_);
567  vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_);
568  n0 /= mag(n0) + VSMALL;
569  n1 /= mag(n1) + VSMALL;
570 
571  if (debug)
572  {
573  Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
574  << " Specified rotation :"
575  << " n0:" << n0 << " n1:" << n1 << endl;
576  }
577 
578  // Rotation (around origin)
579  const tensor reverseT(rotationTensor(n0, -n1));
580 
581  // Rotation
582  forAll(half0Ctrs, faceI)
583  {
584  half0Ctrs[faceI] =
586  (
587  reverseT,
588  half0Ctrs[faceI] - rotationCentre_
589  )
590  + rotationCentre_;
591  anchors0[faceI] =
593  (
594  reverseT,
595  anchors0[faceI] - rotationCentre_
596  )
597  + rotationCentre_;
598  }
599 
600  ppPoints =
602  (
603  reverseT,
604  (pp.points() - rotationCentre_)()
605  )
606  + rotationCentre_;
607 
608  break;
609  }
610 
611  case TRANSLATIONAL:
612  {
613  // Transform 0 points.
614 
615  if (debug)
616  {
617  Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
618  << "Specified translation : " << separationVector_
619  << endl;
620  }
621 
622  half0Ctrs += separationVector_;
623  anchors0 += separationVector_;
624  ppPoints = pp.points() + separationVector_;
625 
626  break;
627  }
628 
629  default:
630  {
631  // Assumes that cyclic is planar. This is also the initial
632  // condition for patches without faces.
633 
634  // Determine the face with max area on both halves. These
635  // two faces are used to determine the transformation tensors
636  label max0I = findMaxArea(pp.points(), half0Faces);
637  vector n0 = half0Faces[max0I].normal(pp.points());
638  n0 /= mag(n0) + VSMALL;
639 
640  label max1I = findMaxArea(pp.points(), half1Faces);
641  vector n1 = half1Faces[max1I].normal(pp.points());
642  n1 /= mag(n1) + VSMALL;
643 
644  if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol)
645  {
646  if (debug)
647  {
648  Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
649  << " Detected rotation :"
650  << " n0:" << n0 << " n1:" << n1 << endl;
651  }
652 
653  // Rotation (around origin)
654  const tensor reverseT(rotationTensor(n0, -n1));
655 
656  // Rotation
657  forAll(half0Ctrs, faceI)
658  {
659  half0Ctrs[faceI] = Foam::transform
660  (
661  reverseT,
662  half0Ctrs[faceI]
663  );
664  anchors0[faceI] = Foam::transform
665  (
666  reverseT,
667  anchors0[faceI]
668  );
669  }
670  ppPoints = Foam::transform(reverseT, pp.points());
671  }
672  else
673  {
674  // Parallel translation. Get average of all used points.
675 
676  primitiveFacePatch half0(half0Faces, pp.points());
677  const pointField& half0Pts = half0.localPoints();
678  const point ctr0(sum(half0Pts)/half0Pts.size());
679 
680  primitiveFacePatch half1(half1Faces, pp.points());
681  const pointField& half1Pts = half1.localPoints();
682  const point ctr1(sum(half1Pts)/half1Pts.size());
683 
684  if (debug)
685  {
686  Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
687  << " Detected translation :"
688  << " n0:" << n0 << " n1:" << n1
689  << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl;
690  }
691 
692  half0Ctrs += ctr1 - ctr0;
693  anchors0 += ctr1 - ctr0;
694  ppPoints = pp.points() + ctr1 - ctr0;
695  }
696  break;
697  }
698  }
699 
700 
701  // Calculate typical distance per face
702  tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs);
703 }
704 
705 
706 // Calculates faceMap and rotation. Returns true if all ok.
707 bool Foam::cyclicPolyPatch::matchAnchors
708 (
709  const bool report,
710  const primitivePatch& pp,
711  const labelList& half0ToPatch,
712  const pointField& anchors0,
713 
714  const labelList& half1ToPatch,
715  const faceList& half1Faces,
716  const labelList& from1To0,
717 
718  const scalarField& tols,
719 
720  labelList& faceMap,
721  labelList& rotation
722 ) const
723 {
724  // Set faceMap such that half0 faces get first and corresponding half1
725  // faces last.
726 
727  forAll(half0ToPatch, half0FaceI)
728  {
729  // Label in original patch
730  label patchFaceI = half0ToPatch[half0FaceI];
731 
732  faceMap[patchFaceI] = half0FaceI;
733 
734  // No rotation
735  rotation[patchFaceI] = 0;
736  }
737 
738  bool fullMatch = true;
739 
740  forAll(from1To0, half1FaceI)
741  {
742  label patchFaceI = half1ToPatch[half1FaceI];
743 
744  // This face has to match the corresponding one on half0.
745  label half0FaceI = from1To0[half1FaceI];
746 
747  label newFaceI = half0FaceI + pp.size()/2;
748 
749  faceMap[patchFaceI] = newFaceI;
750 
751  // Rotate patchFaceI such that its f[0] aligns with that of
752  // the corresponding face
753  // (which after shuffling will be at position half0FaceI)
754 
755  const point& wantedAnchor = anchors0[half0FaceI];
756 
757  rotation[newFaceI] = getRotation
758  (
759  pp.points(),
760  half1Faces[half1FaceI],
761  wantedAnchor,
762  tols[half1FaceI]
763  );
764 
765  if (rotation[newFaceI] == -1)
766  {
767  fullMatch = false;
768 
769  if (report)
770  {
771  const face& f = half1Faces[half1FaceI];
773  (
774  "cyclicPolyPatch::matchAnchors(..)"
775  ) << "Patch:" << name() << " : "
776  << "Cannot find point on face " << f
777  << " with vertices:"
778  << UIndirectList<point>(pp.points(), f)()
779  << " that matches point " << wantedAnchor
780  << " when matching the halves of cyclic patch " << name()
781  << endl
782  << "Continuing with incorrect face ordering from now on!"
783  << endl;
784  }
785  }
786  }
787  return fullMatch;
788 }
789 
790 
791 Foam::label Foam::cyclicPolyPatch::getConsistentRotationFace
792 (
793  const pointField& faceCentres
794 ) const
795 {
796  // Determine a face furthest away from the axis
797 
798  const scalarField magRadSqr =
799  magSqr((faceCentres - rotationCentre_) ^ rotationAxis_);
800 
801  label rotFace = findMax(magRadSqr);
802 
803  if (debug)
804  {
805  Info<< "getConsistentRotationFace(const pointField&)" << nl
806  << " rotFace = " << rotFace << nl
807  << " point = " << faceCentres[rotFace] << nl
808  << " distance = " << Foam::sqrt(magRadSqr[rotFace])
809  << endl;
810  }
811 
812  return rotFace;
813 }
814 
815 
816 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
817 
819 (
820  const word& name,
821  const label size,
822  const label start,
823  const label index,
824  const polyBoundaryMesh& bm
825 )
826 :
827  coupledPolyPatch(name, size, start, index, bm),
828  coupledPointsPtr_(NULL),
829  coupledEdgesPtr_(NULL),
830  featureCos_(0.9),
831  transform_(UNKNOWN),
832  rotationAxis_(vector::zero),
833  rotationCentre_(point::zero),
834  separationVector_(vector::zero)
835 {}
836 
837 
839 (
840  const word& name,
841  const dictionary& dict,
842  const label index,
843  const polyBoundaryMesh& bm
844 )
845 :
846  coupledPolyPatch(name, dict, index, bm),
847  coupledPointsPtr_(NULL),
848  coupledEdgesPtr_(NULL),
849  featureCos_(0.9),
850  transform_(UNKNOWN),
851  rotationAxis_(vector::zero),
852  rotationCentre_(point::zero),
853  separationVector_(vector::zero)
854 {
855  if (dict.found("neighbourPatch"))
856  {
858  (
859  "cyclicPolyPatch::cyclicPolyPatch\n"
860  "(\n"
861  " const word& name,\n"
862  " const dictionary& dict,\n"
863  " const label index,\n"
864  " const polyBoundaryMesh& bm\n"
865  ")",
866  dict
867  ) << "Found \"neighbourPatch\" entry when reading cyclic patch "
868  << name << endl
869  << "Is this mesh already with split cyclics?" << endl
870  << "If so run a newer version that supports it"
871  << ", if not comment out the \"neighbourPatch\" entry and re-run"
872  << exit(FatalIOError);
873  }
874 
875  dict.readIfPresent("featureCos", featureCos_);
876 
877  if (dict.found("transform"))
878  {
879  transform_ = transformTypeNames.read(dict.lookup("transform"));
880  switch (transform_)
881  {
882  case ROTATIONAL:
883  {
884  dict.lookup("rotationAxis") >> rotationAxis_;
885  dict.lookup("rotationCentre") >> rotationCentre_;
886  break;
887  }
888  case TRANSLATIONAL:
889  {
890  dict.lookup("separationVector") >> separationVector_;
891  break;
892  }
893  default:
894  {
895  // no additional info required
896  }
897  }
898  }
899 }
900 
901 
903 (
904  const cyclicPolyPatch& pp,
905  const polyBoundaryMesh& bm
906 )
907 :
908  coupledPolyPatch(pp, bm),
909  coupledPointsPtr_(NULL),
910  coupledEdgesPtr_(NULL),
911  featureCos_(pp.featureCos_),
912  transform_(pp.transform_),
913  rotationAxis_(pp.rotationAxis_),
914  rotationCentre_(pp.rotationCentre_),
915  separationVector_(pp.separationVector_)
916 {}
917 
918 
920 (
921  const cyclicPolyPatch& pp,
922  const polyBoundaryMesh& bm,
923  const label index,
924  const label newSize,
925  const label newStart
926 )
927 :
928  coupledPolyPatch(pp, bm, index, newSize, newStart),
929  coupledPointsPtr_(NULL),
930  coupledEdgesPtr_(NULL),
931  featureCos_(pp.featureCos_),
932  transform_(pp.transform_),
933  rotationAxis_(pp.rotationAxis_),
934  rotationCentre_(pp.rotationCentre_),
935  separationVector_(pp.separationVector_)
936 {}
937 
938 
939 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
940 
942 {
943  deleteDemandDrivenData(coupledPointsPtr_);
944  deleteDemandDrivenData(coupledEdgesPtr_);
945 }
946 
947 
948 
949 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
950 
952 {
954 }
955 
957 {
959  calcTransforms();
960 }
961 
963 {
965 }
966 
968 {
970  calcTransforms();
971 }
972 
974 {
976 }
977 
979 {
981  deleteDemandDrivenData(coupledPointsPtr_);
982  deleteDemandDrivenData(coupledEdgesPtr_);
983 }
984 
985 
987 {
988  if (!coupledPointsPtr_)
989  {
990  // Look at cyclic patch as two halves, A and B.
991  // Now all we know is that relative face index in halfA is same
992  // as coupled face in halfB and also that the 0th vertex
993  // corresponds.
994 
995  // From halfA point to halfB or -1.
996  labelList coupledPoint(nPoints(), -1);
997 
998  for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
999  {
1000  const face& fA = localFaces()[patchFaceA];
1001 
1002  forAll(fA, indexA)
1003  {
1004  label patchPointA = fA[indexA];
1005 
1006  if (coupledPoint[patchPointA] == -1)
1007  {
1008  const face& fB = localFaces()[patchFaceA + size()/2];
1009 
1010  label indexB = (fB.size() - indexA) % fB.size();
1011 
1012  // Filter out points on wedge axis
1013  if (patchPointA != fB[indexB])
1014  {
1015  coupledPoint[patchPointA] = fB[indexB];
1016  }
1017  }
1018  }
1019  }
1020 
1021  coupledPointsPtr_ = new edgeList(nPoints());
1022  edgeList& connected = *coupledPointsPtr_;
1023 
1024  // Extract coupled points.
1025  label connectedI = 0;
1026 
1027  forAll(coupledPoint, i)
1028  {
1029  if (coupledPoint[i] != -1)
1030  {
1031  connected[connectedI++] = edge(i, coupledPoint[i]);
1032  }
1033  }
1034 
1035  connected.setSize(connectedI);
1036 
1037  if (debug)
1038  {
1039  OFstream str
1040  (
1041  boundaryMesh().mesh().time().path()
1042  /"coupledPoints.obj"
1043  );
1044  label vertI = 0;
1045 
1046  Pout<< "Writing file " << str.name() << " with coordinates of "
1047  << "coupled points" << endl;
1048 
1049  forAll(connected, i)
1050  {
1051  const point& a = points()[meshPoints()[connected[i][0]]];
1052  const point& b = points()[meshPoints()[connected[i][1]]];
1053 
1054  str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
1055  str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
1056  vertI += 2;
1057 
1058  str<< "l " << vertI-1 << ' ' << vertI << nl;
1059  }
1060  }
1061 
1062  // Remove any addressing calculated for the coupled edges calculation
1063  const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
1064  .clearOut();
1065  }
1066  return *coupledPointsPtr_;
1067 }
1068 
1069 
1071 {
1072  if (!coupledEdgesPtr_)
1073  {
1074  // Build map from points on halfA to points on halfB.
1075  const edgeList& pointCouples = coupledPoints();
1076 
1077  Map<label> aToB(2*pointCouples.size());
1078 
1079  forAll(pointCouples, i)
1080  {
1081  const edge& e = pointCouples[i];
1082 
1083  aToB.insert(e[0], e[1]);
1084  }
1085 
1086  // Map from edge on half A to points (in halfB indices)
1087  EdgeMap<label> edgeMap(nEdges());
1088 
1089  for (label patchFaceA = 0; patchFaceA < size()/2; patchFaceA++)
1090  {
1091  const labelList& fEdges = faceEdges()[patchFaceA];
1092 
1093  forAll(fEdges, i)
1094  {
1095  label edgeI = fEdges[i];
1096 
1097  const edge& e = edges()[edgeI];
1098 
1099  // Convert edge end points to corresponding points on halfB
1100  // side.
1101  Map<label>::const_iterator fnd0 = aToB.find(e[0]);
1102  if (fnd0 != aToB.end())
1103  {
1104  Map<label>::const_iterator fnd1 = aToB.find(e[1]);
1105  if (fnd1 != aToB.end())
1106  {
1107  edgeMap.insert(edge(fnd0(), fnd1()), edgeI);
1108  }
1109  }
1110  }
1111  }
1112 
1113  coupledEdgesPtr_ = new edgeList(nEdges()/2);
1114  edgeList& coupledEdges = *coupledEdgesPtr_;
1115  label coupleI = 0;
1116 
1117  for (label patchFaceB = size()/2; patchFaceB < size(); patchFaceB++)
1118  {
1119  const labelList& fEdges = faceEdges()[patchFaceB];
1120 
1121  forAll(fEdges, i)
1122  {
1123  label edgeI = fEdges[i];
1124 
1125  const edge& e = edges()[edgeI];
1126 
1127  // Look up halfA edge from HashTable.
1128  EdgeMap<label>::iterator iter = edgeMap.find(e);
1129 
1130  if (iter != edgeMap.end())
1131  {
1132  label halfAEdgeI = iter();
1133 
1134  // Store correspondence. Filter out edges on wedge axis.
1135  if (halfAEdgeI != edgeI)
1136  {
1137  coupledEdges[coupleI++] = edge(halfAEdgeI, edgeI);
1138  }
1139 
1140  // Remove so we build unique list
1141  edgeMap.erase(iter);
1142  }
1143  }
1144  }
1145  coupledEdges.setSize(coupleI);
1146 
1147 
1148  // Some checks
1149 
1150  forAll(coupledEdges, i)
1151  {
1152  const edge& e = coupledEdges[i];
1153 
1154  if (e[0] == e[1] || e[0] < 0 || e[1] < 0)
1155  {
1156  FatalErrorIn("cyclicPolyPatch::coupledEdges() const")
1157  << "Problem : at position " << i
1158  << " illegal couple:" << e
1159  << abort(FatalError);
1160  }
1161  }
1162 
1163  if (debug)
1164  {
1165  OFstream str
1166  (
1167  boundaryMesh().mesh().time().path()
1168  /"coupledEdges.obj"
1169  );
1170  label vertI = 0;
1171 
1172  Pout<< "Writing file " << str.name() << " with centres of "
1173  << "coupled edges" << endl;
1174 
1175  forAll(coupledEdges, i)
1176  {
1177  const edge& e = coupledEdges[i];
1178 
1179  const point& a = edges()[e[0]].centre(localPoints());
1180  const point& b = edges()[e[1]].centre(localPoints());
1181 
1182  str<< "v " << a.x() << ' ' << a.y() << ' ' << a.z() << nl;
1183  str<< "v " << b.x() << ' ' << b.y() << ' ' << b.z() << nl;
1184  vertI += 2;
1185 
1186  str<< "l " << vertI-1 << ' ' << vertI << nl;
1187  }
1188  }
1189 
1190  // Remove any addressing calculated for the coupled edges calculation
1191  const_cast<primitivePatch&>(static_cast<const primitivePatch&>(*this))
1192  .clearOut();
1193  }
1194  return *coupledEdgesPtr_;
1195 }
1196 
1197 
1199 {}
1200 
1201 
1202 // Return new ordering. Ordering is -faceMap: for every face index
1203 // the new face -rotation:for every new face the clockwise shift
1204 // of the original face. Return false if nothing changes (faceMap
1205 // is identity, rotation is 0)
1208  const primitivePatch& pp,
1209  labelList& faceMap,
1210  labelList& rotation
1211 ) const
1212 {
1213  faceMap.setSize(pp.size());
1214  faceMap = -1;
1215 
1216  rotation.setSize(pp.size());
1217  rotation = 0;
1218 
1219  if (pp.empty())
1220  {
1221  // No faces, nothing to change.
1222  return false;
1223  }
1224 
1225  if (pp.size()&1)
1226  {
1227  FatalErrorIn("cyclicPolyPatch::order(..)")
1228  << "Size of cyclic " << name() << " should be a multiple of 2"
1229  << ". It is " << pp.size() << abort(FatalError);
1230  }
1231 
1232  label halfSize = pp.size()/2;
1233 
1234  // Supplied primitivePatch already with new points.
1235  // Cyclics are limited to one transformation tensor
1236  // currently anyway (i.e. straight plane) so should not be too big a
1237  // problem.
1238 
1239 
1240  // Indices of faces on half0
1241  labelList half0ToPatch;
1242  // Indices of faces on half1
1243  labelList half1ToPatch;
1244 
1245 
1246  // 1. Test if already correctly oriented by starting from trivial ordering.
1247  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1248 
1249  half0ToPatch = identity(halfSize);
1250  half1ToPatch = half0ToPatch + halfSize;
1251 
1252  // Get faces
1253  faceList half0Faces(UIndirectList<face>(pp, half0ToPatch));
1254  faceList half1Faces(UIndirectList<face>(pp, half1ToPatch));
1255 
1256  // Get geometric quantities
1257  pointField half0Ctrs, half1Ctrs, anchors0, ppPoints;
1258  scalarField tols;
1259  getCentresAndAnchors
1260  (
1261  pp,
1262  half0Faces,
1263  half1Faces,
1264 
1265  ppPoints,
1266  half0Ctrs,
1267  half1Ctrs,
1268  anchors0,
1269  tols
1270  );
1271 
1272  if (debug)
1273  {
1274  Pout<< "half0 transformed faceCentres (avg) : "
1275  << gAverage(half0Ctrs) << nl
1276  << "half1 untransformed faceCentres (avg) : "
1277  << gAverage(half1Ctrs) << endl;
1278  }
1279 
1280 
1281  // Geometric match of face centre vectors
1282  labelList from1To0;
1283  bool matchedAll = matchPoints
1284  (
1285  half1Ctrs,
1286  half0Ctrs,
1287  tols,
1288  false,
1289  from1To0
1290  );
1291 
1292  if (debug)
1293  {
1294  Pout<< "cyclicPolyPatch::order : test if already ordered:"
1295  << matchedAll << endl;
1296 
1297  // Dump halves
1298  fileName nm0("match1_"+name()+"_half0_faces.obj");
1299  Pout<< "cyclicPolyPatch::order : Writing half0"
1300  << " faces to OBJ file " << nm0 << endl;
1301  writeOBJ(nm0, half0Faces, ppPoints);
1302 
1303  fileName nm1("match1_"+name()+"_half1_faces.obj");
1304  Pout<< "cyclicPolyPatch::order : Writing half1"
1305  << " faces to OBJ file " << nm1 << endl;
1306  writeOBJ(nm1, half1Faces, pp.points());
1307 
1308  if (matchedAll)
1309  {
1310  OFstream ccStr
1311  (
1312  boundaryMesh().mesh().time().path()
1313  /"match1_"+ name() + "_faceCentres.obj"
1314  );
1315  Pout<< "cyclicPolyPatch::order : "
1316  << "Dumping currently found cyclic match as lines between"
1317  << " corresponding face centres to file " << ccStr.name()
1318  << endl;
1319 
1320  label vertI = 0;
1321  forAll(half1Ctrs, i)
1322  {
1323  // Write edge between c1 and c0
1324  const point& c0 = half0Ctrs[from1To0[i]];
1325  const point& c1 = half1Ctrs[i];
1326  writeOBJ(ccStr, c0, c1, vertI);
1327  }
1328  }
1329  }
1330 
1331 
1332  // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.)
1333  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1334 
1335  if (!matchedAll)
1336  {
1337  label faceI = 0;
1338  for (label i = 0; i < halfSize; i++)
1339  {
1340  half0ToPatch[i] = faceI++;
1341  half1ToPatch[i] = faceI++;
1342  }
1343 
1344  // And redo all matching
1345  half0Faces = UIndirectList<face>(pp, half0ToPatch);
1346  half1Faces = UIndirectList<face>(pp, half1ToPatch);
1347 
1348  getCentresAndAnchors
1349  (
1350  pp,
1351  half0Faces,
1352  half1Faces,
1353 
1354  ppPoints,
1355  half0Ctrs,
1356  half1Ctrs,
1357  anchors0,
1358  tols
1359  );
1360 
1361  // Geometric match of face centre vectors
1362  matchedAll = matchPoints
1363  (
1364  half1Ctrs,
1365  half0Ctrs,
1366  tols,
1367  false,
1368  from1To0
1369  );
1370 
1371  if (debug)
1372  {
1373  Pout<< "cyclicPolyPatch::order : test if pairwise ordered:"
1374  << matchedAll << endl;
1375  // Dump halves
1376  fileName nm0("match2_"+name()+"_half0_faces.obj");
1377  Pout<< "cyclicPolyPatch::order : Writing half0"
1378  << " faces to OBJ file " << nm0 << endl;
1379  writeOBJ(nm0, half0Faces, ppPoints);
1380 
1381  fileName nm1("match2_"+name()+"_half1_faces.obj");
1382  Pout<< "cyclicPolyPatch::order : Writing half1"
1383  << " faces to OBJ file " << nm1 << endl;
1384  writeOBJ(nm1, half1Faces, pp.points());
1385 
1386  OFstream ccStr
1387  (
1388  boundaryMesh().mesh().time().path()
1389  /"match2_"+name()+"_faceCentres.obj"
1390  );
1391  Pout<< "cyclicPolyPatch::order : "
1392  << "Dumping currently found cyclic match as lines between"
1393  << " corresponding face centres to file " << ccStr.name()
1394  << endl;
1395 
1396  // Recalculate untransformed face centres
1397  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1398  label vertI = 0;
1399 
1400  forAll(half1Ctrs, i)
1401  {
1402  if (from1To0[i] != -1)
1403  {
1404  // Write edge between c1 and c0
1405  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1406  const point& c0 = half0Ctrs[from1To0[i]];
1407  const point& c1 = half1Ctrs[i];
1408  writeOBJ(ccStr, c0, c1, vertI);
1409  }
1410  }
1411  }
1412  }
1413 
1414 
1415  // 3. Baffles(coincident faces) converted into cyclics (e.g. jump)
1416  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1417 
1418  if (!matchedAll)
1419  {
1420  label baffleI = 0;
1421 
1422  forAll(pp, faceI)
1423  {
1424  const face& f = pp.localFaces()[faceI];
1425  const labelList& pFaces = pp.pointFaces()[f[0]];
1426 
1427  label matchedFaceI = -1;
1428 
1429  forAll(pFaces, i)
1430  {
1431  label otherFaceI = pFaces[i];
1432 
1433  if (otherFaceI > faceI)
1434  {
1435  const face& otherF = pp.localFaces()[otherFaceI];
1436 
1437  // Note: might pick up two similar oriented faces
1438  // (but that is illegal anyway)
1439  if (f == otherF)
1440  {
1441  matchedFaceI = otherFaceI;
1442  break;
1443  }
1444  }
1445  }
1446 
1447  if (matchedFaceI != -1)
1448  {
1449  half0ToPatch[baffleI] = faceI;
1450  half1ToPatch[baffleI] = matchedFaceI;
1451  baffleI++;
1452  }
1453  }
1454 
1455  if (baffleI == halfSize)
1456  {
1457  // And redo all matching
1458  half0Faces = UIndirectList<face>(pp, half0ToPatch);
1459  half1Faces = UIndirectList<face>(pp, half1ToPatch);
1460 
1461  getCentresAndAnchors
1462  (
1463  pp,
1464  half0Faces,
1465  half1Faces,
1466 
1467  ppPoints,
1468  half0Ctrs,
1469  half1Ctrs,
1470  anchors0,
1471  tols
1472  );
1473 
1474  // Geometric match of face centre vectors
1475  matchedAll = matchPoints
1476  (
1477  half1Ctrs,
1478  half0Ctrs,
1479  tols,
1480  false,
1481  from1To0
1482  );
1483 
1484  if (debug)
1485  {
1486  Pout<< "cyclicPolyPatch::order : test if baffles:"
1487  << matchedAll << endl;
1488  // Dump halves
1489  fileName nm0("match3_"+name()+"_half0_faces.obj");
1490  Pout<< "cyclicPolyPatch::order : Writing half0"
1491  << " faces to OBJ file " << nm0 << endl;
1492  writeOBJ(nm0, half0Faces, ppPoints);
1493 
1494  fileName nm1("match3_"+name()+"_half1_faces.obj");
1495  Pout<< "cyclicPolyPatch::order : Writing half1"
1496  << " faces to OBJ file " << nm1 << endl;
1497  writeOBJ(nm1, half1Faces, pp.points());
1498 
1499  OFstream ccStr
1500  (
1501  boundaryMesh().mesh().time().path()
1502  /"match3_"+ name() + "_faceCentres.obj"
1503  );
1504  Pout<< "cyclicPolyPatch::order : "
1505  << "Dumping currently found cyclic match as lines between"
1506  << " corresponding face centres to file " << ccStr.name()
1507  << endl;
1508 
1509  // Recalculate untransformed face centres
1510  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1511  label vertI = 0;
1512 
1513  forAll(half1Ctrs, i)
1514  {
1515  if (from1To0[i] != -1)
1516  {
1517  // Write edge between c1 and c0
1518  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1519  const point& c0 = half0Ctrs[from1To0[i]];
1520  const point& c1 = half1Ctrs[i];
1521  writeOBJ(ccStr, c0, c1, vertI);
1522  }
1523  }
1524  }
1525  }
1526  }
1527 
1528 
1529  // 4. Automatic geometric ordering
1530  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1531 
1532  if (!matchedAll)
1533  {
1534  // Split faces according to feature angle or topology
1535  bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch);
1536 
1537  if (!okSplit)
1538  {
1539  // Did not split into two equal parts.
1540  return false;
1541  }
1542 
1543  // And redo all matching
1544  half0Faces = UIndirectList<face>(pp, half0ToPatch);
1545  half1Faces = UIndirectList<face>(pp, half1ToPatch);
1546 
1547  getCentresAndAnchors
1548  (
1549  pp,
1550  half0Faces,
1551  half1Faces,
1552 
1553  ppPoints,
1554  half0Ctrs,
1555  half1Ctrs,
1556  anchors0,
1557  tols
1558  );
1559 
1560  // Geometric match of face centre vectors
1561  matchedAll = matchPoints
1562  (
1563  half1Ctrs,
1564  half0Ctrs,
1565  tols,
1566  false,
1567  from1To0
1568  );
1569 
1570  if (debug)
1571  {
1572  Pout<< "cyclicPolyPatch::order : automatic ordering result:"
1573  << matchedAll << endl;
1574  // Dump halves
1575  fileName nm0("match4_"+name()+"_half0_faces.obj");
1576  Pout<< "cyclicPolyPatch::order : Writing half0"
1577  << " faces to OBJ file " << nm0 << endl;
1578  writeOBJ(nm0, half0Faces, ppPoints);
1579 
1580  fileName nm1("match4_"+name()+"_half1_faces.obj");
1581  Pout<< "cyclicPolyPatch::order : Writing half1"
1582  << " faces to OBJ file " << nm1 << endl;
1583  writeOBJ(nm1, half1Faces, pp.points());
1584 
1585  OFstream ccStr
1586  (
1587  boundaryMesh().mesh().time().path()
1588  /"match4_"+ name() + "_faceCentres.obj"
1589  );
1590  Pout<< "cyclicPolyPatch::order : "
1591  << "Dumping currently found cyclic match as lines between"
1592  << " corresponding face centres to file " << ccStr.name()
1593  << endl;
1594 
1595  // Recalculate untransformed face centres
1596  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1597  label vertI = 0;
1598 
1599  forAll(half1Ctrs, i)
1600  {
1601  if (from1To0[i] != -1)
1602  {
1603  // Write edge between c1 and c0
1604  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1605  const point& c0 = half0Ctrs[from1To0[i]];
1606  const point& c1 = half1Ctrs[i];
1607  writeOBJ(ccStr, c0, c1, vertI);
1608  }
1609  }
1610  }
1611  }
1612 
1613 
1614  if (!matchedAll || debug)
1615  {
1616  // Dump halves
1617  fileName nm0(name()+"_half0_faces.obj");
1618  Pout<< "cyclicPolyPatch::order : Writing half0"
1619  << " faces to OBJ file " << nm0 << endl;
1620  writeOBJ(nm0, half0Faces, pp.points());
1621 
1622  fileName nm1(name()+"_half1_faces.obj");
1623  Pout<< "cyclicPolyPatch::order : Writing half1"
1624  << " faces to OBJ file " << nm1 << endl;
1625  writeOBJ(nm1, half1Faces, pp.points());
1626 
1627  OFstream ccStr
1628  (
1629  boundaryMesh().mesh().time().path()
1630  /name() + "_faceCentres.obj"
1631  );
1632  Pout<< "cyclicPolyPatch::order : "
1633  << "Dumping currently found cyclic match as lines between"
1634  << " corresponding face centres to file " << ccStr.name()
1635  << endl;
1636 
1637  // Recalculate untransformed face centres
1638  //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
1639  label vertI = 0;
1640 
1641  forAll(half1Ctrs, i)
1642  {
1643  if (from1To0[i] != -1)
1644  {
1645  // Write edge between c1 and c0
1646  //const point& c0 = rawHalf0Ctrs[from1To0[i]];
1647  const point& c0 = half0Ctrs[from1To0[i]];
1648  const point& c1 = half1Ctrs[i];
1649  writeOBJ(ccStr, c0, c1, vertI);
1650  }
1651  }
1652  }
1653 
1654 
1655  if (!matchedAll)
1656  {
1658  (
1659  "cyclicPolyPatch::order"
1660  "(const primitivePatch&, labelList&, labelList&) const"
1661  ) << "Patch:" << name() << " : "
1662  << "Cannot match vectors to faces on both sides of patch" << endl
1663  << " Perhaps your faces do not match?"
1664  << " The obj files written contain the current match." << endl
1665  << " Continuing with incorrect face ordering from now on!"
1666  << endl;
1667 
1668  return false;
1669  }
1670 
1671 
1672  // Set faceMap such that half0 faces get first and corresponding half1
1673  // faces last.
1674  matchAnchors
1675  (
1676  true, // report if anchor matching error
1677  pp,
1678  half0ToPatch,
1679  anchors0,
1680  half1ToPatch,
1681  half1Faces,
1682  from1To0,
1683  tols,
1684  faceMap,
1685  rotation
1686  );
1687 
1688  // Return false if no change necessary, true otherwise.
1689 
1690  forAll(faceMap, faceI)
1691  {
1692  if (faceMap[faceI] != faceI || rotation[faceI] != 0)
1693  {
1694  return true;
1695  }
1696  }
1697 
1698  return false;
1699 }
1700 
1701 
1703 {
1704  polyPatch::write(os);
1705  os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl;
1706  switch (transform_)
1707  {
1708  case ROTATIONAL:
1709  {
1710  os.writeKeyword("transform") << transformTypeNames[ROTATIONAL]
1711  << token::END_STATEMENT << nl;
1712  os.writeKeyword("rotationAxis") << rotationAxis_
1713  << token::END_STATEMENT << nl;
1714  os.writeKeyword("rotationCentre") << rotationCentre_
1715  << token::END_STATEMENT << nl;
1716  break;
1717  }
1718  case TRANSLATIONAL:
1719  {
1720  os.writeKeyword("transform") << transformTypeNames[TRANSLATIONAL]
1721  << token::END_STATEMENT << nl;
1722  os.writeKeyword("separationVector") << separationVector_
1723  << token::END_STATEMENT << nl;
1724  break;
1725  }
1726  default:
1727  {
1728  // no additional info to write
1729  }
1730  }
1731 }
1732 
1733 
1734 // ************************ vim: set sw=4 sts=4 et: ************************ //