hexCellLooper.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 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 \*---------------------------------------------------------------------------*/
25 
26 #include "hexCellLooper.H"
27 #include "cellFeatures.H"
28 #include "polyMesh.H"
29 #include "cellModeller.H"
30 #include "plane.H"
31 #include "ListOps.H"
32 #include "meshTools.H"
33 #include "OFstream.H"
34 
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41 defineTypeNameAndDebug(hexCellLooper, 0);
42 addToRunTimeSelectionTable(cellLooper, hexCellLooper, word);
43 }
44 
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 // Starting from cut edge start walking.
49 bool Foam::hexCellLooper::walkHex
50 (
51  const label cellI,
52  const label startFaceI,
53  const label startEdgeI,
54 
55  labelList& loop,
56  scalarField& loopWeights
57 ) const
58 {
59  label faceI = startFaceI;
60 
61  label edgeI = startEdgeI;
62 
63  label cutI = 0;
64 
65  do
66  {
67  if (debug & 2)
68  {
69  Pout<< " walkHex : inserting cut onto edge:" << edgeI
70  << " vertices:" << mesh().edges()[edgeI] << endl;
71  }
72 
73  // Store cut through edge. For now cut edges halfway.
74  loop[cutI] = edgeToEVert(edgeI);
75  loopWeights[cutI] = 0.5;
76  cutI++;
77 
78  faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
79 
80  const edge& e = mesh().edges()[edgeI];
81 
82  // Walk two edges further
83  edgeI = meshTools::walkFace(mesh(), faceI, edgeI, e.end(), 2);
84 
85  if (edgeI == startEdgeI)
86  {
87  break;
88  }
89  }
90  while (true);
91 
92  // Checks.
93  if (cutI > 4)
94  {
95  Pout<< "hexCellLooper::walkHex" << "Problem : cell:" << cellI
96  << " collected loop:";
97  writeCuts(Pout, loop, loopWeights);
98  Pout<< "loopWeights:" << loopWeights << endl;
99 
100  return false;
101  }
102  else
103  {
104  return true;
105  }
106 }
107 
108 
109 
110 void Foam::hexCellLooper::makeFace
111 (
112  const labelList& loop,
113  const scalarField& loopWeights,
114 
115  labelList& faceVerts,
116  pointField& facePoints
117 ) const
118 {
119  facePoints.setSize(loop.size());
120  faceVerts.setSize(loop.size());
121 
122  forAll(loop, cutI)
123  {
124  label cut = loop[cutI];
125 
126  if (isEdge(cut))
127  {
128  label edgeI = getEdge(cut);
129 
130  const edge& e = mesh().edges()[edgeI];
131 
132  const point& v0 = mesh().points()[e.start()];
133  const point& v1 = mesh().points()[e.end()];
134 
135  facePoints[cutI] =
136  loopWeights[cutI]*v1 + (1.0-loopWeights[cutI])*v0;
137  }
138  else
139  {
140  label vertI = getVertex(cut);
141 
142  facePoints[cutI] = mesh().points()[vertI];
143  }
144 
145  faceVerts[cutI] = cutI;
146  }
147 }
148 
149 
150 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
151 
152 // Construct from components
153 Foam::hexCellLooper::hexCellLooper(const polyMesh& mesh)
154 :
155  geomCellLooper(mesh),
156  hex_(*(cellModeller::lookup("hex")))
157 {}
158 
159 
160 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
161 
163 {}
164 
165 
166 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
167 
169 (
170  const vector& refDir,
171  const label cellI,
172  const boolList& vertIsCut,
173  const boolList& edgeIsCut,
174  const scalarField& edgeWeight,
175 
176  labelList& loop,
177  scalarField& loopWeights
178 ) const
179 {
180  bool success = false;
181 
182  if (mesh().cellShapes()[cellI].model() == hex_)
183  {
184  // Get starting edge. Note: should be compatible with way refDir is
185  // determined.
186  label edgeI = meshTools::cutDirToEdge(mesh(), cellI, refDir);
187 
188  // Get any face using edge
189  label face0;
190  label face1;
191  meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
192 
193  // Walk circumference of hex, cutting edges only
194  loop.setSize(4);
195  loopWeights.setSize(4);
196 
197  success = walkHex(cellI, face0, edgeI, loop, loopWeights);
198  }
199  else
200  {
201  success = geomCellLooper::cut
202  (
203  refDir,
204  cellI,
205  vertIsCut,
206  edgeIsCut,
207  edgeWeight,
208 
209  loop,
210  loopWeights
211  );
212  }
213 
214  if (debug)
215  {
216  if (loop.empty())
217  {
218  WarningIn("hexCellLooper")
219  << "could not cut cell " << cellI << endl;
220 
221  fileName cutsFile("hexCellLooper_" + name(cellI) + ".obj");
222 
223  Pout<< "hexCellLooper : writing cell to " << cutsFile << endl;
224 
225  OFstream cutsStream(cutsFile);
226 
228  (
229  cutsStream,
230  mesh().cells(),
231  mesh().faces(),
232  mesh().points(),
233  labelList(1, cellI)
234  );
235 
236  return false;
237  }
238 
239  // Check for duplicate cuts.
240  labelHashSet loopSet(loop.size());
241 
242  forAll(loop, elemI)
243  {
244  label elem = loop[elemI];
245 
246  if (loopSet.found(elem))
247  {
248  FatalErrorIn("hexCellLooper::walkHex") << " duplicate cut"
249  << abort(FatalError);
250  }
251  loopSet.insert(elem);
252  }
253 
254 
255  face faceVerts(loop.size());
256  pointField facePoints(loop.size());
257 
258  makeFace(loop, loopWeights, faceVerts, facePoints);
259 
260  if ((faceVerts.mag(facePoints) < SMALL) || (loop.size() < 3))
261  {
262  FatalErrorIn("hexCellLooper::walkHex") << "Face:" << faceVerts
263  << " on points:" << facePoints << endl
264  << UIndirectList<point>(facePoints, faceVerts)()
265  << abort(FatalError);
266  }
267  }
268  return success;
269 }
270 
271 
272 // No shortcuts for cutting with arbitrary plane.
274 (
275  const plane& cutPlane,
276  const label cellI,
277  const boolList& vertIsCut,
278  const boolList& edgeIsCut,
279  const scalarField& edgeWeight,
280 
281  labelList& loop,
282  scalarField& loopWeights
283 ) const
284 {
285  return
287  (
288  cutPlane,
289  cellI,
290  vertIsCut,
291  edgeIsCut,
292  edgeWeight,
293 
294  loop,
295  loopWeights
296  );
297 }
298 
299 
300 // ************************************************************************* //
Output to file stream.
Definition: OFstream.H:81
const cellModel & hex_
Reference to hex cell shape.
Definition: hexCellLooper.H:70
const pointField & points
vector point
Point is a vector.
Definition: point.H:41
const cellShapeList & cellShapes
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
bool empty() const
Return true if the UList is empty (ie, size() is zero).
Definition: UListI.H:313
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:470
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.
label cutDirToEdge(const primitiveMesh &, const label cellI, const vector &cutDir)
Reverse of edgeToCutDir: given direction find edge bundle and.
Definition: meshTools.C:819
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
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
Various functions to operate on Lists.
dynamicFvMesh & mesh
A static collection of cell models, and a means of looking them up.
Definition: cellModeller.H:52
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
Implementation of cellLooper. Does pure geometric cut through cell.
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
#define WarningIn(functionName)
Report a warning using Foam::Warning.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
A List with indirect addressing.
Definition: fvMatrix.H:106
#define forAll(list, i)
Definition: UList.H:421
Macros for easy insertion into run-time selection tables.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label otherFace(const primitiveMesh &, const label cellI, const label faceI, const label edgeI)
Return face on cell using edgeI but not faceI. Throws error.
Definition: meshTools.C:554
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
bool success
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
const cellShapeList & cells
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
A class for handling file names.
Definition: fileName.H:69
const polyMesh & mesh() const
Definition: edgeVertex.H:98
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:106
label walkFace(const primitiveMesh &, const label faceI, const label startEdgeI, const label startVertI, const label nEdges)
Returns label of edge nEdges away from startEdge (in the direction.
Definition: meshTools.C:608
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.
virtual ~hexCellLooper()
Destructor.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53