FreeFOAM The Cross-Platform CFD Toolkit
topoCellLooper.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 \*---------------------------------------------------------------------------*/
25 
26 #include "topoCellLooper.H"
27 #include <meshTools/cellFeatures.H>
28 #include <OpenFOAM/polyMesh.H>
30 #include <OpenFOAM/DynamicList.H>
31 #include <OpenFOAM/ListOps.H>
32 #include <meshTools/meshTools.H>
33 #include <OpenFOAM/hexMatcher.H>
34 
36 
37 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 namespace Foam
41 {
42  defineTypeNameAndDebug(topoCellLooper, 0);
43  addToRunTimeSelectionTable(cellLooper, topoCellLooper, word);
44 }
45 
46 // Angle for polys to be considered splitHexes.
47 const Foam::scalar Foam::topoCellLooper::featureCos =
48  Foam::cos(10.0 * mathematicalConstant::pi/180.0);
49 
50 
51 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
52 
53 // In-memory truncate a list
54 template <class T>
55 void Foam::topoCellLooper::subsetList
56 (
57  const label startI,
58  const label freeI,
59  DynamicList<T>& lst
60 )
61 {
62  if (startI == 0)
63  {
64  // Truncate (setSize decides itself not to do anything if nothing
65  // changed)
66  if (freeI < 0)
67  {
68  FatalErrorIn("topoCellLooper::subsetList")
69  << "startI:" << startI << " freeI:" << freeI
70  << " lst:" << lst << abort(FatalError);
71  }
72  lst.setCapacity(freeI);
73  }
74  else
75  {
76  // Shift elements down
77  label newI = 0;
78  for (label elemI = startI; elemI < freeI; elemI++)
79  {
80  lst[newI++] = lst[elemI];
81  }
82 
83  if ((freeI - startI) < 0)
84  {
85  FatalErrorIn("topoCellLooper::subsetList")
86  << "startI:" << startI << " freeI:" << freeI
87  << " lst:" << lst << abort(FatalError);
88  }
89 
90  lst.setCapacity(freeI - startI);
91  }
92 }
93 
94 
95 // Returns label of edge nFeaturePts away from startEdge (in the direction of
96 // startVertI) and not counting non-featurePoints.
97 //
98 // When stepping to this face it can happen in 3 different ways:
99 //
100 // --|------
101 // |
102 // 1 | 0
103 // |A
104 // |
105 // |
106 // --|------
107 //
108 // A: jumping from face0 to face1 across edge A.
109 // startEdge != -1
110 // startVert == -1
111 //
112 // --|------
113 // |
114 // 1 | 0
115 // +B
116 // |
117 // |
118 // --|------
119 //
120 // B: jumping from face0 to face1 across (non-feature) point B
121 // startEdge == -1
122 // startVert != -1
123 //
124 // --|------
125 // 0 | 1
126 // |C
127 // --+
128 // |
129 // |
130 // --|------
131 //
132 // C: jumping from face0 to face1 across (feature) point C on edge.
133 // startEdge != -1
134 // startVert != -1
135 //
137 (
138  const cellFeatures& features,
139  const label faceI,
140  const label startEdgeI,
141  const label startVertI,
142  const label nFeaturePts,
143 
144  label& edgeI,
145  label& vertI
146 ) const
147 {
148  const labelList& fEdges = mesh().faceEdges()[faceI];
149 
150  edgeI = startEdgeI;
151 
152  vertI = startVertI;
153 
154  // Number of feature points crossed so far
155  label nVisited = 0;
156 
157  if (vertI == -1)
158  {
159  // Started on edge. Go to one of its endpoints.
160  vertI = mesh().edges()[edgeI].start();
161 
162  if (features.isFeatureVertex(faceI, vertI))
163  {
164  nVisited++;
165  }
166  }
167 
168  if ((edgeI == -1) || !meshTools::edgeOnFace(mesh(), faceI, edgeI))
169  {
170  // Either edge is not set or not on current face. Just take one of
171  // the edges on this face as starting edge.
172  edgeI = getFirstVertEdge(faceI, vertI);
173  }
174 
175  // Now we should have starting edge on face and a vertex on that edge.
176 
177  do
178  {
179 
180  edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
181 
182  if (nVisited == nFeaturePts)
183  {
184  break;
185  }
186 
187  vertI = mesh().edges()[edgeI].otherVertex(vertI);
188 
189  if (features.isFeatureVertex(faceI, vertI))
190  {
191  nVisited++;
192  }
193  }
194  while (true);
195 }
196 
197 
198 // Returns list of vertices on 'superEdge' i.e. list of edges connected by
199 // non-feature points. First and last are feature points, ones inbetween are
200 // not.
201 Foam::labelList Foam::topoCellLooper::getSuperEdge
202 (
203  const cellFeatures& features,
204  const label faceI,
205  const label startEdgeI,
206  const label startVertI
207 ) const
208 {
209  const labelList& fEdges = mesh().faceEdges()[faceI];
210 
211  labelList superVerts(fEdges.size());
212  label superVertI = 0;
213 
214 
215  label edgeI = startEdgeI;
216 
217  label vertI = startVertI;
218 
219  superVerts[superVertI++] = vertI;
220 
221  label prevEdgeI = -1;
222 
223  do
224  {
225  vertI = mesh().edges()[edgeI].otherVertex(vertI);
226 
227  superVerts[superVertI++] = vertI;
228 
229  prevEdgeI = edgeI;
230 
231  edgeI = meshTools::otherEdge(mesh(), fEdges, edgeI, vertI);
232  }
233  while (!features.isFeaturePoint(prevEdgeI, edgeI));
234 
235  superVerts.setSize(superVertI);
236 
237  return superVerts;
238 }
239 
240 
241 // Return non-feature edge from cells' edges that is most perpendicular
242 // to refinement direction.
243 Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge
244 (
245  const vector& refDir,
246  const label cellI,
247  const cellFeatures& features
248 ) const
249 {
250  const labelList& cEdges = mesh().cellEdges()[cellI];
251 
252  const point& ctr = mesh().cellCentres()[cellI];
253 
254  label cutEdgeI = -1;
255  scalar maxCos = -GREAT;
256 
257  forAll(cEdges, cEdgeI)
258  {
259  label edgeI = cEdges[cEdgeI];
260 
261  if (!features.isFeatureEdge(edgeI))
262  {
263  const edge& e = mesh().edges()[edgeI];
264 
265  // Get plane spanned by e.start, e.end and cell centre.
266  vector e0 = mesh().points()[e.start()] - ctr;
267  vector e1 = mesh().points()[e.end()] - ctr;
268 
269  vector n = e0 ^ e1;
270  n /= mag(n);
271 
272  scalar cosAngle = mag(refDir & n);
273 
274  if (cosAngle > maxCos)
275  {
276  maxCos = cosAngle;
277 
278  cutEdgeI = edgeI;
279  }
280  }
281  }
282 
283  return cutEdgeI;
284 }
285 
286 
287 // Starts from edge and vertex on edge on face (or neighbouring face)
288 // and steps either to existing vertex (vertI != -1) or to edge (vertI == -1)
289 // by walking point-edge and crossing nFeats featurePoints.
290 void Foam::topoCellLooper::walkAcrossFace
291 (
292  const cellFeatures& features,
293  const label faceI,
294  const label startEdgeI,
295  const label startVertI,
296  const label nFeats,
297 
298  label& edgeI,
299  label& vertI
300 ) const
301 {
302  label oppositeVertI = -1;
303  label oppositeEdgeI = -1;
304 
305  // Go to oppositeEdge and a vertex on that.
306  walkFace
307  (
308  features,
309  faceI,
310  startEdgeI,
311  startVertI,
312  nFeats,
313 
314  oppositeEdgeI,
315  oppositeVertI
316  );
317 
318  // Loop over super edge to find internal points if there are any.
319 
320  labelList superEdge =
321  getSuperEdge
322  (
323  features,
324  faceI,
325  oppositeEdgeI,
326  oppositeVertI
327  );
328 
329  label sz = superEdge.size();
330 
331  if (sz == 2)
332  {
333  // No non-feature point inbetween feature points.
334  // Mark edge.
335 
336  vertI = -1;
337  edgeI = oppositeEdgeI;
338  }
339  else if (sz == 3)
340  {
341  vertI = superEdge[1];
342  edgeI = -1;
343  }
344  else
345  {
346  //Should choose acc. to geometry!
347  label index = sz/2;
348 
349  if (debug)
350  {
351  Pout<< " Don't know what to do. Stepped to non-feature point "
352  << "at index " << index << " in superEdge:" << superEdge
353  << endl;
354  }
355 
356  vertI = superEdge[index];
357  edgeI = -1;
358  }
359 }
360 
361 
362 // Walks cell circumference. Updates face, edge, vertex.
363 //
364 // Position on face is given by:
365 //
366 // vertI == -1, faceI != -1, edgeI != -1
367 // on edge of face. Cross edge to neighbouring face.
368 //
369 // vertI != -1, edgeI != -1, faceI == -1
370 // coming from edge onto vertex vertI. Need to step to one
371 // of the faces not using edgeI.
372 //
373 // vertI != -1, edgeI == -1, faceI != -1
374 // coming from vertex on side of face. Step to one of the faces
375 // using vertI but not faceI
376 void Foam::topoCellLooper::walkSplitHex
377 (
378  const label cellI,
379  const cellFeatures& features,
380  const label fromFaceI,
381  const label fromEdgeI,
382  const label fromVertI,
383 
384  DynamicList<label>& loop,
385  DynamicList<scalar>& loopWeights
386 ) const
387 {
388  // Work vars giving position on cell
389  label faceI = fromFaceI;
390  label edgeI = fromEdgeI;
391  label vertI = fromVertI;
392 
393  do
394  {
395  if (debug)
396  {
397  Pout<< "Entering walk with : cell:" << cellI << " face:" << faceI;
398  if (faceI != -1)
399  {
400  Pout<< " verts:" << mesh().faces()[faceI];
401  }
402  Pout<< " edge:" << edgeI;
403  if (edgeI != -1)
404  {
405  Pout<< " verts:" << mesh().edges()[edgeI];
406  }
407  Pout<< " vert:" << vertI << endl;
408  }
409 
410  label startLoop = -1;
411 
412  if
413  (
414  (vertI != -1)
415  && (
416  (startLoop =
417  findIndex
418  (
419  loop,
420  vertToEVert(vertI)
421  )
422  )
423  != -1
424  )
425  )
426  {
427  // Breaking walk since vertI already cut
428  label firstFree = loop.size();
429 
430  subsetList(startLoop, firstFree, loop);
431  subsetList(startLoop, firstFree, loopWeights);
432 
433  break;
434  }
435  if
436  (
437  (edgeI != -1)
438  && (
439  (startLoop =
440  findIndex
441  (
442  loop,
443  edgeToEVert(edgeI)
444  )
445  )
446  != -1
447  )
448  )
449  {
450  // Breaking walk since edgeI already cut
451  label firstFree = loop.size();
452 
453  subsetList(startLoop, firstFree, loop);
454  subsetList(startLoop, firstFree, loopWeights);
455 
456  break;
457  }
458 
459 
460  if (vertI == -1)
461  {
462  // On edge
463  if (edgeI == -1)
464  {
465  FatalErrorIn("walkSplitHex") << "Illegal edge and vert"
466  << abort(FatalError);
467  }
468 
469  loop.append(edgeToEVert(edgeI));
470  loopWeights.append(0.5);
471 
472  // Cross edge to next face
473  faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
474 
475  if (debug)
476  {
477  Pout<< " stepped across edge " << mesh().edges()[edgeI]
478  << " to face " << faceI << " verts:"
479  << mesh().faces()[faceI] << endl;
480  }
481 
482  label nextEdgeI = -1;
483  label nextVertI = -1;
484 
485  // Walk face along its edges
486  walkAcrossFace
487  (
488  features,
489  faceI,
490  edgeI,
491  vertI,
492  2,
493 
494  nextEdgeI,
495  nextVertI
496  );
497 
498  edgeI = nextEdgeI;
499  vertI = nextVertI;
500  }
501  else
502  {
503  // On vertex.
504 
505  loop.append(vertToEVert(vertI));
506  loopWeights.append(-GREAT);
507 
508  if (edgeI == -1)
509  {
510  // Normal vertex on edge of face. Get edges connected to it
511  // which are not on faceI.
512  labelList nextEdges = getVertEdgesNonFace
513  (
514  cellI,
515  faceI,
516  vertI
517  );
518 
519  if (nextEdges.empty())
520  {
521  // Cross to other face (there is only one since no edges)
522  const labelList& pFaces = mesh().pointFaces()[vertI];
523 
524  forAll(pFaces, pFaceI)
525  {
526  label thisFaceI = pFaces[pFaceI];
527 
528  if
529  (
530  (thisFaceI != faceI)
531  && meshTools::faceOnCell(mesh(), cellI, thisFaceI)
532  )
533  {
534  faceI = thisFaceI;
535  break;
536  }
537  }
538 
539  if (debug)
540  {
541  Pout<< " stepped from non-edge vertex " << vertI
542  << " to face " << faceI << " verts:"
543  << mesh().faces()[faceI]
544  << " since candidate edges:" << nextEdges << endl;
545  }
546 
547  label nextEdgeI = -1;
548  label nextVertI = -1;
549 
550  walkAcrossFace
551  (
552  features,
553  faceI,
554  edgeI,
555  vertI,
556  2, // 2 vertices to cross
557 
558  nextEdgeI,
559  nextVertI
560  );
561 
562  edgeI = nextEdgeI;
563  vertI = nextVertI;
564  }
565  else if (nextEdges.size() == 1)
566  {
567  // One edge. Go along this one.
568  edgeI = nextEdges[0];
569 
570  if (debug)
571  {
572  Pout<< " stepped from non-edge vertex " << vertI
573  << " along edge " << edgeI << " verts:"
574  << mesh().edges()[edgeI]
575  << " out of candidate edges:"
576  << nextEdges << endl;
577  }
578 
579  vertI = mesh().edges()[edgeI].otherVertex(vertI);
580 
581  faceI = -1;
582  }
583  else
584  {
585  // Multiple edges to choose from. Pick any one.
586  // (ideally should be geometric)
587 
588  label index = nextEdges.size()/2;
589 
590  edgeI = nextEdges[index];
591 
592  if (debug)
593  {
594  Pout<< " stepped from non-edge vertex " << vertI
595  << " along edge " << edgeI << " verts:"
596  << mesh().edges()[edgeI]
597  << " out of candidate edges:" << nextEdges << endl;
598  }
599 
600  vertI = mesh().edges()[edgeI].otherVertex(vertI);
601 
602  faceI = -1;
603  }
604  }
605  else
606  {
607  // Get faces connected to startVertI but not startEdgeI
608  labelList nextFaces =
609  getVertFacesNonEdge
610  (
611  cellI,
612  edgeI,
613  vertI
614  );
615 
616  if (nextFaces.size() == 1)
617  {
618  // Only one face to cross.
619  faceI = nextFaces[0];
620 
621  label nextEdgeI = -1;
622  label nextVertI = -1;
623 
624  walkAcrossFace
625  (
626  features,
627  faceI,
628  edgeI,
629  vertI,
630  2, // 2 vertices to cross
631 
632  nextEdgeI,
633  nextVertI
634  );
635 
636  edgeI = nextEdgeI;
637  vertI = nextVertI;
638  }
639  else if (nextFaces.size() == 2)
640  {
641  // Split face. Get edge inbetween.
642  faceI = -1;
643 
644  edgeI =
646  (
647  mesh(),
648  nextFaces[0],
649  nextFaces[1]
650  );
651 
652  vertI = mesh().edges()[edgeI].otherVertex(vertI);
653  }
654  else
655  {
656  FatalErrorIn("walkFromVert") << "Not yet implemented"
657  << "Choosing from more than "
658  << "two candidates:" << nextFaces
659  << " when coming from vertex " << vertI << " on cell "
660  << cellI << abort(FatalError);
661  }
662  }
663  }
664 
665  if (debug)
666  {
667  Pout<< "Walked to : face:" << faceI;
668  if (faceI != -1)
669  {
670  Pout<< " verts:" << mesh().faces()[faceI];
671  }
672  Pout<< " edge:" << edgeI;
673  if (edgeI != -1)
674  {
675  Pout<< " verts:" << mesh().edges()[edgeI];
676  }
677  Pout<< " vert:" << vertI << endl;
678  }
679  }
680  while (true);
681 }
682 
683 
684 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
685 
686 // Construct from components
687 Foam::topoCellLooper::topoCellLooper(const polyMesh& mesh)
688 :
689  hexCellLooper(mesh)
690 {}
691 
692 
693 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
694 
696 {}
697 
698 
699 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
700 
702 (
703  const vector& refDir,
704  const label cellI,
705  const boolList& vertIsCut,
706  const boolList& edgeIsCut,
707  const scalarField& edgeWeight,
708 
709  labelList& loop,
710  scalarField& loopWeights
711 ) const
712 {
713  if (mesh().cellShapes()[cellI].model() == hex_)
714  {
715  // Let parent handle hex case.
716  return
718  (
719  refDir,
720  cellI,
721  vertIsCut,
722  edgeIsCut,
723  edgeWeight,
724  loop,
725  loopWeights
726  );
727  }
728  else
729  {
730  cellFeatures superCell(mesh(), featureCos, cellI);
731 
732  if (hexMatcher().isA(superCell.faces()))
733  {
734  label edgeI =
735  getAlignedNonFeatureEdge
736  (
737  refDir,
738  cellI,
739  superCell
740  );
741 
742  label vertI = -1;
743 
744  label faceI = -1;
745 
746  if (edgeI != -1)
747  {
748  // Found non-feature edge. Start walking from vertex on edge.
749  vertI = mesh().edges()[edgeI].start();
750  }
751  else
752  {
753  // No 'matching' non-feature edge found on cell. Get starting
754  // normal i.e. feature edge.
755  edgeI = getMisAlignedEdge(refDir, cellI);
756 
757  // Get any face using edge
758  label face0;
759  label face1;
760  meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
761 
762  faceI = face0;
763  }
764 
765  label nEstCuts = 2*mesh().cells()[cellI].size();
766 
767  DynamicList<label> localLoop(nEstCuts);
768  DynamicList<scalar> localLoopWeights(nEstCuts);
769 
770  walkSplitHex
771  (
772  cellI,
773  superCell,
774  faceI,
775  edgeI,
776  vertI,
777 
778  localLoop,
779  localLoopWeights
780  );
781 
782  if (localLoop.size() <=2)
783  {
784  return false;
785  }
786  else
787  {
788  loop.transfer(localLoop);
789  loopWeights.transfer(localLoopWeights);
790 
791  return true;
792  }
793  }
794  else
795  {
796  // Let parent handle poly case.
797  return hexCellLooper::cut
798  (
799  refDir,
800  cellI,
801  vertIsCut,
802  edgeIsCut,
803  edgeWeight,
804  loop,
805  loopWeights
806  );
807  }
808  }
809 }
810 
811 
813 (
814  const plane& cutPlane,
815  const label cellI,
816  const boolList& vertIsCut,
817  const boolList& edgeIsCut,
818  const scalarField& edgeWeight,
819 
820  labelList& loop,
821  scalarField& loopWeights
822 ) const
823 {
824  // Let parent handle cut with plane.
825  return
827  (
828  cutPlane,
829  cellI,
830  vertIsCut,
831  edgeIsCut,
832  edgeWeight,
833 
834  loop,
835  loopWeights
836  );
837 }
838 
839 
840 // ************************ vim: set sw=4 sts=4 et: ************************ //