calcPointCells.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-2021 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 Description
25  calculate point cells - ie, the cells attached to each point
26 
27  - remove unused points, adjust pointCells and cellFaces accordingly
28 \*---------------------------------------------------------------------------*/
29 
30 #include "meshReader.H"
31 
32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
33 
34 void Foam::meshReader::calcPointCells() const
35 {
36  static const label UNIT_POINT_CELLS = 12;
37 
38  if (pointCellsPtr_)
39  {
41  << "pointCells already calculated"
42  << abort(FatalError);
43  }
44 
46 
47  pointCellsPtr_ = new labelListList(nPoints);
48  labelListList& ptCells = *pointCellsPtr_;
49 
50  forAll(ptCells, i)
51  {
52  ptCells[i].setSize(UNIT_POINT_CELLS);
53  }
54 
55  // Initialise the list of labels which will hold the count of the
56  // actual number of cells per point during the analysis
57  labelList cellCount(nPoints, 0);
58 
59  // Note. Unlike the standard point-cell algorithm, which asks the cell for
60  // the supporting point labels, we need to work based on the cell faces.
61  // This is because some of the faces do not come from the cell shape.
62  // It is also advantageous to remove duplicates from the point-cell
63  // addressing, because this removes a lot of waste later.
64 
65  faceListList& cFaces = cellFaces();
66 
67  // For each cell
68  forAll(cFaces, celli)
69  {
70  const faceList& faces = cFaces[celli];
71 
72  forAll(faces, i)
73  {
74  // For each vertex
75  const labelList& labels = faces[i];
76 
77  forAll(labels, j)
78  {
79  // Set working point label
80  label curPoint = labels[j];
81  labelList& curPointCells = ptCells[curPoint];
82  label curCount = cellCount[curPoint];
83 
84  // check if the cell has been added before
85  bool found = false;
86 
87  for (label f = 0; f < curCount; f++)
88  {
89  if (curPointCells[f] == celli)
90  {
91  found = true;
92  break;
93  }
94  }
95 
96  if (!found)
97  {
98  // If the list of pointCells is not big enough, double it
99  if (curPointCells.size() <= curCount)
100  {
101  curPointCells.setSize(curPointCells.size()*2);
102  }
103 
104  // Enter the cell label in the point's cell list
105  curPointCells[curCount] = celli;
106 
107  // Increment the cell count for the point addressed
108  cellCount[curPoint]++;
109  }
110  }
111  }
112  }
113 
114  // report and remove unused points
115  // - adjust points, pointCells, and cellFaces accordingly
116  label pointi = 0;
117  labelList oldToNew(nPoints, -1);
118 
119  forAll(ptCells, i)
120  {
121  ptCells[i].setSize(cellCount[i]);
122  if (cellCount[i] > 0)
123  {
124  oldToNew[i] = pointi++;
125  }
126  }
127 
128  // report unused points
129  if (nPoints > pointi)
130  {
131  Info<< "removing " << (nPoints - pointi) << " unused points" << endl;
132 
133  nPoints = pointi;
134 
135  // adjust points and truncate - bend const-ness
136  pointField& adjustedPoints = const_cast<pointField&>(points_);
137 
138  inplaceReorder(oldToNew, adjustedPoints);
139  adjustedPoints.setSize(nPoints);
140 
141  // adjust pointCells and truncate
142  inplaceReorder(oldToNew, ptCells);
143  ptCells.setSize(nPoints);
144 
145  // adjust cellFaces - this could be faster
146  // For each cell
147  forAll(cFaces, celli)
148  {
149  faceList& faces = cFaces[celli];
150 
151  // For each face
152  forAll(faces, i)
153  {
154  inplaceRenumber(oldToNew, faces[i]);
155  }
156  }
157  }
158 }
159 
160 
161 const Foam::labelListList& Foam::meshReader::pointCells() const
162 {
163  if (!pointCellsPtr_)
164  {
165  calcPointCells();
166  }
167 
168  return *pointCellsPtr_;
169 }
170 
171 
172 // ************************************************************************* //
bool found
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
pointField points_
Points supporting the mesh.
Definition: meshReader.H:249
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
label nPoints
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
errorManip< error > abort(error &err)
Definition: errorManip.H:131
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
List< faceList > faceListList
Definition: faceListFwd.H:43
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
error FatalError
List< face > faceList
Definition: faceListFwd.H:41
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
labelList f(nPoints)