cellCuts.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-2023 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 "cellCuts.H"
27 #include "polyMesh.H"
28 #include "Time.H"
29 #include "ListOps.H"
30 #include "cellLooper.H"
31 #include "refineCell.H"
32 #include "meshTools.H"
33 #include "geomCellLooper.H"
34 #include "OFstream.H"
35 #include "plane.H"
36 #include "syncTools.H"
37 #include "dummyTransform.H"
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
44 
45  //- Template specialisation for pTraits<edge> so we can use syncTools
46  // functionality
47  template<>
48  class pTraits<edge>
49  {
50  public:
51 
52  //- Component type
53  typedef edge cmptType;
54  };
55 }
56 
57 
58 // * * * * * * * * * * * * * Private Static Functions * * * * * * * * * * * //
59 
60 Foam::label Foam::cellCuts::findPartIndex
61 (
62  const labelList& elems,
63  const label nElems,
64  const label val
65 )
66 {
67  for (label i = 0; i < nElems; i++)
68  {
69  if (elems[i] == val)
70  {
71  return i;
72  }
73  }
74  return -1;
75 }
76 
77 
78 Foam::boolList Foam::cellCuts::expand
79 (
80  const label size,
81  const labelList& labels
82 )
83 {
84  boolList result(size, false);
85 
86  forAll(labels, labelI)
87  {
88  result[labels[labelI]] = true;
89  }
90  return result;
91 }
92 
93 
94 Foam::scalarField Foam::cellCuts::expand
95 (
96  const label size,
97  const labelList& labels,
98  const scalarField& weights
99 )
100 {
101  scalarField result(size, -great);
102 
103  forAll(labels, labelI)
104  {
105  result[labels[labelI]] = weights[labelI];
106  }
107  return result;
108 }
109 
110 
111 Foam::label Foam::cellCuts::firstUnique
112 (
113  const labelList& lst,
114  const Map<label>& map
115 )
116 {
117  forAll(lst, i)
118  {
119  if (!map.found(lst[i]))
120  {
121  return i;
122  }
123  }
124  return -1;
125 }
126 
127 
128 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
129 
130 void Foam::cellCuts::syncProc()
131 {
132  if (!Pstream::parRun())
133  {
134  return;
135  }
136 
137  syncTools::syncPointList(mesh(), pointIsCut_, orEqOp<bool>(), false);
138  syncTools::syncEdgeList(mesh(), edgeIsCut_, orEqOp<bool>(), false);
139  syncTools::syncEdgeList(mesh(), edgeWeight_, maxEqOp<scalar>(), -great);
140 
141  {
142  const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
143 
144  // Convert faceSplitCut into face-local data: vertex and edge w.r.t.
145  // vertex 0: (since this is same on both sides)
146  //
147  // Sending-side vertex Receiving-side vertex
148  // 0 0
149  // 1 3
150  // 2 2
151  // 3 1
152  //
153  // Sending-side edge Receiving side edge
154  // 0-1 3-0
155  // 1-2 2-3
156  // 2-3 1-2
157  // 3-0 0-1
158  //
159  // Encoding is as index:
160  // 0 : not set
161  // >0 : vertex, vertex is index-1
162  // <0 : edge, edge is -index-1
163 
164  edgeList relCuts(nBnd, edge(0, 0));
165 
166  const polyBoundaryMesh& pbm = mesh().boundaryMesh();
167 
168  forAll(pbm, patchi)
169  {
170  const polyPatch& pp = pbm[patchi];
171 
172  if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
173  {
174  forAll(pp, i)
175  {
176  label facei = pp.start()+i;
177  label bFacei = facei-mesh().nInternalFaces();
178 
179  const Map<edge>::const_iterator iter =
180  faceSplitCut_.find(facei);
181  if (iter != faceSplitCut_.end())
182  {
183  const face& f = mesh().faces()[facei];
184  const labelList& fEdges = mesh().faceEdges()[facei];
185  const edge& cuts = iter();
186 
187  forAll(cuts, i)
188  {
189  if (isEdge(cuts[i]))
190  {
191  label edgei = getEdge(cuts[i]);
192  label index = findIndex(fEdges, edgei);
193  relCuts[bFacei][i] = -index-1;
194  }
195  else
196  {
197  label index = findIndex(f, getVertex(cuts[i]));
198  relCuts[bFacei][i] = index+1;
199  }
200  }
201  }
202  }
203  }
204  }
205 
206  // Exchange
208  (
209  mesh(),
210  relCuts,
211  eqOp<edge>(),
212  dummyTransform()
213  );
214 
215  // Convert relCuts back into mesh based data
216  forAll(pbm, patchi)
217  {
218  const polyPatch& pp = pbm[patchi];
219 
220  if (isA<processorPolyPatch>(pp) || isA<cyclicPolyPatch>(pp))
221  {
222  forAll(pp, i)
223  {
224  label facei = pp.start()+i;
225  label bFacei = facei-mesh().nInternalFaces();
226 
227  const edge& relCut = relCuts[bFacei];
228  if (relCut != edge(0, 0))
229  {
230  const face& f = mesh().faces()[facei];
231  const labelList& fEdges = mesh().faceEdges()[facei];
232 
233  edge absoluteCut(0, 0);
234  forAll(relCut, i)
235  {
236  if (relCut[i] < 0)
237  {
238  label oppFp = -relCut[i]-1;
239  label fp = f.size()-1-oppFp;
240  absoluteCut[i] = edgeToEVert(fEdges[fp]);
241  }
242  else
243  {
244  label oppFp = relCut[i]-1;
245  label fp = f.size()-1-oppFp;
246  absoluteCut[i] = vertToEVert(f[fp]);
247  }
248  }
249 
250  if
251  (
252  !faceSplitCut_.insert(facei, absoluteCut)
253  && faceSplitCut_[facei] != absoluteCut
254  )
255  {
257  << "Cut " << faceSplitCut_[facei]
258  << " on face " << mesh().faceCentres()[facei]
259  << " of coupled patch " << pp.name()
260  << " is not consistent with coupled cut "
261  << absoluteCut
262  << exit(FatalError);
263  }
264  }
265  }
266  }
267  }
268  }
269 }
270 
271 
272 void Foam::cellCuts::writeUncutOBJ
273 (
274  const fileName& dir,
275  const label celli
276 ) const
277 {
278  // Cell edges
279  OFstream cutsStream(dir / "cell_" + name(celli) + ".obj");
280 
281  Pout<< "Writing cell for time " << mesh().time().name()
282  << " to " << cutsStream.name() << nl;
283 
285  (
286  cutsStream,
287  mesh().cells(),
288  mesh().faces(),
289  mesh().points(),
290  labelList(1, celli)
291  );
292 
293  // Loop cutting cell in two
294  OFstream cutStream(dir / "cellCuts_" + name(celli) + ".obj");
295 
296  Pout<< "Writing raw cuts on cell for time " << mesh().time().name()
297  << " to " << cutStream.name() << nl;
298 
299  const labelList& cPoints = mesh().cellPoints()[celli];
300 
301  forAll(cPoints, i)
302  {
303  label pointi = cPoints[i];
304  if (pointIsCut_[pointi])
305  {
306  meshTools::writeOBJ(cutStream, mesh().points()[pointi]);
307  }
308  }
309 
310  const pointField& pts = mesh().points();
311 
312  const labelList& cEdges = mesh().cellEdges()[celli];
313 
314  forAll(cEdges, i)
315  {
316  label edgeI = cEdges[i];
317 
318  if (edgeIsCut_[edgeI])
319  {
320  const edge& e = mesh().edges()[edgeI];
321 
322  const scalar w = edgeWeight_[edgeI];
323 
324  meshTools::writeOBJ(cutStream, w*pts[e[1]] + (1-w)*pts[e[0]]);
325  }
326  }
327 }
328 
329 
330 void Foam::cellCuts::writeOBJ
331 (
332  const fileName& dir,
333  const label celli,
334  const pointField& loopPoints,
335  const labelList& anchors
336 ) const
337 {
338  // Cell edges
339  OFstream cutsStream(dir / "cell_" + name(celli) + ".obj");
340 
341  Pout<< "Writing cell for time " << mesh().time().name()
342  << " to " << cutsStream.name() << nl;
343 
345  (
346  cutsStream,
347  mesh().cells(),
348  mesh().faces(),
349  mesh().points(),
350  labelList(1, celli)
351  );
352 
353 
354  // Loop cutting cell in two
355  OFstream loopStream(dir / "cellLoop_" + name(celli) + ".obj");
356 
357  Pout<< "Writing loop for time " << mesh().time().name()
358  << " to " << loopStream.name() << nl;
359 
360  label vertI = 0;
361 
362  writeOBJ(loopStream, loopPoints, vertI);
363 
364 
365  // Anchors for cell
366  OFstream anchorStream(dir / "anchors_" + name(celli) + ".obj");
367 
368  Pout<< "Writing anchors for time " << mesh().time().name()
369  << " to " << anchorStream.name() << endl;
370 
371  forAll(anchors, i)
372  {
373  meshTools::writeOBJ(anchorStream, mesh().points()[anchors[i]]);
374  }
375 }
376 
377 
378 Foam::label Foam::cellCuts::edgeEdgeToFace
379 (
380  const label celli,
381  const label edgeA,
382  const label edgeB
383 ) const
384 {
385  const labelList& cFaces = mesh().cells()[celli];
386 
387  forAll(cFaces, cFacei)
388  {
389  label facei = cFaces[cFacei];
390 
391  const labelList& fEdges = mesh().faceEdges()[facei];
392 
393  if
394  (
395  findIndex(fEdges, edgeA) != -1
396  && findIndex(fEdges, edgeB) != -1
397  )
398  {
399  return facei;
400  }
401  }
402 
403  // Coming here means the loop is illegal since the two edges
404  // are not shared by a face. We just mark loop as invalid and continue.
405 
407  << "cellCuts : Cannot find face on cell "
408  << celli << " that has both edges " << edgeA << ' ' << edgeB << endl
409  << "faces : " << cFaces << endl
410  << "edgeA : " << mesh().edges()[edgeA] << endl
411  << "edgeB : " << mesh().edges()[edgeB] << endl
412  << "Marking the loop across this cell as invalid" << endl;
413 
414  return -1;
415 }
416 
417 
418 Foam::label Foam::cellCuts::edgeVertexToFace
419 (
420  const label celli,
421  const label edgeI,
422  const label vertI
423 ) const
424 {
425  const labelList& cFaces = mesh().cells()[celli];
426 
427  forAll(cFaces, cFacei)
428  {
429  label facei = cFaces[cFacei];
430 
431  const face& f = mesh().faces()[facei];
432 
433  const labelList& fEdges = mesh().faceEdges()[facei];
434 
435  if
436  (
437  findIndex(fEdges, edgeI) != -1
438  && findIndex(f, vertI) != -1
439  )
440  {
441  return facei;
442  }
443  }
444 
446  << "cellCuts : Cannot find face on cell "
447  << celli << " that has both edge " << edgeI << " and vertex "
448  << vertI << endl
449  << "faces : " << cFaces << endl
450  << "edge : " << mesh().edges()[edgeI] << endl
451  << "Marking the loop across this cell as invalid" << endl;
452 
453  return -1;
454 }
455 
456 
457 Foam::label Foam::cellCuts::vertexVertexToFace
458 (
459  const label celli,
460  const label vertA,
461  const label vertB
462 ) const
463 {
464  const labelList& cFaces = mesh().cells()[celli];
465 
466  forAll(cFaces, cFacei)
467  {
468  label facei = cFaces[cFacei];
469 
470  const face& f = mesh().faces()[facei];
471 
472  if
473  (
474  findIndex(f, vertA) != -1
475  && findIndex(f, vertB) != -1
476  )
477  {
478  return facei;
479  }
480  }
481 
483  << "cellCuts : Cannot find face on cell "
484  << celli << " that has vertex " << vertA << " and vertex "
485  << vertB << endl
486  << "faces : " << cFaces << endl
487  << "Marking the loop across this cell as invalid" << endl;
488 
489  return -1;
490 }
491 
492 
493 void Foam::cellCuts::calcFaceCuts() const
494 {
495  if (faceCutsPtr_.valid())
496  {
498  << "faceCuts already calculated" << abort(FatalError);
499  }
500 
501  const faceList& faces = mesh().faces();
502 
503  faceCutsPtr_.reset(new labelListList(mesh().nFaces()));
504  labelListList& faceCuts = faceCutsPtr_();
505 
506  for (label facei = 0; facei < mesh().nFaces(); facei++)
507  {
508  const face& f = faces[facei];
509 
510  // Big enough storage (possibly all points and all edges cut). Shrink
511  // later on.
512  labelList& cuts = faceCuts[facei];
513 
514  cuts.setSize(2*f.size());
515 
516  label cutI = 0;
517 
518  // Do point-edge-point walk over face and collect all cuts.
519  // Problem is that we want to start from one of the endpoints of a
520  // string of connected cuts; we don't want to start somewhere in the
521  // middle.
522 
523  // Pass1: find first point cut not preceded by a cut.
524  label startFp = -1;
525 
526  forAll(f, fp)
527  {
528  label v0 = f[fp];
529 
530  if (pointIsCut_[v0])
531  {
532  label vMin1 = f[f.rcIndex(fp)];
533 
534  if
535  (
536  !pointIsCut_[vMin1]
537  && !edgeIsCut_[findEdge(facei, v0, vMin1)]
538  )
539  {
540  cuts[cutI++] = vertToEVert(v0);
541  startFp = f.fcIndex(fp);
542  break;
543  }
544  }
545  }
546 
547  // Pass2: first edge cut not preceded by point cut
548  if (startFp == -1)
549  {
550  forAll(f, fp)
551  {
552  label fp1 = f.fcIndex(fp);
553 
554  label v0 = f[fp];
555  label v1 = f[fp1];
556 
557  label edgeI = findEdge(facei, v0, v1);
558 
559  if (edgeIsCut_[edgeI] && !pointIsCut_[v0])
560  {
561  cuts[cutI++] = edgeToEVert(edgeI);
562  startFp = fp1;
563  break;
564  }
565  }
566  }
567 
568  // Pass3: nothing found so far. Either face is not cut at all or all
569  // vertices are cut. Start from 0.
570  if (startFp == -1)
571  {
572  startFp = 0;
573  }
574 
575  // Store all cuts starting from startFp;
576  label fp = startFp;
577 
578  bool allVerticesCut = true;
579 
580  forAll(f, i)
581  {
582  label fp1 = f.fcIndex(fp);
583 
584  // Get the three items: current vertex, next vertex and edge
585  // in between
586  label v0 = f[fp];
587  label v1 = f[fp1];
588  label edgeI = findEdge(facei, v0, v1);
589 
590  if (pointIsCut_[v0])
591  {
592  cuts[cutI++] = vertToEVert(v0);
593  }
594  else
595  {
596  // For check if all vertices have been cut (= illegal)
597  allVerticesCut = false;
598  }
599 
600  if (edgeIsCut_[edgeI])
601  {
602  cuts[cutI++] = edgeToEVert(edgeI);
603  }
604 
605  fp = fp1;
606  }
607 
608 
609  if (allVerticesCut)
610  {
612  << "Face " << facei << " vertices " << f
613  << " has all its vertices cut. Not cutting face." << endl;
614 
615  cutI = 0;
616  }
617 
618  // Remove duplicate starting point
619  if (cutI > 1 && cuts[cutI-1] == cuts[0])
620  {
621  cutI--;
622  }
623  cuts.setSize(cutI);
624  }
625 }
626 
627 
628 Foam::label Foam::cellCuts::findEdge
629 (
630  const label facei,
631  const label v0,
632  const label v1
633 ) const
634 {
635  const edgeList& edges = mesh().edges();
636 
637  const labelList& fEdges = mesh().faceEdges()[facei];
638 
639  forAll(fEdges, i)
640  {
641  const edge& e = edges[fEdges[i]];
642 
643  if
644  (
645  (e[0] == v0 && e[1] == v1)
646  || (e[0] == v1 && e[1] == v0)
647  )
648  {
649  return fEdges[i];
650  }
651  }
652 
653  return -1;
654 }
655 
656 
657 Foam::label Foam::cellCuts::loopFace
658 (
659  const label celli,
660  const labelList& loop
661 ) const
662 {
663  const cell& cFaces = mesh().cells()[celli];
664 
665  forAll(cFaces, cFacei)
666  {
667  label facei = cFaces[cFacei];
668 
669  const labelList& fEdges = mesh().faceEdges()[facei];
670  const face& f = mesh().faces()[facei];
671 
672  bool allOnFace = true;
673 
674  forAll(loop, i)
675  {
676  label cut = loop[i];
677 
678  if (isEdge(cut))
679  {
680  if (findIndex(fEdges, getEdge(cut)) == -1)
681  {
682  // Edge not on face. Skip face.
683  allOnFace = false;
684  break;
685  }
686  }
687  else
688  {
689  if (findIndex(f, getVertex(cut)) == -1)
690  {
691  // Vertex not on face. Skip face.
692  allOnFace = false;
693  break;
694  }
695  }
696  }
697 
698  if (allOnFace)
699  {
700  // Found face where all elements of loop are on the face.
701  return facei;
702  }
703  }
704  return -1;
705 }
706 
707 
708 bool Foam::cellCuts::walkPoint
709 (
710  const label celli,
711  const label startCut,
712 
713  const label exclude0,
714  const label exclude1,
715 
716  const label otherCut,
717 
718  label& nVisited,
719  labelList& visited
720 ) const
721 {
722  label vertI = getVertex(otherCut);
723 
724  const labelList& pFaces = mesh().pointFaces()[vertI];
725 
726  forAll(pFaces, pFacei)
727  {
728  label otherFacei = pFaces[pFacei];
729 
730  if
731  (
732  otherFacei != exclude0
733  && otherFacei != exclude1
734  && meshTools::faceOnCell(mesh(), celli, otherFacei)
735  )
736  {
737  label oldNVisited = nVisited;
738 
739  bool foundLoop =
740  walkCell
741  (
742  celli,
743  startCut,
744  otherFacei,
745  otherCut,
746  nVisited,
747  visited
748  );
749 
750  if (foundLoop)
751  {
752  return true;
753  }
754 
755  // No success. Restore state and continue
756  nVisited = oldNVisited;
757  }
758  }
759  return false;
760 }
761 
762 
763 bool Foam::cellCuts::crossEdge
764 (
765  const label celli,
766  const label startCut,
767  const label facei,
768  const label otherCut,
769 
770  label& nVisited,
771  labelList& visited
772 ) const
773 {
774  // Cross edge to other face
775  label edgeI = getEdge(otherCut);
776 
777  label otherFacei = meshTools::otherFace(mesh(), celli, facei, edgeI);
778 
779  // Store old state
780  label oldNVisited = nVisited;
781 
782  // Recurse to otherCut
783  bool foundLoop =
784  walkCell
785  (
786  celli,
787  startCut,
788  otherFacei,
789  otherCut,
790  nVisited,
791  visited
792  );
793 
794  if (foundLoop)
795  {
796  return true;
797  }
798  else
799  {
800  // No success. Restore state (i.e. backtrack)
801  nVisited = oldNVisited;
802 
803  return false;
804  }
805 }
806 
807 
808 bool Foam::cellCuts::addCut
809 (
810  const label celli,
811  const label cut,
812  label& nVisited,
813  labelList& visited
814 ) const
815 {
816  // Check for duplicate cuts.
817  if (findPartIndex(visited, nVisited, cut) != -1)
818  {
819  // Truncate (copy of) visited for ease of printing.
820  labelList truncVisited(visited);
821  truncVisited.setSize(nVisited);
822 
823  Pout<< "For cell " << celli << " : trying to add duplicate cut " << cut;
824  labelList cuts(1, cut);
825  writeCuts(Pout, cuts, loopWeights(cuts));
826 
827  Pout<< " to path:";
828  writeCuts(Pout, truncVisited, loopWeights(truncVisited));
829  Pout<< endl;
830 
831  return false;
832  }
833  else
834  {
835  visited[nVisited++] = cut;
836 
837  return true;
838  }
839 }
840 
841 
842 bool Foam::cellCuts::walkFace
843 (
844  const label celli,
845  const label startCut,
846  const label facei,
847  const label cut,
848 
849  label& lastCut,
850  label& beforeLastCut,
851  label& nVisited,
852  labelList& visited
853 ) const
854 {
855  const labelList& fCuts = faceCuts()[facei];
856 
857  if (fCuts.size() < 2)
858  {
859  return false;
860  }
861 
862  // Easy case : two cuts.
863  if (fCuts.size() == 2)
864  {
865  if (fCuts[0] == cut)
866  {
867  if (!addCut(celli, cut, nVisited, visited))
868  {
869  return false;
870  }
871 
872  beforeLastCut = cut;
873  lastCut = fCuts[1];
874 
875  return true;
876  }
877  else
878  {
879  if (!addCut(celli, cut, nVisited, visited))
880  {
881  return false;
882  }
883 
884  beforeLastCut = cut;
885  lastCut = fCuts[0];
886 
887  return true;
888  }
889  }
890 
891  // Harder case: more than 2 cuts on face.
892  // Should be start or end of string of cuts. Store all but last cut.
893  if (fCuts[0] == cut)
894  {
895  // Walk forward
896  for (label i = 0; i < fCuts.size()-1; i++)
897  {
898  if (!addCut(celli, fCuts[i], nVisited, visited))
899  {
900  return false;
901  }
902  }
903  beforeLastCut = fCuts[fCuts.size()-2];
904  lastCut = fCuts[fCuts.size()-1];
905  }
906  else if (fCuts[fCuts.size()-1] == cut)
907  {
908  for (label i = fCuts.size()-1; i >= 1; --i)
909  {
910  if (!addCut(celli, fCuts[i], nVisited, visited))
911  {
912  return false;
913  }
914  }
915  beforeLastCut = fCuts[1];
916  lastCut = fCuts[0];
917  }
918  else
919  {
921  << "In middle of cut. cell:" << celli << " face:" << facei
922  << " cuts:" << fCuts << " current cut:" << cut << endl;
923 
924  return false;
925  }
926 
927  return true;
928 }
929 
930 
931 
932 bool Foam::cellCuts::walkCell
933 (
934  const label celli,
935  const label startCut,
936  const label facei,
937  const label cut,
938 
939  label& nVisited,
940  labelList& visited
941 ) const
942 {
943  // Walk across face, storing cuts. Return last two cuts visited.
944  label lastCut = -1;
945  label beforeLastCut = -1;
946 
947 
948  if (debug & 2)
949  {
950  Pout<< "For cell:" << celli << " walked across face " << facei
951  << " from cut ";
952  labelList cuts(1, cut);
953  writeCuts(Pout, cuts, loopWeights(cuts));
954  Pout<< endl;
955  }
956 
957  bool validWalk = walkFace
958  (
959  celli,
960  startCut,
961  facei,
962  cut,
963 
964  lastCut,
965  beforeLastCut,
966  nVisited,
967  visited
968  );
969 
970  if (!validWalk)
971  {
972  return false;
973  }
974 
975  if (debug & 2)
976  {
977  Pout<< " to last cut ";
978  labelList cuts(1, lastCut);
979  writeCuts(Pout, cuts, loopWeights(cuts));
980  Pout<< endl;
981  }
982 
983  // Check if starting point reached.
984  if (lastCut == startCut)
985  {
986  if (nVisited >= 3)
987  {
988  if (debug & 2)
989  {
990  // Truncate (copy of) visited for ease of printing.
991  labelList truncVisited(visited);
992  truncVisited.setSize(nVisited);
993 
994  Pout<< "For cell " << celli << " : found closed path:";
995  writeCuts(Pout, truncVisited, loopWeights(truncVisited));
996  Pout<< " closed by " << lastCut << endl;
997  }
998 
999  return true;
1000  }
1001  else
1002  {
1003  // Loop but too short. Probably folded back on itself.
1004  return false;
1005  }
1006  }
1007 
1008 
1009  // Check last cut and one before that to determine type.
1010 
1011  // From beforeLastCut to lastCut:
1012  // - from edge to edge
1013  // (- always not along existing edge)
1014  // - from edge to vertex
1015  // - not along existing edge
1016  // - along existing edge: not allowed?
1017  // - from vertex to edge
1018  // - not along existing edge
1019  // - along existing edge. Not allowed. See above.
1020  // - from vertex to vertex
1021  // - not along existing edge
1022  // - along existing edge
1023 
1024  if (isEdge(beforeLastCut))
1025  {
1026  if (isEdge(lastCut))
1027  {
1028  // beforeLastCut=edge, lastCut=edge.
1029 
1030  // Cross lastCut (=edge) into face not facei
1031  return crossEdge
1032  (
1033  celli,
1034  startCut,
1035  facei,
1036  lastCut,
1037  nVisited,
1038  visited
1039  );
1040  }
1041  else
1042  {
1043  // beforeLastCut=edge, lastCut=vertex.
1044 
1045  // Go from lastCut to all connected faces (except facei)
1046  return walkPoint
1047  (
1048  celli,
1049  startCut,
1050  facei, // exclude0: don't cross back on current face
1051  -1, // exclude1
1052  lastCut,
1053  nVisited,
1054  visited
1055  );
1056  }
1057  }
1058  else
1059  {
1060  if (isEdge(lastCut))
1061  {
1062  // beforeLastCut=vertex, lastCut=edge.
1063  return crossEdge
1064  (
1065  celli,
1066  startCut,
1067  facei,
1068  lastCut,
1069  nVisited,
1070  visited
1071  );
1072  }
1073  else
1074  {
1075  // beforeLastCut=vertex, lastCut=vertex. Check if along existing
1076  // edge.
1077  label edgeI =
1078  findEdge
1079  (
1080  facei,
1081  getVertex(beforeLastCut),
1082  getVertex(lastCut)
1083  );
1084 
1085  if (edgeI != -1)
1086  {
1087  // Cut along existing edge. So is in fact on two faces.
1088  // Get faces on both sides of the edge to make
1089  // sure we don't fold back on to those.
1090 
1091  label f0, f1;
1092  meshTools::getEdgeFaces(mesh(), celli, edgeI, f0, f1);
1093 
1094  return walkPoint
1095  (
1096  celli,
1097  startCut,
1098  f0,
1099  f1,
1100  lastCut,
1101  nVisited,
1102  visited
1103  );
1104 
1105  }
1106  else
1107  {
1108  // Cross cut across face.
1109  return walkPoint
1110  (
1111  celli,
1112  startCut,
1113  facei, // face to exclude
1114  -1, // face to exclude
1115  lastCut,
1116  nVisited,
1117  visited
1118  );
1119  }
1120  }
1121  }
1122 }
1123 
1124 
1125 void Foam::cellCuts::calcCellLoops(const labelList& cutCells)
1126 {
1127  // Determine for every cut cell the loop (= face) it is cut by. Done by
1128  // starting from a cut edge or cut vertex and walking across faces, from
1129  // cut to cut, until starting cut hit.
1130  // If multiple loops are possible across a cell circumference takes the
1131  // first one found.
1132 
1133  // Calculate cuts per face.
1134  const labelListList& allFaceCuts = faceCuts();
1135 
1136  // Per cell the number of faces with valid cuts. Is used as quick
1137  // rejection to see if cell can be cut.
1138  labelList nCutFaces(mesh().nCells(), 0);
1139 
1140  forAll(allFaceCuts, facei)
1141  {
1142  const labelList& fCuts = allFaceCuts[facei];
1143 
1144  if (fCuts.size() == mesh().faces()[facei].size())
1145  {
1146  // Too many cuts on face. WalkCell would get very upset so disable.
1147  nCutFaces[mesh().faceOwner()[facei]] = labelMin;
1148 
1149  if (mesh().isInternalFace(facei))
1150  {
1151  nCutFaces[mesh().faceNeighbour()[facei]] = labelMin;
1152  }
1153  }
1154  else if (fCuts.size() >= 2)
1155  {
1156  // Could be valid cut. Update count for owner and neighbour.
1157  nCutFaces[mesh().faceOwner()[facei]]++;
1158 
1159  if (mesh().isInternalFace(facei))
1160  {
1161  nCutFaces[mesh().faceNeighbour()[facei]]++;
1162  }
1163  }
1164  }
1165 
1166 
1167  // Stack of visited cuts (nVisited used as stack pointer)
1168  // Size big enough.
1169  labelList visited(mesh().nPoints());
1170 
1171  forAll(cutCells, i)
1172  {
1173  label celli = cutCells[i];
1174 
1175  bool validLoop = false;
1176 
1177  // Quick rejection: has enough faces that are cut?
1178  if (nCutFaces[celli] >= 3)
1179  {
1180  const labelList& cFaces = mesh().cells()[celli];
1181 
1182  if (debug & 2)
1183  {
1184  Pout<< "cell:" << celli << " cut faces:" << endl;
1185  forAll(cFaces, i)
1186  {
1187  label facei = cFaces[i];
1188  const labelList& fCuts = allFaceCuts[facei];
1189 
1190  Pout<< " face:" << facei << " cuts:";
1191  writeCuts(Pout, fCuts, loopWeights(fCuts));
1192  Pout<< endl;
1193  }
1194  }
1195 
1196  label nVisited = 0;
1197 
1198  // Determine the first cut face to start walking from.
1199  forAll(cFaces, cFacei)
1200  {
1201  label facei = cFaces[cFacei];
1202 
1203  const labelList& fCuts = allFaceCuts[facei];
1204 
1205  // Take first or last cut of multiple on face.
1206  // Note that in calcFaceCuts
1207  // we have already made sure this is the start or end of a
1208  // string of cuts.
1209  if (fCuts.size() >= 2)
1210  {
1211  // Try walking from start of fCuts.
1212  nVisited = 0;
1213 
1214  if (debug & 2)
1215  {
1216  Pout<< "cell:" << celli
1217  << " start walk at face:" << facei
1218  << " cut:";
1219  labelList cuts(1, fCuts[0]);
1220  writeCuts(Pout, cuts, loopWeights(cuts));
1221  Pout<< endl;
1222  }
1223 
1224  validLoop =
1225  walkCell
1226  (
1227  celli,
1228  fCuts[0],
1229  facei,
1230  fCuts[0],
1231 
1232  nVisited,
1233  visited
1234  );
1235 
1236  if (validLoop)
1237  {
1238  break;
1239  }
1240 
1241  // No need to try and walk from end of cuts since
1242  // all paths already tried by walkCell.
1243  }
1244  }
1245 
1246  if (validLoop)
1247  {
1248  // Copy nVisited elements out of visited (since visited is
1249  // never truncated for efficiency reasons)
1250 
1251  labelList& loop = cellLoops_[celli];
1252 
1253  loop.setSize(nVisited);
1254 
1255  for (label i = 0; i < nVisited; i++)
1256  {
1257  loop[i] = visited[i];
1258  }
1259  }
1260  else
1261  {
1262  // Invalid loop. Leave cellLoops_[celli] zero size which
1263  // flags this.
1264  Pout<< "calcCellLoops(const labelList&) : did not find valid"
1265  << " loop for cell " << celli << endl;
1266  // Dump cell and cuts on cell.
1267  writeUncutOBJ(".", celli);
1268 
1269  cellLoops_[celli].setSize(0);
1270  }
1271  }
1272  else
1273  {
1274  // Pout<< "calcCellLoops(const labelList&) : did not find valid"
1275  // << " loop for cell " << celli << " since not enough cut faces"
1276  // << endl;
1277  cellLoops_[celli].setSize(0);
1278  }
1279  }
1280 }
1281 
1282 
1283 void Foam::cellCuts::walkEdges
1284 (
1285  const label celli,
1286  const label pointi,
1287  const label status,
1288 
1289  Map<label>& edgeStatus,
1290  Map<label>& pointStatus
1291 ) const
1292 {
1293  if (pointStatus.insert(pointi, status))
1294  {
1295  // First visit to pointi
1296 
1297  const labelList& pEdges = mesh().pointEdges()[pointi];
1298 
1299  forAll(pEdges, pEdgeI)
1300  {
1301  label edgeI = pEdges[pEdgeI];
1302 
1303  if
1304  (
1305  meshTools::edgeOnCell(mesh(), celli, edgeI)
1306  && edgeStatus.insert(edgeI, status)
1307  )
1308  {
1309  // First visit to edgeI so recurse.
1310 
1311  label v2 = mesh().edges()[edgeI].otherVertex(pointi);
1312 
1313  walkEdges(celli, v2, status, edgeStatus, pointStatus);
1314  }
1315  }
1316  }
1317 }
1318 
1319 
1321 (
1322  const labelList& cellPoints,
1323  const labelList& anchorPoints,
1324  const labelList& loop
1325 ) const
1326 {
1327  labelList newElems(cellPoints.size());
1328  label newElemI = 0;
1329 
1330  forAll(cellPoints, i)
1331  {
1332  label pointi = cellPoints[i];
1333 
1334  if
1335  (
1336  findIndex(anchorPoints, pointi) == -1
1337  && findIndex(loop, vertToEVert(pointi)) == -1
1338  )
1339  {
1340  newElems[newElemI++] = pointi;
1341  }
1342  }
1343 
1344  newElems.setSize(newElemI);
1345 
1346  return newElems;
1347 }
1348 
1349 
1350 bool Foam::cellCuts::loopAnchorConsistent
1351 (
1352  const label celli,
1353  const pointField& loopPts,
1354  const labelList& anchorPoints
1355 ) const
1356 {
1357  // Create identity face for ease of calculation of area etc.
1358  const face f(identityMap(loopPts.size()));
1359 
1360  const vector a = f.area(loopPts);
1361  const point ctr = f.centre(loopPts);
1362 
1363  // Get average position of anchor points.
1364  vector avg(Zero);
1365 
1366  forAll(anchorPoints, ptI)
1367  {
1368  avg += mesh().points()[anchorPoints[ptI]];
1369  }
1370  avg /= anchorPoints.size();
1371 
1372  if (((avg - ctr) & a) > 0)
1373  {
1374  return true;
1375  }
1376  else
1377  {
1378  return false;
1379  }
1380 }
1381 
1382 
1383 bool Foam::cellCuts::calcAnchors
1384 (
1385  const label celli,
1386  const labelList& loop,
1387  const pointField& loopPts,
1388 
1389  labelList& anchorPoints
1390 ) const
1391 {
1392  const edgeList& edges = mesh().edges();
1393 
1394  const labelList& cPoints = mesh().cellPoints()[celli];
1395  const labelList& cEdges = mesh().cellEdges()[celli];
1396  const cell& cFaces = mesh().cells()[celli];
1397 
1398  // Points on loop
1399 
1400  // Status of point:
1401  // 0 - on loop
1402  // 1 - point set 1
1403  // 2 - point set 2
1404  Map<label> pointStatus(2*cPoints.size());
1405  Map<label> edgeStatus(2*cEdges.size());
1406 
1407  // Mark loop vertices
1408  forAll(loop, i)
1409  {
1410  label cut = loop[i];
1411 
1412  if (isEdge(cut))
1413  {
1414  edgeStatus.insert(getEdge(cut), 0);
1415  }
1416  else
1417  {
1418  pointStatus.insert(getVertex(cut), 0);
1419  }
1420  }
1421  // Since edges between two cut vertices have not been marked, mark them
1422  // explicitly
1423  forAll(cEdges, i)
1424  {
1425  label edgeI = cEdges[i];
1426  const edge& e = edges[edgeI];
1427 
1428  if (pointStatus.found(e[0]) && pointStatus.found(e[1]))
1429  {
1430  edgeStatus.insert(edgeI, 0);
1431  }
1432  }
1433 
1434 
1435  // Find uncut starting vertex
1436  label uncutIndex = firstUnique(cPoints, pointStatus);
1437 
1438  if (uncutIndex == -1)
1439  {
1441  << "Invalid loop " << loop << " for cell " << celli << endl
1442  << "Can not find point on cell which is not cut by loop."
1443  << endl;
1444 
1445  writeOBJ(".", celli, loopPts, labelList(0));
1446 
1447  return false;
1448  }
1449 
1450  // Walk unset vertices and edges and mark with 1 in pointStatus, edgeStatus
1451  walkEdges(celli, cPoints[uncutIndex], 1, edgeStatus, pointStatus);
1452 
1453  // Find new uncut starting vertex
1454  uncutIndex = firstUnique(cPoints, pointStatus);
1455 
1456  if (uncutIndex == -1)
1457  {
1458  // All vertices either in loop or in anchor. So split is along single
1459  // face.
1461  << "Invalid loop " << loop << " for cell " << celli << endl
1462  << "All vertices of cell are either in loop or in anchor set"
1463  << endl;
1464 
1465  writeOBJ(".", celli, loopPts, labelList(0));
1466 
1467  return false;
1468  }
1469 
1470  // Walk unset vertices and edges and mark with 2. These are the
1471  // pointset 2.
1472  walkEdges(celli, cPoints[uncutIndex], 2, edgeStatus, pointStatus);
1473 
1474  // Collect both sets in lists.
1475  DynamicList<label> connectedPoints(cPoints.size());
1476  DynamicList<label> otherPoints(cPoints.size());
1477 
1478  forAllConstIter(Map<label>, pointStatus, iter)
1479  {
1480  if (iter() == 1)
1481  {
1482  connectedPoints.append(iter.key());
1483  }
1484  else if (iter() == 2)
1485  {
1486  otherPoints.append(iter.key());
1487  }
1488  }
1489  connectedPoints.shrink();
1490  otherPoints.shrink();
1491 
1492  // Check that all points have been used.
1493  uncutIndex = firstUnique(cPoints, pointStatus);
1494 
1495  if (uncutIndex != -1)
1496  {
1498  << "Invalid loop " << loop << " for cell " << celli
1499  << " since it splits the cell into more than two cells" << endl;
1500 
1501  writeOBJ(".", celli, loopPts, connectedPoints);
1502 
1503  return false;
1504  }
1505 
1506 
1507  // Check that both parts (connectedPoints, otherPoints) have enough faces.
1508  labelHashSet connectedFaces(2*cFaces.size());
1509  labelHashSet otherFaces(2*cFaces.size());
1510 
1511  forAllConstIter(Map<label>, pointStatus, iter)
1512  {
1513  label pointi = iter.key();
1514 
1515  const labelList& pFaces = mesh().pointFaces()[pointi];
1516 
1517  if (iter() == 1)
1518  {
1519  forAll(pFaces, pFacei)
1520  {
1521  if (meshTools::faceOnCell(mesh(), celli, pFaces[pFacei]))
1522  {
1523  connectedFaces.insert(pFaces[pFacei]);
1524  }
1525  }
1526  }
1527  else if (iter() == 2)
1528  {
1529  forAll(pFaces, pFacei)
1530  {
1531  if (meshTools::faceOnCell(mesh(), celli, pFaces[pFacei]))
1532  {
1533  otherFaces.insert(pFaces[pFacei]);
1534  }
1535  }
1536  }
1537  }
1538 
1539  if (connectedFaces.size() < 3)
1540  {
1542  << "Invalid loop " << loop << " for cell " << celli
1543  << " since would have too few faces on one side." << nl
1544  << "All faces:" << cFaces << endl;
1545 
1546  writeOBJ(".", celli, loopPts, connectedPoints);
1547 
1548  return false;
1549  }
1550 
1551  if (otherFaces.size() < 3)
1552  {
1554  << "Invalid loop " << loop << " for cell " << celli
1555  << " since would have too few faces on one side." << nl
1556  << "All faces:" << cFaces << endl;
1557 
1558  writeOBJ(".", celli, loopPts, otherPoints);
1559 
1560  return false;
1561  }
1562 
1563 
1564  // Check that faces are split into two regions and not more.
1565  // When walking across the face the only transition of pointStatus is
1566  // from set1 to loop to set2 or back. Not allowed is from set1 to loop to
1567  // set1.
1568  {
1569  forAll(cFaces, i)
1570  {
1571  label facei = cFaces[i];
1572 
1573  const face& f = mesh().faces()[facei];
1574 
1575  bool hasSet1 = false;
1576  bool hasSet2 = false;
1577 
1578  label prevStat = pointStatus[f[0]];
1579 
1580  forAll(f, fp)
1581  {
1582  label v0 = f[fp];
1583  label pStat = pointStatus[v0];
1584 
1585  if (pStat == prevStat)
1586  {
1587  }
1588  else if (pStat == 0)
1589  {
1590  // Loop.
1591  }
1592  else if (pStat == 1)
1593  {
1594  if (hasSet1)
1595  {
1596  // Second occurrence of set1.
1598  << "Invalid loop " << loop << " for cell " << celli
1599  << " since face " << f << " would be split into"
1600  << " more than two faces" << endl;
1601 
1602  writeOBJ(".", celli, loopPts, otherPoints);
1603 
1604  return false;
1605  }
1606 
1607  hasSet1 = true;
1608  }
1609  else if (pStat == 2)
1610  {
1611  if (hasSet2)
1612  {
1613  // Second occurrence of set1.
1615  << "Invalid loop " << loop << " for cell " << celli
1616  << " since face " << f << " would be split into"
1617  << " more than two faces" << endl;
1618 
1619  writeOBJ(".", celli, loopPts, otherPoints);
1620 
1621  return false;
1622  }
1623 
1624  hasSet2 = true;
1625  }
1626  else
1627  {
1629  << abort(FatalError);
1630  }
1631 
1632  prevStat = pStat;
1633 
1634 
1635  label v1 = f.nextLabel(fp);
1636  label edgeI = findEdge(facei, v0, v1);
1637 
1638  label eStat = edgeStatus[edgeI];
1639 
1640  if (eStat == prevStat)
1641  {
1642  }
1643  else if (eStat == 0)
1644  {
1645  // Loop.
1646  }
1647  else if (eStat == 1)
1648  {
1649  if (hasSet1)
1650  {
1651  // Second occurrence of set1.
1653  << "Invalid loop " << loop << " for cell " << celli
1654  << " since face " << f << " would be split into"
1655  << " more than two faces" << endl;
1656 
1657  writeOBJ(".", celli, loopPts, otherPoints);
1658 
1659  return false;
1660  }
1661 
1662  hasSet1 = true;
1663  }
1664  else if (eStat == 2)
1665  {
1666  if (hasSet2)
1667  {
1668  // Second occurrence of set1.
1670  << "Invalid loop " << loop << " for cell " << celli
1671  << " since face " << f << " would be split into"
1672  << " more than two faces" << endl;
1673 
1674  writeOBJ(".", celli, loopPts, otherPoints);
1675 
1676  return false;
1677  }
1678 
1679  hasSet2 = true;
1680  }
1681  prevStat = eStat;
1682  }
1683  }
1684  }
1685 
1686 
1687 
1688 
1689  // Check which one of point sets to use.
1690  bool loopOk = loopAnchorConsistent(celli, loopPts, connectedPoints);
1691 
1692  // if (debug)
1693  {
1694  // Additional check: are non-anchor points on other side?
1695  bool otherLoopOk = loopAnchorConsistent(celli, loopPts, otherPoints);
1696 
1697  if (loopOk == otherLoopOk)
1698  {
1699  // Both sets of points are supposedly on the same side as the
1700  // loop normal. Oops.
1701 
1703  << " For cell:" << celli
1704  << " achorpoints and nonanchorpoints are geometrically"
1705  << " on same side!" << endl
1706  << "cellPoints:" << cPoints << endl
1707  << "loop:" << loop << endl
1708  << "anchors:" << connectedPoints << endl
1709  << "otherPoints:" << otherPoints << endl;
1710 
1711  writeOBJ(".", celli, loopPts, connectedPoints);
1712  }
1713  }
1714 
1715  if (loopOk)
1716  {
1717  // connectedPoints on 'outside' of loop so these are anchor points
1718  anchorPoints.transfer(connectedPoints);
1719  connectedPoints.clear();
1720  }
1721  else
1722  {
1723  anchorPoints.transfer(otherPoints);
1724  otherPoints.clear();
1725  }
1726  return true;
1727 }
1728 
1729 
1730 Foam::pointField Foam::cellCuts::loopPoints
1731 (
1732  const labelList& loop,
1733  const scalarField& loopWeights
1734 ) const
1735 {
1736  pointField loopPts(loop.size());
1737 
1738  forAll(loop, fp)
1739  {
1740  loopPts[fp] = coord(loop[fp], loopWeights[fp]);
1741  }
1742  return loopPts;
1743 }
1744 
1745 
1746 Foam::scalarField Foam::cellCuts::loopWeights(const labelList& loop) const
1747 {
1748  scalarField weights(loop.size());
1749 
1750  forAll(loop, fp)
1751  {
1752  label cut = loop[fp];
1753 
1754  if (isEdge(cut))
1755  {
1756  label edgeI = getEdge(cut);
1757 
1758  weights[fp] = edgeWeight_[edgeI];
1759  }
1760  else
1761  {
1762  weights[fp] = -great;
1763  }
1764  }
1765  return weights;
1766 }
1767 
1768 
1769 bool Foam::cellCuts::validEdgeLoop
1770 (
1771  const labelList& loop,
1772  const scalarField& loopWeights
1773 ) const
1774 {
1775  forAll(loop, fp)
1776  {
1777  label cut = loop[fp];
1778 
1779  if (isEdge(cut))
1780  {
1781  label edgeI = getEdge(cut);
1782 
1783  // Check: cut compatible only if can be snapped to existing one.
1784  if (edgeIsCut_[edgeI])
1785  {
1786  scalar edgeLen = mesh().edges()[edgeI].mag(mesh().points());
1787 
1788  if
1789  (
1790  mag(loopWeights[fp] - edgeWeight_[edgeI])
1791  > geomCellLooper::snapTol()*edgeLen
1792  )
1793  {
1794  // Positions differ too much->would create two cuts.
1795  return false;
1796  }
1797  }
1798  }
1799  }
1800  return true;
1801 }
1802 
1803 
1804 Foam::label Foam::cellCuts::countFaceCuts
1805 (
1806  const label facei,
1807  const labelList& loop
1808 ) const
1809 {
1810  // Includes cuts through vertices and through edges.
1811  // Assumes that if edge is cut both in edgeIsCut and in loop that the
1812  // position of the cut is the same.
1813 
1814  label nCuts = 0;
1815 
1816  // Count cut vertices
1817  const face& f = mesh().faces()[facei];
1818 
1819  forAll(f, fp)
1820  {
1821  label vertI = f[fp];
1822 
1823  // Vertex already cut or mentioned in current loop.
1824  if
1825  (
1826  pointIsCut_[vertI]
1827  || (findIndex(loop, vertToEVert(vertI)) != -1)
1828  )
1829  {
1830  nCuts++;
1831  }
1832  }
1833 
1834  // Count cut edges.
1835  const labelList& fEdges = mesh().faceEdges()[facei];
1836 
1837  forAll(fEdges, fEdgeI)
1838  {
1839  label edgeI = fEdges[fEdgeI];
1840 
1841  // Edge already cut or mentioned in current loop.
1842  if
1843  (
1844  edgeIsCut_[edgeI]
1845  || (findIndex(loop, edgeToEVert(edgeI)) != -1)
1846  )
1847  {
1848  nCuts++;
1849  }
1850  }
1851 
1852  return nCuts;
1853 }
1854 
1855 
1856 bool Foam::cellCuts::conservativeValidLoop
1857 (
1858  const label celli,
1859  const labelList& loop
1860 ) const
1861 {
1862 
1863  if (loop.size() < 2)
1864  {
1865  return false;
1866  }
1867 
1868  forAll(loop, cutI)
1869  {
1870  if (isEdge(loop[cutI]))
1871  {
1872  label edgeI = getEdge(loop[cutI]);
1873 
1874  if (edgeIsCut_[edgeI])
1875  {
1876  // edge compatibility already checked.
1877  }
1878  else
1879  {
1880  // Quick rejection: vertices of edge itself cannot be cut.
1881  const edge& e = mesh().edges()[edgeI];
1882 
1883  if (pointIsCut_[e.start()] || pointIsCut_[e.end()])
1884  {
1885  return false;
1886  }
1887 
1888 
1889  // Check faces using this edge
1890  const labelList& eFaces = mesh().edgeFaces()[edgeI];
1891 
1892  forAll(eFaces, eFacei)
1893  {
1894  label nCuts = countFaceCuts(eFaces[eFacei], loop);
1895 
1896  if (nCuts > 2)
1897  {
1898  return false;
1899  }
1900  }
1901  }
1902  }
1903  else
1904  {
1905  // Vertex cut
1906 
1907  label vertI = getVertex(loop[cutI]);
1908 
1909  if (!pointIsCut_[vertI])
1910  {
1911  // New cut through vertex.
1912 
1913  // Check edges using vertex.
1914  const labelList& pEdges = mesh().pointEdges()[vertI];
1915 
1916  forAll(pEdges, pEdgeI)
1917  {
1918  label edgeI = pEdges[pEdgeI];
1919 
1920  if (edgeIsCut_[edgeI])
1921  {
1922  return false;
1923  }
1924  }
1925 
1926  // Check faces using vertex.
1927  const labelList& pFaces = mesh().pointFaces()[vertI];
1928 
1929  forAll(pFaces, pFacei)
1930  {
1931  label nCuts = countFaceCuts(pFaces[pFacei], loop);
1932 
1933  if (nCuts > 2)
1934  {
1935  return false;
1936  }
1937  }
1938  }
1939  }
1940  }
1941 
1942  // All ok.
1943  return true;
1944 }
1945 
1946 
1947 bool Foam::cellCuts::validLoop
1948 (
1949  const label celli,
1950  const labelList& loop,
1951  const scalarField& loopWeights,
1952 
1953  Map<edge>& newFaceSplitCut,
1954  labelList& anchorPoints
1955 ) const
1956 {
1957  // Determine compatibility of loop with existing cut pattern. Does not use
1958  // derived cut-addressing (faceCuts), only pointIsCut, edgeIsCut.
1959  // Adds any cross-cuts found to newFaceSplitCut and sets cell points on
1960  // one side of the loop in anchorPoints.
1961 
1962  if (loop.size() < 2)
1963  {
1964  return false;
1965  }
1966 
1967  if (debug & 4)
1968  {
1969  // Allow as fallback the 'old' loop checking where only a single
1970  // cut per face is allowed.
1971  if (!conservativeValidLoop(celli, loop))
1972  {
1973  Info << "Invalid conservative loop: " << loop << endl;
1974  return false;
1975  }
1976  }
1977 
1978  forAll(loop, fp)
1979  {
1980  label cut = loop[fp];
1981  label nextCut = loop[(fp+1) % loop.size()];
1982 
1983  // Label (if any) of face cut (so cut not along existing edge)
1984  label meshFacei = -1;
1985 
1986  if (isEdge(cut))
1987  {
1988  label edgeI = getEdge(cut);
1989 
1990  // Look one cut ahead to find if it is along existing edge.
1991 
1992  if (isEdge(nextCut))
1993  {
1994  // From edge to edge -> cross cut
1995  label nextEdgeI = getEdge(nextCut);
1996 
1997  // Find face and mark as to be split.
1998  meshFacei = edgeEdgeToFace(celli, edgeI, nextEdgeI);
1999 
2000  if (meshFacei == -1)
2001  {
2002  // Can't find face using both cut edges.
2003  return false;
2004  }
2005  }
2006  else
2007  {
2008  // From edge to vertex -> cross cut only if no existing edge.
2009 
2010  label nextVertI = getVertex(nextCut);
2011  const edge& e = mesh().edges()[edgeI];
2012 
2013  if (e.start() != nextVertI && e.end() != nextVertI)
2014  {
2015  // New edge. Find face and mark as to be split.
2016  meshFacei = edgeVertexToFace(celli, edgeI, nextVertI);
2017 
2018  if (meshFacei == -1)
2019  {
2020  // Can't find face. Illegal.
2021  return false;
2022  }
2023  }
2024  }
2025  }
2026  else
2027  {
2028  // Cut is vertex.
2029  label vertI = getVertex(cut);
2030 
2031  if (isEdge(nextCut))
2032  {
2033  // From vertex to edge -> cross cut only if no existing edge.
2034  label nextEdgeI = getEdge(nextCut);
2035 
2036  const edge& nextE = mesh().edges()[nextEdgeI];
2037 
2038  if (nextE.start() != vertI && nextE.end() != vertI)
2039  {
2040  // Should be cross cut. Find face.
2041  meshFacei = edgeVertexToFace(celli, nextEdgeI, vertI);
2042 
2043  if (meshFacei == -1)
2044  {
2045  return false;
2046  }
2047  }
2048  }
2049  else
2050  {
2051  // From vertex to vertex -> cross cut only if no existing edge.
2052  label nextVertI = getVertex(nextCut);
2053 
2054  if (meshTools::findEdge(mesh(), vertI, nextVertI) == -1)
2055  {
2056  // New edge. Find face.
2057  meshFacei = vertexVertexToFace(celli, vertI, nextVertI);
2058 
2059  if (meshFacei == -1)
2060  {
2061  return false;
2062  }
2063  }
2064  }
2065  }
2066 
2067  if (meshFacei != -1)
2068  {
2069  // meshFacei is cut across along cut-nextCut (so not along existing
2070  // edge). Check if this is compatible with existing pattern.
2071  edge cutEdge(cut, nextCut);
2072 
2073  Map<edge>::const_iterator iter = faceSplitCut_.find(meshFacei);
2074 
2075  if (iter == faceSplitCut_.end())
2076  {
2077  // Face not yet cut so insert.
2078  newFaceSplitCut.insert(meshFacei, cutEdge);
2079  }
2080  else
2081  {
2082  // Face already cut. Ok if same edge.
2083  if (iter() != cutEdge)
2084  {
2085  return false;
2086  }
2087  }
2088  }
2089  }
2090 
2091  // Is there a face on which all cuts are?
2092  label faceContainingLoop = loopFace(celli, loop);
2093 
2094  if (faceContainingLoop != -1)
2095  {
2097  << "Found loop on cell " << celli << " with all points"
2098  << " on face " << faceContainingLoop << endl;
2099 
2100  // writeOBJ(".", celli, loopPoints(loop, loopWeights), labelList(0));
2101 
2102  return false;
2103  }
2104 
2105  // Calculate anchor points
2106  // Final success is determined by whether anchor points can be determined.
2107  return calcAnchors
2108  (
2109  celli,
2110  loop,
2111  loopPoints(loop, loopWeights),
2112  anchorPoints
2113  );
2114 }
2115 
2116 
2117 void Foam::cellCuts::setFromCellLoops()
2118 {
2119  // 'Uncut' edges/vertices that are not used in loops.
2120  pointIsCut_ = false;
2121 
2122  edgeIsCut_ = false;
2123 
2124  faceSplitCut_.clear();
2125 
2126  forAll(cellLoops_, celli)
2127  {
2128  const labelList& loop = cellLoops_[celli];
2129 
2130  if (loop.size())
2131  {
2132  // Storage for cross-face cuts
2133  Map<edge> faceSplitCuts(loop.size());
2134  // Storage for points on one side of cell.
2135  labelList anchorPoints;
2136 
2137  if
2138  (
2139  !validLoop
2140  (
2141  celli,
2142  loop,
2143  loopWeights(loop),
2144  faceSplitCuts,
2145  anchorPoints
2146  )
2147  )
2148  {
2149  // writeOBJ(".", celli, loopPoints(celli), anchorPoints);
2150 
2152  << "Illegal loop " << loop
2153  << " when recreating cut-addressing"
2154  << " from existing cellLoops for cell " << celli
2155  << endl;
2156 
2157  cellLoops_[celli].setSize(0);
2158  cellAnchorPoints_[celli].setSize(0);
2159  }
2160  else
2161  {
2162  // Copy anchor points.
2163  cellAnchorPoints_[celli].transfer(anchorPoints);
2164 
2165  // Copy faceSplitCuts into overall faceSplit info.
2166  forAllConstIter(Map<edge>, faceSplitCuts, iter)
2167  {
2168  faceSplitCut_.insert(iter.key(), iter());
2169  }
2170 
2171  // Update edgeIsCut, pointIsCut information
2172  forAll(loop, cutI)
2173  {
2174  label cut = loop[cutI];
2175 
2176  if (isEdge(cut))
2177  {
2178  edgeIsCut_[getEdge(cut)] = true;
2179  }
2180  else
2181  {
2182  pointIsCut_[getVertex(cut)] = true;
2183  }
2184  }
2185  }
2186  }
2187  }
2188 
2189  // Reset edge weights
2190  forAll(edgeIsCut_, edgeI)
2191  {
2192  if (!edgeIsCut_[edgeI])
2193  {
2194  // Weight not used. Set to illegal value to make any use fall over.
2195  edgeWeight_[edgeI] = -great;
2196  }
2197  }
2198 }
2199 
2200 
2201 bool Foam::cellCuts::setFromCellLoop
2202 (
2203  const label celli,
2204  const labelList& loop,
2205  const scalarField& loopWeights
2206 )
2207 {
2208  // Update basic cut information from single cellLoop. Returns true if loop
2209  // was valid.
2210 
2211  // Dump loop for debugging.
2212  if (debug)
2213  {
2214  OFstream str("last_cell.obj");
2215 
2216  str<< "# edges of cell " << celli << nl;
2217 
2219  (
2220  str,
2221  mesh().cells(),
2222  mesh().faces(),
2223  mesh().points(),
2224  labelList(1, celli)
2225  );
2226 
2227 
2228  OFstream loopStr("last_loop.obj");
2229 
2230  loopStr<< "# looppoints for cell " << celli << nl;
2231 
2232  pointField pointsOfLoop = loopPoints(loop, loopWeights);
2233 
2234  forAll(pointsOfLoop, i)
2235  {
2236  meshTools::writeOBJ(loopStr, pointsOfLoop[i]);
2237  }
2238 
2239  str << 'l';
2240 
2241  forAll(pointsOfLoop, i)
2242  {
2243  str << ' ' << i + 1;
2244  }
2245  str << ' ' << 1 << nl;
2246  }
2247 
2248  bool okLoop = false;
2249 
2250  if (validEdgeLoop(loop, loopWeights))
2251  {
2252  // Storage for cuts across face
2253  Map<edge> faceSplitCuts(loop.size());
2254  // Storage for points on one side of cell
2255  labelList anchorPoints;
2256 
2257  okLoop =
2258  validLoop(celli, loop, loopWeights, faceSplitCuts, anchorPoints);
2259 
2260  if (okLoop)
2261  {
2262  // Valid loop. Copy cellLoops and anchorPoints
2263  cellLoops_[celli] = loop;
2264  cellAnchorPoints_[celli].transfer(anchorPoints);
2265 
2266  // Copy split cuts
2267  forAllConstIter(Map<edge>, faceSplitCuts, iter)
2268  {
2269  faceSplitCut_.insert(iter.key(), iter());
2270  }
2271 
2272 
2273  // Update edgeIsCut, pointIsCut information
2274  forAll(loop, cutI)
2275  {
2276  label cut = loop[cutI];
2277 
2278  if (isEdge(cut))
2279  {
2280  label edgeI = getEdge(cut);
2281 
2282  edgeIsCut_[edgeI] = true;
2283 
2284  edgeWeight_[edgeI] = loopWeights[cutI];
2285  }
2286  else
2287  {
2288  label vertI = getVertex(cut);
2289 
2290  pointIsCut_[vertI] = true;
2291  }
2292  }
2293  }
2294  }
2295 
2296  return okLoop;
2297 }
2298 
2299 
2300 void Foam::cellCuts::setFromCellLoops
2301 (
2302  const labelList& cellLabels,
2303  const labelListList& cellLoops,
2304  const List<scalarField>& cellLoopWeights
2305 )
2306 {
2307  // 'Uncut' edges/vertices that are not used in loops.
2308  pointIsCut_ = false;
2309 
2310  edgeIsCut_ = false;
2311 
2312  forAll(cellLabels, cellLabelI)
2313  {
2314  label celli = cellLabels[cellLabelI];
2315 
2316  const labelList& loop = cellLoops[cellLabelI];
2317 
2318  if (loop.size())
2319  {
2320  const scalarField& loopWeights = cellLoopWeights[cellLabelI];
2321 
2322  if (setFromCellLoop(celli, loop, loopWeights))
2323  {
2324  // Valid loop. Call above will have updated all already.
2325  }
2326  else
2327  {
2328  // Clear cellLoops
2329  cellLoops_[celli].setSize(0);
2330  }
2331  }
2332  }
2333 }
2334 
2335 
2336 void Foam::cellCuts::setFromCellCutter
2337 (
2338  const cellLooper& cellCutter,
2339  const List<refineCell>& refCells
2340 )
2341 {
2342  // 'Uncut' edges/vertices that are not used in loops.
2343  pointIsCut_ = false;
2344 
2345  edgeIsCut_ = false;
2346 
2347  // storage for loop of cuts (cut vertices and/or cut edges)
2348  labelList cellLoop;
2349  scalarField cellLoopWeights;
2350 
2351  // For debugging purposes
2352  DynamicList<label> invalidCutCells(2);
2353  DynamicList<labelList> invalidCutLoops(2);
2354  DynamicList<scalarField> invalidCutLoopWeights(2);
2355 
2356 
2357  forAll(refCells, refCelli)
2358  {
2359  const refineCell& refCell = refCells[refCelli];
2360 
2361  label celli = refCell.cellNo();
2362 
2363  const vector& refDir = refCell.direction();
2364 
2365 
2366  // Cut cell. Determines cellLoop and cellLoopWeights
2367  bool goodCut =
2368  cellCutter.cut
2369  (
2370  refDir,
2371  celli,
2372 
2373  pointIsCut_,
2374  edgeIsCut_,
2375  edgeWeight_,
2376 
2377  cellLoop,
2378  cellLoopWeights
2379  );
2380 
2381  // Check whether edge refinement is on a per face basis compatible with
2382  // current pattern.
2383  if (goodCut)
2384  {
2385  if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
2386  {
2387  // Valid loop. Will have updated all info already.
2388  }
2389  else
2390  {
2391  cellLoops_[celli].setSize(0);
2392 
2394  << "Found loop on cell " << celli
2395  << " that resulted in an unexpected bad cut." << nl
2396  << " Suggestions:" << nl
2397  << " - Turn on the debug switch for 'cellCuts' to get"
2398  << " geometry files that identify this cell." << nl
2399  << " - Also keep in mind to check the defined"
2400  << " reference directions, as these are most likely the"
2401  << " origin of the problem."
2402  << nl << endl;
2403 
2404  // Discarded by validLoop
2405  if (debug)
2406  {
2407  invalidCutCells.append(celli);
2408  invalidCutLoops.append(cellLoop);
2409  invalidCutLoopWeights.append(cellLoopWeights);
2410  }
2411  }
2412  }
2413  else
2414  {
2415  // Clear cellLoops
2416  cellLoops_[celli].setSize(0);
2417  }
2418  }
2419 
2420  if (debug && invalidCutCells.size())
2421  {
2422  invalidCutCells.shrink();
2423  invalidCutLoops.shrink();
2424  invalidCutLoopWeights.shrink();
2425 
2426  fileName cutsFile("invalidLoopCells.obj");
2427 
2428  Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
2429 
2430  OFstream cutsStream(cutsFile);
2431 
2433  (
2434  cutsStream,
2435  mesh().cells(),
2436  mesh().faces(),
2437  mesh().points(),
2438  invalidCutCells
2439  );
2440 
2441  fileName loopsFile("invalidLoops.obj");
2442 
2443  Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
2444 
2445  OFstream loopsStream(loopsFile);
2446 
2447  label vertI = 0;
2448 
2449  forAll(invalidCutLoops, i)
2450  {
2451  writeOBJ
2452  (
2453  loopsStream,
2454  loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2455  vertI
2456  );
2457  }
2458  }
2459 }
2460 
2461 
2462 void Foam::cellCuts::setFromCellCutter
2463 (
2464  const cellLooper& cellCutter,
2465  const labelList& cellLabels,
2466  const List<plane>& cellCutPlanes
2467 )
2468 {
2469  // 'Uncut' edges/vertices that are not used in loops.
2470  pointIsCut_ = false;
2471 
2472  edgeIsCut_ = false;
2473 
2474  // storage for loop of cuts (cut vertices and/or cut edges)
2475  labelList cellLoop;
2476  scalarField cellLoopWeights;
2477 
2478  // For debugging purposes
2479  DynamicList<label> invalidCutCells(2);
2480  DynamicList<labelList> invalidCutLoops(2);
2481  DynamicList<scalarField> invalidCutLoopWeights(2);
2482 
2483 
2484  forAll(cellLabels, i)
2485  {
2486  label celli = cellLabels[i];
2487 
2488  // Cut cell. Determines cellLoop and cellLoopWeights
2489  bool goodCut =
2490  cellCutter.cut
2491  (
2492  cellCutPlanes[i],
2493  celli,
2494 
2495  pointIsCut_,
2496  edgeIsCut_,
2497  edgeWeight_,
2498 
2499  cellLoop,
2500  cellLoopWeights
2501  );
2502 
2503  // Check whether edge refinement is on a per face basis compatible with
2504  // current pattern.
2505  if (goodCut)
2506  {
2507  if (setFromCellLoop(celli, cellLoop, cellLoopWeights))
2508  {
2509  // Valid loop. Will have updated all info already.
2510  }
2511  else
2512  {
2513  cellLoops_[celli].setSize(0);
2514 
2515  // Discarded by validLoop
2516  if (debug)
2517  {
2518  invalidCutCells.append(celli);
2519  invalidCutLoops.append(cellLoop);
2520  invalidCutLoopWeights.append(cellLoopWeights);
2521  }
2522  }
2523  }
2524  else
2525  {
2526  // Clear cellLoops
2527  cellLoops_[celli].setSize(0);
2528  }
2529  }
2530 
2531  if (debug && invalidCutCells.size())
2532  {
2533  invalidCutCells.shrink();
2534  invalidCutLoops.shrink();
2535  invalidCutLoopWeights.shrink();
2536 
2537  fileName cutsFile("invalidLoopCells.obj");
2538 
2539  Pout<< "cellCuts : writing inValidLoops cells to " << cutsFile << endl;
2540 
2541  OFstream cutsStream(cutsFile);
2542 
2544  (
2545  cutsStream,
2546  mesh().cells(),
2547  mesh().faces(),
2548  mesh().points(),
2549  invalidCutCells
2550  );
2551 
2552  fileName loopsFile("invalidLoops.obj");
2553 
2554  Pout<< "cellCuts : writing inValidLoops loops to " << loopsFile << endl;
2555 
2556  OFstream loopsStream(loopsFile);
2557 
2558  label vertI = 0;
2559 
2560  forAll(invalidCutLoops, i)
2561  {
2562  writeOBJ
2563  (
2564  loopsStream,
2565  loopPoints(invalidCutLoops[i], invalidCutLoopWeights[i]),
2566  vertI
2567  );
2568  }
2569  }
2570 }
2571 
2572 
2573 void Foam::cellCuts::orientPlanesAndLoops()
2574 {
2575  // Determine anchorPoints if not yet done by validLoop.
2576  forAll(cellLoops_, celli)
2577  {
2578  const labelList& loop = cellLoops_[celli];
2579 
2580  if (loop.size() && cellAnchorPoints_[celli].empty())
2581  {
2582  // Leave anchor points empty if illegal loop.
2583  calcAnchors
2584  (
2585  celli,
2586  loop,
2587  loopPoints(celli),
2588  cellAnchorPoints_[celli]
2589  );
2590  }
2591  }
2592 
2593  if (debug & 2)
2594  {
2595  Pout<< "cellAnchorPoints:" << endl;
2596  }
2597  forAll(cellAnchorPoints_, celli)
2598  {
2599  if (cellLoops_[celli].size())
2600  {
2601  if (cellAnchorPoints_[celli].empty())
2602  {
2604  << "No anchor points for cut cell " << celli << endl
2605  << "cellLoop:" << cellLoops_[celli] << abort(FatalError);
2606  }
2607 
2608  if (debug & 2)
2609  {
2610  Pout<< " cell:" << celli << " anchored at "
2611  << cellAnchorPoints_[celli] << endl;
2612  }
2613  }
2614  }
2615 
2616  // Calculate number of valid cellLoops
2617  nLoops_ = 0;
2618 
2619  forAll(cellLoops_, celli)
2620  {
2621  if (cellLoops_[celli].size())
2622  {
2623  nLoops_++;
2624  }
2625  }
2626 }
2627 
2628 
2629 void Foam::cellCuts::calcLoopsAndAddressing(const labelList& cutCells)
2630 {
2631  // Sanity check on weights
2632  forAll(edgeIsCut_, edgeI)
2633  {
2634  if (edgeIsCut_[edgeI])
2635  {
2636  scalar weight = edgeWeight_[edgeI];
2637 
2638  if (weight < 0 || weight > 1)
2639  {
2641  << "Weight out of range [0,1]. Edge " << edgeI
2642  << " verts:" << mesh().edges()[edgeI]
2643  << " weight:" << weight << abort(FatalError);
2644  }
2645  }
2646  else
2647  {
2648  // Weight not used. Set to illegal value to make any use fall over.
2649  edgeWeight_[edgeI] = -great;
2650  }
2651  }
2652 
2653 
2654  // Calculate faces that split cells in two
2655  calcCellLoops(cutCells);
2656 
2657  if (debug & 2)
2658  {
2659  Pout<< "-- cellLoops --" << endl;
2660  forAll(cellLoops_, celli)
2661  {
2662  const labelList& loop = cellLoops_[celli];
2663 
2664  if (loop.size())
2665  {
2666  Pout<< "cell:" << celli << " ";
2667  writeCuts(Pout, loop, loopWeights(loop));
2668  Pout<< endl;
2669  }
2670  }
2671  }
2672 
2673  // Redo basic cut information (pointIsCut, edgeIsCut, faceSplitCut)
2674  // using cellLoop only.
2675  setFromCellLoops();
2676 }
2677 
2678 
2679 void Foam::cellCuts::check() const
2680 {
2681  // Check weights for unsnapped values
2682  forAll(edgeIsCut_, edgeI)
2683  {
2684  if (edgeIsCut_[edgeI])
2685  {
2686  if
2687  (
2688  edgeWeight_[edgeI] < geomCellLooper::snapTol()
2689  || edgeWeight_[edgeI] > (1 - geomCellLooper::snapTol())
2690  )
2691  {
2692  // Should have been snapped.
2694  << "edge:" << edgeI << " vertices:"
2695  << mesh().edges()[edgeI]
2696  << " weight:" << edgeWeight_[edgeI] << " should have been"
2697  << " snapped to one of its endpoints"
2698  << endl;
2699  }
2700  }
2701  else
2702  {
2703  if (edgeWeight_[edgeI] > - 1)
2704  {
2706  << "edge:" << edgeI << " vertices:"
2707  << mesh().edges()[edgeI]
2708  << " weight:" << edgeWeight_[edgeI] << " is not cut but"
2709  << " its weight is not marked invalid"
2710  << abort(FatalError);
2711  }
2712  }
2713  }
2714 
2715  // Check that all elements of cellloop are registered
2716  forAll(cellLoops_, celli)
2717  {
2718  const labelList& loop = cellLoops_[celli];
2719 
2720  forAll(loop, i)
2721  {
2722  label cut = loop[i];
2723 
2724  if
2725  (
2726  (isEdge(cut) && !edgeIsCut_[getEdge(cut)])
2727  || (!isEdge(cut) && !pointIsCut_[getVertex(cut)])
2728  )
2729  {
2730  labelList cuts(1, cut);
2731  writeCuts(Pout, cuts, loopWeights(cuts));
2732 
2734  << "cell:" << celli << " loop:"
2735  << loop
2736  << " cut:" << cut << " is not marked as cut"
2737  << abort(FatalError);
2738  }
2739  }
2740  }
2741 
2742  // Check that no elements of cell loop are anchor point.
2743  forAll(cellLoops_, celli)
2744  {
2745  const labelList& anchors = cellAnchorPoints_[celli];
2746 
2747  const labelList& loop = cellLoops_[celli];
2748 
2749  if (loop.size() && anchors.empty())
2750  {
2752  << "cell:" << celli << " loop:" << loop
2753  << " has no anchor points"
2754  << abort(FatalError);
2755  }
2756 
2757 
2758  forAll(loop, i)
2759  {
2760  label cut = loop[i];
2761 
2762  if
2763  (
2764  !isEdge(cut)
2765  && findIndex(anchors, getVertex(cut)) != -1
2766  )
2767  {
2769  << "cell:" << celli << " loop:" << loop
2770  << " anchor points:" << anchors
2771  << " anchor:" << getVertex(cut) << " is part of loop"
2772  << abort(FatalError);
2773  }
2774  }
2775  }
2776 
2777 
2778  // Check that cut faces have a neighbour that is cut.
2779  boolList nbrCellIsCut;
2780  {
2781  boolList cellIsCut(mesh().nCells(), false);
2782  forAll(cellLoops_, celli)
2783  {
2784  cellIsCut[celli] = cellLoops_[celli].size();
2785  }
2786  syncTools::swapBoundaryCellList(mesh(), cellIsCut, nbrCellIsCut);
2787  }
2788 
2789  forAllConstIter(Map<edge>, faceSplitCut_, iter)
2790  {
2791  label facei = iter.key();
2792 
2793  if (mesh().isInternalFace(facei))
2794  {
2795  label own = mesh().faceOwner()[facei];
2796  label nei = mesh().faceNeighbour()[facei];
2797 
2798  if (cellLoops_[own].empty() && cellLoops_[nei].empty())
2799  {
2801  << "Internal face:" << facei << " cut by " << iter()
2802  << " has owner:" << own
2803  << " and neighbour:" << nei
2804  << " that are both uncut"
2805  << abort(FatalError);
2806  }
2807  }
2808  else
2809  {
2810  label bFacei = facei - mesh().nInternalFaces();
2811 
2812  label own = mesh().faceOwner()[facei];
2813 
2814  if (cellLoops_[own].empty() && !nbrCellIsCut[bFacei])
2815  {
2817  << "Boundary face:" << facei << " cut by " << iter()
2818  << " has owner:" << own
2819  << " that is uncut"
2820  << abort(FatalError);
2821  }
2822  }
2823  }
2824 }
2825 
2826 
2827 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
2828 
2830 (
2831  const polyMesh& mesh,
2832  const labelList& cutCells,
2833  const labelList& meshVerts,
2834  const labelList& meshEdges,
2835  const scalarField& meshEdgeWeights
2836 )
2837 :
2838  edgeVertex(mesh),
2839  pointIsCut_(expand(mesh.nPoints(), meshVerts)),
2840  edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
2841  edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
2842  faceSplitCut_(cutCells.size()),
2843  cellLoops_(mesh.nCells()),
2844  nLoops_(-1),
2845  cellAnchorPoints_(mesh.nCells())
2846 {
2847  if (debug)
2848  {
2849  Pout<< "cellCuts : constructor from cut verts and edges" << endl;
2850  }
2851 
2852  calcLoopsAndAddressing(cutCells);
2853 
2854  // Calculate planes and flip cellLoops if necessary
2855  orientPlanesAndLoops();
2856 
2857  if (debug)
2858  {
2859  check();
2860  }
2861 
2862  clearOut();
2863 
2864  if (debug)
2865  {
2866  Pout<< "cellCuts : leaving constructor from cut verts and edges"
2867  << endl;
2868  }
2869 }
2870 
2871 
2873 (
2874  const polyMesh& mesh,
2875  const labelList& meshVerts,
2876  const labelList& meshEdges,
2877  const scalarField& meshEdgeWeights
2878 )
2879 :
2880  edgeVertex(mesh),
2881  pointIsCut_(expand(mesh.nPoints(), meshVerts)),
2882  edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
2883  edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
2884  faceSplitCut_(mesh.nFaces()/10 + 1),
2885  cellLoops_(mesh.nCells()),
2886  nLoops_(-1),
2887  cellAnchorPoints_(mesh.nCells())
2888 {
2889  // Construct from pattern of cuts. Finds out itself which cells are cut.
2890  // (can go wrong if e.g. all neighbours of cell are refined)
2891 
2892  if (debug)
2893  {
2894  Pout<< "cellCuts : constructor from cellLoops" << endl;
2895  }
2896 
2897  calcLoopsAndAddressing(identityMap(mesh.nCells()));
2898 
2899  // Adds cuts on other side of coupled boundaries
2900  syncProc();
2901 
2902  // Calculate planes and flip cellLoops if necessary
2903  orientPlanesAndLoops();
2904 
2905  if (debug)
2906  {
2907  check();
2908  }
2909 
2910  clearOut();
2911 
2912  if (debug)
2913  {
2914  Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
2915  }
2916 }
2917 
2918 
2920 (
2921  const polyMesh& mesh,
2922  const labelList& cellLabels,
2923  const labelListList& cellLoops,
2924  const List<scalarField>& cellEdgeWeights
2925 )
2926 :
2927  edgeVertex(mesh),
2928  pointIsCut_(mesh.nPoints(), false),
2929  edgeIsCut_(mesh.nEdges(), false),
2930  edgeWeight_(mesh.nEdges(), -great),
2931  faceSplitCut_(cellLabels.size()),
2932  cellLoops_(mesh.nCells()),
2933  nLoops_(-1),
2934  cellAnchorPoints_(mesh.nCells())
2935 {
2936  if (debug)
2937  {
2938  Pout<< "cellCuts : constructor from cellLoops" << endl;
2939  }
2940 
2941  // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
2942  // Makes sure cuts are consistent
2943  setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
2944 
2945  // Adds cuts on other side of coupled boundaries
2946  syncProc();
2947 
2948  // Calculate planes and flip cellLoops if necessary
2949  orientPlanesAndLoops();
2950 
2951  if (debug)
2952  {
2953  check();
2954  }
2955 
2956  clearOut();
2957 
2958  if (debug)
2959  {
2960  Pout<< "cellCuts : leaving constructor from cellLoops" << endl;
2961  }
2962 }
2963 
2964 
2966 (
2967  const polyMesh& mesh,
2968  const cellLooper& cellCutter,
2969  const List<refineCell>& refCells
2970 )
2971 :
2972  edgeVertex(mesh),
2973  pointIsCut_(mesh.nPoints(), false),
2974  edgeIsCut_(mesh.nEdges(), false),
2975  edgeWeight_(mesh.nEdges(), -great),
2976  faceSplitCut_(refCells.size()),
2977  cellLoops_(mesh.nCells()),
2978  nLoops_(-1),
2979  cellAnchorPoints_(mesh.nCells())
2980 {
2981  if (debug)
2982  {
2983  Pout<< "cellCuts : constructor from cellCutter" << endl;
2984  }
2985 
2986  // Update pointIsCut, edgeIsCut, faceSplitCut from cell loops.
2987  // Makes sure cuts are consistent
2988  setFromCellCutter(cellCutter, refCells);
2989 
2990  // Adds cuts on other side of coupled boundaries
2991  syncProc();
2992 
2993  // Calculate planes and flip cellLoops if necessary
2994  orientPlanesAndLoops();
2995 
2996  if (debug)
2997  {
2998  check();
2999  }
3000 
3001  clearOut();
3002 
3003  if (debug)
3004  {
3005  Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
3006  }
3007 }
3008 
3009 
3011 (
3012  const polyMesh& mesh,
3013  const boolList& pointIsCut,
3014  const boolList& edgeIsCut,
3015  const scalarField& edgeWeight,
3016  const Map<edge>& faceSplitCut,
3017  const labelListList& cellLoops,
3018  const label nLoops,
3019  const labelListList& cellAnchorPoints
3020 )
3021 :
3022  edgeVertex(mesh),
3023  pointIsCut_(pointIsCut),
3024  edgeIsCut_(edgeIsCut),
3025  edgeWeight_(edgeWeight),
3026  faceSplitCut_(faceSplitCut),
3027  cellLoops_(cellLoops),
3028  nLoops_(nLoops),
3029  cellAnchorPoints_(cellAnchorPoints)
3030 {
3031  if (debug)
3032  {
3033  Pout<< "cellCuts : constructor from components" << endl;
3034  Pout<< "cellCuts : leaving constructor from components" << endl;
3035  }
3036 }
3037 
3038 
3039 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
3040 
3042 {
3043  clearOut();
3044 }
3045 
3046 
3048 {
3049  faceCutsPtr_.clear();
3050 }
3051 
3052 
3053 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
3054 
3055 Foam::pointField Foam::cellCuts::loopPoints(const label celli) const
3056 {
3057  const labelList& loop = cellLoops_[celli];
3058 
3059  pointField loopPts(loop.size());
3060 
3061  forAll(loop, fp)
3062  {
3063  label cut = loop[fp];
3064 
3065  if (isEdge(cut))
3066  {
3067  label edgeI = getEdge(cut);
3068 
3069  loopPts[fp] = coord(cut, edgeWeight_[edgeI]);
3070  }
3071  else
3072  {
3073  loopPts[fp] = coord(cut, -great);
3074  }
3075  }
3076  return loopPts;
3077 }
3078 
3079 
3080 void Foam::cellCuts::flip(const label celli)
3081 {
3082  labelList& loop = cellLoops_[celli];
3083 
3084  reverse(loop);
3085 
3086  // Reverse anchor point set.
3087  cellAnchorPoints_[celli] =
3088  nonAnchorPoints
3089  (
3090  mesh().cellPoints()[celli],
3091  cellAnchorPoints_[celli],
3092  loop
3093  );
3094 }
3095 
3096 
3098 {
3099  labelList& loop = cellLoops_[celli];
3100 
3101  reverse(loop);
3102 }
3103 
3104 
3105 void Foam::cellCuts::writeOBJ
3106 (
3107  Ostream& os,
3108  const pointField& loopPts,
3109  label& vertI
3110 ) const
3111 {
3112  label startVertI = vertI;
3113 
3114  forAll(loopPts, fp)
3115  {
3116  const point& pt = loopPts[fp];
3117 
3118  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
3119 
3120  vertI++;
3121  }
3122 
3123  os << 'f';
3124  forAll(loopPts, fp)
3125  {
3126  os << ' ' << startVertI + fp + 1;
3127  }
3128  os << endl;
3129 }
3130 
3131 
3132 void Foam::cellCuts::writeOBJ(Ostream& os) const
3133 {
3134  label vertI = 0;
3135 
3136  forAll(cellLoops_, celli)
3137  {
3138  writeOBJ(os, loopPoints(celli), vertI);
3139  }
3140 }
3141 
3142 
3143 void Foam::cellCuts::writeCellOBJ(const fileName& dir, const label celli) const
3144 {
3145  const labelList& anchors = cellAnchorPoints_[celli];
3146 
3147  writeOBJ(dir, celli, loopPoints(celli), anchors);
3148 }
3149 
3150 
3151 // ************************************************************************* //
Various functions to operate on Lists.
#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
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:197
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void setSize(const label)
Reset size of List.
Definition: List.C:281
A HashTable to objects of type <T> with a label key.
Definition: Map.H:52
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:85
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
const Cmpt & z() const
Definition: VectorI.H:87
const Cmpt & y() const
Definition: VectorI.H:81
const Cmpt & x() const
Definition: VectorI.H:75
Description of cuts across cells.
Definition: cellCuts.H:111
~cellCuts()
Destructor.
Definition: cellCuts.C:3041
cellCuts(const polyMesh &mesh, const labelList &cutCells, const labelList &meshVerts, const labelList &meshEdges, const scalarField &meshEdgeWeights)
Construct from cells to cut and pattern of cuts.
Definition: cellCuts.C:2830
void writeCellOBJ(const fileName &dir, const label celli) const
debugging:Write edges of cell and loop
Definition: cellCuts.C:3143
void flipLoopOnly(const label celli)
Flip loop for celli. Does not update anchors. Use with care.
Definition: cellCuts.C:3097
void flip(const label celli)
Flip loop for celli. Updates anchor points as well.
Definition: cellCuts.C:3080
const labelListList & cellLoops() const
For each cut cell the cut along the circumference.
Definition: cellCuts.H:563
labelList nonAnchorPoints(const labelList &cellPoints, const labelList &anchorPoints, const labelList &loop) const
Invert anchor point selection.
Definition: cellCuts.C:1321
void clearOut()
Clear out demand driven storage.
Definition: cellCuts.C:3047
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Definition: cellLooper.H:72
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
Definition: edgeVertex.H:53
const polyMesh & mesh() const
Definition: edgeVertex.H:93
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
A class for handling file names.
Definition: fileName.H:82
word name() const
Return file name (part beyond last /)
Definition: fileName.C:195
static scalar snapTol()
edge cmptType
Component type.
Definition: cellCuts.C:53
Traits class for primitives.
Definition: pTraits.H:53
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
label nCells() const
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh points.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronise values on all mesh edges.
static void syncBoundaryFaceList(const polyMesh &, UList< T > &, const CombineOp &cop, const TransformOp &top, const bool parRun=Pstream::parRun())
Synchronise values on boundary faces only.
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
Dummy transform to be used with syncTools.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label patchi
const pointField & points
label nPoints
const cellShapeList & cells
#define WarningInFunction
Report a warning using Foam::Warning.
List< labelPair > faceCuts(const face &f, const scalarField &pAlphas, const scalar isoAlpha)
Return the cuts for a given face. This returns a list of pairs of.
Definition: cutPoly.C:61
void getEdgeFaces(const primitiveMesh &, const label celli, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
Definition: meshTools.C:458
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
label walkFace(const primitiveMesh &, const label facei, const label startEdgeI, const label startVertI, const label nEdges)
Returns label of edge nEdges away from startEdge (in the direction.
Definition: meshTools.C:587
bool faceOnCell(const primitiveMesh &, const label celli, const label facei)
Is face used by cell.
Definition: meshTools.C:308
bool edgeOnCell(const primitiveMesh &, const label celli, const label edgeI)
Is edge used by cell.
Definition: meshTools.C:285
label otherFace(const primitiveMesh &, const label celli, const label facei, const label edgeI)
Return face on cell using edgeI but not facei. Throws error.
Definition: meshTools.C:536
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static const zero Zero
Definition: zero.H:97
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
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
static const labelSphericalTensor labelI(1)
Identity labelTensor.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:213
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
error FatalError
List< edge > edgeList
Definition: edgeList.H:38
labelList identityMap(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
List< face > faceList
Definition: faceListFwd.H:41
static const char nl
Definition: Ostream.H:266
static const label labelMin
Definition: label.H:61
string expand(const string &s, string::size_type &index, const dictionary &dict, const bool allowEnvVars, const bool allowEmpty)
Definition: stringOps.C:146
labelList f(nPoints)
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