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