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