FreeFOAM The Cross-Platform CFD Toolkit
enrichedPatchCutFaces.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 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 <OpenFOAM/boolList.H>
32 #include <OpenFOAM/DynamicList.H>
33 #include <OpenFOAM/labelPair.H>
34 #include <OpenFOAM/primitiveMesh.H>
35 #include <OpenFOAM/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  {
54  FatalErrorIn("void enrichedPatch::calcCutFaces() const")
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 << " local: " << curLocalFace << " or " << curGlobalFace << endl;
108 // if (faceI < slavePatch_.size())
109 // {
110 // Pout<< "original slave: " << slavePatch_[faceI]
111 // << " local: " << slavePatch_.localFaces()[faceI] << endl;
112 // }
113 // else
114 // {
115 // Pout<< "original master: "
116 // << masterPatch_[faceI - slavePatch_.size()] << " "
117 // << masterPatch_.localFaces()[faceI - slavePatch_.size()]
118 // << endl;
119 // }
120 // {
121 // pointField facePoints = curLocalFace.points(lp);
122 // forAll (curLocalFace, pointI)
123 // {
124 // Pout << "v " << facePoints[pointI].x() << " "
125 // << facePoints[pointI].y() << " "
126 // << facePoints[pointI].z() << endl;
127 // }
128 // }
129 
130  // Track the usage of face edges. When all edges are used, the
131  // face decomposition is complete.
132  // The seed edges include all the edges of the original face + all edges
133  // of other faces that have been used in the construction of the
134  // facet. Edges from other faces can be considered as
135  // internal to the current face if used only once.
136 
137  // Track the edge usage to avoid duplicate faces and reset it to unused
138  boolList usedFaceEdges(curLocalFace.size(), false);
139 
140  SLList<edge> edgeSeeds;
141 
142  // Insert the edges of current face into the seed list.
143  edgeList cfe = curLocalFace.edges();
144  forAll (curLocalFace, edgeI)
145  {
146  edgeSeeds.append(cfe[edgeI]);
147  }
148 
149  // Grab face normal
150  vector normal = curLocalFace.normal(lp);
151  normal /= mag(normal);
152 
153  while (edgeSeeds.size())
154  {
155 // Pout << "edgeSeeds.size(): " << edgeSeeds.size() << endl;
156  const edge curEdge = edgeSeeds.removeHead();
157 
158  // Locate the edge in current face
159  const label curEdgeWhich = curLocalFace.which(curEdge.start());
160 
161  // Check if the edge is in current face and if it has been
162  // used already. If so, skip it
163  if
164  (
165  curEdgeWhich > -1
166  && curLocalFace.nextLabel(curEdgeWhich) == curEdge.end()
167  )
168  {
169  // Edge is in the starting face.
170  // If unused, mark it as used; if used, skip it
171  if (usedFaceEdges[curEdgeWhich]) continue;
172 
173  usedFaceEdges[curEdgeWhich] = true;
174  }
175 
176  // If the edge has already been used twice, skip it
177  if (edgesUsedTwice.found(curEdge)) continue;
178 // Pout << "Trying new edge (" << mp[curEdge.start()] << ", " << mp[curEdge.end()] << ") seed: " << curEdge << " used: " << edgesUsedTwice.found(curEdge) << endl;
179 
180  // Estimate the size of cut face as twice the size of original face
181  DynamicList<label> cutFaceGlobalPoints(2*curLocalFace.size());
182  DynamicList<label> cutFaceLocalPoints(2*curLocalFace.size());
183 
184  // Found unused edge.
185  label prevPointLabel = curEdge.start();
186  cutFaceGlobalPoints.append(mp[prevPointLabel]);
187  cutFaceLocalPoints.append(prevPointLabel);
188 // Pout << "prevPointLabel: " << mp[prevPointLabel] << endl;
189  // Grab current point and append it to the list
190  label curPointLabel = curEdge.end();
191  point curPoint = lp[curPointLabel];
192  cutFaceGlobalPoints.append(mp[curPointLabel]);
193  cutFaceLocalPoints.append(curPointLabel);
194 
195  bool completedCutFace = false;
196 
197  label faceSizeDebug = maxFaceSizeDebug_;
198 
199  do
200  {
201  // Grab the next point options
202 // Pout << "curPointLabel: " << mp[curPointLabel] << endl;
203  const labelList& nextPoints = pp[curPointLabel];
204 // Pout << "nextPoints: " << UIndirectList<label>(mp, nextPoints) << endl;
205  // Get the vector along the edge and the right vector
206  vector ahead = curPoint - lp[prevPointLabel];
207  ahead -= normal*(normal & ahead);
208  ahead /= mag(ahead);
209 
210  vector right = normal ^ ahead;
211  right /= mag(right);
212 // Pout<< "normal: " << normal << " ahead: " << ahead << " right: " << right << endl;
213 
214  scalar atanTurn = -GREAT;
215  label bestAtanPoint = -1;
216 
217  forAll (nextPoints, nextI)
218  {
219  // Exclude the point we are coming from; there will always
220  // be more than one edge, so this is safe
221  if (nextPoints[nextI] != prevPointLabel)
222  {
223 // Pout << "cur point: " << curPoint << " trying for point: " << mp[nextPoints[nextI]] << " " << lp[nextPoints[nextI]];
224  vector newDir = lp[nextPoints[nextI]] - curPoint;
225 // Pout << " newDir: " << newDir << " mag: " << mag(newDir) << flush;
226  newDir -= normal*(normal & newDir);
227  scalar magNewDir = mag(newDir);
228 // Pout << " corrected: " << newDir << " mag: " << mag(newDir) << flush;
229 
230  if (magNewDir < SMALL)
231  {
233  (
234  "void enrichedPatch::"
235  "calcCutFaces() const"
236  ) << "Zero length edge detected. Probable "
237  << "projection error: slave patch probably "
238  << "does not project onto master. "
239  << "Please switch on "
240  << "enriched patch debug for more info"
241  << abort(FatalError);
242  }
243 
244  newDir /= magNewDir;
245 
246  scalar curAtanTurn =
247  atan2(newDir & right, newDir & ahead);
248 
249 // Pout << " atan: " << curAtanTurn << endl;
250 
251  if (curAtanTurn > atanTurn)
252  {
253  bestAtanPoint = nextPoints[nextI];
254  atanTurn = curAtanTurn;
255  }
256  } // end of prev point skip
257  } // end of next point selection
258 // Pout<< " bestAtanPoint: " << bestAtanPoint << " or "
259 // << mp[bestAtanPoint] << endl;
260  // Selected next best point.
261 // Pout << "cutFaceGlobalPoints: " << cutFaceGlobalPoints << endl;
262  // Check if the edge about to be added has been used
263  // in the current face or twice in other faces. If
264  // so, the face is bad.
265  const label whichNextPoint = curLocalFace.which(curPointLabel);
266 
267  if
268  (
269  whichNextPoint > -1
270  && curLocalFace.nextLabel(whichNextPoint) == bestAtanPoint
271  && usedFaceEdges[whichNextPoint]
272  )
273  {
274  // This edge is already used in current face
275  // face cannot be good; start on a new one
276 // Pout << "Double usage in current face, cannot be good" << endl;
277  completedCutFace = true;
278  }
279 
280 
281  if (edgesUsedTwice.found(edge(curPointLabel, bestAtanPoint)))
282  {
283  // This edge is already used -
284  // face cannot be good; start on a new one
285  completedCutFace = true;
286 // Pout << "Double usage elsewhere, cannot be good" << endl;
287  }
288 
289  if (completedCutFace) continue;
290 
291  // Insert the next best point into the list
292  if (mp[bestAtanPoint] == cutFaceGlobalPoints[0])
293  {
294  // Face is completed. Save it and move on to the next
295  // available edge
296  completedCutFace = true;
297 
298  if (debug)
299  {
300  Pout<< " local: " << cutFaceLocalPoints
301  << " one side: " << faceI;
302  }
303 
304  // Append the face
305  face cutFaceGlobal;
306  cutFaceGlobal.transfer(cutFaceGlobalPoints);
307 
308  cf.append(cutFaceGlobal);
309 
310  face cutFaceLocal;
311  cutFaceLocal.transfer(cutFaceLocalPoints);
312 // Pout << "\ncutFaceLocal: " << cutFaceLocal.points(lp) << endl;
313  // Go through all edges of the cut faces.
314  // If the edge corresponds to a starting face edge,
315  // mark the starting face edge as true
316 
317  forAll (cutFaceLocal, cutI)
318  {
319  const edge curCutFaceEdge
320  (
321  cutFaceLocal[cutI],
322  cutFaceLocal.nextLabel(cutI)
323  );
324 
325  // Increment the usage count using two hash sets
326  HashSet<edge, Hash<edge> >::iterator euoIter =
327  edgesUsedOnce.find(curCutFaceEdge);
328 
329  if (euoIter == edgesUsedOnce.end())
330  {
331 // Pout << "Found edge not used before: "<< curCutFaceEdge << endl;
332  edgesUsedOnce.insert(curCutFaceEdge);
333  }
334  else
335  {
336 // Pout << "Found edge used once: " << curCutFaceEdge << endl;
337  edgesUsedOnce.erase(euoIter);
338  edgesUsedTwice.insert(curCutFaceEdge);
339  }
340 
341  const label curCutFaceEdgeWhich = curLocalFace.which(curCutFaceEdge.start());
342 
343  if
344  (
345  curCutFaceEdgeWhich > -1
346  && curLocalFace.nextLabel(curCutFaceEdgeWhich) == curCutFaceEdge.end()
347  )
348  {
349  // Found edge in original face
350 // Pout << "Found edge in orig face: " << curCutFaceEdge << ": " << curCutFaceEdgeWhich << endl;
351  usedFaceEdges[curCutFaceEdgeWhich] = true;
352  }
353  else
354  {
355  // Edge not in original face. Add it to seeds
356 // Pout << "Found new edge seed: " << curCutFaceEdge << endl;
357  edgeSeeds.append(curCutFaceEdge.reverseEdge());
358  }
359  }
360 
361 
362  // Find out what the other side is
363 
364  // Algorithm
365  // If the face is in the slave half of the
366  // enrichedFaces lists, it may be matched against
367  // the master face. It can be recognised by the
368  // fact that all its points belong to one master
369  // face. For this purpose:
370  // 1) Grab the master faces around the point zero
371  // of the cut face and collect all master faces
372  // around it.
373  // 2) Loop through the rest of cut face points and
374  // if the point refers to one of the faces marked
375  // by point zero, increment its count.
376  // 3) When completed, the face whose count is
377  // equal to the number of points in the cut face
378  // is the other side. If this is not the case, there is no
379  // face on the other side.
380 
381  if (faceI < slavePatch_.size())
382  {
383  Map<labelList>::const_iterator mpfAddrIter =
384  masterPointFaceAddr.find(cutFaceGlobal[0]);
385 
386  bool otherSideFound = false;
387 
388  if (mpfAddrIter != masterPointFaceAddr.end())
389  {
390  bool miss = false;
391 
392  // Create the label count using point zero
393  const labelList& masterFacesOfPZero = mpfAddrIter();
394 
395  labelList hits(masterFacesOfPZero.size(), 1);
396 
397  for
398  (
399  label pointI = 1;
400  pointI < cutFaceGlobal.size();
401  pointI++
402  )
403  {
405  mpfAddrPointIter =
406  masterPointFaceAddr.find
407  (
408  cutFaceGlobal[pointI]
409  );
410 
411  if
412  (
413  mpfAddrPointIter
414  == masterPointFaceAddr.end()
415  )
416  {
417  // Point is off the master patch. Skip
418  miss = true;
419  break;
420  }
421 
422  const labelList& curMasterFaces =
423  mpfAddrPointIter();
424 
425  // For every current face, try to find it in the
426  // zero-list
427  forAll (curMasterFaces, i)
428  {
429  forAll (masterFacesOfPZero, j)
430  {
431  if
432  (
433  curMasterFaces[i]
434  == masterFacesOfPZero[j]
435  )
436  {
437  hits[j]++;
438  break;
439  }
440  }
441  }
442  }
443 
444  // If all point are found attempt matching
445  if (!miss)
446  {
447  forAll (hits, pointI)
448  {
449  if (hits[pointI] == cutFaceGlobal.size())
450  {
451  // Found other side.
452  otherSideFound = true;
453 
454  cfMaster.append
455  (
456  masterFacesOfPZero[pointI]
457  );
458 
459  cfSlave.append(faceI);
460 
461  // Reverse the face such that it
462  // points out of the master patch
463  cf[cf.size() - 1] =
464  cf[cf.size() - 1].reverseFace();
465 
466  if (debug)
467  {
468  Pout<< " other side: "
469  << masterFacesOfPZero[pointI]
470  << endl;
471  }
472  } // end of hits
473  } // end of for all hits
474 
475  } // end of not miss
476 
477  if (!otherSideFound || miss)
478  {
479  if (debug)
480  {
481  Pout << " solo slave A" << endl;
482  }
483 
484  cfMaster.append(-1);
485  cfSlave.append(faceI);
486  }
487  }
488  else
489  {
490  // First point not in master patch
491  if (debug)
492  {
493  Pout << " solo slave B" << endl;
494  }
495 
496  cfMaster.append(-1);
497  cfSlave.append(faceI);
498  }
499  }
500  else
501  {
502  if (debug)
503  {
504  Pout << " master side" << endl;
505  }
506 
507  cfMaster.append(faceI - slavePatch_.size());
508  cfSlave.append(-1);
509  }
510  }
511  else
512  {
513  // Face not completed. Prepare for the next point search
514  prevPointLabel = curPointLabel;
515  curPointLabel = bestAtanPoint;
516  curPoint = lp[curPointLabel];
517  cutFaceGlobalPoints.append(mp[curPointLabel]);
518  cutFaceLocalPoints.append(curPointLabel);
519 
520  if (debug || cutFaceGlobalPoints.size() > faceSizeDebug)
521  {
522  faceSizeDebug *= 2;
523 
524  // Check for duplicate points in the face
525  forAll (cutFaceGlobalPoints, checkI)
526  {
527  for
528  (
529  label checkJ = checkI + 1;
530  checkJ < cutFaceGlobalPoints.size();
531  checkJ++
532  )
533  {
534  if
535  (
536  cutFaceGlobalPoints[checkI]
537  == cutFaceGlobalPoints[checkJ]
538  )
539  {
540  // Shrink local points for debugging
541  cutFaceLocalPoints.shrink();
542 
543  face origFace;
544  face origFaceLocal;
545  if (faceI < slavePatch_.size())
546  {
547  origFace = slavePatch_[faceI];
548  origFaceLocal =
549  slavePatch_.localFaces()[faceI];
550  }
551  else
552  {
553  origFace =
554  masterPatch_
555  [faceI - slavePatch_.size()];
556 
557  origFaceLocal =
558  masterPatch_.localFaces()
559  [faceI - slavePatch_.size()];
560  }
561 
563  (
564  "void enrichedPatch::"
565  "calcCutFaces() const"
566  ) << "Duplicate point found in cut face. "
567  << "Error in the face cutting "
568  << "algorithm for global face "
569  << origFace << " local face "
570  << origFaceLocal << nl
571  << "Slave size: " << slavePatch_.size()
572  << " Master size: "
573  << masterPatch_.size()
574  << " index: " << faceI << ".\n"
575  << "Face: " << curGlobalFace << nl
576  << "Cut face: " << cutFaceGlobalPoints
577  << " local: " << cutFaceLocalPoints
578  << nl << "Points: "
579  << face(cutFaceLocalPoints).points(lp)
580  << abort(FatalError);
581  }
582  }
583  }
584  } // end of debug
585  }
586  } while (!completedCutFace);
587  } // end of used edges
588 
589  if (debug)
590  {
591  Pout << " Finished face " << faceI << endl;
592  }
593 
594  } // end of local faces
595 
596  // Re-pack the list into compact storage
597  cutFacesPtr_ = new faceList();
598  cutFacesPtr_->transfer(cf);
599 
600  cutFaceMasterPtr_ = new labelList();
601  cutFaceMasterPtr_->transfer(cfMaster);
602 
603  cutFaceSlavePtr_ = new labelList();
604  cutFaceSlavePtr_->transfer(cfSlave);
605 }
606 
607 
608 void Foam::enrichedPatch::clearCutFaces()
609 {
610  deleteDemandDrivenData(cutFacesPtr_);
611  deleteDemandDrivenData(cutFaceMasterPtr_);
612  deleteDemandDrivenData(cutFaceSlavePtr_);
613 }
614 
615 
616 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
617 
619 {
620  if (!cutFacesPtr_)
621  {
622  calcCutFaces();
623  }
624 
625  return *cutFacesPtr_;
626 }
627 
628 
630 {
631  if (!cutFaceMasterPtr_)
632  {
633  calcCutFaces();
634  }
635 
636  return *cutFaceMasterPtr_;
637 }
638 
639 
641 {
642  if (!cutFaceSlavePtr_)
643  {
644  calcCutFaces();
645  }
646 
647  return *cutFaceSlavePtr_;
648 }
649 
650 
651 // ************************ vim: set sw=4 sts=4 et: ************************ //