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 | Copyright (C) 2011-2012 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 "geomCellLooper.H"
27 #include "polyMesh.H"
28 #include "DynamicList.H"
29 #include "plane.H"
30 #include "meshTools.H"
31 #include "SortableList.H"
32 #include "triSurfaceTools.H"
33 #include "HashSet.H"
34 #include "ListOps.H"
35 #include "transform.H"
36 
38 
39 
40 // Extension factor of edges to make sure we catch intersections through
41 // edge endpoints
42 const Foam::scalar Foam::geomCellLooper::pointEqualTol_ = 1e-3;
43 
44 
45 // Snap cuts through edges onto edge endpoints. Fraction of edge length.
46 Foam::scalar Foam::geomCellLooper::snapTol_ = 0.1;
47 
48 
49 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
50 
51 namespace Foam
52 {
53 
54 defineTypeNameAndDebug(geomCellLooper, 0);
55 addToRunTimeSelectionTable(cellLooper, geomCellLooper, word);
56 
57 }
58 
59 
60 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 
62 Foam::scalar Foam::geomCellLooper::minEdgeLen(const label vertI) const
63 {
64  scalar minLen = GREAT;
65 
66  const labelList& pEdges = mesh().pointEdges()[vertI];
67 
68  forAll(pEdges, pEdgeI)
69  {
70  const edge& e = mesh().edges()[pEdges[pEdgeI]];
71 
72  minLen = min(minLen, e.mag(mesh().points()));
73  }
74  return minLen;
75 }
76 
77 
78 // Cut edge with plane. Return true and set weight to fraction between
79 // edge-start and edge-end
80 bool Foam::geomCellLooper::cutEdge
81 (
82  const plane& cutPlane,
83  const label edgeI,
84  scalar& weight
85 ) const
86 {
87  const pointField& pts = mesh().points();
88 
89  const edge& e = mesh().edges()[edgeI];
90 
91  scalar s = cutPlane.normalIntersect(pts[e.start()], e.vec(pts));
92 
93  if ((s > -pointEqualTol_) && (s < 1 + pointEqualTol_))
94  {
95  weight = s;
96 
97  return true;
98  }
99  else
100  {
101  // Make sure we don't use this value
102  weight = -GREAT;
103 
104  return false;
105  }
106 }
107 
108 
109 // If edge close enough to endpoint snap to endpoint.
110 Foam::label Foam::geomCellLooper::snapToVert
111 (
112  const scalar tol,
113  const label edgeI,
114  const scalar weight
115 ) const
116 {
117  const edge& e = mesh().edges()[edgeI];
118 
119  if (weight < tol)
120  {
121  return e.start();
122  }
123  else if (weight > (1-tol))
124  {
125  return e.end();
126  }
127  else
128  {
129  return -1;
130  }
131 }
132 
133 
134 void Foam::geomCellLooper::getBase(const vector& n, vector& e0, vector& e1)
135  const
136 {
137  // Guess for vector normal to n.
138  vector base(1,0,0);
139 
140  scalar nComp = n & base;
141 
142  if (mag(nComp) > 0.8)
143  {
144  // Was bad guess. Try with different vector.
145 
146  base.x() = 0;
147  base.y() = 1;
148 
149  nComp = n & base;
150 
151  if (mag(nComp) > 0.8)
152  {
153  base.y() = 0;
154  base.z() = 1;
155 
156  nComp = n & base;
157  }
158  }
159 
160 
161  // Use component normal to n as base vector.
162  e0 = base - nComp*n;
163 
164  e0 /= mag(e0) + VSMALL;
165 
166  e1 = n ^ e0;
167 
168  //Pout<< "Coord system:" << endl
169  // << " n : " << n << ' ' << mag(n) << endl
170  // << " e0 : " << e0 << ' ' << mag(e0) << endl
171  // << " e1 : " << e1 << ' ' << mag(e1) << endl
172  // << endl;
173 }
174 
175 
176 // Return true if the cut edge at loop[index] is preceded by cuts through
177 // the edge end points.
178 bool Foam::geomCellLooper::edgeEndsCut
179 (
180  const labelList& loop,
181  const label index
182 ) const
183 {
184  label edgeI = getEdge(loop[index]);
185 
186  const edge& e = mesh().edges()[edgeI];
187 
188  const label prevCut = loop[loop.rcIndex(index)];
189  const label nextCut = loop[loop.fcIndex(index)];
190 
191  if (!isEdge(prevCut) && !isEdge(nextCut))
192  {
193  // Cuts before and after are both vertices. Check if both
194  // the edge endpoints
195  label v0 = getVertex(prevCut);
196  label v1 = getVertex(nextCut);
197 
198  if
199  (
200  (v0 == e[0] && v1 == e[1])
201  || (v0 == e[1] && v1 == e[0])
202  )
203  {
204  return true;
205  }
206  }
207  return false;
208 }
209 
210 
211 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
212 
213 // Construct from components
214 Foam::geomCellLooper::geomCellLooper(const polyMesh& mesh)
215 :
216  cellLooper(mesh)
217 {}
218 
219 
220 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
221 
223 {}
224 
225 
226 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
227 
229 (
230  const vector& refDir,
231  const label cellI,
232  const boolList& vertIsCut,
233  const boolList& edgeIsCut,
234  const scalarField& edgeWeight,
235 
236  labelList& loop,
237  scalarField& loopWeights
238 ) const
239 {
240  // Cut through cell centre normal to refDir.
241  return cut
242  (
243  plane(mesh().cellCentres()[cellI], refDir),
244  cellI,
245  vertIsCut,
246  edgeIsCut,
247  edgeWeight,
248  loop,
249  loopWeights
250  );
251 }
252 
253 
255 (
256  const plane& cutPlane,
257  const label cellI,
258  const boolList&,
259  const boolList&,
260  const scalarField&,
261 
262  labelList& loop,
263  scalarField& loopWeights
264 ) const
265 {
266  const pointField& points = mesh().points();
267  const edgeList& edges = mesh().edges();
268 
269  // Find all cuts through edges.
270  // Special cases:
271  // - edge cut close to endpoint. Snap to endpoint.
272  // - edge endpoint close to plane (but edge not cut). Snap to endpoint.
273  // - both endpoints close to plane. Edge parallel to plane. No need to
274  // cut to edge.
275  // Note: any snap-to-point check must be based on all edges using a point
276  // since otherwise cut through close to point but on neighbouring edge
277  // might not be snapped.
278 
279  // Size overly big.
280  label nEstCuts = 2*mesh().cells()[cellI].size();
281 
282  DynamicList<label> localLoop(nEstCuts);
283  DynamicList<scalar> localLoopWeights(nEstCuts);
284 
285  // Points checked. Used to make sure we don't cut edge and edge endpoints
286  // at the same time.
287  labelHashSet checkedPoints(nEstCuts);
288 
289  const labelList& cellEdges = mesh().cellEdges()[cellI];
290 
291  forAll(cellEdges, i)
292  {
293  label edgeI = cellEdges[i];
294 
295  const edge& e = edges[edgeI];
296 
297  bool useStart = false;
298 
299  bool useEnd = false;
300 
301  //
302  // Check distance of endpoints to cutPlane
303  //
304 
305  if (!checkedPoints.found(e.start()))
306  {
307  checkedPoints.insert(e.start());
308 
309  scalar typStartLen = pointEqualTol_ * minEdgeLen(e.start());
310 
311  // Check distance of startPt to plane.
312  if (cutPlane.distance(points[e.start()]) < typStartLen)
313  {
314  // Use point.
315  localLoop.append(vertToEVert(e.start()));
316  localLoopWeights.append(-GREAT);
317 
318  useStart = true;
319  }
320  }
321  if (!checkedPoints.found(e.end()))
322  {
323  checkedPoints.insert(e.end());
324 
325  scalar typEndLen = pointEqualTol_ * minEdgeLen(e.end());
326 
327  // Check distance of endPt to plane.
328  if (cutPlane.distance(points[e.end()]) < typEndLen)
329  {
330  // Use point.
331  localLoop.append(vertToEVert(e.end()));
332  localLoopWeights.append(-GREAT);
333 
334  useEnd = true;
335  }
336  }
337 
338  //
339  // Check cut of edge
340  //
341 
342  if (!useEnd && !useStart)
343  {
344  // Edge end points not close to plane so edge not close to
345  // plane. Cut edge.
346  scalar cutWeight;
347 
348  if (cutEdge(cutPlane, edgeI, cutWeight))
349  {
350  // Snap edge cut to vertex.
351  label cutVertI = snapToVert(snapTol_, edgeI, cutWeight);
352 
353  if (cutVertI == -1)
354  {
355  // Proper cut of edge
356  localLoop.append(edgeToEVert(edgeI));
357  localLoopWeights.append(cutWeight);
358  }
359  else
360  {
361  // Cut through edge end point. Might be duplicate
362  // since connected edge might also be snapped to same
363  // endpoint so only insert if unique.
364  label cut = vertToEVert(cutVertI);
365 
366  if (findIndex(localLoop, cut) == -1)
367  {
368  localLoop.append(vertToEVert(cutVertI));
369  localLoopWeights.append(-GREAT);
370  }
371  }
372  }
373  }
374  }
375 
376  if (localLoop.size() <= 2)
377  {
378  return false;
379  }
380 
381  localLoop.shrink();
382  localLoopWeights.shrink();
383 
384 
385  // Get points on loop and centre of loop
386  pointField loopPoints(localLoop.size());
387  point ctr(vector::zero);
388 
389  forAll(localLoop, i)
390  {
391  loopPoints[i] = coord(localLoop[i], localLoopWeights[i]);
392  ctr += loopPoints[i];
393  }
394  ctr /= localLoop.size();
395 
396 
397  // Get base vectors of coordinate system normal to refDir
398  vector e0, e1;
399  getBase(cutPlane.normal(), e0, e1);
400 
401 
402  // Get sorted angles from point on loop to centre of loop.
403  SortableList<scalar> sortedAngles(localLoop.size());
404 
405  forAll(sortedAngles, i)
406  {
407  vector toCtr(loopPoints[i] - ctr);
408  toCtr /= mag(toCtr);
409 
410  sortedAngles[i] = pseudoAngle(e0, e1, toCtr);
411  }
412  sortedAngles.sort();
413 
414  loop.setSize(localLoop.size());
415  loopWeights.setSize(loop.size());
416 
417 
418  // Fill loop and loopweights according to sorted angles
419 
420  const labelList& indices = sortedAngles.indices();
421 
422  forAll(indices, i)
423  {
424  loop[i] = localLoop[indices[i]];
425  loopWeights[i] = localLoopWeights[indices[i]];
426  }
427 
428 
429  // Check for cut edges along already cut vertices.
430  bool filterLoop = false;
431 
432  forAll(loop, i)
433  {
434  label cut = loop[i];
435 
436  if (isEdge(cut) && edgeEndsCut(loop, i))
437  {
438  filterLoop = true;
439  }
440  }
441 
442  if (filterLoop)
443  {
444  // Found edges in loop where both end points are also in loop. This
445  // is superfluous information so we can remove it.
446 
447  labelList filteredLoop(loop.size());
448  scalarField filteredLoopWeights(loopWeights.size());
449  label filterI = 0;
450 
451  forAll(loop, i)
452  {
453  label cut = loop[i];
454 
455  if (isEdge(cut) && edgeEndsCut(loop, i))
456  {
457  // Superfluous cut. Edge end points cut and edge itself as well.
458  }
459  else
460  {
461  filteredLoop[filterI] = loop[i];
462  filteredLoopWeights[filterI] = loopWeights[i];
463  filterI++;
464  }
465  }
466  filteredLoop.setSize(filterI);
467  filteredLoopWeights.setSize(filterI);
468 
469  loop.transfer(filteredLoop);
470  loopWeights.transfer(filteredLoopWeights);
471  }
472 
473  if (debug&2)
474  {
475  Pout<< "cell:" << cellI << endl;
476  forAll(loop, i)
477  {
478  Pout<< "At angle:" << sortedAngles[i] << endl
479  << " cut:";
480 
481  writeCut(Pout, loop[i], loopWeights[i]);
482 
483  Pout<< " coord:" << coord(loop[i], loopWeights[i]);
484 
485  Pout<< endl;
486  }
487  }
488 
489  return true;
490 }
491 
492 
493 // ************************************************************************* //
iterator end()
Return an iterator to end traversing the UList.
Definition: UListI.H:237
const pointField & points
static bool isEdge(const primitiveMesh &mesh, const label eVert)
Is eVert an edge?
Definition: edgeVertex.H:107
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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 ))
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & pointEdges() const
const vector & normal() const
Return plane normal.
Definition: plane.C:262
virtual ~geomCellLooper()
Destructor.
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 size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
const cellList & cells() const
Various functions to operate on Lists.
dynamicFvMesh & mesh
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const labelListList & cellEdges() const
Namespace for OpenFOAM.
static label edgeToEVert(const primitiveMesh &mesh, const label edgeI)
Convert edgeI to eVert.
Definition: edgeVertex.H:182
label n
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
A list that is sorted upon construction or when explicitly requested with the sort() method...
Definition: List.H:65
#define forAll(list, i)
Definition: UList.H:421
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:258
Macros for easy insertion into run-time selection tables.
static point coord(const primitiveMesh &, const label cut, const scalar weight)
Return coordinate of cut (uses weight if edgeCut)
Definition: edgeVertex.C:172
Geometric class that creates a 2D plane and can return the intersection point between a line and the ...
Definition: plane.H:60
Abstract base class. Concrete implementations know how to cut a cell (i.e. determine a loop around th...
Definition: cellLooper.H:69
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Ostream & writeCut(Ostream &os, const label cut, const scalar) const
Write cut description to Ostream.
Definition: edgeVertex.C:217
label end() const
Return end vertex label.
Definition: edgeI.H:92
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const Vector zero
Definition: Vector.H:80
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:223
const polyMesh & mesh() const
Definition: edgeVertex.H:98
3D tensor transformation operations.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
label start() const
Return start vertex label.
Definition: edgeI.H:81
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.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static label vertToEVert(const primitiveMesh &mesh, const label vertI)
Convert pointI to eVert.
Definition: edgeVertex.H:164
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
scalar distance(const point &p) const
Return distance from the given point to the plane.
Definition: plane.C:331
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116