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