1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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.
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
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/>.
28 Description of cuts across cells.
30 Description of cut is given as list of vertices and
31 list of edges to be cut (and position on edge).
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)
37 Introduces concept of 'cut' which is either an existing vertex
41 -# list of cut vertices and list of cut edges. Constructs cell
42 circumference walks ('cellLoops').
44 -# list of cell circumference walks. Will filter them so they
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
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.
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.
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.
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
81 \*---------------------------------------------------------------------------*/
86 #include "edgeVertex.H"
87 #include "labelList.H"
89 #include "scalarField.H"
90 #include "pointField.H"
91 #include "DynamicList.H"
94 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
99 // Forward declaration of classes
105 /*---------------------------------------------------------------------------*\
106 Class cellCuts Declaration
107 \*---------------------------------------------------------------------------*/
115 // Per point/edge status
117 //- Is mesh point cut
118 boolList pointIsCut_;
123 //- If edge is cut gives weight (0->start() to 1->end())
124 scalarField edgeWeight_;
129 //- Cuts per existing face (includes those along edge of face)
130 // Cuts in no particular order.
131 mutable labelListList* faceCutsPtr_;
133 //- Per face : cut across edge (so not along existing edge)
134 // (can only be one per face)
135 Map<edge> faceSplitCut_;
138 // Cell-cut addressing
140 //- Loop across cell circumference
141 labelListList cellLoops_;
143 //- Number of valid loops in cellLoops_
146 //- For each cut cell the points on the 'anchor' side of the cell.
147 labelListList cellAnchorPoints_;
150 // Private Static Functions
152 //- Find value in first n elements of list.
153 static label findPartIndex
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);
164 //- Create scalarField with all specified labels set to corresponding
165 // value in scalarField.
166 static scalarField expand
173 //- Returns -1 or index of first element of lst that cannot be found
175 static label firstUnique
177 const labelList& lst,
181 // Private Member Functions
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;
187 //- Debugging: write cell's edges, loop and anchors to directory.
192 const pointField& loopPoints,
193 const labelList& anchors
196 //- Find face on cell using the two edges.
205 //- Find face on cell using an edge and a vertex.
206 label edgeVertexToFace
213 //- Find face using two vertices (guaranteed not to be along edge)
214 label vertexVertexToFace
224 //- Calculate faceCuts in face vertex order.
225 void calcFaceCuts() const;
228 // Loop (cuts on cell circumference) calculation
230 //- Find edge (or -1) on faceI using vertices v0,v1
238 //- Find face on which all cuts are (very rare) or -1.
239 label loopFace(const label cellI, const labelList& loop) const;
241 //- Cross otherCut into next faces (not exclude0, exclude1)
245 const label startCut,
247 const label exclude0,
248 const label exclude1,
250 const label otherCut,
256 //- Cross cut (which is edge on faceI) onto next face
260 const label startCut,
262 const label otherCut,
268 // wrapper around visited[nVisited++] = cut. Checks for duplicate
278 //- Walk across faceI following cuts, starting at cut. Stores cuts
283 const label startCut,
288 label& beforeLastCut,
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
299 const label startCut, // overall starting cut
301 const label prevCut, // current cut
306 //- Determine for every cut cell the face it is cut by.
307 void calcCellLoops(const labelList& cutCells);
312 //- Are there enough faces on anchor side of cellI?
316 const labelList& anchorPoints
319 //- Walk unset edges of single cell from starting point and
320 // marks visited edges and vertices with status.
327 Map<label>& edgeStatus,
328 Map<label>& pointStatus
331 //- Check anchor points on 'outside' of loop
332 bool loopAnchorConsistent
335 const pointField& loopPts,
336 const labelList& anchorPoints
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.
345 const labelList& loop,
346 const pointField& loopPts,
348 labelList& anchorPoints
351 //- Returns coordinates of points on loop with explicitly provided
353 pointField loopPoints
355 const labelList& loop,
356 const scalarField& loopWeights
359 //- Returns weights of loop. Inverse of loopPoints.
360 scalarField loopWeights(const labelList& loop) const;
362 //- Check if cut edges in loop are compatible with ones in
366 const labelList& loop,
367 const scalarField& loopWeights
370 //- Counts number of cuts on face.
374 const labelList& loop
377 //- Determines if loop through cellI consistent with existing
379 bool conservativeValidLoop
382 const labelList& loop
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.
392 const labelList& loop,
393 const scalarField& loopWeights,
394 Map<edge>& newFaceSplitCut,
395 labelList& anchorPoints
398 //- Update basic cut information from cellLoops. Assumes cellLoops_
399 // already set and consistent.
400 void setFromCellLoops();
402 //- Update basic cut information for single cell from cellLoop.
406 const labelList& loop,
407 const scalarField& loopWeights
410 //- Update basic cut information from cellLoops. Checks for
411 // consistency with existing cut pattern.
412 void setFromCellLoops
414 const labelList& cellLabels,
415 const labelListList& cellLoops,
416 const List<scalarField>& cellLoopWeights
419 //- Cut cells and update basic cut information from cellLoops.
420 // Checks each loop for consistency with existing cut pattern.
421 void setFromCellCutter
424 const List<refineCell>& refCells
427 //- Same as above but now cut with prescribed plane.
428 void setFromCellCutter
431 const labelList& cellLabels,
435 //- Set orientation of loops
436 void orientPlanesAndLoops();
438 //- top level driver: adressing calculation and loop detection
439 void calcLoopsAndAddressing(const labelList& cutCells);
441 //- Check various consistencies.
445 //- Disallow default bitwise copy construct
446 cellCuts(const cellCuts&);
448 //- Disallow default bitwise assignment
449 void operator=(const cellCuts&);
454 //- Runtime type information
455 ClassName("cellCuts");
460 //- Construct from cells to cut and pattern of cuts
463 const polyMesh& mesh,
464 const labelList& cutCells,
465 const labelList& meshVerts,
466 const labelList& meshEdges,
467 const scalarField& meshEdgeWeights
470 //- Construct from pattern of cuts. Detect cells to cut.
473 const polyMesh& mesh,
474 const labelList& meshVerts,
475 const labelList& meshEdges,
476 const scalarField& meshEdgeWeights
479 //- Construct from complete cellLoops through specified cells.
480 // Checks for consistency.
481 // Constructs cut-cut addressing and cellAnchorPoints.
484 const polyMesh& mesh,
485 const labelList& cellLabels,
486 const labelListList& cellLoops,
487 const List<scalarField>& cellEdgeWeights
490 //- Construct from list of cells to cut and direction to cut in
491 // (always through cell centre) and celllooper.
494 const polyMesh& mesh,
495 const cellLooper& cellCutter,
496 const List<refineCell>& refCells
499 //- Construct from list of cells to cut and plane to cut with and
500 // celllooper. (constructor above always cuts through cell centre)
503 const polyMesh& mesh,
504 const cellLooper& cellCutter,
505 const labelList& cellLabels,
506 const List<plane>& cutPlanes
509 //- Construct from components
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,
519 const labelListList& cellAnchorPoints
526 //- Clear out demand driven storage
534 //- Is mesh point cut
535 const boolList& pointIsCut() const
541 const boolList& edgeIsCut() const
546 //- If edge is cut gives weight (ratio between start() and end())
547 const scalarField& edgeWeight() const
552 //- Cuts per existing face (includes those along edge of face)
553 // Cuts in no particular order
554 const labelListList& faceCuts() const
560 return *faceCutsPtr_;
563 //- Gives for split face the two cuts that split the face into two.
564 const Map<edge>& faceSplitCut() const
566 return faceSplitCut_;
569 //- For each cut cell the cut along the circumference.
570 const labelListList& cellLoops() const
575 //- Number of valid cell loops
581 //- For each cut cell the points on the 'anchor' side of the cell.
582 const labelListList& cellAnchorPoints() const
584 return cellAnchorPoints_;
590 //- Returns coordinates of points on loop for given cell.
591 // Uses cellLoops_ and edgeWeight_
592 pointField loopPoints(const label cellI) const;
594 //- Invert anchor point selection.
595 labelList nonAnchorPoints
597 const labelList& cellPoints,
598 const labelList& anchorPoints,
599 const labelList& loop
602 //- Flip loop for cellI. Updates anchor points as well.
603 void flip(const label cellI);
605 //- Flip loop for cellI. Does not update anchors. Use with care
606 // (only if you're sure loop orientation is wrong)
607 void flipLoopOnly(const label cellI);
612 //- debugging:Write list of cuts to stream in OBJ format
616 const pointField& loopPoints,
620 //- debugging:Write all of cuts to stream in OBJ format
621 void writeOBJ(Ostream& os) const;
623 //- debugging:Write edges of cell and loop
624 void writeCellOBJ(const fileName& dir, const label cellI) const;
629 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
631 } // End namespace Foam
633 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
637 // ************************************************************************* //