thresholdCellFaces.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-2022 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 "thresholdCellFaces.H"
27 #include "emptyPolyPatch.H"
28 #include "processorPolyPatch.H"
29 #include "polygonTriangulate.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(thresholdCellFaces, 0);
36 }
37 
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  // Create a triangulation engine
92  polygonTriangulate triEngine;
93 
94  // internal faces only
95  for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
96  {
97  int side = 0;
98 
99  // check lowerThreshold
100  if (field[own[facei]] > lowerThreshold)
101  {
102  if (field[nei[facei]] < lowerThreshold)
103  {
104  side = +1;
105  }
106  }
107  else if (field[nei[facei]] > lowerThreshold)
108  {
109  side = -1;
110  }
111 
112  // check upperThreshold
113  if (field[own[facei]] < upperThreshold)
114  {
115  if (field[nei[facei]] > upperThreshold)
116  {
117  side = +1;
118  }
119  }
120  else if (field[nei[facei]] < upperThreshold)
121  {
122  side = -1;
123  }
124 
125 
126  if (side)
127  {
128  const face& f = origFaces[facei];
129 
130  forAll(f, fp)
131  {
132  if (oldToNewPoints[f[fp]] == -1)
133  {
134  oldToNewPoints[f[fp]] = nPoints++;
135  }
136  }
137 
138 
139  label cellId;
140  face surfFace;
141 
142  if (side > 0)
143  {
144  surfFace = f;
145  cellId = own[facei];
146  }
147  else
148  {
149  surfFace = f.reverseFace();
150  cellId = nei[facei];
151  }
152 
153 
154  if (triangulate)
155  {
156  triEngine.triangulate
157  (
158  UIndirectList<point>(origPoints, surfFace)
159  );
160  forAll(triEngine.triPoints(), i)
161  {
162  surfFaces.append(triEngine.triPoints(i, surfFace));
163  surfCells.append(cellId);
164  }
165  }
166  else
167  {
168  surfFaces.append(surfFace);
169  surfCells.append(cellId);
170  }
171  }
172  }
173 
174  surfZones[0].size() = surfFaces.size();
175 
176 
177  // nothing special for processor patches?
178  forAll(bMesh, patchi)
179  {
180  const polyPatch& p = bMesh[patchi];
181  surfZone& zone = surfZones[patchi+1];
182 
183  zone.start() = nFaces;
184 
185  if
186  (
187  isA<emptyPolyPatch>(p)
188  || (Pstream::parRun() && isA<processorPolyPatch>(p))
189  )
190  {
191  continue;
192  }
193 
194  label facei = p.start();
195 
196  // patch faces
197  forAll(p, localFacei)
198  {
199  if
200  (
201  field[own[facei]] > lowerThreshold
202  && field[own[facei]] < upperThreshold
203  )
204  {
205  const face& f = origFaces[facei];
206  forAll(f, fp)
207  {
208  if (oldToNewPoints[f[fp]] == -1)
209  {
210  oldToNewPoints[f[fp]] = nPoints++;
211  }
212  }
213 
214  label cellId = own[facei];
215 
216  if (triangulate)
217  {
218  triEngine.triangulate
219  (
220  UIndirectList<point>(origPoints, f)
221  );
222  forAll(triEngine.triPoints(), i)
223  {
224  surfFaces.append(triEngine.triPoints(i, f));
225  surfCells.append(cellId);
226  }
227  }
228  else
229  {
230  surfFaces.append(f);
231  surfCells.append(cellId);
232  }
233  }
234 
235  ++facei;
236  }
237 
238  zone.size() = surfFaces.size() - zone.start();
239  }
240 
241 
242  surfFaces.shrink();
243  surfCells.shrink();
244 
245  // renumber
246  forAll(surfFaces, facei)
247  {
248  inplaceRenumber(oldToNewPoints, surfFaces[facei]);
249  }
250 
251 
252  pointField surfPoints(nPoints);
253  nPoints = 0;
254  forAll(oldToNewPoints, pointi)
255  {
256  if (oldToNewPoints[pointi] >= 0)
257  {
258  surfPoints[oldToNewPoints[pointi]] = origPoints[pointi];
259  nPoints++;
260  }
261  }
262  surfPoints.setSize(nPoints);
263 
264  this->storedPoints().transfer(surfPoints);
265  this->storedFaces().transfer(surfFaces);
266  this->storedZones().transfer(surfZones);
267 
268  meshCells_.transfer(surfCells);
269 }
270 
271 
272 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
273 
275 (
276  const polyMesh& mesh,
277  const scalarField& field,
278  const scalar lowerThreshold,
279  const scalar upperThreshold,
280  const bool triangulate
281 )
282 :
283  mesh_(mesh)
284 {
285 
286  if (lowerThreshold > upperThreshold)
287  {
289  << lowerThreshold << " > " << upperThreshold << endl;
290  }
291 
292  calculate(field, lowerThreshold, upperThreshold, triangulate);
293 }
294 
295 
296 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
void inplaceRenumber(const labelUList &oldToNew, ListType &)
Inplace renumber the values of a list.
thresholdCellFaces(const polyMesh &, const scalarField &, const scalar lowerThreshold, const scalar upperThreshold, const bool triangulate=false)
Construct from mesh, field and threshold value.
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label calculate(const fvMesh &mesh, const labelHashSet &patchIDs, const scalar minFaceFraction, GeometricField< scalar, PatchField, GeoMesh > &distance)
Calculate distance data from patches.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
List< surfZone > surfZoneList
Definition: surfZoneList.H:45
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
label cellId
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
Namespace for OpenFOAM.