FreeFOAM The Cross-Platform CFD Toolkit
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) 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 \*---------------------------------------------------------------------------*/
25 
26 #include "hexCellLooper.H"
27 #include <meshTools/cellFeatures.H>
28 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/cellModeller.H>
31 #include <OpenFOAM/plane.H>
32 #include <OpenFOAM/ListOps.H>
33 #include <meshTools/meshTools.H>
34 #include <OpenFOAM/OFstream.H>
35 
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 
43 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
44 
45 defineTypeNameAndDebug(hexCellLooper, 0);
46 
47 addToRunTimeSelectionTable(cellLooper, hexCellLooper, word);
48 
49 
50 }
51 
52 
53 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
54 
55 // Starting from cut edge start walking.
56 bool Foam::hexCellLooper::walkHex
57 (
58  const label cellI,
59  const label startFaceI,
60  const label startEdgeI,
61 
62  labelList& loop,
63  scalarField& loopWeights
64 ) const
65 {
66  label faceI = startFaceI;
67 
68  label edgeI = startEdgeI;
69 
70  label cutI = 0;
71 
72  do
73  {
74  if (debug & 2)
75  {
76  Pout<< " walkHex : inserting cut onto edge:" << edgeI
77  << " vertices:" << mesh().edges()[edgeI] << endl;
78  }
79 
80  // Store cut through edge. For now cut edges halfway.
81  loop[cutI] = edgeToEVert(edgeI);
82  loopWeights[cutI] = 0.5;
83  cutI++;
84 
85  faceI = meshTools::otherFace(mesh(), cellI, faceI, edgeI);
86 
87  const edge& e = mesh().edges()[edgeI];
88 
89  // Walk two edges further
90  edgeI = meshTools::walkFace(mesh(), faceI, edgeI, e.end(), 2);
91 
92  if (edgeI == startEdgeI)
93  {
94  break;
95  }
96  }
97  while (true);
98 
99  // Checks.
100  if (cutI > 4)
101  {
102  Pout<< "hexCellLooper::walkHex" << "Problem : cell:" << cellI
103  << " collected loop:";
104  writeCuts(Pout, loop, loopWeights);
105  Pout<< "loopWeights:" << loopWeights << endl;
106 
107  return false;
108  }
109  else
110  {
111  return true;
112  }
113 }
114 
115 
116 
117 void Foam::hexCellLooper::makeFace
118 (
119  const labelList& loop,
120  const scalarField& loopWeights,
121 
122  labelList& faceVerts,
123  pointField& facePoints
124 ) const
125 {
126  facePoints.setSize(loop.size());
127  faceVerts.setSize(loop.size());
128 
129  forAll(loop, cutI)
130  {
131  label cut = loop[cutI];
132 
133  if (isEdge(cut))
134  {
135  label edgeI = getEdge(cut);
136 
137  const edge& e = mesh().edges()[edgeI];
138 
139  const point& v0 = mesh().points()[e.start()];
140  const point& v1 = mesh().points()[e.end()];
141 
142  facePoints[cutI] =
143  loopWeights[cutI]*v1 + (1.0-loopWeights[cutI])*v0;
144  }
145  else
146  {
147  label vertI = getVertex(cut);
148 
149  facePoints[cutI] = mesh().points()[vertI];
150  }
151 
152  faceVerts[cutI] = cutI;
153  }
154 }
155 
156 
157 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
158 
159 // Construct from components
160 Foam::hexCellLooper::hexCellLooper(const polyMesh& mesh)
161 :
162  geomCellLooper(mesh),
163  hex_(*(cellModeller::lookup("hex")))
164 {}
165 
166 
167 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
168 
170 {}
171 
172 
173 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
174 
176 (
177  const vector& refDir,
178  const label cellI,
179  const boolList& vertIsCut,
180  const boolList& edgeIsCut,
181  const scalarField& edgeWeight,
182 
183  labelList& loop,
184  scalarField& loopWeights
185 ) const
186 {
187  bool success = false;
188 
189  if (mesh().cellShapes()[cellI].model() == hex_)
190  {
191  // Get starting edge. Note: should be compatible with way refDir is
192  // determined.
193  label edgeI = meshTools::cutDirToEdge(mesh(), cellI, refDir);
194 
195  // Get any face using edge
196  label face0;
197  label face1;
198  meshTools::getEdgeFaces(mesh(), cellI, edgeI, face0, face1);
199 
200  // Walk circumference of hex, cutting edges only
201  loop.setSize(4);
202  loopWeights.setSize(4);
203 
204  success = walkHex(cellI, face0, edgeI, loop, loopWeights);
205  }
206  else
207  {
208  success = geomCellLooper::cut
209  (
210  refDir,
211  cellI,
212  vertIsCut,
213  edgeIsCut,
214  edgeWeight,
215 
216  loop,
217  loopWeights
218  );
219  }
220 
221  if (debug)
222  {
223  if (loop.empty())
224  {
225  WarningIn("hexCellLooper")
226  << "could not cut cell " << cellI << endl;
227 
228  fileName cutsFile("hexCellLooper_" + name(cellI) + ".obj");
229 
230  Pout<< "hexCellLooper : writing cell to " << cutsFile << endl;
231 
232  OFstream cutsStream(cutsFile);
233 
235  (
236  cutsStream,
237  mesh().cells(),
238  mesh().faces(),
239  mesh().points(),
240  labelList(1, cellI)
241  );
242 
243  return false;
244  }
245 
246  // Check for duplicate cuts.
247  labelHashSet loopSet(loop.size());
248 
249  forAll(loop, elemI)
250  {
251  label elem = loop[elemI];
252 
253  if (loopSet.found(elem))
254  {
255  FatalErrorIn("hexCellLooper::walkHex") << " duplicate cut"
256  << abort(FatalError);
257  }
258  loopSet.insert(elem);
259  }
260 
261 
262  face faceVerts(loop.size());
263  pointField facePoints(loop.size());
264 
265  makeFace(loop, loopWeights, faceVerts, facePoints);
266 
267  if ((faceVerts.mag(facePoints) < SMALL) || (loop.size() < 3))
268  {
269  FatalErrorIn("hexCellLooper::walkHex") << "Face:" << faceVerts
270  << " on points:" << facePoints << endl
271  << UIndirectList<point>(facePoints, faceVerts)()
272  << abort(FatalError);
273  }
274  }
275  return success;
276 }
277 
278 
279 // No shortcuts for cutting with arbitrary plane.
281 (
282  const plane& cutPlane,
283  const label cellI,
284  const boolList& vertIsCut,
285  const boolList& edgeIsCut,
286  const scalarField& edgeWeight,
287 
288  labelList& loop,
289  scalarField& loopWeights
290 ) const
291 {
292  return
294  (
295  cutPlane,
296  cellI,
297  vertIsCut,
298  edgeIsCut,
299  edgeWeight,
300 
301  loop,
302  loopWeights
303  );
304 }
305 
306 
307 // ************************ vim: set sw=4 sts=4 et: ************************ //