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