All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cuttingPlane.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2019 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 "cuttingPlane.H"
27 #include "primitiveMesh.H"
28 #include "linePointRef.H"
29 #include "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 void Foam::cuttingPlane::calcCutCells
43 (
44  const primitiveMesh& mesh,
45  const scalarField& dotProducts,
46  const labelUList& cellIdLabels
47 )
48 {
49  const labelListList& cellEdges = mesh.cellEdges();
50  const edgeList& edges = mesh.edges();
51 
52  label listSize = cellEdges.size();
53  if (notNull(cellIdLabels))
54  {
55  listSize = cellIdLabels.size();
56  }
57 
58  cutCells_.setSize(listSize);
59  label cutcelli(0);
60 
61  // fp: check if the list (cellZone) is not empty.
62  const bool isZoneEmpty
63  (
64  (returnReduce(cellIdLabels.size(), sumOp<label>()) > 0) ? false : true
65  );
66 
67  // Find the cut cells by detecting any cell that uses points with
68  // opposing dotProducts.
69  for (label listI = 0; listI < listSize; ++listI)
70  {
71  label celli = listI;
72 
73  if (notNull(cellIdLabels))
74  {
75  celli = cellIdLabels[listI];
76  }
77  else
78  {
79  // fp: in parallel computation, if the cellZone exists globally
80  // but not locally, the postprocessing must be still be limited to
81  // the crossing plane.
82  if (!isZoneEmpty)
83  {
84  cutCells_.setSize(0);
85  return;
86  }
87  }
88  const labelList& cEdges = cellEdges[celli];
89 
90  label nCutEdges = 0;
91 
92  forAll(cEdges, i)
93  {
94  const edge& e = edges[cEdges[i]];
95 
96  if
97  (
98  (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
99  || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
100  )
101  {
102  nCutEdges++;
103 
104  if (nCutEdges > 2)
105  {
106  cutCells_[cutcelli++] = celli;
107 
108  break;
109  }
110  }
111  }
112  }
113 
114  // Set correct list size
115  cutCells_.setSize(cutcelli);
116 }
117 
118 
119 void Foam::cuttingPlane::intersectEdges
120 (
121  const primitiveMesh& mesh,
122  const scalarField& dotProducts,
123  List<label>& edgePoint
124 )
125 {
126  // Use the dotProducts to find out the cut edges.
127  const edgeList& edges = mesh.edges();
128  const pointField& points = mesh.points();
129 
130  // Per edge -1 or the label of the intersection point
131  edgePoint.setSize(edges.size());
132 
133  DynamicList<point> dynCuttingPoints(4*cutCells_.size());
134 
135  forAll(edges, edgeI)
136  {
137  const edge& e = edges[edgeI];
138 
139  if
140  (
141  (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
142  || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
143  )
144  {
145  // Edge is cut
146  edgePoint[edgeI] = dynCuttingPoints.size();
147 
148  const point& p0 = points[e[0]];
149  const point& p1 = points[e[1]];
150 
151  scalar alpha = lineIntersect(linePointRef(p0, p1));
152 
153  if (alpha < zeroish)
154  {
155  dynCuttingPoints.append(p0);
156  }
157  else if (alpha >= 1.0)
158  {
159  dynCuttingPoints.append(p1);
160  }
161  else
162  {
163  dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
164  }
165  }
166  else
167  {
168  edgePoint[edgeI] = -1;
169  }
170  }
171 
172  this->storedPoints().transfer(dynCuttingPoints);
173 }
174 
175 
176 bool Foam::cuttingPlane::walkCell
177 (
178  const primitiveMesh& mesh,
179  const labelUList& edgePoint,
180  const label celli,
181  const label startEdgeI,
182  DynamicList<label>& faceVerts
183 )
184 {
185  label facei = -1;
186  label edgeI = startEdgeI;
187 
188  label nIter = 0;
189 
190  faceVerts.clear();
191  do
192  {
193  faceVerts.append(edgePoint[edgeI]);
194 
195  // Cross edge to other face
196  facei = meshTools::otherFace(mesh, celli, facei, edgeI);
197 
198  // Find next cut edge on face.
199  const labelList& fEdges = mesh.faceEdges()[facei];
200 
201  label nextEdgeI = -1;
202 
203  // Note: here is where we should check for whether there are more
204  // than 2 intersections with the face (warped/non-convex face).
205  // If so should e.g. decompose the cells on both faces and redo
206  // the calculation.
207 
208  forAll(fEdges, i)
209  {
210  label edge2I = fEdges[i];
211 
212  if (edge2I != edgeI && edgePoint[edge2I] != -1)
213  {
214  nextEdgeI = edge2I;
215  break;
216  }
217  }
218 
219  if (nextEdgeI == -1)
220  {
221  // Did not find another cut edge on facei. Do what?
223  << "Did not find closed walk along surface of cell " << celli
224  << " starting from edge " << startEdgeI
225  << " in " << nIter << " iterations." << nl
226  << "Collected cutPoints so far:" << faceVerts
227  << endl;
228 
229  return false;
230  }
231 
232  edgeI = nextEdgeI;
233 
234  nIter++;
235 
236  if (nIter > 1000)
237  {
239  << "Did not find closed walk along surface of cell " << celli
240  << " starting from edge " << startEdgeI
241  << " in " << nIter << " iterations." << nl
242  << "Collected cutPoints so far:" << faceVerts
243  << endl;
244  return false;
245  }
246 
247  } while (edgeI != startEdgeI);
248 
249 
250  if (faceVerts.size() >= 3)
251  {
252  return true;
253  }
254  else
255  {
257  << "Did not find closed walk along surface of cell " << celli
258  << " starting from edge " << startEdgeI << nl
259  << "Collected cutPoints so far:" << faceVerts
260  << endl;
261 
262  return false;
263  }
264 }
265 
266 
267 void Foam::cuttingPlane::walkCellCuts
268 (
269  const primitiveMesh& mesh,
270  const bool triangulate,
271  const labelUList& edgePoint
272 )
273 {
274  const pointField& cutPoints = this->points();
275 
276  // use dynamic lists to handle triangulation and/or missed cuts
277  DynamicList<face> dynCutFaces(cutCells_.size());
278  DynamicList<label> dynCutCells(cutCells_.size());
279 
280  // scratch space for calculating the face vertices
281  DynamicList<label> faceVerts(10);
282 
283  forAll(cutCells_, i)
284  {
285  label celli = cutCells_[i];
286 
287  // Find the starting edge to walk from.
288  const labelList& cEdges = mesh.cellEdges()[celli];
289 
290  label startEdgeI = -1;
291 
292  forAll(cEdges, cEdgeI)
293  {
294  label edgeI = cEdges[cEdgeI];
295 
296  if (edgePoint[edgeI] != -1)
297  {
298  startEdgeI = edgeI;
299  break;
300  }
301  }
302 
303  // Check for the unexpected ...
304  if (startEdgeI == -1)
305  {
307  << "Cannot find cut edge for cut cell " << celli
308  << abort(FatalError);
309  }
310 
311  // Walk from starting edge around the circumference of the cell.
312  bool okCut = walkCell
313  (
314  mesh,
315  edgePoint,
316  celli,
317  startEdgeI,
318  faceVerts
319  );
320 
321  if (okCut)
322  {
323  face f(faceVerts);
324 
325  // Orient face to point in the same direction as the plane normal
326  if ((f.area(cutPoints) & normal()) < 0)
327  {
328  f.flip();
329  }
330 
331  // the cut faces are usually quite ugly, so optionally triangulate
332  if (triangulate)
333  {
334  label nTri = f.triangles(cutPoints, dynCutFaces);
335  while (nTri--)
336  {
337  dynCutCells.append(celli);
338  }
339  }
340  else
341  {
342  dynCutFaces.append(f);
343  dynCutCells.append(celli);
344  }
345  }
346  }
347 
348  this->storedFaces().transfer(dynCutFaces);
349  cutCells_.transfer(dynCutCells);
350 }
351 
352 
353 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
354 
356 (
357  const primitiveMesh& mesh,
358  const bool triangulate,
359  const labelUList& cellIdLabels
360 )
361 {
363  cutCells_.clear();
364 
365  const scalarField dotProducts((mesh.points() - refPoint()) & normal());
366 
367  // Determine cells that are (probably) cut.
368  calcCutCells(mesh, dotProducts, cellIdLabels);
369 
370  // Determine cutPoints and return list of edge cuts.
371  // per edge -1 or the label of the intersection point
372  labelList edgePoint;
373  intersectEdges(mesh, dotProducts, edgePoint);
374 
375  // Do topological walk around cell to find closed loop.
376  walkCellCuts(mesh, triangulate, edgePoint);
377 }
378 
379 
381 (
382  const labelUList& faceMap
383 )
384 {
385  // recalculate the cells cut
386  if (notNull(faceMap) && faceMap.size())
387  {
388  MeshStorage::remapFaces(faceMap);
389 
390  List<label> newCutCells(faceMap.size());
391  forAll(faceMap, facei)
392  {
393  newCutCells[facei] = cutCells_[faceMap[facei]];
394  }
395  cutCells_.transfer(newCutCells);
396  }
397 }
398 
399 
400 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
401 
403 :
404  plane(pln)
405 {}
406 
407 
409 (
410  const plane& pln,
411  const primitiveMesh& mesh,
412  const bool triangulate,
413  const labelUList& cellIdLabels
414 )
415 :
416  plane(pln)
417 {
418  reCut(mesh, triangulate, cellIdLabels);
419 }
420 
421 
422 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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 reCut(const primitiveMesh &, const bool triangulate, const labelUList &cellIdLabels=labelUList::null())
Recut mesh with existing planeDesc, restricted to a list of cells.
Definition: cuttingPlane.C:356
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual void clear()
Clear all storage.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const =0
Return mesh points.
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
UList< label > labelUList
Definition: UList.H:65
plane(const vector &normalVector)
Construct from normal vector through the origin.
Definition: plane.C:97
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
List< edge > edgeList
Definition: edgeList.H:38
List< face > & storedFaces()
Non-const access to the faces.
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
const point & refPoint() const
Return or return plane base point.
Definition: plane.C:239
pointField & storedPoints()
Non-const access to global points.
virtual void remapFaces(const labelUList &faceMap)
Remap action on triangulation or cleanup.
Definition: cuttingPlane.C:381
const Field< PointType > & points() const
Return reference to global points.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
volScalarField scalarField(fieldObject, mesh)
static const char nl
Definition: Ostream.H:260
labelList f(nPoints)
scalar lineIntersect(const line< Point, PointRef > &l) const
Return the cutting point between the plane and.
Definition: plane.H:186
bool notNull(const T &t)
Return true if t is not a reference to the nullObject of type T.
Definition: nullObjectI.H:52
void setSize(const label)
Reset size of List.
Definition: List.C:281
const vector & normal() const
Return plane normal.
Definition: plane.C:233
virtual void remapFaces(const labelUList &faceMap)
Set new zones from faceMap.
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
cuttingPlane(const plane &)
Construct plane description without cutting.
Definition: cuttingPlane.C:402
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
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:536