FreeFOAM The Cross-Platform CFD Toolkit
cuttingPlane.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 "cuttingPlane.H"
27 #include <OpenFOAM/primitiveMesh.H>
28 #include <OpenFOAM/linePointRef.H>
29 #include <meshTools/meshTools.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 // set values for what is close to zero and what is considered to
34 // be positive (and not just rounding noise)
36 const Foam::scalar zeroish = Foam::SMALL;
37 const Foam::scalar positive = Foam::SMALL * 1E3;
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 // Find cut cells
43 void Foam::cuttingPlane::calcCutCells
44 (
45  const primitiveMesh& mesh,
46  const scalarField& dotProducts,
47  const UList<label>& cellIdLabels
48 )
49 {
50  const labelListList& cellEdges = mesh.cellEdges();
51  const edgeList& edges = mesh.edges();
52 
53  label listSize = cellEdges.size();
54  if (&cellIdLabels)
55  {
56  listSize = cellIdLabels.size();
57  }
58 
59  cutCells_.setSize(listSize);
60  label cutcellI(0);
61 
62  // Find the cut cells by detecting any cell that uses points with
63  // opposing dotProducts.
64  for (label listI = 0; listI < listSize; ++listI)
65  {
66  label cellI = listI;
67  if (&cellIdLabels)
68  {
69  cellI = cellIdLabels[listI];
70  }
71 
72  const labelList& cEdges = cellEdges[cellI];
73 
74  label nCutEdges = 0;
75 
76  forAll(cEdges, i)
77  {
78  const edge& e = edges[cEdges[i]];
79 
80  if
81  (
82  (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
83  || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
84  )
85  {
86  nCutEdges++;
87 
88  if (nCutEdges > 2)
89  {
90  cutCells_[cutcellI++] = cellI;
91 
92  break;
93  }
94  }
95  }
96  }
97 
98  // Set correct list size
99  cutCells_.setSize(cutcellI);
100 }
101 
102 
103 // Determine for each edge the intersection point. Calculates
104 // - cutPoints_ : coordinates of all intersection points
105 // - edgePoint : per edge -1 or the index into cutPoints
106 void Foam::cuttingPlane::intersectEdges
107 (
108  const primitiveMesh& mesh,
109  const scalarField& dotProducts,
110  List<label>& edgePoint
111 )
112 {
113  // Use the dotProducts to find out the cut edges.
114  const edgeList& edges = mesh.edges();
115  const pointField& points = mesh.points();
116 
117  // Per edge -1 or the label of the intersection point
118  edgePoint.setSize(edges.size());
119 
120  DynamicList<point> dynCuttingPoints(4*cutCells_.size());
121 
122  forAll(edges, edgeI)
123  {
124  const edge& e = edges[edgeI];
125 
126  if
127  (
128  (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
129  || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
130  )
131  {
132  // Edge is cut
133  edgePoint[edgeI] = dynCuttingPoints.size();
134 
135  const point& p0 = points[e[0]];
136  const point& p1 = points[e[1]];
137 
138  scalar alpha = lineIntersect(linePointRef(p0, p1));
139 
140  if (alpha < zeroish)
141  {
142  dynCuttingPoints.append(p0);
143  }
144  else if (alpha >= 1.0)
145  {
146  dynCuttingPoints.append(p1);
147  }
148  else
149  {
150  dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
151  }
152  }
153  else
154  {
155  edgePoint[edgeI] = -1;
156  }
157  }
158 
159  this->storedPoints().transfer(dynCuttingPoints);
160 }
161 
162 
163 // Coming from startEdgeI cross the edge to the other face
164 // across to the next cut edge.
165 bool Foam::cuttingPlane::walkCell
166 (
167  const primitiveMesh& mesh,
168  const UList<label>& edgePoint,
169  const label cellI,
170  const label startEdgeI,
171  DynamicList<label>& faceVerts
172 )
173 {
174  label faceI = -1;
175  label edgeI = startEdgeI;
176 
177  label nIter = 0;
178 
179  faceVerts.clear();
180  do
181  {
182  faceVerts.append(edgePoint[edgeI]);
183 
184  // Cross edge to other face
185  faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);
186 
187  // Find next cut edge on face.
188  const labelList& fEdges = mesh.faceEdges()[faceI];
189 
190  label nextEdgeI = -1;
191 
192  //Note: here is where we should check for whether there are more
193  // than 2 intersections with the face (warped/non-convex face).
194  // If so should e.g. decompose the cells on both faces and redo
195  // the calculation.
196 
197  forAll(fEdges, i)
198  {
199  label edge2I = fEdges[i];
200 
201  if (edge2I != edgeI && edgePoint[edge2I] != -1)
202  {
203  nextEdgeI = edge2I;
204  break;
205  }
206  }
207 
208  if (nextEdgeI == -1)
209  {
210  // Did not find another cut edge on faceI. Do what?
211  WarningIn("Foam::cuttingPlane::walkCell")
212  << "Did not find closed walk along surface of cell " << cellI
213  << " starting from edge " << startEdgeI
214  << " in " << nIter << " iterations." << nl
215  << "Collected cutPoints so far:" << faceVerts
216  << endl;
217 
218  return false;
219  }
220 
221  edgeI = nextEdgeI;
222 
223  nIter++;
224 
225  if (nIter > 1000)
226  {
227  WarningIn("Foam::cuttingPlane::walkCell")
228  << "Did not find closed walk along surface of cell " << cellI
229  << " starting from edge " << startEdgeI
230  << " in " << nIter << " iterations." << nl
231  << "Collected cutPoints so far:" << faceVerts
232  << endl;
233  return false;
234  }
235 
236  } while (edgeI != startEdgeI);
237 
238 
239  if (faceVerts.size() >= 3)
240  {
241  return true;
242  }
243  else
244  {
245  WarningIn("Foam::cuttingPlane::walkCell")
246  << "Did not find closed walk along surface of cell " << cellI
247  << " starting from edge " << startEdgeI << nl
248  << "Collected cutPoints so far:" << faceVerts
249  << endl;
250 
251  return false;
252  }
253 }
254 
255 
256 // For every cut cell determine a walk through all? its cuts.
257 void Foam::cuttingPlane::walkCellCuts
258 (
259  const primitiveMesh& mesh,
260  const UList<label>& edgePoint
261 )
262 {
263  const pointField& cutPoints = this->points();
264 
265  // use dynamic lists to handle triangulation and/or missed cuts
266  DynamicList<face> dynCutFaces(cutCells_.size());
267  DynamicList<label> dynCutCells(cutCells_.size());
268 
269  // scratch space for calculating the face vertices
270  DynamicList<label> faceVerts(10);
271 
272  forAll(cutCells_, i)
273  {
274  label cellI = cutCells_[i];
275 
276  // Find the starting edge to walk from.
277  const labelList& cEdges = mesh.cellEdges()[cellI];
278 
279  label startEdgeI = -1;
280 
281  forAll(cEdges, cEdgeI)
282  {
283  label edgeI = cEdges[cEdgeI];
284 
285  if (edgePoint[edgeI] != -1)
286  {
287  startEdgeI = edgeI;
288  break;
289  }
290  }
291 
292  // Check for the unexpected ...
293  if (startEdgeI == -1)
294  {
295  FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
296  << "Cannot find cut edge for cut cell " << cellI
297  << abort(FatalError);
298  }
299 
300  // Walk from starting edge around the circumference of the cell.
301  bool okCut = walkCell
302  (
303  mesh,
304  edgePoint,
305  cellI,
306  startEdgeI,
307  faceVerts
308  );
309 
310  if (okCut)
311  {
312  face f(faceVerts);
313 
314  // Orient face to point in the same direction as the plane normal
315  if ((f.normal(cutPoints) & normal()) < 0)
316  {
317  f = f.reverseFace();
318  }
319 
320  // the cut faces are usually quite ugly, so always triangulate
321  label nTri = f.triangles(cutPoints, dynCutFaces);
322  while (nTri--)
323  {
324  dynCutCells.append(cellI);
325  }
326  }
327  }
328 
329  this->storedFaces().transfer(dynCutFaces);
330  cutCells_.transfer(dynCutCells);
331 }
332 
333 
334 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
335 
336 // Construct without cutting
338 :
339  plane(pln)
340 {}
341 
342 
343 // Construct from plane and mesh reference, restricted to a list of cells
345 (
346  const plane& pln,
347  const primitiveMesh& mesh,
348  const UList<label>& cellIdLabels
349 )
350 :
351  plane(pln)
352 {
353  reCut(mesh, cellIdLabels);
354 }
355 
356 
357 
358 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
359 
360 // recut mesh with existing planeDesc
362 (
363  const primitiveMesh& mesh,
364  const UList<label>& cellIdLabels
365 )
366 {
368  cutCells_.clear();
369 
370  scalarField dotProducts = (mesh.points() - refPoint()) & normal();
371 
372  // Determine cells that are (probably) cut.
373  calcCutCells(mesh, dotProducts, cellIdLabels);
374 
375  // Determine cutPoints and return list of edge cuts.
376  // per edge -1 or the label of the intersection point
377  labelList edgePoint;
378  intersectEdges(mesh, dotProducts, edgePoint);
379 
380  // Do topological walk around cell to find closed loop.
381  walkCellCuts(mesh, edgePoint);
382 }
383 
384 
385 // remap action on triangulation
387 (
388  const UList<label>& faceMap
389 )
390 {
391  // recalculate the cells cut
392  if (&faceMap && faceMap.size())
393  {
394  MeshStorage::remapFaces(faceMap);
395 
396  List<label> newCutCells(faceMap.size());
397  forAll(faceMap, faceI)
398  {
399  newCutCells[faceI] = cutCells_[faceMap[faceI]];
400  }
401  cutCells_.transfer(newCutCells);
402  }
403 }
404 
405 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
406 
408 {
409  // Check for assignment to self
410  if (this == &rhs)
411  {
412  FatalErrorIn ("Foam::cuttingPlane::operator=(const cuttingPlane&)")
413  << "Attempted assignment to self"
414  << abort(FatalError);
415  }
416 
417  static_cast<MeshStorage&>(*this) = rhs;
418  static_cast<plane&>(*this) = rhs;
419  cutCells_ = rhs.cutCells();
420 }
421 
422 
423 // ************************ vim: set sw=4 sts=4 et: ************************ //