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