geomCellLooper.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2019-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "geomCellLooper.H"
30 #include "polyMesh.H"
31 #include "DynamicList.H"
32 #include "plane.H"
33 #include "meshTools.H"
34 #include "SortableList.H"
35 #include "triSurfaceTools.H"
36 #include "HashSet.H"
37 #include "ListOps.H"
38 #include "transform.H"
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 namespace Foam
44 {
45  defineTypeNameAndDebug(geomCellLooper, 0);
46  addToRunTimeSelectionTable(cellLooper, geomCellLooper, word);
47 }
48 
49 
50 // Extension factor of edges to make sure we catch intersections through
51 // edge endpoints
52 const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1e-3;
53 
54 
55 // Snap cuts through edges onto edge endpoints. Fraction of edge length.
56 Foam::scalar Foam::geomCellLooper::snapTol_ = 0.1;
57 
58 
59 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
60 
61 Foam::scalar Foam::geomCellLooper::minEdgeLen(const label vertI) const
62 {
63  scalar minLen = GREAT;
64 
65  const labelList& pEdges = mesh().pointEdges()[vertI];
66 
67  forAll(pEdges, pEdgeI)
68  {
69  const edge& e = mesh().edges()[pEdges[pEdgeI]];
70 
71  minLen = min(minLen, e.mag(mesh().points()));
72  }
73  return minLen;
74 }
75 
76 
77 // Cut edge with plane. Return true and set weight to fraction between
78 // edge-start and edge-end
79 bool Foam::geomCellLooper::cutEdge
80 (
81  const plane& cutPlane,
82  const label edgeI,
83  scalar& weight
84 ) const
85 {
86  const pointField& pts = mesh().points();
87 
88  const edge& e = mesh().edges()[edgeI];
89 
90  scalar s = cutPlane.normalIntersect(pts[e.start()], e.vec(pts));
91 
92  if ((s > -pointEqualTol_) && (s < 1 + pointEqualTol_))
93  {
94  weight = s;
95 
96  return true;
97  }
98  else
99  {
100  // Make sure we don't use this value
101  weight = -GREAT;
102 
103  return false;
104  }
105 }
106 
107 
108 // If edge close enough to endpoint snap to endpoint.
109 Foam::label Foam::geomCellLooper::snapToVert
110 (
111  const scalar tol,
112  const label edgeI,
113  const scalar weight
114 ) const
115 {
116  const edge& e = mesh().edges()[edgeI];
117 
118  if (weight < tol)
119  {
120  return e.start();
121  }
122  else if (weight > (1-tol))
123  {
124  return e.end();
125  }
126 
127  return -1;
128 }
129 
130 
131 void Foam::geomCellLooper::getBase(const vector& n, vector& e0, vector& e1)
132  const
133 {
134  // Guess for vector normal to n.
135  vector base(1,0,0);
136 
137  scalar nComp = n & base;
138 
139  if (mag(nComp) > 0.8)
140  {
141  // Was bad guess. Try with different vector.
142 
143  base.x() = 0;
144  base.y() = 1;
145 
146  nComp = n & base;
147 
148  if (mag(nComp) > 0.8)
149  {
150  base.y() = 0;
151  base.z() = 1;
152 
153  nComp = n & base;
154  }
155  }
156 
157 
158  // Use component normal to n as base vector.
159  e0 = normalised(base - nComp * n);
160 
161  e1 = n ^ e0;
162 
163  //Pout<< "Coord system:" << endl
164  // << " n : " << n << ' ' << mag(n) << endl
165  // << " e0 : " << e0 << ' ' << mag(e0) << endl
166  // << " e1 : " << e1 << ' ' << mag(e1) << endl
167  // << endl;
168 }
169 
170 
171 // Return true if the cut edge at loop[index] is preceded by cuts through
172 // the edge end points.
173 bool Foam::geomCellLooper::edgeEndsCut
174 (
175  const labelList& loop,
176  const label index
177 ) const
178 {
179  label edgeI = getEdge(loop[index]);
180 
181  const edge& e = mesh().edges()[edgeI];
182 
183  const label prevCut = loop[loop.rcIndex(index)];
184  const label nextCut = loop[loop.fcIndex(index)];
185 
186  if (!isEdge(prevCut) && !isEdge(nextCut))
187  {
188  // Cuts before and after are both vertices. Check if both
189  // the edge endpoints
190  label v0 = getVertex(prevCut);
191  label v1 = getVertex(nextCut);
192 
193  if
194  (
195  (v0 == e[0] && v1 == e[1])
196  || (v0 == e[1] && v1 == e[0])
197  )
198  {
199  return true;
200  }
201  }
202  return false;
203 }
204 
205 
206 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
207 
208 Foam::geomCellLooper::geomCellLooper(const polyMesh& mesh)
209 :
211 {}
212 
213 
214 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
215 
217 (
218  const vector& refDir,
219  const label celli,
220  const boolList& vertIsCut,
221  const boolList& edgeIsCut,
222  const scalarField& edgeWeight,
223 
224  labelList& loop,
225  scalarField& loopWeights
226 ) const
227 {
228  // Cut through cell centre normal to refDir.
229  return cut
230  (
231  plane(mesh().cellCentres()[celli], refDir),
232  celli,
233  vertIsCut,
234  edgeIsCut,
235  edgeWeight,
236  loop,
237  loopWeights
238  );
239 }
240 
241 
243 (
244  const plane& cutPlane,
245  const label celli,
246  const boolList&,
247  const boolList&,
248  const scalarField&,
249 
250  labelList& loop,
251  scalarField& loopWeights
252 ) const
253 {
254  const pointField& points = mesh().points();
255  const edgeList& edges = mesh().edges();
256 
257  // Find all cuts through edges.
258  // Special cases:
259  // - edge cut close to endpoint. Snap to endpoint.
260  // - edge endpoint close to plane (but edge not cut). Snap to endpoint.
261  // - both endpoints close to plane. Edge parallel to plane. No need to
262  // cut to edge.
263  // Note: any snap-to-point check must be based on all edges using a point
264  // since otherwise cut through close to point but on neighbouring edge
265  // might not be snapped.
266 
267  // Size overly big.
268  label nEstCuts = 2*mesh().cells()[celli].size();
269 
270  DynamicList<label> localLoop(nEstCuts);
271  DynamicList<scalar> localLoopWeights(nEstCuts);
272 
273  // Points checked. Used to make sure we don't cut edge and edge endpoints
274  // at the same time.
275  labelHashSet checkedPoints(nEstCuts);
276 
277  const labelList& cellEdges = mesh().cellEdges()[celli];
278 
279  forAll(cellEdges, i)
280  {
281  label edgeI = cellEdges[i];
282 
283  const edge& e = edges[edgeI];
284 
285  bool useStart = false;
286 
287  bool useEnd = false;
288 
289  //
290  // Check distance of endpoints to cutPlane
291  //
292 
293  if (checkedPoints.insert(e.start()))
294  {
295  const scalar typLen = pointEqualTol_ * minEdgeLen(e.start());
296 
297  // Check distance to plane.
298  if (cutPlane.distance(points[e.start()]) < typLen)
299  {
300  // Use point.
301  localLoop.append(vertToEVert(e.start()));
302  localLoopWeights.append(-GREAT);
303 
304  useStart = true;
305  }
306  }
307 
308  if (checkedPoints.insert(e.end()))
309  {
310  const scalar typLen = pointEqualTol_ * minEdgeLen(e.end());
311 
312  // Check distance to plane.
313  if (cutPlane.distance(points[e.end()]) < typLen)
314  {
315  // Use point.
316  localLoop.append(vertToEVert(e.end()));
317  localLoopWeights.append(-GREAT);
318 
319  useEnd = true;
320  }
321  }
322 
323  //
324  // Check cut of edge
325  //
326 
327  if (!useEnd && !useStart)
328  {
329  // Edge end points not close to plane so edge not close to
330  // plane. Cut edge.
331  scalar cutWeight;
332 
333  if (cutEdge(cutPlane, edgeI, cutWeight))
334  {
335  // Snap edge cut to vertex.
336  label cutVertI = snapToVert(snapTol_, edgeI, cutWeight);
337 
338  if (cutVertI == -1)
339  {
340  // Proper cut of edge
341  localLoop.append(edgeToEVert(edgeI));
342  localLoopWeights.append(cutWeight);
343  }
344  else
345  {
346  // Cut through edge end point. Might be duplicate
347  // since connected edge might also be snapped to same
348  // endpoint so only insert if unique.
349  label cut = vertToEVert(cutVertI);
350 
351  if (!localLoop.found(cut))
352  {
353  localLoop.append(vertToEVert(cutVertI));
354  localLoopWeights.append(-GREAT);
355  }
356  }
357  }
358  }
359  }
360 
361  if (localLoop.size() <= 2)
362  {
363  return false;
364  }
365 
366  localLoop.shrink();
367  localLoopWeights.shrink();
368 
369 
370  // Get points on loop and centre of loop
371  pointField loopPoints(localLoop.size());
372  point ctr(Zero);
373 
374  forAll(localLoop, i)
375  {
376  loopPoints[i] = coord(localLoop[i], localLoopWeights[i]);
377  ctr += loopPoints[i];
378  }
379  ctr /= localLoop.size();
380 
381 
382  // Get base vectors of coordinate system normal to refDir
383  vector e0, e1;
384  getBase(cutPlane.normal(), e0, e1);
385 
386 
387  // Get sorted angles from point on loop to centre of loop.
388  SortableList<scalar> sortedAngles(localLoop.size());
389 
390  forAll(sortedAngles, i)
391  {
392  const vector toCtr = normalised(loopPoints[i] - ctr);
393 
394  sortedAngles[i] = pseudoAngle(e0, e1, toCtr);
395  }
396  sortedAngles.sort();
397 
398  loop.setSize(localLoop.size());
399  loopWeights.setSize(loop.size());
400 
401 
402  // Fill loop and loopweights according to sorted angles
403 
404  const labelList& indices = sortedAngles.indices();
405 
406  forAll(indices, i)
407  {
408  loop[i] = localLoop[indices[i]];
409  loopWeights[i] = localLoopWeights[indices[i]];
410  }
411 
412 
413  // Check for cut edges along already cut vertices.
414  bool filterLoop = false;
415 
416  forAll(loop, i)
417  {
418  label cut = loop[i];
419 
420  if (isEdge(cut) && edgeEndsCut(loop, i))
421  {
422  filterLoop = true;
423  }
424  }
425 
426  if (filterLoop)
427  {
428  // Found edges in loop where both end points are also in loop. This
429  // is superfluous information so we can remove it.
430 
431  labelList filteredLoop(loop.size());
432  scalarField filteredLoopWeights(loopWeights.size());
433  label filterI = 0;
434 
435  forAll(loop, i)
436  {
437  label cut = loop[i];
438 
439  if (isEdge(cut) && edgeEndsCut(loop, i))
440  {
441  // Superfluous cut. Edge end points cut and edge itself as well.
442  }
443  else
444  {
445  filteredLoop[filterI] = loop[i];
446  filteredLoopWeights[filterI] = loopWeights[i];
447  filterI++;
448  }
449  }
450  filteredLoop.setSize(filterI);
451  filteredLoopWeights.setSize(filterI);
452 
453  loop.transfer(filteredLoop);
454  loopWeights.transfer(filteredLoopWeights);
455  }
456 
457  if (debug&2)
458  {
459  Pout<< "cell:" << celli << endl;
460  forAll(loop, i)
461  {
462  Pout<< "At angle:" << sortedAngles[i] << endl
463  << " cut:";
464 
465  writeCut(Pout, loop[i], loopWeights[i]);
466 
467  Pout<< " coord:" << coord(loop[i], loopWeights[i]);
468 
469  Pout<< endl;
470  }
471  }
472 
473  return true;
474 }
475 
476 
477 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::plane::distance
scalar distance(const point &p) const
Return distance (magnitude) from the given point to the plane.
Definition: planeI.H:75
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::polyMesh::points
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1069
meshTools.H
Foam::plane::normal
const vector & normal() const
The plane unit normal.
Definition: planeI.H:39
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::DynamicList< label >
Foam::geomCellLooper::cut
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.
Definition: geomCellLooper.C:217
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:138
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:63
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::primitiveMesh::edges
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
Definition: primitiveMeshEdges.C:505
geomCellLooper.H
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::primitiveMesh::pointEdges
const labelListList & pointEdges() const
Definition: primitiveMeshEdges.C:516
polyMesh.H
Foam::HashSet< label, Hash< label > >
Foam::DynamicList::shrink
DynamicList< T, SizeMin > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:434
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::plane
Geometric class that creates a 3D plane and can return the intersection point between a line and the ...
Definition: plane.H:89
n
label n
Definition: TABSMDCalcMethod2.H:31
SortableList.H
Foam::Field< scalar >
plane.H
Foam::primitiveMesh::cellEdges
const labelListList & cellEdges() const
Definition: primitiveMeshCellEdges.C:115
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::DynamicList::append
DynamicList< T, SizeMin > & append(const T &val)
Append an element to the end of this list.
Definition: DynamicListI.H:511
Foam::List::transfer
void transfer(List< T > &list)
Definition: List.C:456
Foam::cellLooper
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Definition: cellLooper.H:72
cut
Patchify triangles based on orientation w.r.t other (triangulated or triangulatable) surfaces.
HashSet.H
Foam::SortableList
A list that is sorted upon construction or when explicitly requested with the sort() method.
Definition: List.H:60
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam::pseudoAngle
scalar pseudoAngle(const vector &e0, const vector &e1, const vector &vec)
Estimate angle of vec in coordinate system (e0, e1, e0^e1).
Definition: transform.H:401
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::normalised
VectorSpace< Form, Cmpt, Ncmpts > normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:487
Foam::Vector< scalar >
Foam::List< bool >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry, not overwriting existing entries.
Definition: HashSet.H:191
DynamicList.H
ListOps.H
Various functions to operate on Lists.
transform.H
3D tensor transformation operations.
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
triSurfaceTools.H
Foam::edgeVertex::mesh
const polyMesh & mesh() const
Definition: edgeVertex.H:101