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