All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
enrichedPatchCutFaces.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-2018 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  Calculating cut faces of the enriched patch, together with the addressing
26  into master and slave patch.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include "enrichedPatch.H"
31 #include "boolList.H"
32 #include "DynamicList.H"
33 #include "labelPair.H"
34 #include "primitiveMesh.H"
35 #include "HashSet.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 // If the cut face gets more than this number of points, it will be checked
40 const Foam::label Foam::enrichedPatch::maxFaceSizeDebug_ = 100;
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 // Index of debug signs:
46 // x - skip a point
47 // l - left turn
48 // r - right turn
49 
50 void Foam::enrichedPatch::calcCutFaces() const
51 {
52  if (cutFacesPtr_ || cutFaceMasterPtr_ || cutFaceSlavePtr_)
53  {
55  << "Cut faces addressing already calculated."
56  << abort(FatalError);
57  }
58 
59  const faceList& enFaces = enrichedFaces();
60  const labelList& mp = meshPoints();
61  const faceList& lf = localFaces();
62  const pointField& lp = localPoints();
63  const labelListList& pp = pointPoints();
64  // Pout<< "enFaces: " << enFaces << endl;
65  // Pout<< "lf: " << lf << endl;
66  // Pout<< "lp: " << lp << endl;
67  // Pout<< "pp: " << pp << endl;
68  const Map<labelList>& masterPointFaceAddr = masterPointFaces();
69 
70  // Prepare the storage
71  DynamicList<face> cf(2*lf.size());
72  DynamicList<label> cfMaster(2*lf.size());
73  DynamicList<label> cfSlave(2*lf.size());
74 
75  // Algorithm
76  // Go through all the faces
77  // 1) For each face, start at point zero and grab the edge zero.
78  // Both points are added into the cut face. Mark the first edge
79  // as used and position the current point as the end of the first
80  // edge and previous point as point zero.
81  // 2) Grab all possible points from the current point. Excluding
82  // the previous point find the point giving the biggest right
83  // turn. Add the point to the face and mark the edges as used.
84  // Continue doing this until the face is complete, i.e. the next point
85  // to add is the first point of the face.
86  // 3) Once the facet is completed, register its owner from the face
87  // it has been created from (remember that the first lot of faces
88  // inserted into the enriched faces list are the slaves, to allow
89  // matching of the other side).
90  // 4) If the facet is created from an original slave face, go
91  // through the master patch and try to identify the master face
92  // this facet belongs to. This is recognised by the fact that the
93  // faces consists exclusively of the points of the master face and
94  // the points projected onto the face.
95 
96  // Create a set of edge usage parameters
97  HashSet<edge, Hash<edge>> edgesUsedOnce(pp.size());
98  HashSet<edge, Hash<edge>> edgesUsedTwice
100 
101 
102  forAll(lf, facei)
103  {
104  const face& curLocalFace = lf[facei];
105  const face& curGlobalFace = enFaces[facei];
106 
107  // Pout<< "Doing face " << facei
108  // << " local: " << curLocalFace
109  // << " or " << curGlobalFace
110  // << endl;
111 
112  // if (facei < slavePatch_.size())
113  // {
114  // Pout<< "original slave: " << slavePatch_[facei]
115  // << " local: " << slavePatch_.localFaces()[facei] << endl;
116  // }
117  // else
118  // {
119  // Pout<< "original master: "
120  // << masterPatch_[facei - slavePatch_.size()] << " "
121  // << masterPatch_.localFaces()[facei - slavePatch_.size()]
122  // << endl;
123  // }
124  // {
125  // pointField facePoints = curLocalFace.points(lp);
126  // forAll(curLocalFace, pointi)
127  // {
128  // Pout<< "v " << facePoints[pointi].x() << " "
129  // << facePoints[pointi].y() << " "
130  // << facePoints[pointi].z() << endl;
131  // }
132  // }
133 
134  // Track the usage of face edges. When all edges are used, the
135  // face decomposition is complete.
136  // The seed edges include all the edges of the original face + all edges
137  // of other faces that have been used in the construction of the
138  // facet. Edges from other faces can be considered as
139  // internal to the current face if used only once.
140 
141  // Track the edge usage to avoid duplicate faces and reset it to unused
142  boolList usedFaceEdges(curLocalFace.size(), false);
143 
144  SLList<edge> edgeSeeds;
145 
146  // Insert the edges of current face into the seed list.
147  edgeList cfe = curLocalFace.edges();
148  forAll(curLocalFace, edgeI)
149  {
150  edgeSeeds.append(cfe[edgeI]);
151  }
152 
153  // Grab face normal
154  const vector normal = curLocalFace.normal(lp);
155 
156  while (edgeSeeds.size())
157  {
158  // Pout<< "edgeSeeds.size(): "
159  // << edgeSeeds.size()
160  // << endl;
161 
162  const edge curEdge = edgeSeeds.removeHead();
163 
164  // Locate the edge in current face
165  const label curEdgeWhich = curLocalFace.which(curEdge.start());
166 
167  // Check if the edge is in current face and if it has been
168  // used already. If so, skip it
169  if
170  (
171  curEdgeWhich > -1
172  && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
173  )
174  {
175  // Edge is in the starting face.
176  // If unused, mark it as used; if used, skip it
177  if (usedFaceEdges[curEdgeWhich]) continue;
178 
179  usedFaceEdges[curEdgeWhich] = true;
180  }
181 
182  // If the edge has already been used twice, skip it
183  if (edgesUsedTwice.found(curEdge)) continue;
184 
185  // Pout<< "Trying new edge (" << mp[curEdge.start()]
186  // << ", " << mp[curEdge.end()]
187  // << ") seed: " << curEdge
188  // << " used: " << edgesUsedTwice.found(curEdge)
189  // << endl;
190 
191  // Estimate the size of cut face as twice the size of original face
192  DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
193  DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
194 
195  // Found unused edge.
196  label prevPointLabel = curEdge.start();
197  cutFaceGlobalPoints.append(mp[prevPointLabel]);
198  cutFaceLocalPoints.append(prevPointLabel);
199  // Pout<< "prevPointLabel: " << mp[prevPointLabel] << endl;
200  // Grab current point and append it to the list
201  label curPointLabel = curEdge.end();
202  point curPoint = lp[curPointLabel];
203  cutFaceGlobalPoints.append(mp[curPointLabel]);
204  cutFaceLocalPoints.append(curPointLabel);
205 
206  bool completedCutFace = false;
207 
208  label faceSizeDebug = maxFaceSizeDebug_;
209 
210  do
211  {
212  // Grab the next point options
213 
214  // Pout<< "curPointLabel: " << mp[curPointLabel] << endl;
215 
216  const labelList& nextPoints = pp[curPointLabel];
217 
218  // Pout<< "nextPoints: "
219  // << UIndirectList<label>(mp, nextPoints)
220  // << endl;
221 
222  // Get the vector along the edge and the right vector
223  vector ahead = curPoint - lp[prevPointLabel];
224  ahead -= normal*(normal & ahead);
225  ahead /= mag(ahead);
226 
227  vector right = normal ^ ahead;
228  right /= mag(right);
229 
230  // Pout<< "normal: " << normal
231  // << " ahead: " << ahead
232  // << " right: " << right
233  // << endl;
234 
235  scalar atanTurn = -great;
236  label bestAtanPoint = -1;
237 
238  forAll(nextPoints, nextI)
239  {
240  // Exclude the point we are coming from; there will always
241  // be more than one edge, so this is safe
242  if (nextPoints[nextI] != prevPointLabel)
243  {
244  // Pout<< "cur point: " << curPoint
245  // << " trying for point: "
246  // << mp[nextPoints[nextI]]
247  // << " " << lp[nextPoints[nextI]];
248  vector newDir = lp[nextPoints[nextI]] - curPoint;
249  // Pout<< " newDir: " << newDir
250  // << " mag: " << mag(newDir) << flush;
251  newDir -= normal*(normal & newDir);
252  scalar magNewDir = mag(newDir);
253  // Pout<< " corrected: " << newDir
254  // << " mag: " << mag(newDir) << flush;
255 
256  if (magNewDir < small)
257  {
259  << "projection error: slave patch probably "
260  << "does not project onto master. "
261  << "Please switch on "
262  << "enriched patch debug for more info"
263  << abort(FatalError);
264  }
265 
266  newDir /= magNewDir;
267 
268  scalar curAtanTurn =
269  atan2(newDir & right, newDir & ahead);
270 
271  // Pout<< " atan: " << curAtanTurn << endl;
272 
273  if (curAtanTurn > atanTurn)
274  {
275  bestAtanPoint = nextPoints[nextI];
276  atanTurn = curAtanTurn;
277  }
278  } // end of prev point skip
279  } // end of next point selection
280 
281  // Pout<< " bestAtanPoint: " << bestAtanPoint << " or "
282  // << mp[bestAtanPoint]
283  // << endl;
284 
285  // Selected next best point.
286 
287  // Pout<< "cutFaceGlobalPoints: "
288  // << cutFaceGlobalPoints
289  // << endl;
290 
291  // Check if the edge about to be added has been used
292  // in the current face or twice in other faces. If
293  // so, the face is bad.
294  const label whichNextPoint = curLocalFace.which(curPointLabel);
295 
296  if
297  (
298  whichNextPoint > -1
299  && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
300  && usedFaceEdges[whichNextPoint]
301  )
302  {
303  // This edge is already used in current face
304  // face cannot be good; start on a new one
305 
306  // Pout<< "Double usage in current face, cannot be good"
307  // << endl;
308 
309  completedCutFace = true;
310  }
311 
312 
313  if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
314  {
315  // This edge is already used -
316  // face cannot be good; start on a new one
317  completedCutFace = true;
318 
319  // Pout<< "Double usage elsewhere, cannot be good" << endl;
320  }
321 
322  if (completedCutFace) continue;
323 
324  // Insert the next best point into the list
325  if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
326  {
327  // Face is completed. Save it and move on to the next
328  // available edge
329  completedCutFace = true;
330 
331  if (debug)
332  {
333  Pout<< " local: " << cutFaceLocalPoints
334  << " one side: " << facei;
335  }
336 
337  // Append the face
338  face cutFaceGlobal;
339  cutFaceGlobal.transfer(cutFaceGlobalPoints);
340 
341  cf.append(cutFaceGlobal);
342 
343  face cutFaceLocal;
344  cutFaceLocal.transfer(cutFaceLocalPoints);
345 
346  // Pout<< "\ncutFaceLocal: "
347  // << cutFaceLocal.points(lp)
348  // << endl;
349 
350  // Go through all edges of the cut faces.
351  // If the edge corresponds to a starting face edge,
352  // mark the starting face edge as true
353 
354  forAll(cutFaceLocal, cutI)
355  {
356  const edge curCutFaceEdge
357  (
358  cutFaceLocal[cutI],
359  cutFaceLocal.nextLabel(cutI)
360  );
361 
362  // Increment the usage count using two hash sets
363  HashSet<edge, Hash<edge>>::iterator euoIter =
364  edgesUsedOnce.find(curCutFaceEdge);
365 
366  if (euoIter == edgesUsedOnce.end())
367  {
368  // Pout<< "Found edge not used before: "
369  // << curCutFaceEdge
370  // << endl;
371  edgesUsedOnce.insert(curCutFaceEdge);
372  }
373  else
374  {
375  // Pout<< "Found edge used once: "
376  // << curCutFaceEdge
377  // << endl;
378  edgesUsedOnce.erase(euoIter);
379  edgesUsedTwice.insert(curCutFaceEdge);
380  }
381 
382  const label curCutFaceEdgeWhich = curLocalFace.which
383  (
384  curCutFaceEdge.start()
385  );
386 
387  if
388  (
389  curCutFaceEdgeWhich > -1
390  && curLocalFace.nextLabel(curCutFaceEdgeWhich)
391  == curCutFaceEdge.end()
392  )
393  {
394  // Found edge in original face
395 
396  // Pout<< "Found edge in orig face: "
397  // << curCutFaceEdge << ": "
398  // << curCutFaceEdgeWhich
399  // << endl;
400 
401  usedFaceEdges[curCutFaceEdgeWhich] = true;
402  }
403  else
404  {
405  // Edge not in original face. Add it to seeds
406 
407  // Pout<< "Found new edge seed: "
408  // << curCutFaceEdge
409  // << endl;
410 
411  edgeSeeds.append(curCutFaceEdge.reverseEdge());
412  }
413  }
414 
415  // Find out what the other side is
416 
417  // Algorithm
418 
419  // If the face is in the slave half of the
420  // enrichedFaces lists, it may be matched against
421  // the master face. It can be recognised by the
422  // fact that all its points belong to one master
423  // face. For this purpose:
424  // 1) Grab the master faces around the point zero
425  // of the cut face and collect all master faces
426  // around it.
427  // 2) Loop through the rest of cut face points and
428  // if the point refers to one of the faces marked
429  // by point zero, increment its count.
430  // 3) When completed, the face whose count is
431  // equal to the number of points in the cut face
432  // is the other side. If this is not the case, there is no
433  // face on the other side.
434 
435  if (facei < slavePatch_.size())
436  {
437  Map<labelList>::const_iterator mpfAddrIter =
438  masterPointFaceAddr.find(cutFaceGlobal[0]);
439 
440  bool otherSideFound = false;
441 
442  if (mpfAddrIter != masterPointFaceAddr.end())
443  {
444  bool miss = false;
445 
446  // Create the label count using point zero
447  const labelList& masterFacesOfPZero = mpfAddrIter();
448 
449  labelList hits(masterFacesOfPZero.size(), 1);
450 
451  for
452  (
453  label pointi = 1;
454  pointi < cutFaceGlobal.size();
455  pointi++
456  )
457  {
459  mpfAddrPointIter =
460  masterPointFaceAddr.find
461  (
462  cutFaceGlobal[pointi]
463  );
464 
465  if
466  (
467  mpfAddrPointIter
468  == masterPointFaceAddr.end()
469  )
470  {
471  // Point is off the master patch. Skip
472  miss = true;
473  break;
474  }
475 
476  const labelList& curMasterFaces =
477  mpfAddrPointIter();
478 
479  // For every current face, try to find it in the
480  // zero-list
481  forAll(curMasterFaces, i)
482  {
483  forAll(masterFacesOfPZero, j)
484  {
485  if
486  (
487  curMasterFaces[i]
488  == masterFacesOfPZero[j]
489  )
490  {
491  hits[j]++;
492  break;
493  }
494  }
495  }
496  }
497 
498  // If all point are found attempt matching
499  if (!miss)
500  {
501  forAll(hits, pointi)
502  {
503  if (hits[pointi] == cutFaceGlobal.size())
504  {
505  // Found other side.
506  otherSideFound = true;
507 
508  cfMaster.append
509  (
510  masterFacesOfPZero[pointi]
511  );
512 
513  cfSlave.append(facei);
514 
515  // Reverse the face such that it
516  // points out of the master patch
517  cf.last().flip();
518 
519  if (debug)
520  {
521  Pout<< " other side: "
522  << masterFacesOfPZero[pointi]
523  << endl;
524  }
525  } // end of hits
526  } // end of for all hits
527 
528  } // end of not miss
529 
530  if (!otherSideFound || miss)
531  {
532  if (debug)
533  {
534  Pout<< " solo slave A" << endl;
535  }
536 
537  cfMaster.append(-1);
538  cfSlave.append(facei);
539  }
540  }
541  else
542  {
543  // First point not in master patch
544  if (debug)
545  {
546  Pout<< " solo slave B" << endl;
547  }
548 
549  cfMaster.append(-1);
550  cfSlave.append(facei);
551  }
552  }
553  else
554  {
555  if (debug)
556  {
557  Pout<< " master side" << endl;
558  }
559 
560  cfMaster.append(facei - slavePatch_.size());
561  cfSlave.append(-1);
562  }
563  }
564  else
565  {
566  // Face not completed. Prepare for the next point search
567  prevPointLabel = curPointLabel;
568  curPointLabel = bestAtanPoint;
569  curPoint = lp[curPointLabel];
570  cutFaceGlobalPoints.append(mp[curPointLabel]);
571  cutFaceLocalPoints.append(curPointLabel);
572 
573  if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
574  {
575  faceSizeDebug *= 2;
576 
577  // Check for duplicate points in the face
578  forAll(cutFaceGlobalPoints, checkI)
579  {
580  for
581  (
582  label checkJ = checkI + 1;
583  checkJ < cutFaceGlobalPoints.size();
584  checkJ++
585  )
586  {
587  if
588  (
589  cutFaceGlobalPoints[checkI]
590  == cutFaceGlobalPoints[checkJ]
591  )
592  {
593  // Shrink local points for debugging
594  cutFaceLocalPoints.shrink();
595 
596  face origFace;
597  face origFaceLocal;
598  if (facei < slavePatch_.size())
599  {
600  origFace = slavePatch_[facei];
601  origFaceLocal =
602  slavePatch_.localFaces()[facei];
603  }
604  else
605  {
606  origFace =
607  masterPatch_
608  [facei - slavePatch_.size()];
609 
610  origFaceLocal =
611  masterPatch_.localFaces()
612  [facei - slavePatch_.size()];
613  }
614 
616  << "Duplicate point found in cut face. "
617  << "Error in the face cutting "
618  << "algorithm for global face "
619  << origFace << " local face "
620  << origFaceLocal << nl
621  << "Slave size: " << slavePatch_.size()
622  << " Master size: "
623  << masterPatch_.size()
624  << " index: " << facei << ".\n"
625  << "Face: " << curGlobalFace << nl
626  << "Cut face: " << cutFaceGlobalPoints
627  << " local: " << cutFaceLocalPoints
628  << nl << "Points: "
629  << face(cutFaceLocalPoints).points(lp)
630  << abort(FatalError);
631  }
632  }
633  }
634  } // end of debug
635  }
636  } while (!completedCutFace);
637  } // end of used edges
638 
639  if (debug)
640  {
641  Pout<< " Finished face " << facei << endl;
642  }
643 
644  } // end of local faces
645 
646  // Re-pack the list into compact storage
647  cutFacesPtr_ = new faceList();
648  cutFacesPtr_->transfer(cf);
649 
650  cutFaceMasterPtr_ = new labelList();
651  cutFaceMasterPtr_->transfer(cfMaster);
652 
653  cutFaceSlavePtr_ = new labelList();
654  cutFaceSlavePtr_->transfer(cfSlave);
655 }
656 
657 
658 void Foam::enrichedPatch::clearCutFaces()
659 {
660  deleteDemandDrivenData(cutFacesPtr_);
661  deleteDemandDrivenData(cutFaceMasterPtr_);
662  deleteDemandDrivenData(cutFaceSlavePtr_);
663 }
664 
665 
666 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
667 
669 {
670  if (!cutFacesPtr_)
671  {
672  calcCutFaces();
673  }
674 
675  return *cutFacesPtr_;
676 }
677 
678 
680 {
681  if (!cutFaceMasterPtr_)
682  {
683  calcCutFaces();
684  }
685 
686  return *cutFaceMasterPtr_;
687 }
688 
689 
691 {
692  if (!cutFaceSlavePtr_)
693  {
694  calcCutFaces();
695  }
696 
697  return *cutFaceSlavePtr_;
698 }
699 
700 
701 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#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
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const labelList & cutFaceMaster() const
Return cut face master list.
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const faceList & cutFaces() const
Return list of cut faces.
HashTable< T, label, Hash< label > >::const_iterator const_iterator
Definition: Map.H:59
static const unsigned edgesPerPoint_
Estimated number of edges per point.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const Map< labelList > & masterPointFaces() const
Master point face addressing.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
const faceList & enrichedFaces() const
Return enriched faces.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< edge > edgeList
Definition: edgeList.H:38
const labelList & cutFaceSlave() const
Return cut face slave list.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const labelListList & pointPoints() const
Return point-point addressing.
static const char nl
Definition: Ostream.H:260
dimensionedScalar atan2(const dimensionedScalar &x, const dimensionedScalar &y)
const pointField & localPoints() const
Return local points.
vector point
Point is a vector.
Definition: point.H:41
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const List< FaceType > & localFaces() const
Return patch faces addressing into local point list.
dimensioned< scalar > mag(const dimensioned< Type > &)
const faceList & localFaces() const
Return local faces.
const labelList & meshPoints() const
Return mesh points.
void deleteDemandDrivenData(DataPtr &dataPtr)
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342