FreeFOAM The Cross-Platform CFD Toolkit
thresholdCellFaces.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 "thresholdCellFaces.H"
27 
28 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/DynamicList.H>
30 
33 
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 void Foam::thresholdCellFaces::calculate
42 (
43  const scalarField& field,
44  const scalar lowerThreshold,
45  const scalar upperThreshold,
46  const bool triangulate
47 )
48 {
49  const labelList& own = mesh_.faceOwner();
50  const labelList& nei = mesh_.faceNeighbour();
51 
52  const faceList& origFaces = mesh_.faces();
53  const pointField& origPoints = mesh_.points();
54 
55  const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
56 
57 
58  surfZoneList surfZones(bMesh.size()+1);
59 
60  surfZones[0] = surfZone
61  (
62  "internalMesh",
63  0, // size
64  0, // start
65  0 // index
66  );
67 
68  forAll(bMesh, patchI)
69  {
70  surfZones[patchI+1] = surfZone
71  (
72  bMesh[patchI].name(),
73  0, // size
74  0, // start
75  patchI+1 // index
76  );
77  }
78 
79 
80  label nFaces = 0;
81  label nPoints = 0;
82 
83 
84  meshCells_.clear();
85 
86  DynamicList<face> surfFaces(0.5 * mesh_.nFaces());
87  DynamicList<label> surfCells(surfFaces.size());
88 
89  labelList oldToNewPoints(origPoints.size(), -1);
90 
91 
92  // internal faces only
93  for (label faceI = 0; faceI < mesh_.nInternalFaces(); ++faceI)
94  {
95  int side = 0;
96 
97  // check lowerThreshold
98  if (field[own[faceI]] > lowerThreshold)
99  {
100  if (field[nei[faceI]] < lowerThreshold)
101  {
102  side = +1;
103  }
104  }
105  else if (field[nei[faceI]] > lowerThreshold)
106  {
107  side = -1;
108  }
109 
110  // check upperThreshold
111  if (field[own[faceI]] < upperThreshold)
112  {
113  if (field[nei[faceI]] > upperThreshold)
114  {
115  side = +1;
116  }
117  }
118  else if (field[nei[faceI]] < upperThreshold)
119  {
120  side = -1;
121  }
122 
123 
124  if (side)
125  {
126  const face& f = origFaces[faceI];
127 
128  forAll(f, fp)
129  {
130  if (oldToNewPoints[f[fp]] == -1)
131  {
132  oldToNewPoints[f[fp]] = nPoints++;
133  }
134  }
135 
136 
137  label cellId;
138  face surfFace;
139 
140  if (side > 0)
141  {
142  surfFace = f;
143  cellId = own[faceI];
144  }
145  else
146  {
147  surfFace = f.reverseFace();
148  cellId = nei[faceI];
149  }
150 
151 
152  if (triangulate)
153  {
154  label count = surfFace.triangles(origPoints, surfFaces);
155  while (count-- > 0)
156  {
157  surfCells.append(cellId);
158  }
159  }
160  else
161  {
162  surfFaces.append(surfFace);
163  surfCells.append(cellId);
164  }
165  }
166  }
167 
168  surfZones[0].size() = surfFaces.size();
169 
170 
171  // nothing special for processor patches?
172  forAll(bMesh, patchI)
173  {
174  const polyPatch& p = bMesh[patchI];
175  surfZone& zone = surfZones[patchI+1];
176 
177  zone.start() = nFaces;
178 
179  if
180  (
181  isA<emptyPolyPatch>(p)
182  || (Pstream::parRun() && isA<processorPolyPatch>(p))
183  )
184  {
185  continue;
186  }
187 
188  label faceI = p.start();
189 
190  // patch faces
191  forAll(p, localFaceI)
192  {
193  if
194  (
195  field[own[faceI]] > lowerThreshold
196  && field[own[faceI]] < upperThreshold
197  )
198  {
199  const face& f = origFaces[faceI];
200  forAll(f, fp)
201  {
202  if (oldToNewPoints[f[fp]] == -1)
203  {
204  oldToNewPoints[f[fp]] = nPoints++;
205  }
206  }
207 
208  label cellId = own[faceI];
209 
210  if (triangulate)
211  {
212  label count = f.triangles(origPoints, surfFaces);
213  while (count-- > 0)
214  {
215  surfCells.append(cellId);
216  }
217  }
218  else
219  {
220  surfFaces.append(f);
221  surfCells.append(cellId);
222  }
223  }
224 
225  ++faceI;
226  }
227 
228  zone.size() = surfFaces.size() - zone.start();
229  }
230 
231 
232  surfFaces.shrink();
233  surfCells.shrink();
234 
235  // renumber
236  forAll(surfFaces, faceI)
237  {
238  inplaceRenumber(oldToNewPoints, surfFaces[faceI]);
239  }
240 
241 
242  pointField surfPoints(nPoints);
243  nPoints = 0;
244  forAll(oldToNewPoints, pointI)
245  {
246  if (oldToNewPoints[pointI] >= 0)
247  {
248  surfPoints[oldToNewPoints[pointI]] = origPoints[pointI];
249  nPoints++;
250  }
251  }
252  surfPoints.setSize(nPoints);
253 
254  this->storedPoints().transfer(surfPoints);
255  this->storedFaces().transfer(surfFaces);
256  this->storedZones().transfer(surfZones);
257 
258  meshCells_.transfer(surfCells);
259 }
260 
261 
262 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
263 
265 (
266  const polyMesh& mesh,
267  const scalarField& field,
268  const scalar lowerThreshold,
269  const scalar upperThreshold,
270  const bool triangulate
271 )
272 :
273  mesh_(mesh)
274 {
275 
276  if (lowerThreshold > upperThreshold)
277  {
278  WarningIn("thresholdCellFaces::thresholdCellFaces(...)")
279  << "lower > upper limit! "
280  << lowerThreshold << " > " << upperThreshold << endl;
281  }
282 
283  calculate(field, lowerThreshold, upperThreshold, triangulate);
284 }
285 
286 
287 // ************************ vim: set sw=4 sts=4 et: ************************ //