splitCells.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 Application
25  splitCells
26 
27 Description
28  Utility to split cells with flat faces.
29 
30  Uses a geometric cut with a plane dividing the edge angle into two so
31  might produce funny cells. For hexes it will use by default a cut from
32  edge onto opposite edge (i.e. purely topological).
33 
34  Options:
35  - split cells from cellSet only
36  - use geometric cut for hexes as well
37 
38  The angle is the angle between two faces sharing an edge as seen from
39  inside each cell. So a cube will have all angles 90. If you want
40  to split cells with cell centre outside use e.g. angle 200
41 
42 \*---------------------------------------------------------------------------*/
43 
44 #include "argList.H"
45 #include "Time.H"
46 #include "polyTopoChange.H"
47 #include "polyTopoChanger.H"
48 #include "mapPolyMesh.H"
49 #include "polyMesh.H"
50 #include "cellCuts.H"
51 #include "cellSet.H"
52 #include "cellModeller.H"
53 #include "meshCutter.H"
54 #include "unitConversion.H"
55 #include "geomCellLooper.H"
56 #include "plane.H"
57 #include "edgeVertex.H"
58 #include "meshTools.H"
59 #include "ListOps.H"
60 
61 using namespace Foam;
62 
63 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
64 
65 
66 labelList pack(const boolList& lst)
67 {
68  labelList packedLst(lst.size());
69  label packedI = 0;
70 
71  forAll(lst, i)
72  {
73  if (lst[i])
74  {
75  packedLst[packedI++] = i;
76  }
77  }
78  packedLst.setSize(packedI);
79 
80  return packedLst;
81 }
82 
83 
84 scalarField pack(const boolList& lst, const scalarField& elems)
85 {
86  scalarField packedElems(lst.size());
87  label packedI = 0;
88 
89  forAll(lst, i)
90  {
91  if (lst[i])
92  {
93  packedElems[packedI++] = elems[i];
94  }
95  }
96  packedElems.setSize(packedI);
97 
98  return packedElems;
99 }
100 
101 
102 // Given sin and cos of max angle between normals calculate whether f0 and f1
103 // on celli make larger angle. Uses sinAngle only for quadrant detection.
104 bool largerAngle
105 (
106  const primitiveMesh& mesh,
107  const scalar cosAngle,
108  const scalar sinAngle,
109 
110  const label celli,
111  const label f0, // face label
112  const label f1,
113 
114  const vector& n0, // normal at f0
115  const vector& n1
116 )
117 {
118  const labelList& own = mesh.faceOwner();
119 
120  bool sameFaceOrder = !((own[f0] == celli) ^ (own[f1] == celli));
121 
122  // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
123  // gives -1.
124  scalar normalCosAngle = n0 & n1;
125 
126  if (sameFaceOrder)
127  {
128  normalCosAngle = -normalCosAngle;
129  }
130 
131 
132  // Get cos between faceCentre and normal vector to determine in
133  // which quadrant angle is. (Is correct for unwarped faces only!)
134  // Correct for non-outwards pointing normal.
135  vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
136  c1c0 /= mag(c1c0) + VSMALL;
137 
138  scalar fcCosAngle = n0 & c1c0;
139 
140  if (own[f0] != celli)
141  {
142  fcCosAngle = -fcCosAngle;
143  }
144 
145  if (sinAngle < 0.0)
146  {
147  // Looking for concave angles (quadrant 3 or 4)
148  if (fcCosAngle <= 0)
149  {
150  // Angle is convex so smaller.
151  return false;
152  }
153  else
154  {
155  if (normalCosAngle < cosAngle)
156  {
157  return false;
158  }
159  else
160  {
161  return true;
162  }
163  }
164  }
165  else
166  {
167  // Looking for convex angles (quadrant 1 or 2)
168  if (fcCosAngle > 0)
169  {
170  // Concave angle
171  return true;
172  }
173  else
174  {
175  // Convex. Check cos of normal vectors.
176  if (normalCosAngle > cosAngle)
177  {
178  return false;
179  }
180  else
181  {
182  return true;
183  }
184  }
185  }
186 }
187 
188 
189 // Split hex (and hex only) along edgeI creating two prisms
190 bool splitHex
191 (
192  const polyMesh& mesh,
193  const label celli,
194  const label edgeI,
195 
196  DynamicList<label>& cutCells,
197  DynamicList<labelList>& cellLoops,
198  DynamicList<scalarField>& cellEdgeWeights
199 )
200 {
201  // cut handling functions
202  edgeVertex ev(mesh);
203 
204  const edgeList& edges = mesh.edges();
205  const faceList& faces = mesh.faces();
206 
207  const edge& e = edges[edgeI];
208 
209  // Get faces on the side, i.e. faces not using edge but still using one of
210  // the edge endpoints.
211 
212  label leftI = -1;
213  label rightI = -1;
214  label leftFp = -1;
215  label rightFp = -1;
216 
217  const cell& cFaces = mesh.cells()[celli];
218 
219  forAll(cFaces, i)
220  {
221  label facei = cFaces[i];
222 
223  const face& f = faces[facei];
224 
225  label fp0 = findIndex(f, e[0]);
226  label fp1 = findIndex(f, e[1]);
227 
228  if (fp0 == -1)
229  {
230  if (fp1 != -1)
231  {
232  // Face uses e[1] but not e[0]
233  rightI = facei;
234  rightFp = fp1;
235 
236  if (leftI != -1)
237  {
238  // Have both faces so exit
239  break;
240  }
241  }
242  }
243  else
244  {
245  if (fp1 != -1)
246  {
247  // Face uses both e[1] and e[0]
248  }
249  else
250  {
251  leftI = facei;
252  leftFp = fp0;
253 
254  if (rightI != -1)
255  {
256  break;
257  }
258  }
259  }
260  }
261 
262  if (leftI == -1 || rightI == -1)
263  {
265  << " rightI:" << rightI << abort(FatalError);
266  }
267 
268  // Walk two vertices further on faces.
269 
270  const face& leftF = faces[leftI];
271 
272  label leftV = leftF[(leftFp + 2) % leftF.size()];
273 
274  const face& rightF = faces[rightI];
275 
276  label rightV = rightF[(rightFp + 2) % rightF.size()];
277 
278  labelList loop(4);
279  loop[0] = ev.vertToEVert(e[0]);
280  loop[1] = ev.vertToEVert(leftV);
281  loop[2] = ev.vertToEVert(rightV);
282  loop[3] = ev.vertToEVert(e[1]);
283 
284  scalarField loopWeights(4);
285  loopWeights[0] = -GREAT;
286  loopWeights[1] = -GREAT;
287  loopWeights[2] = -GREAT;
288  loopWeights[3] = -GREAT;
289 
290  cutCells.append(celli);
291  cellLoops.append(loop);
292  cellEdgeWeights.append(loopWeights);
293 
294  return true;
295 }
296 
297 
298 // Split celli along edgeI with a plane along halfNorm direction.
299 bool splitCell
300 (
301  const polyMesh& mesh,
302  const geomCellLooper& cellCutter,
303 
304  const label celli,
305  const label edgeI,
306  const vector& halfNorm,
307 
308  const boolList& vertIsCut,
309  const boolList& edgeIsCut,
310  const scalarField& edgeWeight,
311 
312  DynamicList<label>& cutCells,
313  DynamicList<labelList>& cellLoops,
314  DynamicList<scalarField>& cellEdgeWeights
315 )
316 {
317  const edge& e = mesh.edges()[edgeI];
318 
319  vector eVec = e.vec(mesh.points());
320  eVec /= mag(eVec);
321 
322  vector planeN = eVec ^ halfNorm;
323 
324  // Slightly tilt plane to make it not cut edges exactly
325  // halfway on fully regular meshes (since we want cuts
326  // to be snapped to vertices)
327  planeN += 0.01*halfNorm;
328 
329  planeN /= mag(planeN);
330 
331  // Define plane through edge
332  plane cutPlane(mesh.points()[e.start()], planeN);
333 
334  labelList loop;
335  scalarField loopWeights;
336 
337  if
338  (
339  cellCutter.cut
340  (
341  cutPlane,
342  celli,
343  vertIsCut,
344  edgeIsCut,
345  edgeWeight,
346  loop,
347  loopWeights
348  )
349  )
350  {
351  // Did manage to cut cell. Copy into overall list.
352  cutCells.append(celli);
353  cellLoops.append(loop);
354  cellEdgeWeights.append(loopWeights);
355 
356  return true;
357  }
358  else
359  {
360  return false;
361  }
362 }
363 
364 
365 // Collects cuts for all cells in cellSet
366 void collectCuts
367 (
368  const polyMesh& mesh,
369  const geomCellLooper& cellCutter,
370  const bool geometry,
371  const scalar minCos,
372  const scalar minSin,
373  const cellSet& cellsToCut,
374 
375  DynamicList<label>& cutCells,
376  DynamicList<labelList>& cellLoops,
377  DynamicList<scalarField>& cellEdgeWeights
378 )
379 {
380  // Get data from mesh
381  const cellShapeList& cellShapes = mesh.cellShapes();
382  const labelList& own = mesh.faceOwner();
383  const labelListList& cellEdges = mesh.cellEdges();
384  const vectorField& faceAreas = mesh.faceAreas();
385 
386  // Hex shape
387  const cellModel& hex = *(cellModeller::lookup("hex"));
388 
389  // cut handling functions
390  edgeVertex ev(mesh);
391 
392 
393  // Cut information per mesh entity
394  boolList vertIsCut(mesh.nPoints(), false);
395  boolList edgeIsCut(mesh.nEdges(), false);
396  scalarField edgeWeight(mesh.nEdges(), -GREAT);
397 
398  forAllConstIter(cellSet, cellsToCut, iter)
399  {
400  const label celli = iter.key();
401  const labelList& cEdges = cellEdges[celli];
402 
403  forAll(cEdges, i)
404  {
405  label edgeI = cEdges[i];
406 
407  label f0, f1;
408  meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1);
409 
410  vector n0 = faceAreas[f0];
411  n0 /= mag(n0);
412 
413  vector n1 = faceAreas[f1];
414  n1 /= mag(n1);
415 
416  if
417  (
418  largerAngle
419  (
420  mesh,
421  minCos,
422  minSin,
423 
424  celli,
425  f0,
426  f1,
427  n0,
428  n1
429  )
430  )
431  {
432  bool splitOk = false;
433 
434  if (!geometry && cellShapes[celli].model() == hex)
435  {
436  splitOk =
437  splitHex
438  (
439  mesh,
440  celli,
441  edgeI,
442 
443  cutCells,
444  cellLoops,
445  cellEdgeWeights
446  );
447  }
448  else
449  {
450  vector halfNorm;
451 
452  if ((own[f0] == celli) ^ (own[f1] == celli))
453  {
454  // Opposite owner orientation
455  halfNorm = 0.5*(n0 - n1);
456  }
457  else
458  {
459  // Faces have same owner or same neighbour so
460  // normals point in same direction
461  halfNorm = 0.5*(n0 + n1);
462  }
463 
464  splitOk =
465  splitCell
466  (
467  mesh,
468  cellCutter,
469  celli,
470  edgeI,
471  halfNorm,
472 
473  vertIsCut,
474  edgeIsCut,
475  edgeWeight,
476 
477  cutCells,
478  cellLoops,
479  cellEdgeWeights
480  );
481  }
482 
483 
484  if (splitOk)
485  {
486  // Update cell/edge/vertex wise info.
487  label index = cellLoops.size() - 1;
488  const labelList& loop = cellLoops[index];
489  const scalarField& loopWeights = cellEdgeWeights[index];
490 
491  forAll(loop, i)
492  {
493  label cut = loop[i];
494 
495  if (ev.isEdge(cut))
496  {
497  edgeIsCut[ev.getEdge(cut)] = true;
498  edgeWeight[ev.getEdge(cut)] = loopWeights[i];
499  }
500  else
501  {
502  vertIsCut[ev.getVertex(cut)] = true;
503  }
504  }
505 
506  // Stop checking edges for this cell.
507  break;
508  }
509  }
510  }
511  }
512 
513  cutCells.shrink();
514  cellLoops.shrink();
515  cellEdgeWeights.shrink();
516 }
517 
518 
519 
520 int main(int argc, char *argv[])
521 {
523  (
524  "split cells with flat faces"
525  );
526  #include "addOverwriteOption.H"
528  argList::validArgs.append("edgeAngle [0..360]");
529 
531  (
532  "set",
533  "name",
534  "split cells from specified cellSet only"
535  );
537  (
538  "geometry",
539  "use geometric cut for hexes as well"
540  );
542  (
543  "tol",
544  "scalar", "edge snap tolerance (default 0.2)"
545  );
546 
547  #include "setRootCase.H"
548  #include "createTime.H"
549  runTime.functionObjects().off();
550  #include "createPolyMesh.H"
551  const word oldInstance = mesh.pointsInstance();
552 
553  const scalar featureAngle = args.argRead<scalar>(1);
554  const scalar minCos = Foam::cos(degToRad(featureAngle));
555  const scalar minSin = Foam::sin(degToRad(featureAngle));
556 
557  const bool readSet = args.optionFound("set");
558  const bool geometry = args.optionFound("geometry");
559  const bool overwrite = args.optionFound("overwrite");
560 
561  const scalar edgeTol = args.optionLookupOrDefault("tol", 0.2);
562 
563  Info<< "Trying to split cells with internal angles > feature angle\n" << nl
564  << "featureAngle : " << featureAngle << nl
565  << "edge snapping tol : " << edgeTol << nl;
566  if (readSet)
567  {
568  Info<< "candidate cells : cellSet " << args["set"] << nl;
569  }
570  else
571  {
572  Info<< "candidate cells : all cells" << nl;
573  }
574  if (geometry)
575  {
576  Info<< "hex cuts : geometric; using edge tolerance" << nl;
577  }
578  else
579  {
580  Info<< "hex cuts : topological; cut to opposite edge" << nl;
581  }
582  Info<< endl;
583 
584 
585  // Cell circumference cutter
586  geomCellLooper cellCutter(mesh);
587  // Snap all edge cuts close to endpoints to vertices.
589 
590  // Candidate cells to cut
591  cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
592 
593  if (readSet)
594  {
595  // Read cells to cut from cellSet
596  cellSet cells(mesh, args["set"]);
597 
598  cellsToCut = cells;
599  }
600 
601  while (true)
602  {
603  if (!readSet)
604  {
605  // Try all cells for cutting
606  for (label celli = 0; celli < mesh.nCells(); celli++)
607  {
608  cellsToCut.insert(celli);
609  }
610  }
611 
612 
613  // Cut information per cut cell
614  DynamicList<label> cutCells(mesh.nCells()/10 + 10);
615  DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
616  DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
617 
618  collectCuts
619  (
620  mesh,
621  cellCutter,
622  geometry,
623  minCos,
624  minSin,
625  cellsToCut,
626 
627  cutCells,
628  cellLoops,
629  cellEdgeWeights
630  );
631 
632  cellSet cutSet(mesh, "cutSet", cutCells.size());
633  forAll(cutCells, i)
634  {
635  cutSet.insert(cutCells[i]);
636  }
637 
638  // Gets cuts across cells from cuts through edges.
639  Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
640  << cutSet.instance()/cutSet.local()/cutSet.name()
641  << endl << endl;
642  cutSet.write();
643 
644  // Analyze cuts for clashes.
645  cellCuts cuts
646  (
647  mesh,
648  cutCells, // cells candidate for cutting
649  cellLoops,
650  cellEdgeWeights
651  );
652 
653  Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
654 
655  if (cuts.nLoops() == 0)
656  {
657  break;
658  }
659 
660  // Remove cut cells from cellsToCut (Note:only relevant if -readSet)
661  forAll(cuts.cellLoops(), celli)
662  {
663  if (cuts.cellLoops()[celli].size())
664  {
665  //Info<< "Removing cut cell " << celli << " from wishlist"
666  // << endl;
667  cellsToCut.erase(celli);
668  }
669  }
670 
671  // At least some cells are cut.
672  polyTopoChange meshMod(mesh);
673 
674  // Cutting engine
675  meshCutter cutter(mesh);
676 
677  // Insert mesh refinement into polyTopoChange.
678  cutter.setRefinement(cuts, meshMod);
679 
680  // Do all changes
681  Info<< "Morphing ..." << endl;
682 
683  if (!overwrite)
684  {
685  runTime++;
686  }
687 
688  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
689 
690  if (morphMap().hasMotionPoints())
691  {
692  mesh.movePoints(morphMap().preMotionPoints());
693  }
694 
695  // Update stored labels on meshCutter
696  cutter.updateMesh(morphMap());
697 
698  // Update cellSet
699  cellsToCut.updateMesh(morphMap());
700 
701  Info<< "Remaining:" << cellsToCut.size() << endl;
702 
703  // Write resulting mesh
704  if (overwrite)
705  {
706  mesh.setInstance(oldInstance);
707  }
708 
709  Info<< "Writing refined morphMesh to time " << runTime.timeName()
710  << endl;
711 
712  mesh.write();
713  }
714 
715  Info<< "End\n" << endl;
716 
717  return 0;
718 }
719 
720 
721 // ************************************************************************* //
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: polyMesh.C:1048
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
virtual bool cut(const vector &refDir, const label celli, const boolList &vertIsCut, const boolList &edgeIsCut, const scalarField &edgeWeight, labelList &loop, scalarField &loopWeights) const
Create cut along circumference of celli. Gets current mesh cuts.
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 void setSnapTol(const scalar tol)
static const cellModel * lookup(const word &)
Look up a model by name and return a pointer to the model or NULL.
Definition: cellModeller.C:88
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
const vectorField & faceAreas() const
Unit conversion functions.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual void updateMesh(const mapPolyMesh &morphMap)
Update any stored data for new labels.
Definition: cellSet.C:225
static void noParallel()
Remove the parallel options.
Definition: argList.C:146
static SLList< string > validArgs
A list of valid (mandatory) arguments.
Definition: argList.H:154
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const vectorField & faceCentres() const
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
Description of cuts across cells.
Definition: cellCuts.H:108
const cellList & cells() const
scalar f1
Definition: createFields.H:28
T optionLookupOrDefault(const word &opt, const T &deflt) const
Read a value from the named option if present.
Definition: argListI.H:237
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
Various functions to operate on Lists.
bool erase(const iterator &)
Erase a hashedEntry specified by given iterator.
Definition: HashTable.C:367
const cellShapeList & cellShapes() const
Return cell shapes.
dimensionedScalar cos(const dimensionedScalar &ds)
const cellShapeList & cells
bool optionFound(const word &opt) const
Return true if the named option is found.
Definition: argListI.H:108
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Description of cell after splitting. Contains cellLabel and pointers to cells it it split in...
Definition: splitCell.H:51
A class for handling words, derived from string.
Definition: word.H:59
Cuts (splits) cells.
Definition: meshCutter.H:134
static void addOption(const word &opt, const string &param="", const string &usage="")
Add to an option to validOptions with usage information.
Definition: argList.C:93
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:292
Combines edge or vertex in single label. Used to specify cuts across cell circumference.
Definition: edgeVertex.H:52
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensionedScalar sin(const dimensionedScalar &ds)
label start() const
Return start vertex label.
Definition: edgeI.H:81
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:240
static const char nl
Definition: Ostream.H:262
const labelListList & cellEdges() const
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
label nEdges() const
T argRead(const label index) const
Read a value from the argument at index.
Definition: argListI.H:177
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
void setSize(const label)
Reset size of List.
Definition: List.C:295
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:772
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
A collection of cell labels.
Definition: cellSet.H:48
Direct mesh changes based on v1.3 polyTopoChange syntax.
Maps a geometry to a set of cell primitives, which enables geometric cell data to be calculated witho...
Definition: cellModel.H:64
vector vec(const pointField &) const
Return the vector (end - start)
Definition: edgeI.H:157
virtual bool write() const
Write using setting from DB.
label nPoints() const
messageStream Info
virtual const labelList & faceOwner() const =0
Face face-owner addresing.
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
virtual Ostream & write(const token &)=0
Write next token to stream.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void addBoolOption(const word &opt, const string &usage="")
Add to a bool option to validOptions with usage information.
Definition: argList.C:83
static void addNote(const string &)
Add extra notes for the usage information.
Definition: argList.C:124
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
Foam::argList args(argc, argv)
void getEdgeFaces(const primitiveMesh &, const label celli, const label edgeI, label &face0, label &face1)
Get faces on cell using edgeI. Throws error if no two found.
Definition: meshTools.C:458
Namespace for OpenFOAM.
Implementation of cellLooper. Does pure geometric cut through cell.