meshCutter.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-2022 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 "meshCutter.H"
27 #include "polyMesh.H"
28 #include "polyTopoChange.H"
29 #include "cellCuts.H"
30 #include "polyTopoChangeMap.H"
31 #include "meshTools.H"
32 #include "polyModifyFace.H"
33 #include "polyAddPoint.H"
34 #include "polyAddFace.H"
35 #include "polyAddCell.H"
36 #include "syncTools.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(meshCutter, 0);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * //
47 
48 bool Foam::meshCutter::uses(const labelList& elems1, const labelList& elems2)
49 {
50  forAll(elems1, elemI)
51  {
52  if (findIndex(elems2, elems1[elemI]) != -1)
53  {
54  return true;
55  }
56  }
57  return false;
58 }
59 
60 
61 bool Foam::meshCutter::isIn
62 (
63  const edge& twoCuts,
64  const labelList& cuts
65 )
66 {
67  label index = findIndex(cuts, twoCuts[0]);
68 
69  if (index == -1)
70  {
71  return false;
72  }
73 
74  return
75  (
76  cuts[cuts.fcIndex(index)] == twoCuts[1]
77  || cuts[cuts.rcIndex(index)] == twoCuts[1]
78  );
79 }
80 
81 
82 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
83 
84 Foam::label Foam::meshCutter::findCutCell
85 (
86  const cellCuts& cuts,
87  const labelList& cellLabels
88 ) const
89 {
90  forAll(cellLabels, labelI)
91  {
92  label celli = cellLabels[labelI];
93 
94  if (cuts.cellLoops()[celli].size())
95  {
96  return celli;
97  }
98  }
99  return -1;
100 }
101 
102 
103 Foam::label Foam::meshCutter::findInternalFacePoint
104 (
105  const labelList& pointLabels
106 ) const
107 {
108  forAll(pointLabels, labelI)
109  {
110  label pointi = pointLabels[labelI];
111 
112  const labelList& pFaces = mesh().pointFaces()[pointi];
113 
114  forAll(pFaces, pFacei)
115  {
116  label facei = pFaces[pFacei];
117 
118  if (mesh().isInternalFace(facei))
119  {
120  return pointi;
121  }
122  }
123  }
124 
125  if (pointLabels.empty())
126  {
128  << "Empty pointLabels" << abort(FatalError);
129  }
130 
131  return -1;
132 }
133 
134 
135 void Foam::meshCutter::faceCells
136 (
137  const cellCuts& cuts,
138  const label facei,
139  label& own,
140  label& nei
141 ) const
142 {
143  const labelListList& anchorPts = cuts.cellAnchorPoints();
144  const labelListList& cellLoops = cuts.cellLoops();
145 
146  const face& f = mesh().faces()[facei];
147 
148  own = mesh().faceOwner()[facei];
149 
150  if (cellLoops[own].size() && uses(f, anchorPts[own]))
151  {
152  own = addedCells_[own];
153  }
154 
155  nei = -1;
156 
157  if (mesh().isInternalFace(facei))
158  {
159  nei = mesh().faceNeighbour()[facei];
160 
161  if (cellLoops[nei].size() && uses(f, anchorPts[nei]))
162  {
163  nei = addedCells_[nei];
164  }
165  }
166 }
167 
168 
169 void Foam::meshCutter::getFaceInfo
170 (
171  const label facei,
172  label& patchID,
173  label& zoneID,
174  label& zoneFlip
175 ) const
176 {
177  patchID = -1;
178 
179  if (!mesh().isInternalFace(facei))
180  {
181  patchID = mesh().boundaryMesh().whichPatch(facei);
182  }
183 
184  zoneID = mesh().faceZones().whichZone(facei);
185 
186  zoneFlip = false;
187 
188  if (zoneID >= 0)
189  {
190  const faceZone& fZone = mesh().faceZones()[zoneID];
191 
192  zoneFlip = fZone.flipMap()[fZone.whichFace(facei)];
193  }
194 }
195 
196 
197 void Foam::meshCutter::addFace
198 (
199  polyTopoChange& meshMod,
200  const label facei,
201  const face& newFace,
202  const label own,
203  const label nei
204 )
205 {
206  label patchID, zoneID, zoneFlip;
207 
208  getFaceInfo(facei, patchID, zoneID, zoneFlip);
209 
210  if ((nei == -1) || (own < nei))
211  {
212  // Ordering ok.
213  if (debug & 2)
214  {
215  Pout<< "Adding face " << newFace
216  << " with new owner:" << own
217  << " with new neighbour:" << nei
218  << " patchID:" << patchID
219  << " zoneID:" << zoneID
220  << " zoneFlip:" << zoneFlip
221  << endl;
222  }
223 
224  meshMod.setAction
225  (
226  polyAddFace
227  (
228  newFace, // face
229  own, // owner
230  nei, // neighbour
231  -1, // master point
232  -1, // master edge
233  facei, // master face for addition
234  false, // flux flip
235  patchID, // patch for face
236  zoneID, // zone for face
237  zoneFlip // face zone flip
238  )
239  );
240  }
241  else
242  {
243  // Reverse owner/neighbour
244  if (debug & 2)
245  {
246  Pout<< "Adding (reversed) face " << newFace.reverseFace()
247  << " with new owner:" << nei
248  << " with new neighbour:" << own
249  << " patchID:" << patchID
250  << " zoneID:" << zoneID
251  << " zoneFlip:" << zoneFlip
252  << endl;
253  }
254 
255  meshMod.setAction
256  (
257  polyAddFace
258  (
259  newFace.reverseFace(), // face
260  nei, // owner
261  own, // neighbour
262  -1, // master point
263  -1, // master edge
264  facei, // master face for addition
265  false, // flux flip
266  patchID, // patch for face
267  zoneID, // zone for face
268  zoneFlip // face zone flip
269  )
270  );
271  }
272 }
273 
274 
275 void Foam::meshCutter::modFace
276 (
277  polyTopoChange& meshMod,
278  const label facei,
279  const face& newFace,
280  const label own,
281  const label nei
282 )
283 {
284  label patchID, zoneID, zoneFlip;
285 
286  getFaceInfo(facei, patchID, zoneID, zoneFlip);
287 
288  if
289  (
290  (own != mesh().faceOwner()[facei])
291  || (
292  mesh().isInternalFace(facei)
293  && (nei != mesh().faceNeighbour()[facei])
294  )
295  || (newFace != mesh().faces()[facei])
296  )
297  {
298  if (debug & 2)
299  {
300  Pout<< "Modifying face " << facei
301  << " old vertices:" << mesh().faces()[facei]
302  << " new vertices:" << newFace
303  << " new owner:" << own
304  << " new neighbour:" << nei
305  << " new zoneID:" << zoneID
306  << " new zoneFlip:" << zoneFlip
307  << endl;
308  }
309 
310  if ((nei == -1) || (own < nei))
311  {
312  meshMod.setAction
313  (
314  polyModifyFace
315  (
316  newFace, // modified face
317  facei, // label of face being modified
318  own, // owner
319  nei, // neighbour
320  false, // face flip
321  patchID, // patch for face
322  false, // remove from zone
323  zoneID, // zone for face
324  zoneFlip // face flip in zone
325  )
326  );
327  }
328  else
329  {
330  meshMod.setAction
331  (
332  polyModifyFace
333  (
334  newFace.reverseFace(), // modified face
335  facei, // label of face being modified
336  nei, // owner
337  own, // neighbour
338  false, // face flip
339  patchID, // patch for face
340  false, // remove from zone
341  zoneID, // zone for face
342  zoneFlip // face flip in zone
343  )
344  );
345  }
346  }
347 }
348 
349 
350 void Foam::meshCutter::copyFace
351 (
352  const face& f,
353  const label startFp,
354  const label endFp,
355  face& newFace
356 ) const
357 {
358  label fp = startFp;
359 
360  label newFp = 0;
361 
362  while (fp != endFp)
363  {
364  newFace[newFp++] = f[fp];
365 
366  fp = (fp + 1) % f.size();
367  }
368  newFace[newFp] = f[fp];
369 }
370 
371 
372 void Foam::meshCutter::splitFace
373 (
374  const face& f,
375  const label v0,
376  const label v1,
377 
378  face& f0,
379  face& f1
380 ) const
381 {
382  // Check if we find any new vertex which is part of the splitEdge.
383  label startFp = findIndex(f, v0);
384 
385  if (startFp == -1)
386  {
388  << "Cannot find vertex (new numbering) " << v0
389  << " on face " << f
390  << abort(FatalError);
391  }
392 
393  label endFp = findIndex(f, v1);
394 
395  if (endFp == -1)
396  {
398  << "Cannot find vertex (new numbering) " << v1
399  << " on face " << f
400  << abort(FatalError);
401  }
402 
403 
404  f0.setSize((endFp + 1 + f.size() - startFp) % f.size());
405  f1.setSize(f.size() - f0.size() + 2);
406 
407  copyFace(f, startFp, endFp, f0);
408  copyFace(f, endFp, startFp, f1);
409 }
410 
411 
412 Foam::face Foam::meshCutter::addEdgeCutsToFace(const label facei) const
413 {
414  const face& f = mesh().faces()[facei];
415 
416  face newFace(2 * f.size());
417 
418  label newFp = 0;
419 
420  forAll(f, fp)
421  {
422  // Duplicate face vertex .
423  newFace[newFp++] = f[fp];
424 
425  // Check if edge has been cut.
426  label fp1 = f.fcIndex(fp);
427 
428  HashTable<label, edge, Hash<edge>>::const_iterator fnd =
429  addedPoints_.find(edge(f[fp], f[fp1]));
430 
431  if (fnd != addedPoints_.end())
432  {
433  // edge has been cut. Introduce new vertex.
434  newFace[newFp++] = fnd();
435  }
436  }
437 
438  newFace.setSize(newFp);
439 
440  return newFace;
441 }
442 
443 
444 Foam::face Foam::meshCutter::loopToFace
445 (
446  const label celli,
447  const labelList& loop
448 ) const
449 {
450  face newFace(2*loop.size());
451 
452  label newFacei = 0;
453 
454  forAll(loop, fp)
455  {
456  label cut = loop[fp];
457 
458  if (isEdge(cut))
459  {
460  label edgeI = getEdge(cut);
461 
462  const edge& e = mesh().edges()[edgeI];
463 
464  label vertI = addedPoints_[e];
465 
466  newFace[newFacei++] = vertI;
467  }
468  else
469  {
470  // cut is vertex.
471  label vertI = getVertex(cut);
472 
473  newFace[newFacei++] = vertI;
474 
475  label nextCut = loop[loop.fcIndex(fp)];
476 
477  if (!isEdge(nextCut))
478  {
479  // From vertex to vertex -> cross cut only if no existing edge.
480  label nextVertI = getVertex(nextCut);
481 
482  label edgeI = meshTools::findEdge(mesh(), vertI, nextVertI);
483 
484  if (edgeI != -1)
485  {
486  // Existing edge. Insert split-edge point if any.
487  HashTable<label, edge, Hash<edge>>::const_iterator fnd =
488  addedPoints_.find(mesh().edges()[edgeI]);
489 
490  if (fnd != addedPoints_.end())
491  {
492  newFace[newFacei++] = fnd();
493  }
494  }
495  }
496  }
497  }
498  newFace.setSize(newFacei);
499 
500  return newFace;
501 }
502 
503 
504 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
505 
507 :
508  edgeVertex(mesh),
509  addedCells_(),
510  addedFaces_(),
511  addedPoints_()
512 
513 {}
514 
515 
516 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
517 
519 {}
520 
521 
522 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
523 
525 (
526  const cellCuts& cuts,
527  polyTopoChange& meshMod
528 )
529 {
530  // Clear and size maps here since mesh size will change.
531  addedCells_.clear();
532  addedCells_.resize(cuts.nLoops());
533 
534  addedFaces_.clear();
535  addedFaces_.resize(cuts.nLoops());
536 
537  addedPoints_.clear();
538  addedPoints_.resize(cuts.nLoops());
539 
540  if (returnReduce(cuts.nLoops(), sumOp<label>()) == 0)
541  {
542  return;
543  }
544 
545  const labelListList& anchorPts = cuts.cellAnchorPoints();
546  const labelListList& cellLoops = cuts.cellLoops();
547 
548  if (debug)
549  {
550  // Check that any edge is cut only if any cell using it is cut
551  boolList edgeOnCutCell(mesh().nEdges(), false);
552  forAll(cuts.cellLoops(), celli)
553  {
554  if (cuts.cellLoops()[celli].size())
555  {
556  const labelList& cEdges = mesh().cellEdges(celli);
557  forAll(cEdges, i)
558  {
559  edgeOnCutCell[cEdges[i]] = true;
560  }
561  }
562  }
563  syncTools::syncEdgeList(mesh(), edgeOnCutCell, orEqOp<bool>(), false);
564 
565  forAll(cuts.edgeIsCut(), edgeI)
566  {
567  if (cuts.edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
568  {
569  const edge& e = mesh().edges()[edgeI];
570 
572  << "Problem: cut edge but none of the cells using"
573  << " it is cut\n"
574  << "edge:" << edgeI << " verts:" << e
575  << " at:" << e.line(mesh().points())
576  << endl; // abort(FatalError);
577  }
578  }
579  }
580 
581 
582  //
583  // Add new points along cut edges.
584  //
585 
586  forAll(cuts.edgeIsCut(), edgeI)
587  {
588  if (cuts.edgeIsCut()[edgeI])
589  {
590  const edge& e = mesh().edges()[edgeI];
591 
592  // One of the edge end points should be master point of nbCelli.
593  label masterPointi = e.start();
594 
595  const point& v0 = mesh().points()[e.start()];
596  const point& v1 = mesh().points()[e.end()];
597 
598  scalar weight = cuts.edgeWeight()[edgeI];
599 
600  point newPt = weight*v1 + (1.0-weight)*v0;
601 
602  label addedPointi =
603  meshMod.setAction
604  (
606  (
607  newPt, // point
608  masterPointi, // master point
609  -1, // zone for point
610  true // supports a cell
611  )
612  );
613 
614  // Store on (hash of) edge.
615  addedPoints_.insert(e, addedPointi);
616 
617  if (debug & 2)
618  {
619  Pout<< "Added point " << addedPointi
620  << " to vertex "
621  << masterPointi << " of edge " << edgeI
622  << " vertices " << e << endl;
623  }
624  }
625  }
626 
627  //
628  // Add cells (on 'anchor' side of cell)
629  //
630 
631  forAll(cellLoops, celli)
632  {
633  if (cellLoops[celli].size())
634  {
635  // Add a cell to the existing cell
636  label addedCelli =
637  meshMod.setAction
638  (
640  (
641  -1, // master point
642  -1, // master edge
643  -1, // master face
644  celli, // master cell
645  mesh().cellZones().whichZone(celli) // zone for cell
646  )
647  );
648 
649  addedCells_.insert(celli, addedCelli);
650 
651  if (debug & 2)
652  {
653  Pout<< "Added cell " << addedCells_[celli] << " to cell "
654  << celli << endl;
655  }
656  }
657  }
658 
659 
660  //
661  // For all cut cells add an internal face
662  //
663 
664  forAll(cellLoops, celli)
665  {
666  const labelList& loop = cellLoops[celli];
667 
668  if (loop.size())
669  {
670  // Convert loop (=list of cuts) into proper face.
671  // Orientation should already be ok. (done by cellCuts)
672  //
673  face newFace(loopToFace(celli, loop));
674 
675  // Pick any anchor point on cell
676  label masterPointi = findInternalFacePoint(anchorPts[celli]);
677 
678  label addedFacei =
679  meshMod.setAction
680  (
682  (
683  newFace, // face
684  celli, // owner
685  addedCells_[celli], // neighbour
686  masterPointi, // master point
687  -1, // master edge
688  -1, // master face for addition
689  false, // flux flip
690  -1, // patch for face
691  -1, // zone for face
692  false // face zone flip
693  )
694  );
695 
696  addedFaces_.insert(celli, addedFacei);
697 
698  if (debug & 2)
699  {
700  // Gets edgeweights of loop
701  scalarField weights(loop.size());
702  forAll(loop, i)
703  {
704  label cut = loop[i];
705 
706  weights[i] =
707  (
708  isEdge(cut)
709  ? cuts.edgeWeight()[getEdge(cut)]
710  : -great
711  );
712  }
713 
714  Pout<< "Added splitting face " << newFace << " index:"
715  << addedFacei
716  << " to owner " << celli
717  << " neighbour " << addedCells_[celli]
718  << " from Loop:";
719  writeCuts(Pout, loop, weights);
720  Pout<< endl;
721  }
722  }
723  }
724 
725 
726  //
727  // Modify faces on the outside and create new ones
728  // (in effect split old faces into two)
729  //
730 
731  // Maintain whether face has been updated (for -split edges
732  // -new owner/neighbour)
733  boolList faceUptodate(mesh().nFaces(), false);
734 
735  const Map<edge>& faceSplitCuts = cuts.faceSplitCut();
736 
737  forAllConstIter(Map<edge>, faceSplitCuts, iter)
738  {
739  label facei = iter.key();
740 
741  // Renumber face to include split edges.
742  face newFace(addEdgeCutsToFace(facei));
743 
744  // Edge splitting the face. Convert cuts to new vertex numbering.
745  const edge& splitEdge = iter();
746 
747  label cut0 = splitEdge[0];
748 
749  label v0;
750  if (isEdge(cut0))
751  {
752  label edgeI = getEdge(cut0);
753  v0 = addedPoints_[mesh().edges()[edgeI]];
754  }
755  else
756  {
757  v0 = getVertex(cut0);
758  }
759 
760  label cut1 = splitEdge[1];
761  label v1;
762  if (isEdge(cut1))
763  {
764  label edgeI = getEdge(cut1);
765  v1 = addedPoints_[mesh().edges()[edgeI]];
766  }
767  else
768  {
769  v1 = getVertex(cut1);
770  }
771 
772  // Split face along the elements of the splitEdge.
773  face f0, f1;
774  splitFace(newFace, v0, v1, f0, f1);
775 
776  label own = mesh().faceOwner()[facei];
777 
778  label nei = -1;
779 
780  if (mesh().isInternalFace(facei))
781  {
782  nei = mesh().faceNeighbour()[facei];
783  }
784 
785  if (debug & 2)
786  {
787  Pout<< "Split face " << mesh().faces()[facei]
788  << " own:" << own << " nei:" << nei
789  << " into f0:" << f0
790  << " and f1:" << f1 << endl;
791  }
792 
793  // Check which face uses anchorPoints (connects to addedCell)
794  // and which one doesn't (connects to original cell)
795 
796  // Bit tricky. We have to know whether this faceSplit splits owner/
797  // neighbour or both. Even if cell is cut we have to make sure this is
798  // the one that cuts it (this face cut might not be the one splitting
799  // the cell)
800 
801  const face& f = mesh().faces()[facei];
802 
803  label f0Owner = -1;
804  label f1Owner = -1;
805 
806  if (cellLoops[own].empty())
807  {
808  f0Owner = own;
809  f1Owner = own;
810  }
811  else if (isIn(splitEdge, cellLoops[own]))
812  {
813  // Owner is cut by this splitCut. See which of f0, f1 gets
814  // owner, which gets addedCells_[owner]
815  if (uses(f0, anchorPts[own]))
816  {
817  f0Owner = addedCells_[own];
818  f1Owner = own;
819  }
820  else
821  {
822  f0Owner = own;
823  f1Owner = addedCells_[own];
824  }
825  }
826  else
827  {
828  // Owner not cut by this splitCut. Check on original face whether
829  // use anchorPts.
830  if (uses(f, anchorPts[own]))
831  {
832  label newCelli = addedCells_[own];
833  f0Owner = newCelli;
834  f1Owner = newCelli;
835  }
836  else
837  {
838  f0Owner = own;
839  f1Owner = own;
840  }
841  }
842 
843 
844  label f0Neighbour = -1;
845  label f1Neighbour = -1;
846 
847  if (nei != -1)
848  {
849  if (cellLoops[nei].empty())
850  {
851  f0Neighbour = nei;
852  f1Neighbour = nei;
853  }
854  else if (isIn(splitEdge, cellLoops[nei]))
855  {
856  // Neighbour is cut by this splitCut. See which of f0, f1
857  // gets which neighbour/addedCells_[neighbour]
858  if (uses(f0, anchorPts[nei]))
859  {
860  f0Neighbour = addedCells_[nei];
861  f1Neighbour = nei;
862  }
863  else
864  {
865  f0Neighbour = nei;
866  f1Neighbour = addedCells_[nei];
867  }
868  }
869  else
870  {
871  // neighbour not cut by this splitCut. Check on original face
872  // whether use anchorPts.
873  if (uses(f, anchorPts[nei]))
874  {
875  label newCelli = addedCells_[nei];
876  f0Neighbour = newCelli;
877  f1Neighbour = newCelli;
878  }
879  else
880  {
881  f0Neighbour = nei;
882  f1Neighbour = nei;
883  }
884  }
885  }
886 
887  // f0 is the added face, f1 the modified one
888  addFace(meshMod, facei, f0, f0Owner, f0Neighbour);
889 
890  modFace(meshMod, facei, f1, f1Owner, f1Neighbour);
891 
892  faceUptodate[facei] = true;
893  }
894 
895 
896  //
897  // Faces that have not been split but just appended to. Are guaranteed
898  // to be reachable from an edgeCut.
899  //
900 
901  const boolList& edgeIsCut = cuts.edgeIsCut();
902 
903  forAll(edgeIsCut, edgeI)
904  {
905  if (edgeIsCut[edgeI])
906  {
907  const labelList& eFaces = mesh().edgeFaces()[edgeI];
908 
909  forAll(eFaces, i)
910  {
911  label facei = eFaces[i];
912 
913  if (!faceUptodate[facei])
914  {
915  // Renumber face to include split edges.
916  face newFace(addEdgeCutsToFace(facei));
917 
918  if (debug & 2)
919  {
920  Pout<< "Added edge cuts to face " << facei
921  << " f:" << mesh().faces()[facei]
922  << " newFace:" << newFace << endl;
923  }
924 
925  // Get (new or original) owner and neighbour of facei
926  label own, nei;
927  faceCells(cuts, facei, own, nei);
928 
929  modFace(meshMod, facei, newFace, own, nei);
930 
931  faceUptodate[facei] = true;
932  }
933  }
934  }
935  }
936 
937 
938  //
939  // Correct any original faces on split cell for new neighbour/owner
940  //
941 
942  forAll(cellLoops, celli)
943  {
944  if (cellLoops[celli].size())
945  {
946  const labelList& cllFaces = mesh().cells()[celli];
947 
948  forAll(cllFaces, cllFacei)
949  {
950  label facei = cllFaces[cllFacei];
951 
952  if (!faceUptodate[facei])
953  {
954  // Update face with new owner/neighbour (if any)
955  const face& f = mesh().faces()[facei];
956 
957  if (debug && (f != addEdgeCutsToFace(facei)))
958  {
960  << "Problem: edges added to face which does "
961  << " not use a marked cut" << endl
962  << "facei:" << facei << endl
963  << "face:" << f << endl
964  << "newFace:" << addEdgeCutsToFace(facei)
965  << abort(FatalError);
966  }
967 
968  // Get (new or original) owner and neighbour of facei
969  label own, nei;
970  faceCells(cuts, facei, own, nei);
971 
972  modFace
973  (
974  meshMod,
975  facei,
976  f,
977  own,
978  nei
979  );
980 
981  faceUptodate[facei] = true;
982  }
983  }
984  }
985  }
986 
987  if (debug)
988  {
989  Pout<< "meshCutter:" << nl
990  << " cells split:" << addedCells_.size() << nl
991  << " faces added:" << addedFaces_.size() << nl
992  << " points added on edges:" << addedPoints_.size() << nl
993  << endl;
994  }
995 }
996 
997 
999 {
1000  // Update stored labels for mesh change.
1001 
1002  {
1003  // Create copy since new label might (temporarily) clash with existing
1004  // key.
1005  Map<label> newAddedCells(addedCells_.size());
1006 
1007  forAllConstIter(Map<label>, addedCells_, iter)
1008  {
1009  label celli = iter.key();
1010  label newCelli = map.reverseCellMap()[celli];
1011 
1012  label addedCelli = iter();
1013 
1014  label newAddedCelli = map.reverseCellMap()[addedCelli];
1015 
1016  if (newCelli >= 0 && newAddedCelli >= 0)
1017  {
1018  if
1019  (
1020  (debug & 2)
1021  && (newCelli != celli || newAddedCelli != addedCelli)
1022  )
1023  {
1024  Pout<< "meshCutter::topoChange :"
1025  << " updating addedCell for cell " << celli
1026  << " from " << addedCelli
1027  << " to " << newAddedCelli << endl;
1028  }
1029  newAddedCells.insert(newCelli, newAddedCelli);
1030  }
1031  }
1032 
1033  // Copy
1034  addedCells_.transfer(newAddedCells);
1035  }
1036 
1037  {
1038  Map<label> newAddedFaces(addedFaces_.size());
1039 
1040  forAllConstIter(Map<label>, addedFaces_, iter)
1041  {
1042  label celli = iter.key();
1043  label newCelli = map.reverseCellMap()[celli];
1044 
1045  label addedFacei = iter();
1046 
1047  label newAddedFacei = map.reverseFaceMap()[addedFacei];
1048 
1049  if ((newCelli >= 0) && (newAddedFacei >= 0))
1050  {
1051  if
1052  (
1053  (debug & 2)
1054  && (newCelli != celli || newAddedFacei != addedFacei)
1055  )
1056  {
1057  Pout<< "meshCutter::topoChange :"
1058  << " updating addedFace for cell " << celli
1059  << " from " << addedFacei
1060  << " to " << newAddedFacei
1061  << endl;
1062  }
1063  newAddedFaces.insert(newCelli, newAddedFacei);
1064  }
1065  }
1066 
1067  // Copy
1068  addedFaces_.transfer(newAddedFaces);
1069  }
1070 
1071  {
1072  HashTable<label, edge, Hash<edge>> newAddedPoints(addedPoints_.size());
1073 
1074  for
1075  (
1076  HashTable<label, edge, Hash<edge>>::const_iterator iter =
1077  addedPoints_.begin();
1078  iter != addedPoints_.end();
1079  ++iter
1080  )
1081  {
1082  const edge& e = iter.key();
1083 
1084  label newStart = map.reversePointMap()[e.start()];
1085 
1086  label newEnd = map.reversePointMap()[e.end()];
1087 
1088  label addedPointi = iter();
1089 
1090  label newAddedPointi = map.reversePointMap()[addedPointi];
1091 
1092  if ((newStart >= 0) && (newEnd >= 0) && (newAddedPointi >= 0))
1093  {
1094  edge newE = edge(newStart, newEnd);
1095 
1096  if
1097  (
1098  (debug & 2)
1099  && (e != newE || newAddedPointi != addedPointi)
1100  )
1101  {
1102  Pout<< "meshCutter::topoChange :"
1103  << " updating addedPoints for edge " << e
1104  << " from " << addedPointi
1105  << " to " << newAddedPointi
1106  << endl;
1107  }
1108 
1109  newAddedPoints.insert(newE, newAddedPointi);
1110  }
1111  }
1112 
1113  // Copy
1114  addedPoints_.transfer(newAddedPoints);
1115  }
1116 }
1117 
1118 
1119 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
const labelListList & cellEdges() const
const boolList & edgeIsCut() const
Is edge cut.
Definition: cellCuts.H:534
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
meshCutter(const polyMesh &mesh)
Construct from mesh.
Definition: meshCutter.C:506
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:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
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
Holds information (coordinate and normal) regarding nearest wall point.
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
const cellList & cells() const
void transfer(HashTable< T, Key, Hash > &)
Transfer the contents of the argument table into this table.
Definition: HashTable.C:513
fvMesh & mesh
Description of cuts across cells.
Definition: cellCuts.H:108
scalar f1
Definition: createFields.H:15
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
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:1211
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
Class containing data for cell addition.
Definition: polyAddCell.H:46
const labelList & reverseFaceMap() const
Reverse face map.
linePointRef line(const pointField &) const
Return edge line.
Definition: edgeI.H:187
static const labelSphericalTensor labelI(1)
Identity labelTensor.
const pointField & points
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
void clear()
Clear all entries from table.
Definition: HashTable.C:468
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
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
const labelList & reversePointMap() const
Reverse point map.
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1224
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
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
Definition: meshCutter.C:525
const polyMesh & mesh() const
Definition: edgeVertex.H:93
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
const labelList & reverseCellMap() const
Reverse cell map.
static bool isEdge(const primitiveMesh &mesh, const label eVert)
Is eVert an edge?
Definition: edgeVertex.H:102
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:540
label nLoops() const
Number of valid cell loops.
Definition: cellCuts.H:569
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: MeshZones.C:221
const Map< edge > & faceSplitCut() const
Gives for split face the two cuts that split the face into two.
Definition: cellCuts.H:557
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
Definition: cellCuts.H:563
#define WarningInFunction
Report a warning using Foam::Warning.
label end() const
Return end vertex label.
Definition: edgeI.H:92
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
~meshCutter()
Destructor.
Definition: meshCutter.C:518
const meshFaceZones & faceZones() const
Return face zones.
Definition: polyMesh.H:495
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:52
const labelListList & pointFaces() const
void resize(const label newSize)
Resize the hash table for efficiency.
Definition: HashTable.C:432
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void topoChange(const polyTopoChangeMap &)
Force recalculation of locally stored data on topological change.
Definition: meshCutter.C:998
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
const labelListList & edgeFaces() const
label start() const
Return start vertex label.
Definition: edgeI.H:81
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:575
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.