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