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