polyMeshTools.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2012-2016 OpenFOAM Foundation
9  Copyright (C) 2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "polyMeshTools.H"
30 #include "syncTools.H"
31 #include "pyramidPointFaceRef.H"
32 #include "primitiveMeshTools.H"
33 
34 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
35 
37 (
38  const polyMesh& mesh,
39  const vectorField& areas,
40  const vectorField& cc
41 )
42 {
43  const labelList& own = mesh.faceOwner();
44  const labelList& nei = mesh.faceNeighbour();
45  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
46 
47  tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
48  scalarField& ortho = tortho.ref();
49 
50  // Internal faces
51  forAll(nei, facei)
52  {
53  ortho[facei] = primitiveMeshTools::faceOrthogonality
54  (
55  cc[own[facei]],
56  cc[nei[facei]],
57  areas[facei]
58  );
59  }
60 
61 
62  // Coupled faces
63 
64  pointField neighbourCc;
65  syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc);
66 
67  forAll(pbm, patchi)
68  {
69  const polyPatch& pp = pbm[patchi];
70  if (pp.coupled())
71  {
72  forAll(pp, i)
73  {
74  label facei = pp.start() + i;
75  label bFacei = facei - mesh.nInternalFaces();
76 
77  ortho[facei] = primitiveMeshTools::faceOrthogonality
78  (
79  cc[own[facei]],
80  neighbourCc[bFacei],
81  areas[facei]
82  );
83  }
84  }
85  }
86 
87  return tortho;
88 }
89 
90 
92 (
93  const polyMesh& mesh,
94  const pointField& p,
95  const vectorField& fCtrs,
96  const vectorField& fAreas,
97  const vectorField& cellCtrs
98 )
99 {
100  const labelList& own = mesh.faceOwner();
101  const labelList& nei = mesh.faceNeighbour();
102  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
103 
104  tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
105  scalarField& skew = tskew.ref();
106 
107  forAll(nei, facei)
108  {
109  skew[facei] = primitiveMeshTools::faceSkewness
110  (
111  mesh,
112  p,
113  fCtrs,
114  fAreas,
115 
116  facei,
117  cellCtrs[own[facei]],
118  cellCtrs[nei[facei]]
119  );
120  }
121 
122 
123  // Boundary faces: consider them to have only skewness error.
124  // (i.e. treat as if mirror cell on other side)
125 
126  pointField neighbourCc;
127  syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc);
128 
129  forAll(pbm, patchi)
130  {
131  const polyPatch& pp = pbm[patchi];
132  if (pp.coupled())
133  {
134  forAll(pp, i)
135  {
136  label facei = pp.start() + i;
137  label bFacei = facei - mesh.nInternalFaces();
138 
139  skew[facei] = primitiveMeshTools::faceSkewness
140  (
141  mesh,
142  p,
143  fCtrs,
144  fAreas,
145 
146  facei,
147  cellCtrs[own[facei]],
148  neighbourCc[bFacei]
149  );
150  }
151  }
152  else
153  {
154  forAll(pp, i)
155  {
156  label facei = pp.start() + i;
157 
158  skew[facei] = primitiveMeshTools::boundaryFaceSkewness
159  (
160  mesh,
161  p,
162  fCtrs,
163  fAreas,
164 
165  facei,
166  cellCtrs[own[facei]]
167  );
168  }
169  }
170  }
171 
172  return tskew;
173 }
174 
175 
177 (
178  const polyMesh& mesh,
179  const vectorField& fCtrs,
180  const vectorField& fAreas,
181  const vectorField& cellCtrs
182 )
183 {
184  const labelList& own = mesh.faceOwner();
185  const labelList& nei = mesh.faceNeighbour();
186  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
187 
188  tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
189  scalarField& weight = tweight.ref();
190 
191  // Internal faces
192  forAll(nei, facei)
193  {
194  const point& fc = fCtrs[facei];
195  const vector& fa = fAreas[facei];
196 
197  scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
198  scalar dNei = mag(fa & (cellCtrs[nei[facei]]-fc));
199 
200  weight[facei] = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
201  }
202 
203 
204  // Coupled faces
205 
206  pointField neiCc;
207  syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
208 
209  forAll(pbm, patchi)
210  {
211  const polyPatch& pp = pbm[patchi];
212  if (pp.coupled())
213  {
214  forAll(pp, i)
215  {
216  label facei = pp.start() + i;
217  label bFacei = facei - mesh.nInternalFaces();
218 
219  const point& fc = fCtrs[facei];
220  const vector& fa = fAreas[facei];
221 
222  scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
223  scalar dNei = mag(fa & (neiCc[bFacei]-fc));
224 
225  weight[facei] = min(dNei,dOwn)/(dNei+dOwn+VSMALL);
226  }
227  }
228  }
229 
230  return tweight;
231 }
232 
233 
235 (
236  const polyMesh& mesh,
237  const scalarField& vol
238 )
239 {
240  const labelList& own = mesh.faceOwner();
241  const labelList& nei = mesh.faceNeighbour();
242  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
243 
244  tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
245  scalarField& ratio = tratio.ref();
246 
247  // Internal faces
248  forAll(nei, facei)
249  {
250  scalar volOwn = vol[own[facei]];
251  scalar volNei = vol[nei[facei]];
252 
253  ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+VSMALL);
254  }
255 
256 
257  // Coupled faces
258 
259  scalarField neiVol;
260  syncTools::swapBoundaryCellList(mesh, vol, neiVol);
261 
262  forAll(pbm, patchi)
263  {
264  const polyPatch& pp = pbm[patchi];
265  if (pp.coupled())
266  {
267  forAll(pp, i)
268  {
269  label facei = pp.start() + i;
270  label bFacei = facei - mesh.nInternalFaces();
271 
272  scalar volOwn = vol[own[facei]];
273  scalar volNei = neiVol[bFacei];
274 
275  ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+VSMALL);
276  }
277  }
278  }
279 
280  return tratio;
281 }
282 
283 
285 (
286  const polyMesh::readUpdateState& state0,
287  const polyMesh::readUpdateState& state1
288 )
289 {
290  if
291  (
292  (
293  state0 == polyMesh::UNCHANGED
294  && state1 != polyMesh::UNCHANGED
295  )
296  || (
297  state0 == polyMesh::POINTS_MOVED
298  && (state1 != polyMesh::UNCHANGED && state1 != polyMesh::POINTS_MOVED)
299  )
300  || (
301  state0 == polyMesh::TOPO_CHANGE
302  && state1 == polyMesh::TOPO_PATCH_CHANGE
303  )
304  )
305  {
306  return state1;
307  }
308  else
309  {
310  return state0;
311  }
312 }
313 
314 
315 // ************************************************************************* //
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
p
volScalarField & p
Definition: createFieldRefs.H:8
Foam::skew
dimensionedTensor skew(const dimensionedTensor &dt)
Definition: dimensionedTensor.C:138
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Definition: polyBoundaryMesh.H:63
primitiveMeshTools.H
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::polyPatch::coupled
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:377
Foam::polyMeshTools::volRatio
static tmp< scalarField > volRatio(const polyMesh &mesh, const scalarField &vol)
Generate volume ratio field.
Definition: polyMeshTools.C:235
Foam::polyMeshTools::faceSkewness
static tmp< scalarField > faceSkewness(const polyMesh &mesh, const pointField &points, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate skewness field.
Definition: polyMeshTools.C:92
syncTools.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:77
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
Foam::Field< vector >
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:68
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::polyPatch::start
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:361
Foam::polyMesh::readUpdateState
readUpdateState
Enumeration defining the state of the mesh after a read update.
Definition: polyMesh.H:90
polyMeshTools.H
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::polyMeshTools::combine
static polyMesh::readUpdateState combine(const polyMesh::readUpdateState &state0, const polyMesh::readUpdateState &state1)
Combine readUpdateState. topo change trumps geom-only.
Definition: polyMeshTools.C:285
Foam::polyMeshTools::faceWeights
static tmp< scalarField > faceWeights(const polyMesh &mesh, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate interpolation factors field.
Definition: polyMeshTools.C:177
pyramidPointFaceRef.H
Foam::polyMeshTools::faceOrthogonality
static tmp< scalarField > faceOrthogonality(const polyMesh &mesh, const vectorField &fAreas, const vectorField &cellCtrs)
Generate orthogonality field. (1 for fully orthogonal, < 1 for.
Definition: polyMeshTools.C:37