FreeFOAM The Cross-Platform CFD Toolkit
meshCutAndRemove.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 "meshCutAndRemove.H"
27 #include <OpenFOAM/polyMesh.H>
34 #include <dynamicMesh/cellCuts.H>
35 #include <OpenFOAM/mapPolyMesh.H>
36 #include <meshTools/meshTools.H>
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
41 
42 // * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * //
43 
44 // Returns -1 or index in elems1 of first shared element.
45 Foam::label Foam::meshCutAndRemove::firstCommon
46 (
47  const labelList& elems1,
48  const labelList& elems2
49 )
50 {
51  forAll(elems1, elemI)
52  {
53  label index1 = findIndex(elems2, elems1[elemI]);
54 
55  if (index1 != -1)
56  {
57  return index1;
58  }
59  }
60  return -1;
61 }
62 
63 
64 // Check if twoCuts at two consecutive position in cuts.
65 bool Foam::meshCutAndRemove::isIn
66 (
67  const edge& twoCuts,
68  const labelList& cuts
69 )
70 {
71  label index = findIndex(cuts, twoCuts[0]);
72 
73  if (index == -1)
74  {
75  return false;
76  }
77 
78  return
79  (
80  cuts[cuts.fcIndex(index)] == twoCuts[1]
81  || cuts[cuts.rcIndex(index)] == twoCuts[1]
82  );
83 }
84 
85 
86 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
87 
88 // Returns the cell in cellLabels that is cut. Or -1.
89 Foam::label Foam::meshCutAndRemove::findCutCell
90 (
91  const cellCuts& cuts,
92  const labelList& cellLabels
93 ) const
94 {
95  forAll(cellLabels, labelI)
96  {
97  label cellI = cellLabels[labelI];
98 
99  if (cuts.cellLoops()[cellI].size())
100  {
101  return cellI;
102  }
103  }
104  return -1;
105 }
106 
107 
108 //- Returns first pointI in pointLabels that uses an internal
109 // face. Used to find point to inflate cell/face from (has to be
110 // connected to internal face). Returns -1 (so inflate from nothing) if
111 // none found.
112 Foam::label Foam::meshCutAndRemove::findInternalFacePoint
113 (
114  const labelList& pointLabels
115 ) const
116 {
117  forAll(pointLabels, labelI)
118  {
119  label pointI = pointLabels[labelI];
120 
121  const labelList& pFaces = mesh().pointFaces()[pointI];
122 
123  forAll(pFaces, pFaceI)
124  {
125  label faceI = pFaces[pFaceI];
126 
127  if (mesh().isInternalFace(faceI))
128  {
129  return pointI;
130  }
131  }
132  }
133 
134  if (pointLabels.empty())
135  {
136  FatalErrorIn("meshCutAndRemove::findInternalFacePoint(const labelList&)")
137  << "Empty pointLabels" << abort(FatalError);
138  }
139 
140  return -1;
141 }
142 
143 
144 // Find point on face that is part of original mesh and that is point connected
145 // to the patch
146 Foam::label Foam::meshCutAndRemove::findPatchFacePoint
147 (
148  const face& f,
149  const label exposedPatchI
150 ) const
151 {
152  const labelListList& pointFaces = mesh().pointFaces();
153  const polyBoundaryMesh& patches = mesh().boundaryMesh();
154 
155  forAll(f, fp)
156  {
157  label pointI = f[fp];
158 
159  if (pointI < mesh().nPoints())
160  {
161  const labelList& pFaces = pointFaces[pointI];
162 
163  forAll(pFaces, i)
164  {
165  if (patches.whichPatch(pFaces[i]) == exposedPatchI)
166  {
167  return pointI;
168  }
169  }
170  }
171  }
172  return -1;
173 }
174 
175 
176 // Get new owner and neighbour of face. Checks anchor points to see if
177 // cells have been removed.
178 void Foam::meshCutAndRemove::faceCells
179 (
180  const cellCuts& cuts,
181  const label exposedPatchI,
182  const label faceI,
183  label& own,
184  label& nei,
185  label& patchID
186 ) const
187 {
188  const labelListList& anchorPts = cuts.cellAnchorPoints();
189  const labelListList& cellLoops = cuts.cellLoops();
190 
191  const face& f = mesh().faces()[faceI];
192 
193  own = mesh().faceOwner()[faceI];
194 
195  if (cellLoops[own].size() && firstCommon(f, anchorPts[own]) == -1)
196  {
197  // owner has been split and this is the removed part.
198  own = -1;
199  }
200 
201  nei = -1;
202 
203  if (mesh().isInternalFace(faceI))
204  {
205  nei = mesh().faceNeighbour()[faceI];
206 
207  if (cellLoops[nei].size() && firstCommon(f, anchorPts[nei]) == -1)
208  {
209  nei = -1;
210  }
211  }
212 
213  patchID = mesh().boundaryMesh().whichPatch(faceI);
214 
215  if (patchID == -1 && (own == -1 || nei == -1))
216  {
217  // Face was internal but becomes external
218  patchID = exposedPatchI;
219  }
220 }
221 
222 
223 void Foam::meshCutAndRemove::getZoneInfo
224 (
225  const label faceI,
226  label& zoneID,
227  bool& zoneFlip
228 ) const
229 {
230  zoneID = mesh().faceZones().whichZone(faceI);
231 
232  zoneFlip = false;
233 
234  if (zoneID >= 0)
235  {
236  const faceZone& fZone = mesh().faceZones()[zoneID];
237 
238  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
239  }
240 }
241 
242 
243 // Adds a face from point.
244 void Foam::meshCutAndRemove::addFace
245 (
246  polyTopoChange& meshMod,
247  const label faceI,
248  const label masterPointI,
249  const face& newFace,
250  const label own,
251  const label nei,
252  const label patchID
253 )
254 {
255  label zoneID;
256  bool zoneFlip;
257 
258  getZoneInfo(faceI, zoneID, zoneFlip);
259 
260  if ((nei == -1) || (own != -1 && own < nei))
261  {
262  // Ordering ok.
263  if (debug & 2)
264  {
265  Pout<< "Adding face " << newFace
266  << " with new owner:" << own
267  << " with new neighbour:" << nei
268  << " patchID:" << patchID
269  << " anchor:" << masterPointI
270  << " zoneID:" << zoneID
271  << " zoneFlip:" << zoneFlip
272  << endl;
273  }
274 
275  meshMod.setAction
276  (
277  polyAddFace
278  (
279  newFace, // face
280  own, // owner
281  nei, // neighbour
282  masterPointI, // master point
283  -1, // master edge
284  -1, // master face for addition
285  false, // flux flip
286  patchID, // patch for face
287  zoneID, // zone for face
288  zoneFlip // face zone flip
289  )
290  );
291  }
292  else
293  {
294  // Reverse owner/neighbour
295  if (debug & 2)
296  {
297  Pout<< "Adding (reversed) face " << newFace.reverseFace()
298  << " with new owner:" << nei
299  << " with new neighbour:" << own
300  << " patchID:" << patchID
301  << " anchor:" << masterPointI
302  << " zoneID:" << zoneID
303  << " zoneFlip:" << zoneFlip
304  << endl;
305  }
306 
307  meshMod.setAction
308  (
309  polyAddFace
310  (
311  newFace.reverseFace(), // face
312  nei, // owner
313  own, // neighbour
314  masterPointI, // master point
315  -1, // master edge
316  -1, // master face for addition
317  false, // flux flip
318  patchID, // patch for face
319  zoneID, // zone for face
320  zoneFlip // face zone flip
321  )
322  );
323  }
324 }
325 
326 
327 // Modifies existing faceI for either new owner/neighbour or new face points.
328 void Foam::meshCutAndRemove::modFace
329 (
330  polyTopoChange& meshMod,
331  const label faceI,
332  const face& newFace,
333  const label own,
334  const label nei,
335  const label patchID
336 )
337 {
338  label zoneID;
339  bool zoneFlip;
340 
341  getZoneInfo(faceI, zoneID, zoneFlip);
342 
343  if
344  (
345  (own != mesh().faceOwner()[faceI])
346  || (
347  mesh().isInternalFace(faceI)
348  && (nei != mesh().faceNeighbour()[faceI])
349  )
350  || (newFace != mesh().faces()[faceI])
351  )
352  {
353  if (debug & 2)
354  {
355  Pout<< "Modifying face " << faceI
356  << " old vertices:" << mesh().faces()[faceI]
357  << " new vertices:" << newFace
358  << " new owner:" << own
359  << " new neighbour:" << nei
360  << " new patch:" << patchID
361  << " new zoneID:" << zoneID
362  << " new zoneFlip:" << zoneFlip
363  << endl;
364  }
365 
366  if ((nei == -1) || (own != -1 && own < nei))
367  {
368  meshMod.setAction
369  (
370  polyModifyFace
371  (
372  newFace, // modified face
373  faceI, // label of face being modified
374  own, // owner
375  nei, // neighbour
376  false, // face flip
377  patchID, // patch for face
378  false, // remove from zone
379  zoneID, // zone for face
380  zoneFlip // face flip in zone
381  )
382  );
383  }
384  else
385  {
386  meshMod.setAction
387  (
388  polyModifyFace
389  (
390  newFace.reverseFace(), // modified face
391  faceI, // label of face being modified
392  nei, // owner
393  own, // neighbour
394  false, // face flip
395  patchID, // patch for face
396  false, // remove from zone
397  zoneID, // zone for face
398  zoneFlip // face flip in zone
399  )
400  );
401  }
402  }
403 }
404 
405 
406 // Copies face starting from startFp up to and including endFp.
407 void Foam::meshCutAndRemove::copyFace
408 (
409  const face& f,
410  const label startFp,
411  const label endFp,
412  face& newFace
413 ) const
414 {
415  label fp = startFp;
416 
417  label newFp = 0;
418 
419  while (fp != endFp)
420  {
421  newFace[newFp++] = f[fp];
422 
423  fp = (fp + 1) % f.size();
424  }
425  newFace[newFp] = f[fp];
426 }
427 
428 
429 // Actually split face in two along splitEdge v0, v1 (the two vertices in new
430 // vertex numbering). Generates faces in same ordering
431 // as original face. Replaces cutEdges by the points introduced on them
432 // (addedPoints_).
433 void Foam::meshCutAndRemove::splitFace
434 (
435  const face& f,
436  const label v0,
437  const label v1,
438 
439  face& f0,
440  face& f1
441 ) const
442 {
443  // Check if we find any new vertex which is part of the splitEdge.
444  label startFp = findIndex(f, v0);
445 
446  if (startFp == -1)
447  {
449  (
450  "meshCutAndRemove::splitFace"
451  ", const face&, const label, const label, face&, face&)"
452  ) << "Cannot find vertex (new numbering) " << v0
453  << " on face " << f
454  << abort(FatalError);
455  }
456 
457  label endFp = findIndex(f, v1);
458 
459  if (endFp == -1)
460  {
462  (
463  "meshCutAndRemove::splitFace("
464  ", const face&, const label, const label, face&, face&)"
465  ) << "Cannot find vertex (new numbering) " << v1
466  << " on face " << f
467  << abort(FatalError);
468  }
469 
470 
471  f0.setSize((endFp + 1 + f.size() - startFp) % f.size());
472  f1.setSize(f.size() - f0.size() + 2);
473 
474  copyFace(f, startFp, endFp, f0);
475  copyFace(f, endFp, startFp, f1);
476 }
477 
478 
479 // Adds additional vertices (from edge cutting) to face. Used for faces which
480 // are not split but still might use edge that has been cut.
481 Foam::face Foam::meshCutAndRemove::addEdgeCutsToFace(const label faceI) const
482 {
483  const face& f = mesh().faces()[faceI];
484 
485  face newFace(2 * f.size());
486 
487  label newFp = 0;
488 
489  forAll(f, fp)
490  {
491  // Duplicate face vertex.
492  newFace[newFp++] = f[fp];
493 
494  // Check if edge has been cut.
495  label fp1 = f.fcIndex(fp);
496 
497  HashTable<label, edge, Hash<edge> >::const_iterator fnd =
498  addedPoints_.find(edge(f[fp], f[fp1]));
499 
500  if (fnd != addedPoints_.end())
501  {
502  // edge has been cut. Introduce new vertex.
503  newFace[newFp++] = fnd();
504  }
505  }
506 
507  newFace.setSize(newFp);
508 
509  return newFace;
510 }
511 
512 
513 // Walk loop (loop of cuts) across circumference of cellI. Returns face in
514 // new vertices.
515 // Note: tricky bit is that it can use existing edges which have been split.
516 Foam::face Foam::meshCutAndRemove::loopToFace
517 (
518  const label cellI,
519  const labelList& loop
520 ) const
521 {
522  face newFace(2*loop.size());
523 
524  label newFaceI = 0;
525 
526  forAll(loop, fp)
527  {
528  label cut = loop[fp];
529 
530  if (isEdge(cut))
531  {
532  label edgeI = getEdge(cut);
533 
534  const edge& e = mesh().edges()[edgeI];
535 
536  label vertI = addedPoints_[e];
537 
538  newFace[newFaceI++] = vertI;
539  }
540  else
541  {
542  // cut is vertex.
543  label vertI = getVertex(cut);
544 
545  newFace[newFaceI++] = vertI;
546 
547  label nextCut = loop[loop.fcIndex(fp)];
548 
549  if (!isEdge(nextCut))
550  {
551  // From vertex to vertex -> cross cut only if no existing edge.
552  label nextVertI = getVertex(nextCut);
553 
554  label edgeI = meshTools::findEdge(mesh(), vertI, nextVertI);
555 
556  if (edgeI != -1)
557  {
558  // Existing edge. Insert split-edge point if any.
559  HashTable<label, edge, Hash<edge> >::const_iterator fnd =
560  addedPoints_.find(mesh().edges()[edgeI]);
561 
562  if (fnd != addedPoints_.end())
563  {
564  newFace[newFaceI++] = fnd();
565  }
566  }
567  }
568  }
569  }
570  newFace.setSize(newFaceI);
571 
572  return newFace;
573 }
574 
575 
576 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
577 
578 // Construct from components
579 Foam::meshCutAndRemove::meshCutAndRemove(const polyMesh& mesh)
580 :
581  edgeVertex(mesh),
582  addedFaces_(),
583  addedPoints_()
584 {}
585 
586 
587 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
588 
590 (
591  const label exposedPatchI,
592  const cellCuts& cuts,
593  const labelList& cutPatch,
594  polyTopoChange& meshMod
595 )
596 {
597  // Clear and size maps here since mesh size will change.
598  addedFaces_.clear();
599  addedFaces_.resize(cuts.nLoops());
600 
601  addedPoints_.clear();
602  addedPoints_.resize(cuts.nLoops());
603 
604  if (cuts.nLoops() == 0)
605  {
606  return;
607  }
608 
609  const labelListList& anchorPts = cuts.cellAnchorPoints();
610  const labelListList& cellLoops = cuts.cellLoops();
611  const polyBoundaryMesh& patches = mesh().boundaryMesh();
612 
613  if (exposedPatchI < 0 || exposedPatchI >= patches.size())
614  {
616  (
617  "meshCutAndRemove::setRefinement("
618  ", const label, const cellCuts&, const labelList&"
619  ", polyTopoChange&)"
620  ) << "Illegal exposed patch " << exposedPatchI
621  << abort(FatalError);
622  }
623 
624 
625  //
626  // Add new points along cut edges.
627  //
628 
629  forAll(cuts.edgeIsCut(), edgeI)
630  {
631  if (cuts.edgeIsCut()[edgeI])
632  {
633  const edge& e = mesh().edges()[edgeI];
634 
635  // Check if there is any cell using this edge.
636  if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
637  {
639  (
640  "meshCutAndRemove::setRefinement("
641  ", const label, const cellCuts&, const labelList&"
642  ", polyTopoChange&)"
643  ) << "Problem: cut edge but none of the cells using it is\n"
644  << "edge:" << edgeI << " verts:" << e
645  << abort(FatalError);
646  }
647 
648  // One of the edge end points should be master point of nbCellI.
649  label masterPointI = e.start();
650 
651  const point& v0 = mesh().points()[e.start()];
652  const point& v1 = mesh().points()[e.end()];
653 
654  scalar weight = cuts.edgeWeight()[edgeI];
655 
656  point newPt = weight*v1 + (1.0-weight)*v0;
657 
658  label addedPointI =
659  meshMod.setAction
660  (
662  (
663  newPt, // point
664  masterPointI, // master point
665  -1, // zone for point
666  true // supports a cell
667  )
668  );
669 
670  // Store on (hash of) edge.
671  addedPoints_.insert(e, addedPointI);
672 
673  if (debug & 2)
674  {
675  Pout<< "Added point " << addedPointI
676  << " to vertex "
677  << masterPointI << " of edge " << edgeI
678  << " vertices " << e << endl;
679  }
680  }
681  }
682 
683 
684  //
685  // Remove all points that will not be used anymore
686  //
687  {
688  boolList usedPoint(mesh().nPoints(), false);
689 
690  forAll(cellLoops, cellI)
691  {
692  const labelList& loop = cellLoops[cellI];
693 
694  if (loop.size())
695  {
696  // Cell is cut. Uses only anchor points and loop itself.
697  forAll(loop, fp)
698  {
699  label cut = loop[fp];
700 
701  if (!isEdge(cut))
702  {
703  usedPoint[getVertex(cut)] = true;
704  }
705  }
706 
707  const labelList& anchors = anchorPts[cellI];
708 
709  forAll(anchors, i)
710  {
711  usedPoint[anchors[i]] = true;
712  }
713  }
714  else
715  {
716  // Cell is not cut so use all its points
717  const labelList& cPoints = mesh().cellPoints()[cellI];
718 
719  forAll(cPoints, i)
720  {
721  usedPoint[cPoints[i]] = true;
722  }
723  }
724  }
725 
726 
727  // Check
728  const Map<edge>& faceSplitCut = cuts.faceSplitCut();
729 
730  forAllConstIter(Map<edge>, faceSplitCut, iter)
731  {
732  const edge& fCut = iter();
733 
734  forAll(fCut, i)
735  {
736  label cut = fCut[i];
737 
738  if (!isEdge(cut))
739  {
740  label pointI = getVertex(cut);
741 
742  if (!usedPoint[pointI])
743  {
745  (
746  "meshCutAndRemove::setRefinement("
747  ", const label, const cellCuts&, const labelList&"
748  ", polyTopoChange&)"
749  ) << "Problem: faceSplitCut not used by any loop"
750  << " or cell anchor point"
751  << "face:" << iter.key() << " point:" << pointI
752  << " coord:" << mesh().points()[pointI]
753  << abort(FatalError);
754  }
755  }
756  }
757  }
758 
759  forAll(cuts.pointIsCut(), pointI)
760  {
761  if (cuts.pointIsCut()[pointI])
762  {
763  if (!usedPoint[pointI])
764  {
766  (
767  "meshCutAndRemove::setRefinement("
768  ", const label, const cellCuts&, const labelList&"
769  ", polyTopoChange&)"
770  ) << "Problem: point is marked as cut but"
771  << " not used by any loop"
772  << " or cell anchor point"
773  << "point:" << pointI
774  << " coord:" << mesh().points()[pointI]
775  << abort(FatalError);
776  }
777  }
778  }
779 
780 
781  // Remove unused points.
782  forAll(usedPoint, pointI)
783  {
784  if (!usedPoint[pointI])
785  {
786  meshMod.setAction(polyRemovePoint(pointI));
787 
788  if (debug & 2)
789  {
790  Pout<< "Removing unused point " << pointI << endl;
791  }
792  }
793  }
794  }
795 
796 
797  //
798  // For all cut cells add an internal or external face
799  //
800 
801  forAll(cellLoops, cellI)
802  {
803  const labelList& loop = cellLoops[cellI];
804 
805  if (loop.size())
806  {
807  if (cutPatch[cellI] < 0 || cutPatch[cellI] >= patches.size())
808  {
810  (
811  "meshCutAndRemove::setRefinement("
812  ", const label, const cellCuts&, const labelList&"
813  ", polyTopoChange&)"
814  ) << "Illegal patch " << cutPatch[cellI]
815  << " provided for cut cell " << cellI
816  << abort(FatalError);
817  }
818 
819  //
820  // Convert loop (=list of cuts) into proper face.
821  // cellCuts sets orientation is towards anchor side so reverse.
822  //
823  face newFace(loopToFace(cellI, loop));
824 
825  reverse(newFace);
826 
827  // Pick any anchor point on cell
828  label masterPointI = findPatchFacePoint(newFace, exposedPatchI);
829 
830  label addedFaceI =
831  meshMod.setAction
832  (
834  (
835  newFace, // face
836  cellI, // owner
837  -1, // neighbour
838  masterPointI, // master point
839  -1, // master edge
840  -1, // master face for addition
841  false, // flux flip
842  cutPatch[cellI], // patch for face
843  -1, // zone for face
844  false // face zone flip
845  )
846  );
847 
848  addedFaces_.insert(cellI, addedFaceI);
849 
850  if (debug & 2)
851  {
852  Pout<< "Added splitting face " << newFace << " index:"
853  << addedFaceI << " from masterPoint:" << masterPointI
854  << " to owner " << cellI << " with anchors:"
855  << anchorPts[cellI]
856  << " from Loop:";
857 
858  // Gets edgeweights of loop
859  scalarField weights(loop.size());
860  forAll(loop, i)
861  {
862  label cut = loop[i];
863 
864  weights[i] =
865  (
866  isEdge(cut)
867  ? cuts.edgeWeight()[getEdge(cut)]
868  : -GREAT
869  );
870  }
871  writeCuts(Pout, loop, weights);
872  Pout<< endl;
873  }
874  }
875  }
876 
877 
878  //
879  // Modify faces to use only anchorpoints and loop points
880  // (so throw away part without anchorpoints)
881  //
882 
883 
884  // Maintain whether face has been updated (for -split edges
885  // -new owner/neighbour)
886  boolList faceUptodate(mesh().nFaces(), false);
887 
888  const Map<edge>& faceSplitCuts = cuts.faceSplitCut();
889 
890  for
891  (
892  Map<edge>::const_iterator iter = faceSplitCuts.begin();
893  iter != faceSplitCuts.end();
894  ++iter
895  )
896  {
897  label faceI = iter.key();
898 
899  // Renumber face to include split edges.
900  face newFace(addEdgeCutsToFace(faceI));
901 
902  // Edge splitting the face. Convert edge to new vertex numbering.
903  const edge& splitEdge = iter();
904 
905  label cut0 = splitEdge[0];
906 
907  label v0;
908  if (isEdge(cut0))
909  {
910  label edgeI = getEdge(cut0);
911  v0 = addedPoints_[mesh().edges()[edgeI]];
912  }
913  else
914  {
915  v0 = getVertex(cut0);
916  }
917 
918  label cut1 = splitEdge[1];
919  label v1;
920  if (isEdge(cut1))
921  {
922  label edgeI = getEdge(cut1);
923  v1 = addedPoints_[mesh().edges()[edgeI]];
924  }
925  else
926  {
927  v1 = getVertex(cut1);
928  }
929 
930  // Split face along the elements of the splitEdge.
931  face f0, f1;
932  splitFace(newFace, v0, v1, f0, f1);
933 
934  label own = mesh().faceOwner()[faceI];
935 
936  label nei = -1;
937 
938  if (mesh().isInternalFace(faceI))
939  {
940  nei = mesh().faceNeighbour()[faceI];
941  }
942 
943  if (debug & 2)
944  {
945  Pout<< "Split face " << mesh().faces()[faceI]
946  << " own:" << own << " nei:" << nei
947  << " into f0:" << f0
948  << " and f1:" << f1 << endl;
949  }
950 
951 
952  // Check which cell using face uses anchorPoints (so is kept)
953  // and which one doesn't (gets removed)
954 
955  // Bit tricky. We have to know whether this faceSplit splits owner/
956  // neighbour or both. Even if cell is cut we have to make sure this is
957  // the one that cuts it (this face cut might not be the one splitting
958  // the cell)
959  // The face f gets split into two parts, f0 and f1.
960  // Each of these can have a different owner and or neighbour.
961 
962  const face& f = mesh().faces()[faceI];
963 
964  label f0Own = -1;
965  label f1Own = -1;
966 
967  if (cellLoops[own].empty())
968  {
969  // Owner side is not split so keep both halves.
970  f0Own = own;
971  f1Own = own;
972  }
973  else if (isIn(splitEdge, cellLoops[own]))
974  {
975  // Owner is cut by this splitCut. See which of f0, f1 gets
976  // preserved and becomes owner, and which gets removed.
977  if (firstCommon(f0, anchorPts[own]) != -1)
978  {
979  // f0 preserved so f1 gets deleted
980  f0Own = own;
981  f1Own = -1;
982  }
983  else
984  {
985  f0Own = -1;
986  f1Own = own;
987  }
988  }
989  else
990  {
991  // Owner not cut by this splitCut but by another.
992  // Check on original face whether
993  // use anchorPts.
994  if (firstCommon(f, anchorPts[own]) != -1)
995  {
996  // both f0 and f1 owner side preserved
997  f0Own = own;
998  f1Own = own;
999  }
1000  else
1001  {
1002  // both f0 and f1 owner side removed
1003  f0Own = -1;
1004  f1Own = -1;
1005  }
1006  }
1007 
1008 
1009  label f0Nei = -1;
1010  label f1Nei = -1;
1011 
1012  if (nei != -1)
1013  {
1014  if (cellLoops[nei].empty())
1015  {
1016  f0Nei = nei;
1017  f1Nei = nei;
1018  }
1019  else if (isIn(splitEdge, cellLoops[nei]))
1020  {
1021  // Neighbour is cut by this splitCut. So anchor part of it
1022  // gets kept, non-anchor bit gets removed. See which of f0, f1
1023  // connects to which part.
1024 
1025  if (firstCommon(f0, anchorPts[nei]) != -1)
1026  {
1027  f0Nei = nei;
1028  f1Nei = -1;
1029  }
1030  else
1031  {
1032  f0Nei = -1;
1033  f1Nei = nei;
1034  }
1035  }
1036  else
1037  {
1038  // neighbour not cut by this splitCut. Check on original face
1039  // whether use anchorPts.
1040 
1041  if (firstCommon(f, anchorPts[nei]) != -1)
1042  {
1043  f0Nei = nei;
1044  f1Nei = nei;
1045  }
1046  else
1047  {
1048  // both f0 and f1 on neighbour side removed
1049  f0Nei = -1;
1050  f1Nei = -1;
1051  }
1052  }
1053  }
1054 
1055 
1056  if (debug & 2)
1057  {
1058  Pout<< "f0 own:" << f0Own << " nei:" << f0Nei
1059  << " f1 own:" << f1Own << " nei:" << f1Nei
1060  << endl;
1061  }
1062 
1063 
1064  // If faces were internal but now become external set a patch.
1065  // If they were external already keep the patch.
1066  label patchID = patches.whichPatch(faceI);
1067 
1068  if (patchID == -1)
1069  {
1070  patchID = exposedPatchI;
1071  }
1072 
1073 
1074  // Do as much as possible by modifying faceI. Delay any remove
1075  // face. Keep track of whether faceI has been used.
1076 
1077  bool modifiedFaceI = false;
1078 
1079  if (f0Own == -1)
1080  {
1081  if (f0Nei != -1)
1082  {
1083  // f0 becomes external face (note:modFace will reverse face)
1084  modFace(meshMod, faceI, f0, f0Own, f0Nei, patchID);
1085  modifiedFaceI = true;
1086  }
1087  }
1088  else
1089  {
1090  if (f0Nei == -1)
1091  {
1092  // f0 becomes external face
1093  modFace(meshMod, faceI, f0, f0Own, f0Nei, patchID);
1094  modifiedFaceI = true;
1095  }
1096  else
1097  {
1098  // f0 stays internal face.
1099  modFace(meshMod, faceI, f0, f0Own, f0Nei, -1);
1100  modifiedFaceI = true;
1101  }
1102  }
1103 
1104 
1105  // f1 is added face (if at all)
1106 
1107  if (f1Own == -1)
1108  {
1109  if (f1Nei == -1)
1110  {
1111  // f1 not needed.
1112  }
1113  else
1114  {
1115  // f1 becomes external face (note:modFace will reverse face)
1116  if (!modifiedFaceI)
1117  {
1118  modFace(meshMod, faceI, f1, f1Own, f1Nei, patchID);
1119  modifiedFaceI = true;
1120  }
1121  else
1122  {
1123  label masterPointI = findPatchFacePoint(f1, patchID);
1124 
1125  addFace
1126  (
1127  meshMod,
1128  faceI, // face for zone info
1129  masterPointI, // inflation point
1130  f1, // vertices of face
1131  f1Own,
1132  f1Nei,
1133  patchID // patch for new face
1134  );
1135  }
1136  }
1137  }
1138  else
1139  {
1140  if (f1Nei == -1)
1141  {
1142  // f1 becomes external face
1143  if (!modifiedFaceI)
1144  {
1145  modFace(meshMod, faceI, f1, f1Own, f1Nei, patchID);
1146  modifiedFaceI = true;
1147  }
1148  else
1149  {
1150  label masterPointI = findPatchFacePoint(f1, patchID);
1151 
1152  addFace
1153  (
1154  meshMod,
1155  faceI,
1156  masterPointI,
1157  f1,
1158  f1Own,
1159  f1Nei,
1160  patchID
1161  );
1162  }
1163  }
1164  else
1165  {
1166  // f1 is internal face.
1167  if (!modifiedFaceI)
1168  {
1169  modFace(meshMod, faceI, f1, f1Own, f1Nei, -1);
1170  modifiedFaceI = true;
1171  }
1172  else
1173  {
1174  label masterPointI = findPatchFacePoint(f1, -1);
1175 
1176  addFace(meshMod, faceI, masterPointI, f1, f1Own, f1Nei, -1);
1177  }
1178  }
1179  }
1180 
1181  if (f0Own == -1 && f0Nei == -1 && !modifiedFaceI)
1182  {
1183  meshMod.setAction(polyRemoveFace(faceI));
1184 
1185  if (debug & 2)
1186  {
1187  Pout<< "Removed face " << faceI << endl;
1188  }
1189  }
1190 
1191  faceUptodate[faceI] = true;
1192  }
1193 
1194 
1195  //
1196  // Faces that have not been split but just appended to. Are guaranteed
1197  // to be reachable from an edgeCut.
1198  //
1199 
1200  const boolList& edgeIsCut = cuts.edgeIsCut();
1201 
1202  forAll(edgeIsCut, edgeI)
1203  {
1204  if (edgeIsCut[edgeI])
1205  {
1206  const labelList& eFaces = mesh().edgeFaces()[edgeI];
1207 
1208  forAll(eFaces, i)
1209  {
1210  label faceI = eFaces[i];
1211 
1212  if (!faceUptodate[faceI])
1213  {
1214  // So the face has not been split itself (i.e. its owner
1215  // or neighbour have not been split) so it only
1216  // borders by edge a cell which has been split.
1217 
1218  // Get (new or original) owner and neighbour of faceI
1219  label own, nei, patchID;
1220  faceCells(cuts, exposedPatchI, faceI, own, nei, patchID);
1221 
1222 
1223  if (own == -1 && nei == -1)
1224  {
1225  meshMod.setAction(polyRemoveFace(faceI));
1226 
1227  if (debug & 2)
1228  {
1229  Pout<< "Removed face " << faceI << endl;
1230  }
1231  }
1232  else
1233  {
1234  // Renumber face to include split edges.
1235  face newFace(addEdgeCutsToFace(faceI));
1236 
1237  if (debug & 2)
1238  {
1239  Pout<< "Added edge cuts to face " << faceI
1240  << " f:" << mesh().faces()[faceI]
1241  << " newFace:" << newFace << endl;
1242  }
1243 
1244  modFace
1245  (
1246  meshMod,
1247  faceI,
1248  newFace,
1249  own,
1250  nei,
1251  patchID
1252  );
1253  }
1254 
1255  faceUptodate[faceI] = true;
1256  }
1257  }
1258  }
1259  }
1260 
1261 
1262  //
1263  // Remove any faces on the non-anchor side of a split cell.
1264  // Note: could loop through all cut cells only and check their faces but
1265  // looping over all faces is cleaner and probably faster for dense
1266  // cut patterns.
1267 
1268  const faceList& faces = mesh().faces();
1269 
1270  forAll(faces, faceI)
1271  {
1272  if (!faceUptodate[faceI])
1273  {
1274  // Get (new or original) owner and neighbour of faceI
1275  label own, nei, patchID;
1276  faceCells(cuts, exposedPatchI, faceI, own, nei, patchID);
1277 
1278  if (own == -1 && nei == -1)
1279  {
1280  meshMod.setAction(polyRemoveFace(faceI));
1281 
1282  if (debug & 2)
1283  {
1284  Pout<< "Removed face " << faceI << endl;
1285  }
1286  }
1287  else
1288  {
1289  modFace(meshMod, faceI, faces[faceI], own, nei, patchID);
1290  }
1291 
1292  faceUptodate[faceI] = true;
1293  }
1294  }
1295 
1296  if (debug)
1297  {
1298  Pout<< "meshCutAndRemove:" << nl
1299  << " cells split:" << cuts.nLoops() << nl
1300  << " faces added:" << addedFaces_.size() << nl
1301  << " points added on edges:" << addedPoints_.size() << nl
1302  << endl;
1303  }
1304 }
1305 
1306 
1308 {
1309  // Update stored labels for mesh change.
1310  {
1311  Map<label> newAddedFaces(addedFaces_.size());
1312 
1313  for
1314  (
1315  Map<label>::const_iterator iter = addedFaces_.begin();
1316  iter != addedFaces_.end();
1317  ++iter
1318  )
1319  {
1320  label cellI = iter.key();
1321 
1322  label newCellI = map.reverseCellMap()[cellI];
1323 
1324  label addedFaceI = iter();
1325 
1326  label newAddedFaceI = map.reverseFaceMap()[addedFaceI];
1327 
1328  if ((newCellI >= 0) && (newAddedFaceI >= 0))
1329  {
1330  if
1331  (
1332  (debug & 2)
1333  && (newCellI != cellI || newAddedFaceI != addedFaceI)
1334  )
1335  {
1336  Pout<< "meshCutAndRemove::updateMesh :"
1337  << " updating addedFace for cell " << cellI
1338  << " from " << addedFaceI
1339  << " to " << newAddedFaceI
1340  << endl;
1341  }
1342  newAddedFaces.insert(newCellI, newAddedFaceI);
1343  }
1344  }
1345 
1346  // Copy
1347  addedFaces_.transfer(newAddedFaces);
1348  }
1349 
1350  {
1351  HashTable<label, edge, Hash<edge> > newAddedPoints(addedPoints_.size());
1352 
1353  for
1354  (
1355  HashTable<label, edge, Hash<edge> >::const_iterator iter =
1356  addedPoints_.begin();
1357  iter != addedPoints_.end();
1358  ++iter
1359  )
1360  {
1361  const edge& e = iter.key();
1362 
1363  label newStart = map.reversePointMap()[e.start()];
1364 
1365  label newEnd = map.reversePointMap()[e.end()];
1366 
1367  label addedPointI = iter();
1368 
1369  label newAddedPointI = map.reversePointMap()[addedPointI];
1370 
1371  if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointI >= 0))
1372  {
1373  edge newE = edge(newStart, newEnd);
1374 
1375  if
1376  (
1377  (debug & 2)
1378  && (e != newE || newAddedPointI != addedPointI)
1379  )
1380  {
1381  Pout<< "meshCutAndRemove::updateMesh :"
1382  << " updating addedPoints for edge " << e
1383  << " from " << addedPointI
1384  << " to " << newAddedPointI
1385  << endl;
1386  }
1387 
1388  newAddedPoints.insert(newE, newAddedPointI);
1389  }
1390  }
1391 
1392  // Copy
1393  addedPoints_.transfer(newAddedPoints);
1394  }
1395 }
1396 
1397 
1398 // ************************ vim: set sw=4 sts=4 et: ************************ //