createCoupleMatches.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 Description
25  Create coupled match faces and add them to the cells
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "starMesh.H"
30 #include "boolList.H"
31 #include "pointHit.H"
32 #include "IOmanip.H"
33 #include "boundBox.H"
34 #include "Map.H"
35 #include "unitConversion.H"
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 void Foam::starMesh::createCoupleMatches()
40 {
41  // Loop through all couples and create intersection faces. Add all points
42  // of intersection faces to the couple points lists. The numbering of
43  // the list is set such that the list can be appended to the
44  // existing points list
45 
46  // Estimate the number of cells affected by couple matches
47  const label cellMapSize = min
48  (
49  cellShapes_.size()/10,
50  couples_.size()*2
51  );
52 
53  // Store newly created faces for each cell
54  Map<SLList<face>> cellAddedFaces(cellMapSize);
55 
56  Map<SLList<label>> cellRemovedFaces(cellMapSize);
57 
58  // In order to remove often allocation, remember the number of live points.
59  // If you run out of space in point creation, increase it by the number of
60  // couples (good scale) and resize at the end;
61  label nLivePoints = points_.size();
62 
63  const label infoJump = max(1000, couples_.size()/20);
64 
65  forAll(couples_, coupleI)
66  {
67  if (coupleI % infoJump == 0)
68  {
69  Info<< "Doing couple " << coupleI << ". STAR couple ID: "
70  << couples_[coupleI].coupleID() << endl;
71  }
72 
73  // Initialise cell edges for master and slave cells
74  const coupledFacePair& fp = couples_[coupleI];
75  const face& masterFace = cellFaces_[fp.masterCell()][fp.masterFace()];
76  const face& slaveFace = cellFaces_[fp.slaveCell()][fp.slaveFace()];
77 
78  #ifdef DEBUG_COUPLE
79  Info<< "coupleI: " << coupleI << endl
80  << "masterFace: " << masterFace << endl
81  << "master points: " << masterFace.points(points_) << endl
82  << "slaveFace: " << slaveFace << endl
83  << "slave points: " << slaveFace.points(points_)
84  << endl << endl;
85  #endif
86 
87  // check the angle of face area vectors
88  scalar faceAreaAngle =
89  mag
90  (
91  -(masterFace.area(points_) & slaveFace.area(points_))/
92  (masterFace.mag(points_)*slaveFace.mag(points_) + vSmall)
93  );
94 
95  if (faceAreaAngle < 0.94)
96  {
97  Info<< "Couple direction mismatch in the couple match "
98  << coupleI << ". STAR couple ID: "
99  << couples_[coupleI].coupleID() << endl
100  << "The angle between face normals is "
101  << radToDeg(Foam::acos(faceAreaAngle))
102  << " deg." << endl
103  << "master cell: " << fp.masterCell()
104  << " STAR number: " << starCellID_[fp.masterCell()]
105  << " type: " << cellShapes_[fp.masterCell()].model().name()
106  << " face: " << fp.masterFace() << endl
107  << "slave cell : " << fp.slaveCell()
108  << " STAR number: " << starCellID_[fp.slaveCell()]
109  << " type: " << cellShapes_[fp.slaveCell()].model().name()
110  << " face: " << fp.slaveFace() << endl;
111  }
112 
113  // Deal with integral patches
114  if (fp.integralMatch())
115  {
116  // Master face is replaced by a set of slave faces
117 
118  Map<SLList<label>>::iterator crfIter =
119  cellRemovedFaces.find(fp.masterCell());
120 
121  if (crfIter == cellRemovedFaces.end())
122  {
123  cellRemovedFaces.insert
124  (
125  fp.masterCell(),
126  SLList<label>(fp.masterFace())
127  );
128  }
129  else
130  {
131  crfIter().append(fp.masterFace());
132  }
133 
134  Map<SLList<face>>::iterator cafIter =
135  cellAddedFaces.find(fp.masterCell());
136  if (cafIter == cellAddedFaces.end())
137  {
138  cellAddedFaces.insert
139  (
140  fp.masterCell(),
141  SLList<face>(slaveFace.reverseFace())
142  );
143  }
144  else
145  {
146  cafIter().append(slaveFace.reverseFace());
147  }
148  }
149  else
150  {
151  // Create cut faces, which replace both master and slave faces
152 
153  // Store newly created points
154  SLList<point> coupleFacePoints;
155 
156  // Master data
157  edgeList masterEdges = masterFace.edges();
158  List<SLList<label>> masterEdgePoints(masterEdges.size());
159 
160  // Slave data
161  edgeList slaveEdges = slaveFace.edges();
162  List<SLList<label>> slaveEdgePoints(slaveEdges.size());
163 
164  // Find common plane
165  vector n = masterFace.area(points_);
166  n /= mag(n) + vSmall;
167 
168  // Loop through all edges of the master face. For every edge,
169  // intersect it with all edges of the cutting face.
170  forAll(masterEdges, masterEdgeI)
171  {
172  const edge& curMasterEdge = masterEdges[masterEdgeI];
173 
174  point P = points_[curMasterEdge.start()];
175 
176  // get d and return it into plane
177  vector d = curMasterEdge.vec(points_);
178  d -= n*(n & d);
179 
180  #ifdef DEBUG_COUPLE_INTERSECTION
181  Info<< "curMasterEdge: " << curMasterEdge << endl
182  << "P: " << P << endl << "d: " << d << endl;
183  #endif
184 
185  // go through all slave edges and try to get an intersection.
186  // The point is created along the original master edge rather
187  // than its corrected direction.
188  forAll(slaveEdges, slaveEdgeI)
189  {
190  const edge& curSlaveEdge = slaveEdges[slaveEdgeI];
191 
192  point S = points_[curSlaveEdge.start()];
193 
194  // get e and return it into plane
195  vector e = curSlaveEdge.vec(points_);
196  e -= n*(n & e);
197  scalar det = -(e & (n ^ d));
198 
199  #ifdef DEBUG_COUPLE_INTERSECTION
200  Info<< "curSlaveEdge: " << curSlaveEdge << endl
201  << "S: " << S << endl
202  << "e: " << e << endl;
203  #endif
204 
205  if (mag(det) > small)
206  {
207  // non-singular matrix. Look for intersection
208  scalar beta = ((S - P) & (n ^ d))/det;
209 
210  #ifdef DEBUG_COUPLE_INTERSECTION
211  Info<< " beta: " << beta << endl;
212  #endif
213 
214  if (beta > -smallMergeTol_ && beta < 1 + smallMergeTol_)
215  {
216  // slave intersection OK. Try master intersection
217  scalar alpha =
218  (((S - P) & d) + beta*(d & e))/magSqr(d);
219 
220  #ifdef DEBUG_COUPLE_INTERSECTION
221  Info<< " alpha: " << alpha << endl;
222  #endif
223 
224  if
225  (
226  alpha > -smallMergeTol_
227  && alpha < 1 + smallMergeTol_
228  )
229  {
230  // intersection of non-parallel edges
231  #ifdef DEBUG_COUPLE_INTERSECTION
232  Info<< "intersection of non-parallel edges"
233  << endl;
234  #endif
235 
236 
237  // check for insertion of start-end
238  // points in the middle of the other
239  // edge
240  if (alpha < smallMergeTol_)
241  {
242  // inserting the start of master edge
243  if
244  (
245  beta > smallMergeTol_
246  && beta < 1 - smallMergeTol_
247  )
248  {
249  slaveEdgePoints[slaveEdgeI].append
250  (
251  curMasterEdge.start()
252  );
253  }
254  }
255  else if (alpha > 1 - smallMergeTol_)
256  {
257  // inserting the end of master edge
258  if
259  (
260  beta > smallMergeTol_
261  && beta < 1 - smallMergeTol_
262  )
263  {
264  slaveEdgePoints[slaveEdgeI].append
265  (
266  curMasterEdge.end()
267  );
268  }
269  }
270  else if (beta < smallMergeTol_)
271  {
272  // inserting the start of the slave edge
273  if
274  (
275  alpha > smallMergeTol_
276  && alpha < 1 - smallMergeTol_
277  )
278  {
279  masterEdgePoints[masterEdgeI].append
280  (
281  curSlaveEdge.start()
282  );
283  }
284  }
285  else if (beta > 1 - smallMergeTol_)
286  {
287  // inserting the start of the slave edge
288  if
289  (
290  alpha > smallMergeTol_
291  && alpha < 1 - smallMergeTol_
292  )
293  {
294  masterEdgePoints[masterEdgeI].append
295  (
296  curSlaveEdge.end()
297  );
298  }
299  }
300  else
301  {
302  masterEdgePoints[masterEdgeI].append
303  (
304  nLivePoints + coupleFacePoints.size()
305  );
306 
307  slaveEdgePoints[slaveEdgeI].append
308  (
309  nLivePoints + coupleFacePoints.size()
310  );
311 
312  #ifdef DEBUG_COUPLE_INTERSECTION
313  Info<< "regular intersection. "
314  << "Adding point: "
315  << coupleFacePoints.size()
316  << " which is "
317  << P + alpha*curMasterEdge.vec(points_)
318  << endl;
319  #endif
320 
321  // A new point is created. Warning:
322  // using original edge for accuracy.
323  //
324  coupleFacePoints.append
325  (P + alpha*curMasterEdge.vec(points_));
326  }
327  }
328  }
329  }
330  else
331  {
332  // Add special cases, for intersection of two
333  // parallel line Warning. Here, typically, no new
334  // points will be created. Either one or two of
335  // the slave edge points need to be added to the
336  // master edge and vice versa. The problem is that
337  // no symmetry exists, i.e. both operations needs
338  // to be done separately for both master and slave
339  // side.
340 
341  // Master side
342  // check if the first or second point of slave edge is
343  // on the master edge
344  vector ps = S - P;
345 
346  bool colinear = false;
347 
348  if (mag(ps) < small)
349  {
350  // colinear because P and S are the same point
351  colinear = true;
352  }
353  else if
354  (
355  (ps & d)/(mag(ps)*mag(d)) > 1.0 - smallMergeTol_
356  )
357  {
358  // colinear because ps and d are parallel
359  colinear = true;
360  }
361 
362  if (colinear)
363  {
364  scalar alpha1 = (ps & d)/magSqr(d);
365 
366  if
367  (
368  alpha1 > -smallMergeTol_
369  && alpha1 < 1 + smallMergeTol_
370  )
371  {
372  #ifdef DEBUG_COUPLE_INTERSECTION
373  Info<< "adding irregular master "
374  << "intersection1: "
375  << points_[slaveEdges[slaveEdgeI].start()]
376  << endl;
377  #endif
378 
379  masterEdgePoints[masterEdgeI].append
380  (
381  slaveEdges[slaveEdgeI].start()
382  );
383  }
384 
385  scalar alpha2 = ((ps + e) & d)/magSqr(d);
386 
387  if
388  (
389  alpha2 > -smallMergeTol_
390  && alpha2 < 1 + smallMergeTol_
391  )
392  {
393  #ifdef DEBUG_COUPLE_INTERSECTION
394  Info<< "adding irregular master "
395  << "intersection2: "
396  << points_[slaveEdges[slaveEdgeI].end()]
397  << endl;
398  #endif
399 
400  masterEdgePoints[masterEdgeI].append
401  (
402  slaveEdges[slaveEdgeI].end()
403  );
404  }
405 
406  // Slave side
407  // check if the first or second point of
408  // master edge is on the slave edge
409 
410  vector sp = P - S;
411 
412  scalar beta1 = (sp & e)/magSqr(e);
413 
414  #ifdef DEBUG_COUPLE_INTERSECTION
415  Info<< "P: " << P << " S: " << S << " d: " << d
416  << " e: " << e << " sp: " << sp
417  << " beta1: " << beta1 << endl;
418  #endif
419 
420  if
421  (
422  beta1 > -smallMergeTol_
423  && beta1 < 1 + smallMergeTol_
424  )
425  {
426  #ifdef DEBUG_COUPLE_INTERSECTION
427  Info<< "adding irregular slave "
428  << "intersection1: "
429  << points_[masterEdges[masterEdgeI].start()]
430  << endl;
431  #endif
432 
433  slaveEdgePoints[slaveEdgeI].append
434  (
435  masterEdges[masterEdgeI].start()
436  );
437  }
438 
439  scalar beta2 = ((sp + d) & e)/magSqr(e);
440 
441  if
442  (
443  beta2 > -smallMergeTol_
444  && beta2 < 1 + smallMergeTol_
445  )
446  {
447  #ifdef DEBUG_COUPLE_INTERSECTION
448  Info<< "adding irregular slave "
449  << "intersection2: "
450  << points_[masterEdges[masterEdgeI].end()]
451  << endl;
452  #endif
453 
454  slaveEdgePoints[slaveEdgeI].append
455  (
456  masterEdges[masterEdgeI].end()
457  );
458  }
459  } // end of colinear
460  } // end of singular intersection
461  } // end of slave edges
462  } // end of master edges
463 
464  #ifdef DEBUG_COUPLE_INTERSECTION
465  Info<< "additional slave edge points: " << endl;
466  forAll(slaveEdgePoints, edgeI)
467  {
468  Info<< "edge: " << edgeI << ": " << slaveEdgePoints[edgeI]
469  << endl;
470  }
471  #endif
472 
473  // Add new points
474  if (nLivePoints + coupleFacePoints.size() >= points_.size())
475  {
476  // increase the size of the points list
477  Info<< "Resizing points list" << endl;
478  points_.setSize(points_.size() + couples_.size());
479  }
480 
481  for
482  (
483  SLList<point>::iterator coupleFacePointsIter =
484  coupleFacePoints.begin();
485  coupleFacePointsIter != coupleFacePoints.end();
486  ++coupleFacePointsIter
487  )
488  {
489  points_[nLivePoints] = coupleFacePointsIter();
490  nLivePoints++;
491  }
492 
493  // edge intersection finished
494 
495  // Creating new master side
496 
497  // count the number of additional points for face
498  label nAdditionalMasterPoints = 0;
499 
500  forAll(masterEdgePoints, edgeI)
501  {
502  nAdditionalMasterPoints += masterEdgePoints[edgeI].size();
503  }
504 
505  face tmpMasterFace
506  (
507  masterFace.size()
508  + nAdditionalMasterPoints
509  );
510  label nTmpMasterLabels = 0;
511 
512  #ifdef DEBUG_COUPLE_INTERSECTION
513  Info<< "masterFace: " << masterFace << endl
514  << "nAdditionalMasterPoints: " << nAdditionalMasterPoints
515  << endl;
516  #endif
517 
518  forAll(masterEdges, masterEdgeI)
519  {
520  // Insert the starting point of the edge
521  tmpMasterFace[nTmpMasterLabels] =
522  masterEdges[masterEdgeI].start();
523  nTmpMasterLabels++;
524 
525  // get reference to added points of current edge
526  const SLList<label>& curMEdgePoints =
527  masterEdgePoints[masterEdgeI];
528 
529  // create a markup list of points that have been used
530  boolList usedMasterPoint(curMEdgePoints.size(), false);
531 
532  vector edgeVector = masterEdges[masterEdgeI].vec(points_);
533 
534  #ifdef DEBUG_FACE_ORDERING
535  Info<< "edgeVector: " << edgeVector << endl
536  << "curMEdgePoints.size(): " << curMEdgePoints.size()
537  << endl;
538  #endif
539 
540  // renormalise
541  edgeVector /= magSqr(edgeVector);
542 
543  point edgeStartPoint =
544  points_[masterEdges[masterEdgeI].start()];
545 
546  // loop until the next label to add is -1
547  for (;;)
548  {
549  label nextPointLabel = -1;
550  label usedI = -1;
551  scalar minAlpha = great;
552 
553  label i = 0;
554 
555  for
556  (
557  SLList<label>::const_iterator curMEdgePointsIter =
558  curMEdgePoints.begin();
559  curMEdgePointsIter != curMEdgePoints.end();
560  ++curMEdgePointsIter
561  )
562  {
563  if (!usedMasterPoint[i])
564  {
565  scalar alpha =
566  edgeVector
567  & (
568  points_[curMEdgePointsIter()]
569  - edgeStartPoint
570  );
571 
572  #ifdef DEBUG_FACE_ORDERING
573  Info<< " edgeStartPoint: " << edgeStartPoint
574  << " edgeEndPoint: "
575  << points_[masterEdges[masterEdgeI].end()]
576  << " other point: "
577  << points_[curMEdgePointsIter()]
578  << " alpha: " << alpha << endl;
579  #endif
580 
581  if (alpha < minAlpha)
582  {
583  minAlpha = alpha;
584  usedI = i;
585  nextPointLabel = curMEdgePointsIter();
586  }
587  }
588 
589  #ifdef DEBUG_FACE_ORDERING
590  Info<< "nextPointLabel: " << nextPointLabel << endl;
591  #endif
592 
593  i++;
594  }
595 
596  if (nextPointLabel > -1)
597  {
598  #ifdef DEBUG_FACE_ORDERING
599  Info<< "added nextPointLabel: " << nextPointLabel
600  << " nTmpMasterLabels: " << nTmpMasterLabels
601  << " to place " << nTmpMasterLabels << endl;
602  #endif
603 
604  usedMasterPoint[usedI] = true;
605  // add the next point
606  tmpMasterFace[nTmpMasterLabels] =
607  nextPointLabel;
608  nTmpMasterLabels++;
609  }
610  else
611  {
612  break;
613  }
614  }
615  }
616 
617  // reset the size of master
618  tmpMasterFace.setSize(nTmpMasterLabels);
619 
620  #ifdef DEBUG_FACE_ORDERING
621  Info<< "tmpMasterFace: " << tmpMasterFace << endl;
622  #endif
623 
624  // Eliminate all zero-length edges
625  face newMasterFace(labelList(tmpMasterFace.size(), labelMax));
626 
627  // insert first point by hand. Careful: the first one is
628  // used for comparison to allow the edge collapse across
629  // point zero.
630  //
631  newMasterFace[0] = tmpMasterFace[0];
632  label nMaster = 0;
633 
634  edgeList mstEdgesToCollapse = tmpMasterFace.edges();
635 
636  scalar masterTol =
637  cpMergePointTol_*boundBox(tmpMasterFace.points(points_)).mag();
638 
639  forAll(mstEdgesToCollapse, edgeI)
640  {
641  #ifdef DEBUG_FACE_ORDERING
642  Info<< "edgeI: " << edgeI << " curEdge: "
643  << mstEdgesToCollapse[edgeI] << endl
644  << "master edge " << edgeI << ", "
645  << mstEdgesToCollapse[edgeI].mag(points_) << endl;
646  #endif
647 
648  // Edge merge tolerance = masterTol
649  if (mstEdgesToCollapse[edgeI].mag(points_) < masterTol)
650  {
651  newMasterFace[nMaster] =
652  min
653  (
654  newMasterFace[nMaster],
655  mstEdgesToCollapse[edgeI].end()
656  );
657 
658  #ifdef DEBUG_FACE_ORDERING
659  Info<< "Collapsed: nMaster: " << nMaster
660  << " label: " << newMasterFace[nMaster] << endl;
661  #endif
662 
663  }
664  else
665  {
666  nMaster++;
667 
668  if (edgeI < mstEdgesToCollapse.size() - 1)
669  {
670  // last edge does not add the point
671  #ifdef DEBUG_FACE_ORDERING
672  Info<< "Added: nMaster: " << nMaster
673  << " label: " << mstEdgesToCollapse[edgeI].end()
674  << endl;
675  #endif
676 
677  newMasterFace[nMaster] =
678  mstEdgesToCollapse[edgeI].end();
679  }
680  }
681  }
682 
683  newMasterFace.setSize(nMaster);
684 
685  #ifdef DEBUG_COUPLE
686  Info<< "newMasterFace: " << newMasterFace << endl
687  << "points: " << newMasterFace.points(points_) << endl;
688  #endif
689 
690  // Creating new slave side
691 
692  // count the number of additional points for face
693  label nAdditionalSlavePoints = 0;
694 
695  forAll(slaveEdgePoints, edgeI)
696  {
697  nAdditionalSlavePoints += slaveEdgePoints[edgeI].size();
698  }
699 
700  face tmpSlaveFace
701  (
702  slaveFace.size()
703  + nAdditionalSlavePoints
704  );
705  label nTmpSlaveLabels = 0;
706 
707  #ifdef DEBUG_COUPLE_INTERSECTION
708  Info<< "slaveFace: " << slaveFace << endl
709  << "nAdditionalSlavePoints: " << nAdditionalSlavePoints << endl;
710  #endif
711 
712  forAll(slaveEdges, slaveEdgeI)
713  {
714  // Insert the starting point of the edge
715  tmpSlaveFace[nTmpSlaveLabels] =
716  slaveEdges[slaveEdgeI].start();
717  nTmpSlaveLabels++;
718 
719  // get reference to added points of current edge
720  const SLList<label>& curSEdgePoints =
721  slaveEdgePoints[slaveEdgeI];
722 
723  // create a markup list of points that have been used
724  boolList usedSlavePoint(curSEdgePoints.size(), false);
725 
726  vector edgeVector = slaveEdges[slaveEdgeI].vec(points_);
727 
728  #ifdef DEBUG_FACE_ORDERING
729  Info<< "curSEdgePoints.size(): "
730  << curSEdgePoints.size() << endl
731  << "edgeVector: " << edgeVector << endl;
732  #endif
733 
734  // renormalise
735  edgeVector /= magSqr(edgeVector);
736 
737  point edgeStartPoint =
738  points_[slaveEdges[slaveEdgeI].start()];
739 
740  // loop until the next label to add is -1
741  for (;;)
742  {
743  label nextPointLabel = -1;
744  label usedI = -1;
745  scalar minAlpha = great;
746 
747  label i = 0;
748 
749  for
750  (
751  SLList<label>::const_iterator curSEdgePointsIter =
752  curSEdgePoints.begin();
753  curSEdgePointsIter != curSEdgePoints.end();
754  ++curSEdgePointsIter
755  )
756  {
757  if (!usedSlavePoint[i])
758  {
759  scalar alpha =
760  edgeVector
761  & (
762  points_[curSEdgePointsIter()]
763  - edgeStartPoint
764  );
765 
766  #ifdef DEBUG_FACE_ORDERING
767  Info<< " edgeStartPoint: " << edgeStartPoint
768  << " edgeEndPoint: "
769  << points_[slaveEdges[slaveEdgeI].end()]
770  << " other point: "
771  << points_[curSEdgePointsIter()]
772  << " alpha: " << alpha << endl;
773  #endif
774 
775  if (alpha < minAlpha)
776  {
777  minAlpha = alpha;
778  usedI = i;
779  nextPointLabel = curSEdgePointsIter();
780  }
781  }
782 
783  #ifdef DEBUG_FACE_ORDERING
784  Info<< "nextPointLabel: " << nextPointLabel << endl;
785  #endif
786 
787  i++;
788  }
789 
790  if (nextPointLabel > -1)
791  {
792  #ifdef DEBUG_FACE_ORDERING
793  Info<< "added nextPointLabel: " << nextPointLabel
794  << " nTmpSlaveLabels: " << nTmpSlaveLabels
795  << " to place " << nTmpSlaveLabels << endl;
796  #endif
797 
798  usedSlavePoint[usedI] = true;
799  // add the next point
800  tmpSlaveFace[nTmpSlaveLabels] =
801  nextPointLabel;
802  nTmpSlaveLabels++;
803  }
804  else
805  {
806  break;
807  }
808  }
809  }
810 
811  // reset the size of slave
812  tmpSlaveFace.setSize(nTmpSlaveLabels);
813 
814  #ifdef DEBUG_FACE_ORDERING
815  Info<< "tmpSlaveFace: " << tmpSlaveFace << endl;
816  #endif
817 
818  // Eliminate all zero-length edges
819  face newSlaveFace(labelList(tmpSlaveFace.size(), labelMax));
820 
821  // insert first point by hand. Careful: the first one is
822  // used for comparison to allow the edge collapse across
823  // point zero.
824  //
825  newSlaveFace[0] = tmpSlaveFace[0];
826  label nSlave = 0;
827 
828  edgeList slvEdgesToCollapse = tmpSlaveFace.edges();
829 
830  scalar slaveTol =
831  cpMergePointTol_*boundBox(tmpSlaveFace.points(points_)).mag();
832 
833  forAll(slvEdgesToCollapse, edgeI)
834  {
835  #ifdef DEBUG_FACE_ORDERING
836  Info<< "slave edge length: " << edgeI << ", "
837  << slvEdgesToCollapse[edgeI].mag(points_)<< endl;
838  #endif
839 
840  // edge merge tolerance = slaveTol
841  if (slvEdgesToCollapse[edgeI].mag(points_) < slaveTol)
842  {
843  newSlaveFace[nSlave] =
844  min
845  (
846  newSlaveFace[nSlave],
847  slvEdgesToCollapse[edgeI].end()
848  );
849  }
850  else
851  {
852  nSlave++;
853  if (edgeI < slvEdgesToCollapse.size() - 1)
854  {
855  // last edge does not add the point
856  newSlaveFace[nSlave] = slvEdgesToCollapse[edgeI].end();
857  }
858  }
859  }
860 
861  newSlaveFace.setSize(nSlave);
862 
863  #ifdef DEBUG_COUPLE
864  Info<< "newSlaveFace: " << newSlaveFace << endl
865  << "points: " << newSlaveFace.points(points_) << endl << endl;
866  #endif
867 
868  // Create the intersection face
869 
870  // Algorithm:
871  // Loop through
872  // points of the master and try to find one which falls
873  // within the slave. If not found, look through all
874  // edges of the slave and find one which falls within the
875  // master. This point will be the starting location for
876  // the cut face.
877 
878  edgeList newMasterEdges = newMasterFace.edges();
879  edgeList newSlaveEdges = newSlaveFace.edges();
880 
881  #ifdef DEBUG_RIGHT_HAND_WALK
882  Info<< "newMasterEdges: " << newMasterEdges << endl
883  << "newSlaveEdges: " << newSlaveEdges << endl;
884  #endif
885 
886  edge startEdge(-1, -1);
887 
888  // Remember where the start edge was found:
889  // 0 for not found
890  // 1 for master
891  // 2 for slave
892  label startEdgeFound = 0;
893 
894  vector masterProjDir = -newMasterFace.area(points_);
895 
896  forAll(newSlaveEdges, edgeI)
897  {
898  // Take the slave edge points and project into the master.
899  // In order to create a good intersection, move the
900  // point away from the master in the direction of its
901  // normal.
902  point pointStart = points_[newSlaveEdges[edgeI].start()];
903 
904  point pointEnd = points_[newSlaveEdges[edgeI].end()];
905 
906  if
907  (
908  newMasterFace.ray
909  (
910  pointStart,
911  masterProjDir,
912  points_,
914  ).hit()
915  && newMasterFace.ray
916  (
917  pointEnd,
918  masterProjDir,
919  points_,
921  ).hit()
922  )
923  {
924  startEdge = newSlaveEdges[edgeI];
925  startEdgeFound = 2;
926 
927  #ifdef DEBUG_RIGHT_HAND_WALK
928  Info<< "slave edge found" << endl;
929  #endif
930 
931  break;
932  }
933  }
934 
935  if (startEdgeFound == 0)
936  {
937  vector slaveProjDir = -newSlaveFace.area(points_);
938 
939  forAll(newMasterEdges, edgeI)
940  {
941  // Take the edge master points and project into the slave.
942  // In order to create a good intersection, move the
943  // point away from the slave in the direction of its
944  // normal.
945  point pointStart = points_[newMasterEdges[edgeI].start()];
946 
947  point pointEnd = points_[newMasterEdges[edgeI].end()];
948 
949  if
950  (
951  newSlaveFace.ray
952  (
953  pointStart,
954  slaveProjDir,
955  points_,
957  ).hit()
958  && newSlaveFace.ray
959  (
960  pointEnd,
961  slaveProjDir,
962  points_,
964  ).hit()
965  )
966  {
967  startEdge = newMasterEdges[edgeI];
968  startEdgeFound = 1;
969 
970  #ifdef DEBUG_RIGHT_HAND_WALK
971  Info<< "master edge found" << endl;
972  #endif
973 
974  break;
975  }
976  }
977  }
978 
979  // create the intersected face using right-hand walk rule
980  face intersectedFace
981  (
982  labelList(newMasterFace.size() + newSlaveFace.size(), -1)
983  );
984 
985  if (startEdgeFound > 0)
986  {
987  #ifdef DEBUG_RIGHT_HAND_WALK
988  Info<< "start edge: " << startEdge << endl;
989  #endif
990 
991  // Loop through both faces and add all edges
992  // containing the current point and add them to the
993  // list of edges to consider. Make sure all edges are
994  // added such that the current point is their start.
995  // Loop through all edges to consider and find the one
996  // which produces the biggest right-hand-turn. This
997  // is the next edge to be added to the face. If its
998  // end is the same as the starting point, the face is
999  // complete; resize it to the number of active points
1000  // and exit.
1001 
1002  vector planeNormal = newMasterFace.area(points_);
1003  planeNormal /= mag(planeNormal) + vSmall;
1004 
1005  #ifdef DEBUG_RIGHT_HAND_WALK
1006  Info<< "planeNormal: " << planeNormal << endl;
1007  #endif
1008 
1009  // Do a check to control the right-hand turn. This is
1010  // based on the triple product of the edge start
1011  // vector to face centre, the edge vector and the
1012  // plane normal. If the triple product is negative,
1013  // the edge needs to be reversed to allow the
1014  // right-hand-turn rule to work.
1015 
1016  vector faceCentre;
1017 
1018  if (startEdgeFound == 1)
1019  {
1020  faceCentre = newMasterFace.centre(points_);
1021  }
1022  else
1023  {
1024  faceCentre = newSlaveFace.centre(points_);
1025  }
1026 
1027  scalar tripleProduct =
1028  (
1029  (faceCentre - points_[startEdge.start()])
1030  ^ startEdge.vec(points_)
1031  ) & planeNormal;
1032 
1033  if (tripleProduct < 0)
1034  {
1035  #ifdef DEBUG_RIGHT_HAND_WALK
1036  Info<< "Turning edge for right-hand turn rule" << endl;
1037  #endif
1038  startEdge.flip();
1039  }
1040 
1041  // prepare the loop for the right-hand walk
1042  intersectedFace[0] = startEdge.start();
1043  intersectedFace[1] = startEdge.end();
1044  label nIntFacePoints = 2;
1045 
1046  edge curEdge = startEdge;
1047 
1048  bool completedFace = false;
1049 
1050  do
1051  {
1052  SLList<edge> edgesToConsider;
1053 
1054  // collect master edges
1055  forAll(newMasterEdges, edgeI)
1056  {
1057  const edge& cme = newMasterEdges[edgeI];
1058 
1059  if (cme != curEdge)
1060  {
1061  if (cme.start() == curEdge.end())
1062  {
1063  edgesToConsider.append(cme);
1064  }
1065  else if (cme.end() == curEdge.end())
1066  {
1067  edgesToConsider.append(cme.reverseEdge());
1068  }
1069  // otherwise, it does not have the current point
1070  }
1071  }
1072 
1073  // collect slave edges
1074  forAll(newSlaveEdges, edgeI)
1075  {
1076  const edge& cse = newSlaveEdges[edgeI];
1077 
1078  if (cse != curEdge)
1079  {
1080  if (cse.start() == curEdge.end())
1081  {
1082  edgesToConsider.append(cse);
1083  }
1084  else if (cse.end() == curEdge.end())
1085  {
1086  edgesToConsider.append(cse.reverseEdge());
1087  }
1088  // otherwise, it does not have the current point
1089  }
1090  }
1091 
1092  #ifdef DEBUG_RIGHT_HAND_WALK
1093  Info<< "number of edges to consider: "
1094  << edgesToConsider.size() << endl
1095  << "edges to consider: " << edgesToConsider << endl;
1096  #endif
1097 
1098  if (edgesToConsider.empty())
1099  {
1101  << setprecision(12)
1102  << "void starMesh::createCoupleMatches() : "
1103  << endl << "error in face intersection: "
1104  << "no edges to consider for closing the loop"
1105  << coupleI << ". STAR couple ID: "
1106  << couples_[coupleI].coupleID() << endl
1107  << "Cut Master Face: " << newMasterFace << endl
1108  << "points: " << newMasterFace.points(points_)
1109  << endl
1110  << "Cut Slave Face: " << newSlaveFace << endl
1111  << "points: " << newSlaveFace.points(points_)
1112  << endl << "intersected face: "
1113  << intersectedFace
1114  << abort(FatalError);
1115  }
1116 
1117  // vector along the edge
1118  vector ahead = curEdge.vec(points_);
1119  ahead -= planeNormal*(planeNormal & ahead);
1120  ahead /= mag(ahead) + vSmall;
1121 
1122  // vector pointing right
1123  vector right = ahead ^ planeNormal;
1124  right /= mag(right) + vSmall;
1125 
1126  // first edge taken for reference
1127  edge nextEdge = edgesToConsider.first();
1128  vector nextEdgeVec = nextEdge.vec(points_);
1129  nextEdgeVec -= planeNormal*(planeNormal & nextEdgeVec);
1130  nextEdgeVec /= mag(nextEdgeVec) + vSmall;
1131 
1132  scalar rightTurn = nextEdgeVec & right;
1133  scalar goStraight = nextEdgeVec & ahead;
1134 
1135  #ifdef DEBUG_RIGHT_HAND_WALK
1136  Info<< "rightTurn: " << rightTurn
1137  << " goStraight: " << goStraight << endl;
1138  #endif
1139 
1140  for
1141  (
1142  SLList<edge>::iterator etcIter =
1143  edgesToConsider.begin();
1144  etcIter != edgesToConsider.end();
1145  ++etcIter
1146  )
1147  {
1148  // right-hand walk rule
1149  vector newDir = etcIter().vec(points_);
1150  newDir -= planeNormal*(planeNormal & newDir);
1151  newDir /= mag(newDir) + vSmall;
1152 
1153  scalar curRightTurn = newDir & right;
1154  scalar curGoStraight = newDir & ahead;
1155 
1156  #ifdef DEBUG_RIGHT_HAND_WALK
1157  Info<< "curRightTurn: " << curRightTurn
1158  << " curGoStraight: " << curGoStraight << endl;
1159  #endif
1160 
1161  if (rightTurn < 0) // old edge turning left
1162  {
1163  if (curRightTurn < 0) // new edge turning left
1164  {
1165  // both go left. Grab one with greater ahead
1166  if (curGoStraight > goStraight)
1167  {
1168  #ifdef DEBUG_RIGHT_HAND_WALK
1169  Info<< "a" << endl;
1170  #endif
1171 
1172  // Good edge, turning left less than before
1173  nextEdge = etcIter();
1174  rightTurn = curRightTurn;
1175  goStraight = curGoStraight;
1176  }
1177  }
1178  else // new edge turning right
1179  {
1180  #ifdef DEBUG_RIGHT_HAND_WALK
1181  Info<< "b" << endl;
1182  #endif
1183 
1184  // good edge, turning right
1185  nextEdge = etcIter();
1186  rightTurn = curRightTurn;
1187  goStraight = curGoStraight;
1188  }
1189  }
1190  else // old edge turning right
1191  {
1192  // new edge turning left rejected
1193  if (curRightTurn >= 0) // new edge turning right
1194  {
1195  // grab one with smaller ahead
1196  if (curGoStraight < goStraight)
1197  {
1198  #ifdef DEBUG_RIGHT_HAND_WALK
1199  Info<< "c" << endl;
1200  #endif
1201 
1202  // good edge, turning right more than before
1203  nextEdge = etcIter();
1204  rightTurn = curRightTurn;
1205  goStraight = curGoStraight;
1206  }
1207  }
1208  }
1209  }
1210 
1211  // check if the loop is completed
1212  if (nextEdge.end() == intersectedFace[0])
1213  {
1214  // loop is completed. No point to add
1215  completedFace = true;
1216  }
1217  else
1218  {
1219  // Check if there is room for more points
1220  if (nIntFacePoints >= intersectedFace.size())
1221  {
1223  << setprecision(12)
1224  << "void starMesh::createCoupleMatches() : "
1225  << endl << "error in intersected face: "
1226  << "lost thread for intersection of couple "
1227  << coupleI << ". STAR couple ID: "
1228  << couples_[coupleI].coupleID() << endl
1229  << "Cut Master Face: " << newMasterFace << endl
1230  << "points: " << newMasterFace.points(points_)
1231  << endl
1232  << "Cut Slave Face: " << newSlaveFace << endl
1233  << "points: " << newSlaveFace.points(points_)
1234  << endl << "intersected face: "
1235  << intersectedFace
1236  << abort(FatalError);
1237  }
1238 
1239  // insert the point
1240  intersectedFace[nIntFacePoints] = nextEdge.end();
1241  nIntFacePoints++;
1242 
1243  // grab the current point and the current edge
1244  curEdge = nextEdge;
1245 
1246  #ifdef DEBUG_RIGHT_HAND_WALK
1247  Info<< "inserted point " << nextEdge.end() << endl
1248  << "curEdge: " << curEdge << endl;
1249  #endif
1250  }
1251  }
1252  while (!completedFace);
1253 
1254  // resize the face
1255  intersectedFace.setSize(nIntFacePoints);
1256 
1257  #ifdef DEBUG_COUPLE
1258  Info<< "intersectedFace: " << intersectedFace << endl;
1259  #endif
1260 
1261  // check the intersection face for duplicate points
1262  forAll(intersectedFace, checkI)
1263  {
1264  for
1265  (
1266  label checkJ = checkI + 1;
1267  checkJ < intersectedFace.size();
1268  checkJ++
1269  )
1270  {
1271  if (intersectedFace[checkI] == intersectedFace[checkJ])
1272  {
1274  << setprecision(12)
1275  << "void starMesh::createCoupleMatches() : "
1276  << endl << "error in intersected face: "
1277  << "duplicate point in intersected face "
1278  << "for couple no " << coupleI
1279  << ". STAR couple ID: "
1280  << couples_[coupleI].coupleID() << endl
1281  << "Duplicate point label: "
1282  << intersectedFace[checkI] << endl
1283  << "Cut Master Face: " << newMasterFace << endl
1284  << "points: " << newMasterFace.points(points_)
1285  << endl
1286  << "Cut Slave Face: " << newSlaveFace << endl
1287  << "points: " << newSlaveFace.points(points_)
1288  << endl << "intersected face: "
1289  << intersectedFace
1290  << abort(FatalError);
1291  }
1292  }
1293  }
1294  }
1295  else
1296  {
1298  << setprecision(12)
1299  << "void starMesh::createCoupleMatches() : " << endl
1300  << "could not find start edge for intersection of couple "
1301  << coupleI << ". STAR couple ID: "
1302  << couples_[coupleI].coupleID() << endl
1303  << "Cut Master Face: " << newMasterFace << endl
1304  << "points: " << newMasterFace.points(points_) << endl
1305  << "Cut Slave Face: " << newSlaveFace << endl
1306  << "points: " << newSlaveFace.points(points_)
1307  << abort(FatalError);
1308  }
1309 
1310  // Project all points of the intersected face
1311  // onto the master face to ensure closedness
1312  vector pointProjectionNormal = -masterFace.area(points_);
1313 
1314  forAll(intersectedFace, intPointi)
1315  {
1316  #ifdef DEBUG_COUPLE_PROJECTION
1317  Info<< "Proj: old point: "
1318  << points_[intersectedFace[intPointi]] << endl;
1319  #endif
1320 
1321  pointHit projHit =
1322  masterFace.ray
1323  (
1324  points_[intersectedFace[intPointi]],
1325  pointProjectionNormal,
1326  points_,
1328  );
1329 
1330  if (projHit.hit())
1331  {
1332  points_[intersectedFace[intPointi]] =
1333  projHit.hitPoint();
1334 
1335  #ifdef DEBUG_COUPLE_PROJECTION
1336  Info<< " new point: "
1337  << points_[intersectedFace[intPointi]] << endl;
1338  #endif
1339  }
1340  }
1341 
1342  // Check the direction of the intersection face
1343  if
1344  (
1345  (
1346  masterFace.area(points_)
1347  & intersectedFace.area(points_)
1348  ) < vSmall
1349  )
1350  {
1351  intersectedFace.flip();
1352  }
1353 
1354  // Add the new face to both master and slave
1355 
1356  // Master face is replaced by a set of slave faces
1357  Map<SLList<label>>::iterator crfMasterIter =
1358  cellRemovedFaces.find(fp.masterCell());
1359 
1360  if (crfMasterIter == cellRemovedFaces.end())
1361  {
1362  cellRemovedFaces.insert
1363  (
1364  fp.masterCell(),
1365  SLList<label>(fp.masterFace())
1366  );
1367  }
1368  else
1369  {
1370  crfMasterIter().append(fp.masterFace());
1371  }
1372 
1373  Map<SLList<label>>::iterator crfSlaveIter =
1374  cellRemovedFaces.find(fp.slaveCell());
1375 
1376  if (crfSlaveIter == cellRemovedFaces.end())
1377  {
1378  cellRemovedFaces.insert
1379  (
1380  fp.slaveCell(),
1381  SLList<label>(fp.slaveFace())
1382  );
1383  }
1384  else
1385  {
1386  crfSlaveIter().append(fp.slaveFace());
1387  }
1388 
1389  Map<SLList<face>>::iterator cafMasterIter =
1390  cellAddedFaces.find(fp.masterCell());
1391  if (cafMasterIter == cellAddedFaces.end())
1392  {
1393  cellAddedFaces.insert
1394  (
1395  fp.masterCell(),
1396  SLList<face>(intersectedFace)
1397  );
1398  }
1399  else
1400  {
1401  cafMasterIter().append(intersectedFace);
1402  }
1403 
1404  Map<SLList<face>>::iterator cafSlaveIter =
1405  cellAddedFaces.find(fp.slaveCell());
1406  if (cafSlaveIter == cellAddedFaces.end())
1407  {
1408  cellAddedFaces.insert
1409  (
1410  fp.slaveCell(),
1411  SLList<face>(intersectedFace.reverseFace())
1412  );
1413  }
1414  else
1415  {
1416  cafSlaveIter().append(intersectedFace.reverseFace());
1417  }
1418  } // end of arbitrary match
1419  }
1420 
1421  if (couples_.size())
1422  {
1423  // Loop through all cells and reset faces for removal to zero size
1424  const labelList crfToc = cellRemovedFaces.toc();
1425 
1426  forAll(crfToc, celli)
1427  {
1428  const label curCell = crfToc[celli];
1429 
1430  const SLList<label>& curRemovedFaces = cellRemovedFaces[curCell];
1431 
1432  for
1433  (
1434  SLList<label>::const_iterator curRemovedFacesIter =
1435  curRemovedFaces.begin();
1436  curRemovedFacesIter != curRemovedFaces.end();
1437  ++curRemovedFacesIter
1438  )
1439  {
1440  cellFaces_[curCell][curRemovedFacesIter()].setSize(0);
1441  }
1442 
1443  if (curRemovedFaces.size())
1444  {
1445  // reset the shape pointer to unknown
1446  cellShapes_[curCell] = cellShape(*unknownPtr_, labelList(0));
1447  }
1448  }
1449 
1450  const labelList cafToc = cellAddedFaces.toc();
1451 
1452  // Insert the new faces into the list
1453  forAll(cafToc, celli)
1454  {
1455  const label curCell = cafToc[celli];
1456 
1457  const SLList<face>& curAddedFaces = cellAddedFaces[curCell];
1458 
1459  faceList oldFaces = cellFaces_[curCell];
1460 
1461  faceList& newFaces = cellFaces_[curCell];
1462 
1463  newFaces.setSize(oldFaces.size() + curAddedFaces.size());
1464  label nNewFaces = 0;
1465 
1466  // copy original faces that have not been removed
1467  forAll(oldFaces, facei)
1468  {
1469  if (oldFaces[facei].size())
1470  {
1471  newFaces[nNewFaces] = oldFaces[facei];
1472  nNewFaces++;
1473  }
1474  }
1475 
1476  // add new faces
1477  for
1478  (
1479  SLList<face>::const_iterator curAddedFacesIter =
1480  curAddedFaces.begin();
1481  curAddedFacesIter != curAddedFaces.end();
1482  ++curAddedFacesIter
1483  )
1484  {
1485  newFaces[nNewFaces] = curAddedFacesIter();
1486  nNewFaces++;
1487  }
1488 
1489  // reset the size of the face list
1490  newFaces.setSize(nNewFaces);
1491 
1492  if (curAddedFaces.size())
1493  {
1494  // reset the shape pointer to unknown
1495  cellShapes_[curCell] = cellShape(*unknownPtr_, labelList(0));
1496  }
1497  }
1498 
1499  // Resize the point list to the number of created points
1500  points_.setSize(nLivePoints);
1501 
1502  // Finished
1503  Info<< "Finished doing couples" << endl;
1504  }
1505 }
1506 
1507 
1508 // ************************************************************************* //
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
scalar radToDeg(const scalar rad)
Conversion from radians to degrees.
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
friend class iterator
Definition: LList.H:82
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:116
Omanip< int > setprecision(const int i)
Definition: IOmanip.H:205
List< edge > edgeList
Definition: edgeList.H:38
static const label labelMax
Definition: label.H:62
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Istream and Ostream manipulators taking arguments.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
friend class const_iterator
Definition: LList.H:85
void setSize(const label)
Reset size of List.
Definition: List.C:281
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
PointHit< point > pointHit
Definition: pointHit.H:41
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))