primitiveMeshEdges.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 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 "primitiveMesh.H"
27 #include "DynamicList.H"
28 #include "demandDrivenData.H"
29 #include "SortableList.H"
30 #include "ListOps.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 // Returns edgeI between two points.
35 Foam::label Foam::primitiveMesh::getEdge
36 (
37  List<DynamicList<label> >& pe,
38  DynamicList<edge>& es,
39 
40  const label pointI,
41  const label nextPointI
42 )
43 {
44  // Find connection between pointI and nextPointI
45  forAll(pe[pointI], ppI)
46  {
47  label eI = pe[pointI][ppI];
48 
49  const edge& e = es[eI];
50 
51  if (e.start() == nextPointI || e.end() == nextPointI)
52  {
53  return eI;
54  }
55  }
56 
57  // Make new edge.
58  label edgeI = es.size();
59  pe[pointI].append(edgeI);
60  pe[nextPointI].append(edgeI);
61  if (pointI < nextPointI)
62  {
63  es.append(edge(pointI, nextPointI));
64  }
65  else
66  {
67  es.append(edge(nextPointI, pointI));
68  }
69  return edgeI;
70 }
71 
72 
73 void Foam::primitiveMesh::calcEdges(const bool doFaceEdges) const
74 {
75  if (debug)
76  {
77  Pout<< "primitiveMesh::calcEdges(const bool) : "
78  << "calculating edges, pointEdges and optionally faceEdges"
79  << endl;
80  }
81 
82  // It is an error to attempt to recalculate edges
83  // if the pointer is already set
84  if ((edgesPtr_ || pePtr_) || (doFaceEdges && fePtr_))
85  {
86  FatalErrorIn("primitiveMesh::calcEdges(const bool) const")
87  << "edges or pointEdges or faceEdges already calculated"
88  << abort(FatalError);
89  }
90  else
91  {
92  // ALGORITHM:
93  // Go through the faces list. Search pointEdges for existing edge.
94  // If not found create edge and add to pointEdges.
95  // In second loop sort edges in order of increasing neighbouring
96  // point.
97  // This algorithm replaces the one using pointFaces which used more
98  // allocations but less memory and was on practical cases
99  // quite a bit slower.
100 
101  const faceList& fcs = faces();
102 
103  // Size up lists
104  // ~~~~~~~~~~~~~
105 
106  // Estimate pointEdges storage
107  List<DynamicList<label> > pe(nPoints());
108  forAll(pe, pointI)
109  {
110  pe[pointI].setCapacity(primitiveMesh::edgesPerPoint_);
111  }
112 
113  // Estimate edges storage
114  DynamicList<edge> es(pe.size()*primitiveMesh::edgesPerPoint_/2);
115 
116  // Estimate faceEdges storage
117  if (doFaceEdges)
118  {
119  fePtr_ = new labelListList(fcs.size());
120  labelListList& faceEdges = *fePtr_;
121  forAll(fcs, faceI)
122  {
123  faceEdges[faceI].setSize(fcs[faceI].size());
124  }
125  }
126 
127 
128  // Find consecutive face points in edge list
129  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
130 
131  // Edges on boundary faces
132  label nExtEdges = 0;
133  // Edges using no boundary point
134  nInternal0Edges_ = 0;
135  // Edges using 1 boundary point
136  label nInt1Edges = 0;
137  // Edges using two boundary points but not on boundary face:
138  // edges.size()-nExtEdges-nInternal0Edges_-nInt1Edges
139 
140  // Ordering is different if points are ordered.
141  if (nInternalPoints_ == -1)
142  {
143  // No ordering. No distinction between types.
144  forAll(fcs, faceI)
145  {
146  const face& f = fcs[faceI];
147 
148  forAll(f, fp)
149  {
150  label pointI = f[fp];
151  label nextPointI = f[f.fcIndex(fp)];
152 
153  label edgeI = getEdge(pe, es, pointI, nextPointI);
154 
155  if (doFaceEdges)
156  {
157  (*fePtr_)[faceI][fp] = edgeI;
158  }
159  }
160  }
161  // Assume all edges internal
162  nExtEdges = 0;
163  nInternal0Edges_ = es.size();
164  nInt1Edges = 0;
165  }
166  else
167  {
168  // 1. Do external faces first. This creates external edges.
169  for (label faceI = nInternalFaces_; faceI < fcs.size(); faceI++)
170  {
171  const face& f = fcs[faceI];
172 
173  forAll(f, fp)
174  {
175  label pointI = f[fp];
176  label nextPointI = f[f.fcIndex(fp)];
177 
178  label oldNEdges = es.size();
179  label edgeI = getEdge(pe, es, pointI, nextPointI);
180 
181  if (es.size() > oldNEdges)
182  {
183  nExtEdges++;
184  }
185  if (doFaceEdges)
186  {
187  (*fePtr_)[faceI][fp] = edgeI;
188  }
189  }
190  }
191 
192  // 2. Do internal faces. This creates internal edges.
193  for (label faceI = 0; faceI < nInternalFaces_; faceI++)
194  {
195  const face& f = fcs[faceI];
196 
197  forAll(f, fp)
198  {
199  label pointI = f[fp];
200  label nextPointI = f[f.fcIndex(fp)];
201 
202  label oldNEdges = es.size();
203  label edgeI = getEdge(pe, es, pointI, nextPointI);
204 
205  if (es.size() > oldNEdges)
206  {
207  if (pointI < nInternalPoints_)
208  {
209  if (nextPointI < nInternalPoints_)
210  {
211  nInternal0Edges_++;
212  }
213  else
214  {
215  nInt1Edges++;
216  }
217  }
218  else
219  {
220  if (nextPointI < nInternalPoints_)
221  {
222  nInt1Edges++;
223  }
224  else
225  {
226  // Internal edge with two points on boundary
227  }
228  }
229  }
230  if (doFaceEdges)
231  {
232  (*fePtr_)[faceI][fp] = edgeI;
233  }
234  }
235  }
236  }
237 
238 
239  // For unsorted meshes the edges will be in order of occurrence inside
240  // the faces. For sorted meshes the first nExtEdges will be the external
241  // edges.
242 
243  if (nInternalPoints_ != -1)
244  {
245  nInternalEdges_ = es.size()-nExtEdges;
246  nInternal1Edges_ = nInternal0Edges_+nInt1Edges;
247 
248  //Pout<< "Edge overview:" << nl
249  // << " total number of edges : " << es.size()
250  // << nl
251  // << " boundary edges : " << nExtEdges
252  // << nl
253  // << " edges using no boundary point : "
254  // << nInternal0Edges_
255  // << nl
256  // << " edges using one boundary points : " << nInt1Edges
257  // << nl
258  // << " edges using two boundary points : "
259  // << es.size()-nExtEdges-nInternal0Edges_-nInt1Edges << endl;
260  }
261 
262 
263  // Like faces sort edges in order of increasing neigbouring point.
264  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
265  // Automatically if points are sorted into internal and external points
266  // the edges will be sorted into
267  // - internal point to internal point
268  // - internal point to external point
269  // - external point to external point
270  // The problem is that the last one mixes external edges with internal
271  // edges with two points on the boundary.
272 
273 
274  // Map to sort into new upper-triangular order
275  labelList oldToNew(es.size());
276  if (debug)
277  {
278  oldToNew = -1;
279  }
280 
281  // start of edges with 0 boundary points
282  label internal0EdgeI = 0;
283 
284  // start of edges with 1 boundary points
285  label internal1EdgeI = nInternal0Edges_;
286 
287  // start of edges with 2 boundary points
288  label internal2EdgeI = nInternal1Edges_;
289 
290  // start of external edges
291  label externalEdgeI = nInternalEdges_;
292 
293 
294  // To sort neighbouring points in increasing order. Defined outside
295  // for optimisation reasons: if all connectivity size same will need
296  // no reallocations
297  SortableList<label> nbrPoints(primitiveMesh::edgesPerPoint_);
298 
299  forAll(pe, pointI)
300  {
301  const DynamicList<label>& pEdges = pe[pointI];
302 
303  nbrPoints.setSize(pEdges.size());
304 
305  forAll(pEdges, i)
306  {
307  const edge& e = es[pEdges[i]];
308 
309  label nbrPointI = e.otherVertex(pointI);
310 
311  if (nbrPointI < pointI)
312  {
313  nbrPoints[i] = -1;
314  }
315  else
316  {
317  nbrPoints[i] = nbrPointI;
318  }
319  }
320  nbrPoints.sort();
321 
322 
323  if (nInternalPoints_ == -1)
324  {
325  // Sort all upper-triangular
326  forAll(nbrPoints, i)
327  {
328  if (nbrPoints[i] != -1)
329  {
330  label edgeI = pEdges[nbrPoints.indices()[i]];
331 
332  oldToNew[edgeI] = internal0EdgeI++;
333  }
334  }
335  }
336  else
337  {
338  if (pointI < nInternalPoints_)
339  {
340  forAll(nbrPoints, i)
341  {
342  label nbrPointI = nbrPoints[i];
343 
344  label edgeI = pEdges[nbrPoints.indices()[i]];
345 
346  if (nbrPointI != -1)
347  {
348  if (edgeI < nExtEdges)
349  {
350  // External edge
351  oldToNew[edgeI] = externalEdgeI++;
352  }
353  else if (nbrPointI < nInternalPoints_)
354  {
355  // Both points inside
356  oldToNew[edgeI] = internal0EdgeI++;
357  }
358  else
359  {
360  // One points inside, one outside
361  oldToNew[edgeI] = internal1EdgeI++;
362  }
363  }
364  }
365  }
366  else
367  {
368  forAll(nbrPoints, i)
369  {
370  label nbrPointI = nbrPoints[i];
371 
372  label edgeI = pEdges[nbrPoints.indices()[i]];
373 
374  if (nbrPointI != -1)
375  {
376  if (edgeI < nExtEdges)
377  {
378  // External edge
379  oldToNew[edgeI] = externalEdgeI++;
380  }
381  else if (nbrPointI < nInternalPoints_)
382  {
383  // Not possible!
384  FatalErrorIn("primitiveMesh::calcEdges(..)")
385  << abort(FatalError);
386  }
387  else
388  {
389  // Both points outside
390  oldToNew[edgeI] = internal2EdgeI++;
391  }
392  }
393  }
394  }
395  }
396  }
397 
398 
399  if (debug)
400  {
401  label edgeI = findIndex(oldToNew, -1);
402 
403  if (edgeI != -1)
404  {
405  const edge& e = es[edgeI];
406 
407  FatalErrorIn("primitiveMesh::calcEdges(..)")
408  << "Did not sort edge " << edgeI << " points:" << e
409  << " coords:" << points()[e[0]] << points()[e[1]]
410  << endl
411  << "Current buckets:" << endl
412  << " internal0EdgeI:" << internal0EdgeI << endl
413  << " internal1EdgeI:" << internal1EdgeI << endl
414  << " internal2EdgeI:" << internal2EdgeI << endl
415  << " externalEdgeI:" << externalEdgeI << endl
416  << abort(FatalError);
417  }
418  }
419 
420 
421 
422  // Renumber and transfer.
423 
424  // Edges
425  edgesPtr_ = new edgeList(es.size());
426  edgeList& edges = *edgesPtr_;
427  forAll(es, edgeI)
428  {
429  edges[oldToNew[edgeI]] = es[edgeI];
430  }
431 
432  // pointEdges
433  pePtr_ = new labelListList(nPoints());
434  labelListList& pointEdges = *pePtr_;
435  forAll(pe, pointI)
436  {
437  DynamicList<label>& pEdges = pe[pointI];
438  pEdges.shrink();
439  inplaceRenumber(oldToNew, pEdges);
440  pointEdges[pointI].transfer(pEdges);
441  Foam::sort(pointEdges[pointI]);
442  }
443 
444  // faceEdges
445  if (doFaceEdges)
446  {
447  labelListList& faceEdges = *fePtr_;
448  forAll(faceEdges, faceI)
449  {
450  inplaceRenumber(oldToNew, faceEdges[faceI]);
451  }
452  }
453  }
454 }
455 
456 
457 Foam::label Foam::primitiveMesh::findFirstCommonElementFromSortedLists
458 (
459  const labelList& list1,
460  const labelList& list2
461 )
462 {
463  label result = -1;
464 
465  labelList::const_iterator iter1 = list1.begin();
466  labelList::const_iterator iter2 = list2.begin();
467 
468  while (iter1 != list1.end() && iter2 != list2.end())
469  {
470  if (*iter1 < *iter2)
471  {
472  ++iter1;
473  }
474  else if (*iter1 > *iter2)
475  {
476  ++iter2;
477  }
478  else
479  {
480  result = *iter1;
481  break;
482  }
483  }
484  if (result == -1)
485  {
487  (
488  "primitiveMesh::findFirstCommonElementFromSortedLists"
489  "(const labelList&, const labelList&)"
490  ) << "No common elements in lists " << list1 << " and " << list2
491  << abort(FatalError);
492  }
493  return result;
494 }
495 
496 
497 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
498 
500 {
501  if (!edgesPtr_)
502  {
503  //calcEdges(true);
504  calcEdges(false);
505  }
506 
507  return *edgesPtr_;
508 }
509 
511 {
512  if (!pePtr_)
513  {
514  //calcEdges(true);
515  calcEdges(false);
516  }
517 
518  return *pePtr_;
519 }
520 
521 
523 {
524  if (!fePtr_)
525  {
526  if (debug)
527  {
528  Pout<< "primitiveMesh::faceEdges() : "
529  << "calculating faceEdges" << endl;
530  }
531 
532  //calcEdges(true);
533  const faceList& fcs = faces();
534  const labelListList& pe = pointEdges();
535  const edgeList& es = edges();
536 
537  fePtr_ = new labelListList(fcs.size());
538  labelListList& faceEdges = *fePtr_;
539 
540  forAll(fcs, faceI)
541  {
542  const face& f = fcs[faceI];
543 
544  labelList& fEdges = faceEdges[faceI];
545  fEdges.setSize(f.size());
546 
547  forAll(f, fp)
548  {
549  label pointI = f[fp];
550  label nextPointI = f[f.fcIndex(fp)];
551 
552  // Find edge between pointI, nextPontI
553  const labelList& pEdges = pe[pointI];
554 
555  forAll(pEdges, i)
556  {
557  label edgeI = pEdges[i];
558 
559  if (es[edgeI].otherVertex(pointI) == nextPointI)
560  {
561  fEdges[fp] = edgeI;
562  break;
563  }
564  }
565  }
566  }
567  }
568 
569  return *fePtr_;
570 }
571 
572 
573 void Foam::primitiveMesh::clearOutEdges()
574 {
575  deleteDemandDrivenData(edgesPtr_);
576  deleteDemandDrivenData(pePtr_);
577  deleteDemandDrivenData(fePtr_);
578  labels_.clear();
579  labelSet_.clear();
580 }
581 
582 
584 (
585  const label faceI,
586  DynamicList<label>& storage
587 ) const
588 {
589  if (hasFaceEdges())
590  {
591  return faceEdges()[faceI];
592  }
593  else
594  {
595  const labelListList& pointEs = pointEdges();
596  const face& f = faces()[faceI];
597 
598  storage.clear();
599  if (f.size() > storage.capacity())
600  {
601  storage.setCapacity(f.size());
602  }
603 
604  forAll(f, fp)
605  {
606  storage.append
607  (
608  findFirstCommonElementFromSortedLists
609  (
610  pointEs[f[fp]],
611  pointEs[f.nextLabel(fp)]
612  )
613  );
614  }
615 
616  return storage;
617  }
618 }
619 
620 
622 {
623  return faceEdges(faceI, labels_);
624 }
625 
626 
628 (
629  const label cellI,
630  DynamicList<label>& storage
631 ) const
632 {
633  if (hasCellEdges())
634  {
635  return cellEdges()[cellI];
636  }
637  else
638  {
639  const labelList& cFaces = cells()[cellI];
640 
641  labelSet_.clear();
642 
643  forAll(cFaces, i)
644  {
645  const labelList& fe = faceEdges(cFaces[i]);
646 
647  forAll(fe, feI)
648  {
649  labelSet_.insert(fe[feI]);
650  }
651  }
652 
653  storage.clear();
654 
655  if (labelSet_.size() > storage.capacity())
656  {
657  storage.setCapacity(labelSet_.size());
658  }
659 
660  forAllConstIter(labelHashSet, labelSet_, iter)
661  {
662  storage.append(iter.key());
663  }
664 
665  return storage;
666  }
667 }
668 
669 
671 {
672  return cellEdges(cellI, labels_);
673 }
674 
675 
676 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
677 
678 // ************************************************************************* //
virtual const pointField & points() const =0
Return mesh points.
void setCapacity(const label)
Alter the size of the underlying storage.
Definition: DynamicListI.H:109
void sort(UList< T > &)
Definition: UList.C:107
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
const labelListList & pointEdges() const
bool hasCellEdges() const
void deleteDemandDrivenData(DataPtr &dataPtr)
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:242
label capacity() const
Size of the underlying storage.
Definition: DynamicListI.H:100
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
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const cellList & cells() const
Various functions to operate on Lists.
List< face > faceList
Definition: faceListFwd.H:43
virtual const faceList & faces() const =0
Return faces.
label nextLabel(const label i) const
Next vertex on face.
Definition: faceI.H:117
const labelListList & faceEdges() const
const labelListList & cellEdges() const
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
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
void clear()
Clear all entries from table.
Definition: HashTable.C:473
#define forAll(list, i)
Definition: UList.H:421
Template functions to aid in the implementation of demand driven data.
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const unsigned edgesPerPoint_
Estimated number of edges per point.
const T * const_iterator
Random access iterator for traversing UList.
Definition: UList.H:262
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
bool hasFaceEdges() const
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
label nPoints() const
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116