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