FreeFOAM The Cross-Platform CFD Toolkit
boundaryCutter.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 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "boundaryCutter.H"
29 #include <OpenFOAM/polyMesh.H>
36 #include <OpenFOAM/mapPolyMesh.H>
37 #include <meshTools/meshTools.H>
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43 
44 defineTypeNameAndDebug(boundaryCutter, 0);
45 
46 }
47 
48 
49 // * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * //
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 void Foam::boundaryCutter::getFaceInfo
55 (
56  const label faceI,
57  label& patchID,
58  label& zoneID,
59  label& zoneFlip
60 ) const
61 {
62  patchID = -1;
63 
64  if (!mesh_.isInternalFace(faceI))
65  {
66  patchID = mesh_.boundaryMesh().whichPatch(faceI);
67  }
68 
69  zoneID = mesh_.faceZones().whichZone(faceI);
70 
71  zoneFlip = false;
72 
73  if (zoneID >= 0)
74  {
75  const faceZone& fZone = mesh_.faceZones()[zoneID];
76 
77  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
78  }
79 }
80 
81 
82 // Adds additional vertices (from edge cutting) to face. Used for faces which
83 // are not split but still might use edge that has been cut.
84 Foam::face Foam::boundaryCutter::addEdgeCutsToFace
85 (
86  const label faceI,
87  const Map<labelList>& edgeToAddedPoints
88 ) const
89 {
90  const edgeList& edges = mesh_.edges();
91  const face& f = mesh_.faces()[faceI];
92  const labelList& fEdges = mesh_.faceEdges()[faceI];
93 
94  // Storage for face
95  DynamicList<label> newFace(2 * f.size());
96 
97  forAll(f, fp)
98  {
99  // Duplicate face vertex .
100  newFace.append(f[fp]);
101 
102  // Check if edge has been cut.
103  label v1 = f.nextLabel(fp);
104 
105  label edgeI = meshTools::findEdge(edges, fEdges, f[fp], v1);
106 
107  Map<labelList>::const_iterator fnd = edgeToAddedPoints.find(edgeI);
108 
109  if (fnd != edgeToAddedPoints.end())
110  {
111  // edge has been cut. Introduce new vertices. Check order.
112  const labelList& addedPoints = fnd();
113 
114  if (edges[edgeI].start() == f[fp])
115  {
116  // Introduce in same order.
117  forAll(addedPoints, i)
118  {
119  newFace.append(addedPoints[i]);
120  }
121  }
122  else
123  {
124  // Introduce in opposite order.
125  forAllReverse(addedPoints, i)
126  {
127  newFace.append(addedPoints[i]);
128  }
129  }
130  }
131  }
132 
133  face returnFace;
134  returnFace.transfer(newFace);
135 
136  if (debug)
137  {
138  Pout<< "addEdgeCutsToFace:" << nl
139  << " from : " << f << nl
140  << " to : " << returnFace << endl;
141  }
142 
143  return returnFace;
144 }
145 
146 
147 void Foam::boundaryCutter::addFace
148 (
149  const label faceI,
150  const face& newFace,
151 
152  bool& modifiedFace, // have we already 'used' faceI
153  polyTopoChange& meshMod
154 ) const
155 {
156  // Information about old face
157  label patchID, zoneID, zoneFlip;
158  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
159  label own = mesh_.faceOwner()[faceI];
160  label masterPoint = mesh_.faces()[faceI][0];
161 
162  if (!modifiedFace)
163  {
164  meshMod.setAction
165  (
166  polyModifyFace
167  (
168  newFace, // face
169  faceI,
170  own, // owner
171  -1, // neighbour
172  false, // flux flip
173  patchID, // patch for face
174  false, // remove from zone
175  zoneID, // zone for face
176  zoneFlip // face zone flip
177  )
178  );
179 
180  modifiedFace = true;
181  }
182  else
183  {
184  meshMod.setAction
185  (
186  polyAddFace
187  (
188  newFace, // face
189  own, // owner
190  -1, // neighbour
191  masterPoint, // master point
192  -1, // master edge
193  -1, // master face for addition
194  false, // flux flip
195  patchID, // patch for face
196  zoneID, // zone for face
197  zoneFlip // face zone flip
198  )
199  );
200  }
201 }
202 
203 
204 
205 // Splits a face using the cut edges and modified points
206 bool Foam::boundaryCutter::splitFace
207 (
208  const label faceI,
209  const Map<point>& pointToPos,
210  const Map<labelList>& edgeToAddedPoints,
211  polyTopoChange& meshMod
212 ) const
213 {
214  const edgeList& edges = mesh_.edges();
215  const face& f = mesh_.faces()[faceI];
216  const labelList& fEdges = mesh_.faceEdges()[faceI];
217 
218  // Count number of split edges and total number of splits.
219  label nSplitEdges = 0;
220  label nModPoints = 0;
221  label nTotalSplits = 0;
222 
223  forAll(f, fp)
224  {
225  if (pointToPos.found(f[fp]))
226  {
227  nModPoints++;
228  nTotalSplits++;
229  }
230 
231  // Check if edge has been cut.
232  label nextV = f.nextLabel(fp);
233 
234  label edgeI = meshTools::findEdge(edges, fEdges, f[fp], nextV);
235 
236  Map<labelList>::const_iterator fnd = edgeToAddedPoints.find(edgeI);
237 
238  if (fnd != edgeToAddedPoints.end())
239  {
240  nSplitEdges++;
241  nTotalSplits += fnd().size();
242  }
243  }
244 
245  if (debug)
246  {
247  Pout<< "Face:" << faceI
248  << " nModPoints:" << nModPoints
249  << " nSplitEdges:" << nSplitEdges
250  << " nTotalSplits:" << nTotalSplits << endl;
251  }
252 
253  if (nSplitEdges == 0 && nModPoints == 0)
254  {
255  FatalErrorIn("boundaryCutter::splitFace") << "Problem : face:" << faceI
256  << " nSplitEdges:" << nSplitEdges
257  << " nTotalSplits:" << nTotalSplits
258  << abort(FatalError);
259  return false;
260  }
261  else if (nSplitEdges + nModPoints == 1)
262  {
263  // single or multiple cuts on a single edge or single modified point
264  // Dont cut and let caller handle this.
265  Warning << "Face " << faceI << " has only one edge cut " << endl;
266  return false;
267  }
268  else
269  {
270  // So guaranteed to have two edges cut or points modified. Split face:
271  // - find starting cut
272  // - walk to next cut. Make face
273  // - loop until face done.
274 
275  // Information about old face
276  label patchID, zoneID, zoneFlip;
277  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
278 
279  // Get face with new points on cut edges for ease of looping
280  face extendedFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
281 
282  // Find first added point. This is the starting vertex for splitting.
283  label startFp = -1;
284 
285  forAll(extendedFace, fp)
286  {
287  if (extendedFace[fp] >= mesh_.nPoints())
288  {
289  startFp = fp;
290  break;
291  }
292  }
293 
294  if (startFp == -1)
295  {
296  // No added point. Maybe there is a modified point?
297  forAll(extendedFace, fp)
298  {
299  if (pointToPos.found(extendedFace[fp]))
300  {
301  startFp = fp;
302  break;
303  }
304  }
305  }
306 
307  if (startFp == -1)
308  {
309  FatalErrorIn("boundaryCutter::splitFace")
310  << "Problem" << abort(FatalError);
311  }
312 
313  // Have we already modified existing face (first face gets done
314  // as modification; all following ones as polyAddFace)
315  bool modifiedFace = false;
316 
317  // Example face:
318  // +--+
319  // / |
320  // / |
321  // + +
322  // \ |
323  // \ |
324  // +--+
325  //
326  // Needs to get split into:
327  // - three 'side' faces a,b,c
328  // - one middle face d
329  // +--+
330  // /|\A|
331  // / | \|
332  // + C|D +
333  // \ | /|
334  // \|/B|
335  // +--+
336 
337 
338  // Storage for new face
339  DynamicList<label> newFace(extendedFace.size());
340 
341  label fp = startFp;
342 
343  forAll(extendedFace, i)
344  {
345  label pointI = extendedFace[fp];
346 
347  newFace.append(pointI);
348 
349  if
350  (
351  newFace.size() > 2
352  && (
353  pointI >= mesh_.nPoints()
354  || pointToPos.found(pointI)
355  )
356  )
357  {
358  // Enough vertices to create a face from.
359  face tmpFace;
360  tmpFace.transfer(newFace);
361 
362  // Add face tmpFace
363  addFace(faceI, tmpFace, modifiedFace, meshMod);
364 
365  // Starting point is also the starting point for the new face
366  newFace.append(extendedFace[startFp]);
367  newFace.append(extendedFace[fp]);
368  }
369 
370  fp = (fp+1) % extendedFace.size();
371  }
372 
373  // Check final face.
374  if (newFace.size() > 2)
375  {
376  // Enough vertices to create a face from.
377  face tmpFace;
378  tmpFace.transfer(newFace);
379 
380  // Add face tmpFace
381  addFace(faceI, tmpFace, modifiedFace, meshMod);
382  }
383 
384  // Split something
385  return true;
386  }
387 }
388 
389 
390 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
391 
392 // Construct from components
393 Foam::boundaryCutter::boundaryCutter(const polyMesh& mesh)
394 :
395  mesh_(mesh),
396  edgeAddedPoints_(),
397  faceAddedPoint_()
398 {}
399 
400 
401 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
402 
404 {}
405 
406 
407 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
408 
410 (
411  const Map<point>& pointToPos,
412  const Map<List<point> >& edgeToCuts,
413  const Map<labelPair>& faceToSplit,
414  const Map<point>& faceToFeaturePoint,
415  polyTopoChange& meshMod
416 )
417 {
418  // Clear and size maps here since mesh size will change.
419  edgeAddedPoints_.clear();
420 
421  faceAddedPoint_.clear();
422  faceAddedPoint_.resize(faceToFeaturePoint.size());
423 
424 
425  //
426  // Points that just need to be moved
427  // Note: could just as well be handled outside of setRefinement.
428  //
429 
430  forAllConstIter(Map<point>, pointToPos, iter)
431  {
432  meshMod.setAction
433  (
435  (
436  iter.key(), // point
437  iter(), // position
438  false, // no zone
439  -1, // zone for point
440  true // supports a cell
441  )
442  );
443  }
444 
445 
446  //
447  // Add new points along cut edges.
448  //
449 
450  // Map from edge label to sorted list of points
451  Map<labelList> edgeToAddedPoints(edgeToCuts.size());
452 
453  forAllConstIter(Map<List<point> >, edgeToCuts, iter)
454  {
455  label edgeI = iter.key();
456 
457  const edge& e = mesh_.edges()[edgeI];
458 
459  // Sorted (from start to end) list of cuts on edge
460  const List<point>& cuts = iter();
461 
462  forAll(cuts, cutI)
463  {
464  // point on feature to move to
465  const point& featurePoint = cuts[cutI];
466 
467  label addedPointI =
468  meshMod.setAction
469  (
471  (
472  featurePoint, // point
473  e.start(), // master point
474  -1, // zone for point
475  true // supports a cell
476  )
477  );
478 
479  Map<labelList>::iterator fnd = edgeToAddedPoints.find(edgeI);
480 
481  if (fnd != edgeToAddedPoints.end())
482  {
483  labelList& addedPoints = fnd();
484 
485  label sz = addedPoints.size();
486  addedPoints.setSize(sz+1);
487  addedPoints[sz] = addedPointI;
488  }
489  else
490  {
491  edgeToAddedPoints.insert(edgeI, labelList(1, addedPointI));
492  }
493 
494  if (debug)
495  {
496  Pout<< "Added point " << addedPointI << " for edge " << edgeI
497  << " with cuts:" << edgeToAddedPoints[edgeI] << endl;
498  }
499  }
500  }
501 
502 
503  //
504  // Introduce feature points.
505  //
506 
507  forAllConstIter(Map<point>, faceToFeaturePoint, iter)
508  {
509  label faceI = iter.key();
510 
511  const face& f = mesh_.faces()[faceI];
512 
513  if (faceToSplit.found(faceI))
514  {
515  FatalErrorIn("boundaryCutter::setRefinement")
516  << "Face " << faceI << " vertices " << f
517  << " is both marked for face-centre decomposition and"
518  << " diagonal splitting."
519  << abort(FatalError);
520  }
521 
522  if (mesh_.isInternalFace(faceI))
523  {
524  FatalErrorIn("boundaryCutter::setRefinement")
525  << "Face " << faceI << " vertices " << f
526  << " is not an external face. Cannot split it"
527  << abort(FatalError);
528  }
529 
530  label addedPointI =
531  meshMod.setAction
532  (
534  (
535  iter(), // point
536  f[0], // master point
537  -1, // zone for point
538  true // supports a cell
539  )
540  );
541  faceAddedPoint_.insert(faceI, addedPointI);
542 
543  if (debug)
544  {
545  Pout<< "Added point " << addedPointI << " for feature point "
546  << iter() << " on face " << faceI << " with centre "
547  << mesh_.faceCentres()[faceI] << endl;
548  }
549  }
550 
551 
552  //
553  // Split or retriangulate faces
554  //
555 
556 
557  // Maintain whether face has been updated (for -split edges
558  // -new owner/neighbour)
559  boolList faceUptodate(mesh_.nFaces(), false);
560 
561 
562  // Triangulate faces containing feature points
563  forAllConstIter(Map<label>, faceAddedPoint_, iter)
564  {
565  label faceI = iter.key();
566 
567  // Get face with new points on cut edges.
568  face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
569 
570  label addedPointI = iter();
571 
572  // Information about old face
573  label patchID, zoneID, zoneFlip;
574  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
575  label own = mesh_.faceOwner()[faceI];
576  label masterPoint = mesh_.faces()[faceI][0];
577 
578  // Triangulate face around mid point
579 
580  face tri(3);
581 
582  forAll(newFace, fp)
583  {
584  label nextV = newFace.nextLabel(fp);
585 
586  tri[0] = newFace[fp];
587  tri[1] = nextV;
588  tri[2] = addedPointI;
589 
590  if (fp == 0)
591  {
592  // Modify the existing face.
593  meshMod.setAction
594  (
596  (
597  tri, // face
598  faceI,
599  own, // owner
600  -1, // neighbour
601  false, // flux flip
602  patchID, // patch for face
603  false, // remove from zone
604  zoneID, // zone for face
605  zoneFlip // face zone flip
606  )
607  );
608  }
609  else
610  {
611  // Add additional faces
612  meshMod.setAction
613  (
615  (
616  tri, // face
617  own, // owner
618  -1, // neighbour
619  masterPoint, // master point
620  -1, // master edge
621  -1, // master face for addition
622  false, // flux flip
623  patchID, // patch for face
624  zoneID, // zone for face
625  zoneFlip // face zone flip
626  )
627  );
628  }
629  }
630 
631  faceUptodate[faceI] = true;
632  }
633 
634 
635  // Diagonally split faces
636  forAllConstIter(Map<labelPair>, faceToSplit, iter)
637  {
638  label faceI = iter.key();
639 
640  const face& f = mesh_.faces()[faceI];
641 
642  if (faceAddedPoint_.found(faceI))
643  {
644  FatalErrorIn("boundaryCutter::setRefinement")
645  << "Face " << faceI << " vertices " << f
646  << " is both marked for face-centre decomposition and"
647  << " diagonal splitting."
648  << abort(FatalError);
649  }
650 
651 
652  // Get face with new points on cut edges.
653  face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
654 
655  // Information about old face
656  label patchID, zoneID, zoneFlip;
657  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
658  label own = mesh_.faceOwner()[faceI];
659  label masterPoint = mesh_.faces()[faceI][0];
660 
661  // Split face from one side of diagonal to other.
662  const labelPair& diag = iter();
663 
664  label fp0 = findIndex(newFace, f[diag[0]]);
665  label fp1 = findIndex(newFace, f[diag[1]]);
666 
667  if (fp0 == -1 || fp1 == -1 || fp0 == fp1)
668  {
669  FatalErrorIn("boundaryCutter::setRefinement")
670  << "Problem : Face " << faceI << " vertices " << f
671  << " newFace:" << newFace << " diagonal:" << f[diag[0]]
672  << ' ' << f[diag[1]]
673  << abort(FatalError);
674  }
675 
676  // Replace existing face by newFace from fp0 to fp1 and add new one
677  // from fp1 to fp0.
678 
679  DynamicList<label> newVerts(newFace.size());
680 
681  // Get vertices from fp0 to (and including) fp1
682  label fp = fp0;
683 
684  do
685  {
686  newVerts.append(newFace[fp]);
687 
688  fp = (fp == newFace.size()-1 ? 0 : fp+1);
689 
690  } while (fp != fp1);
691 
692  newVerts.append(newFace[fp1]);
693 
694 
695  // Modify the existing face.
696  meshMod.setAction
697  (
699  (
700  face(newVerts.shrink()), // face
701  faceI,
702  own, // owner
703  -1, // neighbour
704  false, // flux flip
705  patchID, // patch for face
706  false, // remove from zone
707  zoneID, // zone for face
708  zoneFlip // face zone flip
709  )
710  );
711 
712 
713  newVerts.clear();
714 
715  // Get vertices from fp1 to (and including) fp0
716 
717  do
718  {
719  newVerts.append(newFace[fp]);
720 
721  fp = (fp == newFace.size()-1 ? 0 : fp+1);
722 
723  } while (fp != fp0);
724 
725  newVerts.append(newFace[fp0]);
726 
727  // Add additional face
728  meshMod.setAction
729  (
731  (
732  face(newVerts.shrink()), // face
733  own, // owner
734  -1, // neighbour
735  masterPoint, // master point
736  -1, // master edge
737  -1, // master face for addition
738  false, // flux flip
739  patchID, // patch for face
740  zoneID, // zone for face
741  zoneFlip // face zone flip
742  )
743  );
744 
745  faceUptodate[faceI] = true;
746  }
747 
748 
749  // Split external faces without feature point but using cut edges.
750  // Does right handed walk but not really.
751  forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
752  {
753  label edgeI = iter.key();
754 
755  const labelList& eFaces = mesh_.edgeFaces()[edgeI];
756 
757  forAll(eFaces, i)
758  {
759  label faceI = eFaces[i];
760 
761  if (!faceUptodate[faceI] && !mesh_.isInternalFace(faceI))
762  {
763  // Is external face so split
764  if (splitFace(faceI, pointToPos, edgeToAddedPoints, meshMod))
765  {
766  // Successfull split
767  faceUptodate[faceI] = true;
768  }
769  }
770  }
771  }
772 
773 
774  // Add cut edges (but don't split) any other faces using any cut edge.
775  // These can be external faces where splitFace hasn't cut them or
776  // internal faces.
777  forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
778  {
779  label edgeI = iter.key();
780 
781  const labelList& eFaces = mesh_.edgeFaces()[edgeI];
782 
783  forAll(eFaces, i)
784  {
785  label faceI = eFaces[i];
786 
787  if (!faceUptodate[faceI])
788  {
789  // Renumber face to include split edges.
790  face newFace(addEdgeCutsToFace(faceI, edgeToAddedPoints));
791 
792  label own = mesh_.faceOwner()[faceI];
793 
794  label nei = -1;
795 
796  if (mesh_.isInternalFace(faceI))
797  {
798  nei = mesh_.faceNeighbour()[faceI];
799  }
800 
801  label patchID, zoneID, zoneFlip;
802  getFaceInfo(faceI, patchID, zoneID, zoneFlip);
803 
804  meshMod.setAction
805  (
807  (
808  newFace, // modified face
809  faceI, // label of face being modified
810  own, // owner
811  nei, // neighbour
812  false, // face flip
813  patchID, // patch for face
814  false, // remove from zone
815  zoneID, // zone for face
816  zoneFlip // face flip in zone
817  )
818  );
819 
820  faceUptodate[faceI] = true;
821  }
822  }
823  }
824 
825  // Convert edge to points storage from edge labels (not preserved)
826  // to point labels
827  edgeAddedPoints_.resize(edgeToCuts.size());
828 
829  forAllConstIter(Map<labelList>, edgeToAddedPoints, iter)
830  {
831  edgeAddedPoints_.insert(mesh_.edges()[iter.key()], iter());
832  }
833 }
834 
835 
837 {
838  // Update stored labels for mesh change.
839 
840  //
841  // Do faceToAddedPoint
842  //
843 
844  {
845  // Create copy since we're deleting entries.
846  Map<label> newAddedPoints(faceAddedPoint_.size());
847 
848  forAllConstIter(Map<label>, faceAddedPoint_, iter)
849  {
850  label oldFaceI = iter.key();
851 
852  label newFaceI = morphMap.reverseFaceMap()[oldFaceI];
853 
854  label oldPointI = iter();
855 
856  label newPointI = morphMap.reversePointMap()[oldPointI];
857 
858  if (newFaceI >= 0 && newPointI >= 0)
859  {
860  newAddedPoints.insert(newFaceI, newPointI);
861  }
862  }
863 
864  // Copy
865  faceAddedPoint_.transfer(newAddedPoints);
866  }
867 
868 
869  //
870  // Do edgeToAddedPoints
871  //
872 
873 
874  {
875  // Create copy since we're deleting entries
877  newEdgeAddedPoints(edgeAddedPoints_.size());
878 
879  for
880  (
881  HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
882  edgeAddedPoints_.begin();
883  iter != edgeAddedPoints_.end();
884  ++iter
885  )
886  {
887  const edge& e = iter.key();
888 
889  label newStart = morphMap.reversePointMap()[e.start()];
890 
891  label newEnd = morphMap.reversePointMap()[e.end()];
892 
893  if (newStart >= 0 && newEnd >= 0)
894  {
895  const labelList& addedPoints = iter();
896 
897  labelList newAddedPoints(addedPoints.size());
898  label newI = 0;
899 
900  forAll(addedPoints, i)
901  {
902  label newAddedPointI =
903  morphMap.reversePointMap()[addedPoints[i]];
904 
905  if (newAddedPointI >= 0)
906  {
907  newAddedPoints[newI++] = newAddedPointI;
908  }
909  }
910  if (newI > 0)
911  {
912  newAddedPoints.setSize(newI);
913 
914  edge newE = edge(newStart, newEnd);
915 
916  newEdgeAddedPoints.insert(newE, newAddedPoints);
917  }
918  }
919  }
920 
921  // Copy
922  edgeAddedPoints_.transfer(newEdgeAddedPoints);
923  }
924 }
925 
926 
927 // ************************ vim: set sw=4 sts=4 et: ************************ //