FreeFOAM The Cross-Platform CFD Toolkit
cellCuts.H
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 Class
25  Foam::cellCuts
26 
27 Description
28  Description of cuts across cells.
29 
30  Description of cut is given as list of vertices and
31  list of edges to be cut (and position on edge).
32 
33  Does some checking of correctness/non-overlapping of cuts. 2x2x2
34  refinement has to be done in three passes since cuts can not overlap
35  (would make addressing too complicated)
36 
37  Introduces concept of 'cut' which is either an existing vertex
38  or a edge.
39 
40  Input can either be
41  -# list of cut vertices and list of cut edges. Constructs cell
42  circumference walks ('cellLoops').
43 
44  -# list of cell circumference walks. Will filter them so they
45  don't overlap.
46 
47  -# cellWalker and list of cells to refine (refineCell). Constructs
48  cellLoops and does B. cellWalker is class which can cut a single
49  cell using a plane through the cell centre and in a certain normal
50  direction
51 
52  CellCuts constructed from cellLoops (B, C) can have multiple cut-edges
53  and/or cut-point as long as there is per face only one (or none) cut across
54  a face, i.e. a face can only be split into two faces.
55 
56  The information available after construction:
57  - pointIsCut, edgeIsCut.
58  - faceSplitCut : the cross-cut of a face.
59  - cellLoops : the loop which cuts across a cell
60  - cellAnchorPoints : per cell the vertices on one side which are
61  considered the anchor points.
62 
63  AnchorPoints: connected loops have to be oriented in the same way to
64  be able to grow new internal faces out of the same bottom faces.
65  (limitation of the mapping procedure). The loop is cellLoop is oriented
66  such that the normal of it points towards the anchorPoints.
67  This is the only routine which uses geometric information.
68 
69 
70  Limitations:
71  - cut description is very precise. Hard to get right.
72  - do anchor points need to be unique per cell? Very inefficient.
73  - is orientation guaranteed to be correct? Uses geometric info so can go
74  wrong on highly distorted meshes.
75  - is memory inefficient. Probably need to use Maps instead of
76  labelLists.
77 
78 SourceFiles
79  cellCuts.C
80 
81 \*---------------------------------------------------------------------------*/
82 
83 #ifndef cellCuts_H
84 #define cellCuts_H
85 
86 #include <dynamicMesh/edgeVertex.H>
87 #include <OpenFOAM/labelList.H>
88 #include <OpenFOAM/boolList.H>
89 #include <OpenFOAM/scalarField.H>
90 #include <OpenFOAM/pointField.H>
91 #include <OpenFOAM/DynamicList.H>
92 #include <OpenFOAM/typeInfo.H>
93 
94 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
95 
96 namespace Foam
97 {
98 
99 // Forward declaration of classes
100 class polyMesh;
101 class cellLooper;
102 class refineCell;
103 class plane;
104 
105 /*---------------------------------------------------------------------------*\
106  Class cellCuts Declaration
107 \*---------------------------------------------------------------------------*/
108 
109 class cellCuts
110 :
111  public edgeVertex
112 {
113  // Private data
114 
115  // Per point/edge status
116 
117  //- Is mesh point cut
118  boolList pointIsCut_;
119 
120  //- Is edge cut
121  boolList edgeIsCut_;
122 
123  //- If edge is cut gives weight (0->start() to 1->end())
124  scalarField edgeWeight_;
125 
126 
127  // Cut addressing
128 
129  //- Cuts per existing face (includes those along edge of face)
130  // Cuts in no particular order.
131  mutable labelListList* faceCutsPtr_;
132 
133  //- Per face : cut across edge (so not along existing edge)
134  // (can only be one per face)
135  Map<edge> faceSplitCut_;
136 
137 
138  // Cell-cut addressing
139 
140  //- Loop across cell circumference
141  labelListList cellLoops_;
142 
143  //- Number of valid loops in cellLoops_
144  label nLoops_;
145 
146  //- For each cut cell the points on the 'anchor' side of the cell.
147  labelListList cellAnchorPoints_;
148 
149 
150  // Private Static Functions
151 
152  //- Find value in first n elements of list.
153  static label findPartIndex
154  (
155  const labelList&,
156  const label n,
157  const label val
158  );
159 
160  //- Create boolList with all labels specified set to true
161  // (and rest to false)
162  static boolList expand(const label size, const labelList& labels);
163 
164  //- Create scalarField with all specified labels set to corresponding
165  // value in scalarField.
166  static scalarField expand
167  (
168  const label,
169  const labelList&,
170  const scalarField&
171  );
172 
173  //- Returns -1 or index of first element of lst that cannot be found
174  // in map.
175  static label firstUnique
176  (
177  const labelList& lst,
178  const Map<label>&
179  );
180 
181  // Private Member Functions
182 
183  //- Debugging: write cell's edges and any cut vertices and edges
184  // (so no cell loop determined yet)
185  void writeUncutOBJ(const fileName&, const label cellI) const;
186 
187  //- Debugging: write cell's edges, loop and anchors to directory.
188  void writeOBJ
189  (
190  const fileName& dir,
191  const label cellI,
192  const pointField& loopPoints,
193  const labelList& anchors
194  ) const;
195 
196  //- Find face on cell using the two edges.
197  label edgeEdgeToFace
198  (
199  const label cellI,
200  const label edgeA,
201  const label edgeB
202  ) const;
203 
204 
205  //- Find face on cell using an edge and a vertex.
206  label edgeVertexToFace
207  (
208  const label cellI,
209  const label edgeI,
210  const label vertI
211  ) const;
212 
213  //- Find face using two vertices (guaranteed not to be along edge)
214  label vertexVertexToFace
215  (
216  const label cellI,
217  const label vertA,
218  const label vertB
219  ) const;
220 
221 
222  // Cut addressing
223 
224  //- Calculate faceCuts in face vertex order.
225  void calcFaceCuts() const;
226 
227 
228  // Loop (cuts on cell circumference) calculation
229 
230  //- Find edge (or -1) on faceI using vertices v0,v1
231  label findEdge
232  (
233  const label faceI,
234  const label v0,
235  const label v1
236  ) const;
237 
238  //- Find face on which all cuts are (very rare) or -1.
239  label loopFace(const label cellI, const labelList& loop) const;
240 
241  //- Cross otherCut into next faces (not exclude0, exclude1)
242  bool walkPoint
243  (
244  const label cellI,
245  const label startCut,
246 
247  const label exclude0,
248  const label exclude1,
249 
250  const label otherCut,
251 
252  label& nVisited,
253  labelList& visited
254  ) const;
255 
256  //- Cross cut (which is edge on faceI) onto next face
257  bool crossEdge
258  (
259  const label cellI,
260  const label startCut,
261  const label faceI,
262  const label otherCut,
263 
264  label& nVisited,
265  labelList& visited
266  ) const;
267 
268  // wrapper around visited[nVisited++] = cut. Checks for duplicate
269  // cuts.
270  bool addCut
271  (
272  const label cellI,
273  const label cut,
274  label& nVisited,
275  labelList& visited
276  ) const;
277 
278  //- Walk across faceI following cuts, starting at cut. Stores cuts
279  // visited
280  bool walkFace
281  (
282  const label cellI,
283  const label startCut,
284  const label faceI,
285  const label cut,
286 
287  label& lastCut,
288  label& beforeLastCut,
289  label& nVisited,
290  labelList& visited
291  ) const;
292 
293  //- Walk across cuts (cut edges or cut vertices) of cell. Stops when
294  // hit cut already visited. Returns true when loop of 3 or more
295  // vertices found.
296  bool walkCell
297  (
298  const label cellI,
299  const label startCut, // overall starting cut
300  const label faceI,
301  const label prevCut, // current cut
302  label& nVisited,
303  labelList& visited
304  ) const;
305 
306  //- Determine for every cut cell the face it is cut by.
307  void calcCellLoops(const labelList& cutCells);
308 
309 
310  // Cell anchoring
311 
312  //- Are there enough faces on anchor side of cellI?
313  bool checkFaces
314  (
315  const label cellI,
316  const labelList& anchorPoints
317  ) const;
318 
319  //- Walk unset edges of single cell from starting point and
320  // marks visited edges and vertices with status.
321  void walkEdges
322  (
323  const label cellI,
324  const label pointI,
325  const label status,
326 
327  Map<label>& edgeStatus,
328  Map<label>& pointStatus
329  ) const;
330 
331  //- Check anchor points on 'outside' of loop
332  bool loopAnchorConsistent
333  (
334  const label cellI,
335  const pointField& loopPts,
336  const labelList& anchorPoints
337  ) const;
338 
339  //- Determines set of anchor points given a loop. The loop should
340  // split the cell into two. Returns true if a valid set of anchor
341  // points determined, false otherwise.
342  bool calcAnchors
343  (
344  const label cellI,
345  const labelList& loop,
346  const pointField& loopPts,
347 
348  labelList& anchorPoints
349  ) const;
350 
351  //- Returns coordinates of points on loop with explicitly provided
352  // weights.
353  pointField loopPoints
354  (
355  const labelList& loop,
356  const scalarField& loopWeights
357  ) const;
358 
359  //- Returns weights of loop. Inverse of loopPoints.
360  scalarField loopWeights(const labelList& loop) const;
361 
362  //- Check if cut edges in loop are compatible with ones in
363  // edgeIsCut_
364  bool validEdgeLoop
365  (
366  const labelList& loop,
367  const scalarField& loopWeights
368  ) const;
369 
370  //- Counts number of cuts on face.
371  label countFaceCuts
372  (
373  const label faceI,
374  const labelList& loop
375  ) const;
376 
377  //- Determines if loop through cellI consistent with existing
378  // pattern.
379  bool conservativeValidLoop
380  (
381  const label cellI,
382  const labelList& loop
383  ) const;
384 
385  //- Check if loop is compatible with existing cut pattern in
386  // pointIsCut, edgeIsCut, faceSplitCut.
387  // Calculates and returns for current cell the cut faces and the
388  // points on one side of the loop.
389  bool validLoop
390  (
391  const label cellI,
392  const labelList& loop,
393  const scalarField& loopWeights,
394  Map<edge>& newFaceSplitCut,
395  labelList& anchorPoints
396  ) const;
397 
398  //- Update basic cut information from cellLoops. Assumes cellLoops_
399  // already set and consistent.
400  void setFromCellLoops();
401 
402  //- Update basic cut information for single cell from cellLoop.
403  bool setFromCellLoop
404  (
405  const label cellI,
406  const labelList& loop,
407  const scalarField& loopWeights
408  );
409 
410  //- Update basic cut information from cellLoops. Checks for
411  // consistency with existing cut pattern.
412  void setFromCellLoops
413  (
414  const labelList& cellLabels,
415  const labelListList& cellLoops,
416  const List<scalarField>& cellLoopWeights
417  );
418 
419  //- Cut cells and update basic cut information from cellLoops.
420  // Checks each loop for consistency with existing cut pattern.
421  void setFromCellCutter
422  (
423  const cellLooper&,
424  const List<refineCell>& refCells
425  );
426 
427  //- Same as above but now cut with prescribed plane.
428  void setFromCellCutter
429  (
430  const cellLooper&,
431  const labelList& cellLabels,
432  const List<plane>&
433  );
434 
435  //- Set orientation of loops
436  void orientPlanesAndLoops();
437 
438  //- top level driver: adressing calculation and loop detection
439  void calcLoopsAndAddressing(const labelList& cutCells);
440 
441  //- Check various consistencies.
442  void check() const;
443 
444 
445  //- Disallow default bitwise copy construct
446  cellCuts(const cellCuts&);
447 
448  //- Disallow default bitwise assignment
449  void operator=(const cellCuts&);
450 
451 
452 public:
453 
454  //- Runtime type information
455  ClassName("cellCuts");
456 
457 
458  // Constructors
459 
460  //- Construct from cells to cut and pattern of cuts
461  cellCuts
462  (
463  const polyMesh& mesh,
464  const labelList& cutCells,
465  const labelList& meshVerts,
466  const labelList& meshEdges,
467  const scalarField& meshEdgeWeights
468  );
469 
470  //- Construct from pattern of cuts. Detect cells to cut.
471  cellCuts
472  (
473  const polyMesh& mesh,
474  const labelList& meshVerts,
475  const labelList& meshEdges,
476  const scalarField& meshEdgeWeights
477  );
478 
479  //- Construct from complete cellLoops through specified cells.
480  // Checks for consistency.
481  // Constructs cut-cut addressing and cellAnchorPoints.
482  cellCuts
483  (
484  const polyMesh& mesh,
485  const labelList& cellLabels,
486  const labelListList& cellLoops,
487  const List<scalarField>& cellEdgeWeights
488  );
489 
490  //- Construct from list of cells to cut and direction to cut in
491  // (always through cell centre) and celllooper.
492  cellCuts
493  (
494  const polyMesh& mesh,
495  const cellLooper& cellCutter,
496  const List<refineCell>& refCells
497  );
498 
499  //- Construct from list of cells to cut and plane to cut with and
500  // celllooper. (constructor above always cuts through cell centre)
501  cellCuts
502  (
503  const polyMesh& mesh,
504  const cellLooper& cellCutter,
505  const labelList& cellLabels,
506  const List<plane>& cutPlanes
507  );
508 
509  //- Construct from components
510  cellCuts
511  (
512  const polyMesh& mesh,
513  const boolList& pointIsCut,
514  const boolList& edgeIsCut,
515  const scalarField& edgeWeight,
516  const Map<edge>& faceSplitCut,
517  const labelListList& cellLoops,
518  const label nLoops,
520  );
521 
522 
523  // Destructor
524 
525  ~cellCuts();
526 
527  //- Clear out demand driven storage
528  void clearOut();
529 
530 
531  // Member Functions
532 
533  // Access
534 
535  //- Is mesh point cut
536  const boolList& pointIsCut() const
537  {
538  return pointIsCut_;
539  }
540 
541  //- Is edge cut
542  const boolList& edgeIsCut() const
543  {
544  return edgeIsCut_;
545  }
546 
547  //- If edge is cut gives weight (ratio between start() and end())
548  const scalarField& edgeWeight() const
549  {
550  return edgeWeight_;
551  }
552 
553  //- Cuts per existing face (includes those along edge of face)
554  // Cuts in no particular order
555  const labelListList& faceCuts() const
556  {
557  if (!faceCutsPtr_)
558  {
559  calcFaceCuts();
560  }
561  return *faceCutsPtr_;
562  }
563 
564  //- Gives for split face the two cuts that split the face into two.
565  const Map<edge>& faceSplitCut() const
566  {
567  return faceSplitCut_;
568  }
569 
570  //- For each cut cell the cut along the circumference.
571  const labelListList& cellLoops() const
572  {
573  return cellLoops_;
574  }
575 
576  //- Number of valid cell loops
577  label nLoops() const
578  {
579  return nLoops_;
580  }
581 
582  //- For each cut cell the points on the 'anchor' side of the cell.
584  {
585  return cellAnchorPoints_;
586  }
587 
588 
589  // Other
590 
591  //- Returns coordinates of points on loop for given cell.
592  // Uses cellLoops_ and edgeWeight_
593  pointField loopPoints(const label cellI) const;
594 
595  //- Invert anchor point selection.
597  (
598  const labelList& cellPoints,
599  const labelList& anchorPoints,
600  const labelList& loop
601  ) const;
602 
603  //- Flip loop for cellI. Updates anchor points as well.
604  void flip(const label cellI);
605 
606  //- Flip loop for cellI. Does not update anchors. Use with care
607  // (only if you're sure loop orientation is wrong)
608  void flipLoopOnly(const label cellI);
609 
610 
611  // Write
612 
613  //- debugging:Write list of cuts to stream in OBJ format
614  void writeOBJ
615  (
616  Ostream& os,
617  const pointField& loopPoints,
618  label& vertI
619  ) const;
620 
621  //- debugging:Write all of cuts to stream in OBJ format
622  void writeOBJ(Ostream& os) const;
623 
624  //- debugging:Write edges of cell and loop
625  void writeCellOBJ(const fileName& dir, const label cellI) const;
626 
627 };
628 
629 
630 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
631 
632 } // End namespace Foam
633 
634 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
635 
636 #endif
637 
638 // ************************ vim: set sw=4 sts=4 et: ************************ //