FreeFOAM The Cross-Platform CFD Toolkit
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) 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 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 Usage
43 
44  - splitCells [OPTIONS] <edge angle [0..360]>
45 
46  @param <edge angle [0..360]> \n
47  @todo Detailed description of argument.
48 
49  @param -geometry \n
50  Geometry based splitting for hex cells.
51 
52  @param -tol <tolerance>\n
53  Edge snap tolerance.
54 
55  @param -set <cellSet name>\n
56  Cell set to split.
57 
58  @param -overwrite \n
59  Overwrite existing data.
60 
61  @param -case <dir>\n
62  Case directory.
63 
64  @param -help \n
65  Display help message.
66 
67  @param -doc \n
68  Display Doxygen API documentation page for this application.
69 
70  @param -srcDoc \n
71  Display Doxygen source documentation page for this application.
72 
73 \*---------------------------------------------------------------------------*/
74 
75 #include <OpenFOAM/argList.H>
76 #include <OpenFOAM/Time.H>
79 #include <OpenFOAM/mapPolyMesh.H>
80 #include <OpenFOAM/polyMesh.H>
81 #include <dynamicMesh/cellCuts.H>
82 #include <meshTools/cellSet.H>
83 #include <OpenFOAM/cellModeller.H>
84 #include <dynamicMesh/meshCutter.H>
87 #include <OpenFOAM/plane.H>
88 #include <dynamicMesh/edgeVertex.H>
89 #include <meshTools/meshTools.H>
90 #include <OpenFOAM/ListOps.H>
91 
92 using namespace Foam;
93 
94 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
95 
96 
97 labelList pack(const boolList& lst)
98 {
99  labelList packedLst(lst.size());
100  label packedI = 0;
101 
102  forAll(lst, i)
103  {
104  if (lst[i])
105  {
106  packedLst[packedI++] = i;
107  }
108  }
109  packedLst.setSize(packedI);
110 
111  return packedLst;
112 }
113 
114 
115 scalarField pack(const boolList& lst, const scalarField& elems)
116 {
117  scalarField packedElems(lst.size());
118  label packedI = 0;
119 
120  forAll(lst, i)
121  {
122  if (lst[i])
123  {
124  packedElems[packedI++] = elems[i];
125  }
126  }
127  packedElems.setSize(packedI);
128 
129  return packedElems;
130 }
131 
132 
133 // Given sin and cos of max angle between normals calculate whether f0 and f1
134 // on cellI make larger angle. Uses sinAngle only for quadrant detection.
135 bool largerAngle
136 (
137  const primitiveMesh& mesh,
138  const scalar cosAngle,
139  const scalar sinAngle,
140 
141  const label cellI,
142  const label f0, // face label
143  const label f1,
144 
145  const vector& n0, // normal at f0
146  const vector& n1
147 )
148 {
149  const labelList& own = mesh.faceOwner();
150 
151  bool sameFaceOrder = !((own[f0] == cellI) ^ (own[f1] == cellI));
152 
153  // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
154  // gives -1.
155  scalar normalCosAngle = n0 & n1;
156 
157  if (sameFaceOrder)
158  {
159  normalCosAngle = -normalCosAngle;
160  }
161 
162 
163  // Get cos between faceCentre and normal vector to determine in
164  // which quadrant angle is. (Is correct for unwarped faces only!)
165  // Correct for non-outwards pointing normal.
166  vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
167  c1c0 /= mag(c1c0) + VSMALL;
168 
169  scalar fcCosAngle = n0 & c1c0;
170 
171  if (own[f0] != cellI)
172  {
173  fcCosAngle = -fcCosAngle;
174  }
175 
176  if (sinAngle < 0.0)
177  {
178  // Looking for concave angles (quadrant 3 or 4)
179  if (fcCosAngle <= 0)
180  {
181  // Angle is convex so smaller.
182  return false;
183  }
184  else
185  {
186  if (normalCosAngle < cosAngle)
187  {
188  return false;
189  }
190  else
191  {
192  return true;
193  }
194  }
195  }
196  else
197  {
198  // Looking for convex angles (quadrant 1 or 2)
199  if (fcCosAngle > 0)
200  {
201  // Concave angle
202  return true;
203  }
204  else
205  {
206  // Convex. Check cos of normal vectors.
207  if (normalCosAngle > cosAngle)
208  {
209  return false;
210  }
211  else
212  {
213  return true;
214  }
215  }
216  }
217 }
218 
219 
220 // Split hex (and hex only) along edgeI creating two prisms
221 bool splitHex
222 (
223  const polyMesh& mesh,
224  const label cellI,
225  const label edgeI,
226 
227  DynamicList<label>& cutCells,
228  DynamicList<labelList>& cellLoops,
229  DynamicList<scalarField>& cellEdgeWeights
230 )
231 {
232  // cut handling functions
233  edgeVertex ev(mesh);
234 
235  const edgeList& edges = mesh.edges();
236  const faceList& faces = mesh.faces();
237 
238  const edge& e = edges[edgeI];
239 
240  // Get faces on the side, i.e. faces not using edge but still using one of
241  // the edge endpoints.
242 
243  label leftI = -1;
244  label rightI = -1;
245  label leftFp = -1;
246  label rightFp = -1;
247 
248  const cell& cFaces = mesh.cells()[cellI];
249 
250  forAll(cFaces, i)
251  {
252  label faceI = cFaces[i];
253 
254  const face& f = faces[faceI];
255 
256  label fp0 = findIndex(f, e[0]);
257  label fp1 = findIndex(f, e[1]);
258 
259  if (fp0 == -1)
260  {
261  if (fp1 != -1)
262  {
263  // Face uses e[1] but not e[0]
264  rightI = faceI;
265  rightFp = fp1;
266 
267  if (leftI != -1)
268  {
269  // Have both faces so exit
270  break;
271  }
272  }
273  }
274  else
275  {
276  if (fp1 != -1)
277  {
278  // Face uses both e[1] and e[0]
279  }
280  else
281  {
282  leftI = faceI;
283  leftFp = fp0;
284 
285  if (rightI != -1)
286  {
287  break;
288  }
289  }
290  }
291  }
292 
293  if (leftI == -1 || rightI == -1)
294  {
295  FatalErrorIn("splitHex") << "Problem : leftI:" << leftI
296  << " rightI:" << rightI << abort(FatalError);
297  }
298 
299  // Walk two vertices further on faces.
300 
301  const face& leftF = faces[leftI];
302 
303  label leftV = leftF[(leftFp + 2) % leftF.size()];
304 
305  const face& rightF = faces[rightI];
306 
307  label rightV = rightF[(rightFp + 2) % rightF.size()];
308 
309  labelList loop(4);
310  loop[0] = ev.vertToEVert(e[0]);
311  loop[1] = ev.vertToEVert(leftV);
312  loop[2] = ev.vertToEVert(rightV);
313  loop[3] = ev.vertToEVert(e[1]);
314 
315  scalarField loopWeights(4);
316  loopWeights[0] = -GREAT;
317  loopWeights[1] = -GREAT;
318  loopWeights[2] = -GREAT;
319  loopWeights[3] = -GREAT;
320 
321  cutCells.append(cellI);
322  cellLoops.append(loop);
323  cellEdgeWeights.append(loopWeights);
324 
325  return true;
326 }
327 
328 
329 // Split cellI along edgeI with a plane along halfNorm direction.
330 bool splitCell
331 (
332  const polyMesh& mesh,
333  const geomCellLooper& cellCutter,
334 
335  const label cellI,
336  const label edgeI,
337  const vector& halfNorm,
338 
339  const boolList& vertIsCut,
340  const boolList& edgeIsCut,
341  const scalarField& edgeWeight,
342 
343  DynamicList<label>& cutCells,
344  DynamicList<labelList>& cellLoops,
345  DynamicList<scalarField>& cellEdgeWeights
346 )
347 {
348  const edge& e = mesh.edges()[edgeI];
349 
350  vector eVec = e.vec(mesh.points());
351  eVec /= mag(eVec);
352 
353  vector planeN = eVec ^ halfNorm;
354 
355  // Slightly tilt plane to make it not cut edges exactly
356  // halfway on fully regular meshes (since we want cuts
357  // to be snapped to vertices)
358  planeN += 0.01*halfNorm;
359 
360  planeN /= mag(planeN);
361 
362  // Define plane through edge
363  plane cutPlane(mesh.points()[e.start()], planeN);
364 
365  labelList loop;
366  scalarField loopWeights;
367 
368  if
369  (
370  cellCutter.cut
371  (
372  cutPlane,
373  cellI,
374  vertIsCut,
375  edgeIsCut,
376  edgeWeight,
377  loop,
378  loopWeights
379  )
380  )
381  {
382  // Did manage to cut cell. Copy into overall list.
383  cutCells.append(cellI);
384  cellLoops.append(loop);
385  cellEdgeWeights.append(loopWeights);
386 
387  return true;
388  }
389  else
390  {
391  return false;
392  }
393 }
394 
395 
396 // Collects cuts for all cells in cellSet
397 void collectCuts
398 (
399  const polyMesh& mesh,
400  const geomCellLooper& cellCutter,
401  const bool geometry,
402  const scalar minCos,
403  const scalar minSin,
404  const cellSet& cellsToCut,
405 
406  DynamicList<label>& cutCells,
407  DynamicList<labelList>& cellLoops,
408  DynamicList<scalarField>& cellEdgeWeights
409 )
410 {
411  // Get data from mesh
412  const cellShapeList& cellShapes = mesh.cellShapes();
413  const labelList& own = mesh.faceOwner();
414  const labelListList& cellEdges = mesh.cellEdges();
415  const vectorField& faceAreas = mesh.faceAreas();
416 
417  // Hex shape
418  const cellModel& hex = *(cellModeller::lookup("hex"));
419 
420  // cut handling functions
421  edgeVertex ev(mesh);
422 
423 
424  // Cut information per mesh entity
425  boolList vertIsCut(mesh.nPoints(), false);
426  boolList edgeIsCut(mesh.nEdges(), false);
427  scalarField edgeWeight(mesh.nEdges(), -GREAT);
428 
429  for
430  (
431  cellSet::const_iterator iter = cellsToCut.begin();
432  iter != cellsToCut.end();
433  ++iter
434  )
435  {
436  label cellI = iter.key();
437 
438  const labelList& cEdges = cellEdges[cellI];
439 
440  forAll(cEdges, i)
441  {
442  label edgeI = cEdges[i];
443 
444  label f0, f1;
445  meshTools::getEdgeFaces(mesh, cellI, edgeI, f0, f1);
446 
447  vector n0 = faceAreas[f0];
448  n0 /= mag(n0);
449 
450  vector n1 = faceAreas[f1];
451  n1 /= mag(n1);
452 
453  if
454  (
455  largerAngle
456  (
457  mesh,
458  minCos,
459  minSin,
460 
461  cellI,
462  f0,
463  f1,
464  n0,
465  n1
466  )
467  )
468  {
469  bool splitOk = false;
470 
471  if (!geometry && cellShapes[cellI].model() == hex)
472  {
473  splitOk =
474  splitHex
475  (
476  mesh,
477  cellI,
478  edgeI,
479 
480  cutCells,
481  cellLoops,
482  cellEdgeWeights
483  );
484  }
485  else
486  {
487  vector halfNorm;
488 
489  if ((own[f0] == cellI) ^ (own[f1] == cellI))
490  {
491  // Opposite owner orientation
492  halfNorm = 0.5*(n0 - n1);
493  }
494  else
495  {
496  // Faces have same owner or same neighbour so
497  // normals point in same direction
498  halfNorm = 0.5*(n0 + n1);
499  }
500 
501  splitOk =
502  splitCell
503  (
504  mesh,
505  cellCutter,
506  cellI,
507  edgeI,
508  halfNorm,
509 
510  vertIsCut,
511  edgeIsCut,
512  edgeWeight,
513 
514  cutCells,
515  cellLoops,
516  cellEdgeWeights
517  );
518  }
519 
520 
521  if (splitOk)
522  {
523  // Update cell/edge/vertex wise info.
524  label index = cellLoops.size() - 1;
525  const labelList& loop = cellLoops[index];
526  const scalarField& loopWeights = cellEdgeWeights[index];
527 
528  forAll(loop, i)
529  {
530  label cut = loop[i];
531 
532  if (ev.isEdge(cut))
533  {
534  edgeIsCut[ev.getEdge(cut)] = true;
535  edgeWeight[ev.getEdge(cut)] = loopWeights[i];
536  }
537  else
538  {
539  vertIsCut[ev.getVertex(cut)] = true;
540  }
541  }
542 
543  // Stop checking edges for this cell.
544  break;
545  }
546  }
547  }
548  }
549 
550  cutCells.shrink();
551  cellLoops.shrink();
552  cellEdgeWeights.shrink();
553 }
554 
555 
556 // Main program:
557 
558 int main(int argc, char *argv[])
559 {
561  argList::validOptions.insert("set", "cellSet name");
562  argList::validOptions.insert("geometry", "");
563  argList::validOptions.insert("tol", "edge snap tolerance");
564  argList::validOptions.insert("overwrite", "");
565  argList::validArgs.append("edge angle [0..360]");
566 
567 # include <OpenFOAM/setRootCase.H>
568 # include <OpenFOAM/createTime.H>
569  runTime.functionObjects().off();
570 # include <OpenFOAM/createPolyMesh.H>
571  const word oldInstance = mesh.pointsInstance();
572 
573  scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
574 
575  scalar radAngle = featureAngle * mathematicalConstant::pi/180.0;
576  scalar minCos = Foam::cos(radAngle);
577  scalar minSin = Foam::sin(radAngle);
578 
579  bool readSet = args.optionFound("set");
580  bool geometry = args.optionFound("geometry");
581  bool overwrite = args.optionFound("overwrite");
582 
583  scalar edgeTol = 0.2;
584  args.optionReadIfPresent("tol", edgeTol);
585 
586  Info<< "Trying to split cells with internal angles > feature angle\n" << nl
587  << "featureAngle : " << featureAngle << nl
588  << "edge snapping tol : " << edgeTol << nl;
589  if (readSet)
590  {
591  Info<< "candidate cells : cellSet " << args.option("set") << nl;
592  }
593  else
594  {
595  Info<< "candidate cells : all cells" << nl;
596  }
597  if (geometry)
598  {
599  Info<< "hex cuts : geometric; using edge tolerance" << nl;
600  }
601  else
602  {
603  Info<< "hex cuts : topological; cut to opposite edge" << nl;
604  }
605  Info<< endl;
606 
607 
608  // Cell circumference cutter
609  geomCellLooper cellCutter(mesh);
610  // Snap all edge cuts close to endpoints to vertices.
612 
613  // Candidate cells to cut
614  cellSet cellsToCut(mesh, "cellsToCut", mesh.nCells()/100);
615 
616  if (readSet)
617  {
618  // Read cells to cut from cellSet
619  cellSet cells(mesh, args.option("set"));
620 
621  cellsToCut = cells;
622  }
623 
624  while (true)
625  {
626  if (!readSet)
627  {
628  // Try all cells for cutting
629  for (label cellI = 0; cellI < mesh.nCells(); cellI++)
630  {
631  cellsToCut.insert(cellI);
632  }
633  }
634 
635 
636  // Cut information per cut cell
637  DynamicList<label> cutCells(mesh.nCells()/10 + 10);
638  DynamicList<labelList> cellLoops(mesh.nCells()/10 + 10);
639  DynamicList<scalarField> cellEdgeWeights(mesh.nCells()/10 + 10);
640 
641  collectCuts
642  (
643  mesh,
644  cellCutter,
645  geometry,
646  minCos,
647  minSin,
648  cellsToCut,
649 
650  cutCells,
651  cellLoops,
652  cellEdgeWeights
653  );
654 
655  cellSet cutSet(mesh, "cutSet", cutCells.size());
656  forAll(cutCells, i)
657  {
658  cutSet.insert(cutCells[i]);
659  }
660 
661  // Gets cuts across cells from cuts through edges.
662  Info<< "Writing " << cutSet.size() << " cells to cut to cellSet "
663  << cutSet.instance()/cutSet.local()/cutSet.name()
664  << endl << endl;
665  cutSet.write();
666 
667  // Analyze cuts for clashes.
668  cellCuts cuts
669  (
670  mesh,
671  cutCells, // cells candidate for cutting
672  cellLoops,
673  cellEdgeWeights
674  );
675 
676  Info<< "Actually cut cells:" << cuts.nLoops() << nl << endl;
677 
678  if (cuts.nLoops() == 0)
679  {
680  break;
681  }
682 
683  // Remove cut cells from cellsToCut (Note:only relevant if -readSet)
684  forAll(cuts.cellLoops(), cellI)
685  {
686  if (cuts.cellLoops()[cellI].size())
687  {
688  //Info<< "Removing cut cell " << cellI << " from wishlist"
689  // << endl;
690  cellsToCut.erase(cellI);
691  }
692  }
693 
694  // At least some cells are cut.
695  polyTopoChange meshMod(mesh);
696 
697  // Cutting engine
698  meshCutter cutter(mesh);
699 
700  // Insert mesh refinement into polyTopoChange.
701  cutter.setRefinement(cuts, meshMod);
702 
703  // Do all changes
704  Info<< "Morphing ..." << endl;
705 
706  if (!overwrite)
707  {
708  runTime++;
709  }
710 
711  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);
712 
713  if (morphMap().hasMotionPoints())
714  {
715  mesh.movePoints(morphMap().preMotionPoints());
716  }
717 
718  // Update stored labels on meshCutter
719  cutter.updateMesh(morphMap());
720 
721  // Update cellSet
722  cellsToCut.updateMesh(morphMap());
723 
724  Info<< "Remaining:" << cellsToCut.size() << endl;
725 
726  // Write resulting mesh
727  if (overwrite)
728  {
729  mesh.setInstance(oldInstance);
730  }
731 
732  Info<< "Writing refined morphMesh to time " << runTime.timeName()
733  << endl;
734 
735  mesh.write();
736  }
737 
738  Info << "End\n" << endl;
739 
740  return 0;
741 }
742 
743 
744 // ************************ vim: set sw=4 sts=4 et: ************************ //