FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/DynamicList.H>
29 #include <OpenFOAM/SortableList.H>
30 #include <OpenFOAM/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 
621 const Foam::labelList& Foam::primitiveMesh::faceEdges(const label faceI) const
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 
670 const Foam::labelList& Foam::primitiveMesh::cellEdges(const label cellI) const
671 {
672  return cellEdges(cellI, labels_);
673 }
674 
675 
676 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
677 
678 // ************************ vim: set sw=4 sts=4 et: ************************ //