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