polyMeshTools.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) 2012-2018 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 "polyMeshTools.H"
27 #include "syncTools.H"
28 #include "pyramidPointFaceRef.H"
29 #include "primitiveMeshTools.H"
30 #include "polyMeshTools.H"
31 
32 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
33 
35 (
36  const polyMesh& mesh,
37  const vectorField& areas,
38  const vectorField& cc
39 )
40 {
41  const labelList& own = mesh.faceOwner();
42  const labelList& nei = mesh.faceNeighbour();
43  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
44 
45  tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
46  scalarField& ortho = tortho.ref();
47 
48  // Internal faces
49  forAll(nei, facei)
50  {
51  ortho[facei] = primitiveMeshTools::faceOrthogonality
52  (
53  cc[own[facei]],
54  cc[nei[facei]],
55  areas[facei]
56  );
57  }
58 
59 
60  // Coupled faces
61 
62  pointField neighbourCc;
63  syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc);
64 
65  forAll(pbm, patchi)
66  {
67  const polyPatch& pp = pbm[patchi];
68  if (pp.coupled())
69  {
70  forAll(pp, i)
71  {
72  label facei = pp.start() + i;
73  label bFacei = facei - mesh.nInternalFaces();
74 
75  ortho[facei] = primitiveMeshTools::faceOrthogonality
76  (
77  cc[own[facei]],
78  neighbourCc[bFacei],
79  areas[facei]
80  );
81  }
82  }
83  }
84 
85  return tortho;
86 }
87 
88 
90 (
91  const polyMesh& mesh,
92  const pointField& p,
93  const vectorField& fCtrs,
94  const vectorField& fAreas,
95  const vectorField& cellCtrs
96 )
97 {
98  const labelList& own = mesh.faceOwner();
99  const labelList& nei = mesh.faceNeighbour();
100  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
101 
102  tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
103  scalarField& skew = tskew.ref();
104 
105  forAll(nei, facei)
106  {
107  skew[facei] = primitiveMeshTools::faceSkewness
108  (
109  mesh,
110  p,
111  fCtrs,
112  fAreas,
113 
114  facei,
115  cellCtrs[own[facei]],
116  cellCtrs[nei[facei]]
117  );
118  }
119 
120 
121  // Boundary faces: consider them to have only skewness error.
122  // (i.e. treat as if mirror cell on other side)
123 
124  pointField neighbourCc;
125  syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc);
126 
127  forAll(pbm, patchi)
128  {
129  const polyPatch& pp = pbm[patchi];
130  if (pp.coupled())
131  {
132  forAll(pp, i)
133  {
134  label facei = pp.start() + i;
135  label bFacei = facei - mesh.nInternalFaces();
136 
137  skew[facei] = primitiveMeshTools::faceSkewness
138  (
139  mesh,
140  p,
141  fCtrs,
142  fAreas,
143 
144  facei,
145  cellCtrs[own[facei]],
146  neighbourCc[bFacei]
147  );
148  }
149  }
150  else
151  {
152  forAll(pp, i)
153  {
154  label facei = pp.start() + i;
155 
156  skew[facei] = primitiveMeshTools::boundaryFaceSkewness
157  (
158  mesh,
159  p,
160  fCtrs,
161  fAreas,
162 
163  facei,
164  cellCtrs[own[facei]]
165  );
166  }
167  }
168  }
169 
170  return tskew;
171 }
172 
173 
175 (
176  const polyMesh& mesh,
177  const vectorField& fCtrs,
178  const vectorField& fAreas,
179  const vectorField& cellCtrs
180 )
181 {
182  const labelList& own = mesh.faceOwner();
183  const labelList& nei = mesh.faceNeighbour();
184  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
185 
186  tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
187  scalarField& weight = tweight.ref();
188 
189  // Internal faces
190  forAll(nei, facei)
191  {
192  const point& fc = fCtrs[facei];
193  const vector& fa = fAreas[facei];
194 
195  scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
196  scalar dNei = mag(fa & (cellCtrs[nei[facei]]-fc));
197 
198  weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
199  }
200 
201 
202  // Coupled faces
203 
204  pointField neiCc;
205  syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
206 
207  forAll(pbm, patchi)
208  {
209  const polyPatch& pp = pbm[patchi];
210  if (pp.coupled())
211  {
212  forAll(pp, i)
213  {
214  label facei = pp.start() + i;
215  label bFacei = facei - mesh.nInternalFaces();
216 
217  const point& fc = fCtrs[facei];
218  const vector& fa = fAreas[facei];
219 
220  scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
221  scalar dNei = mag(fa & (neiCc[bFacei]-fc));
222 
223  weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
224  }
225  }
226  }
227 
228  return tweight;
229 }
230 
231 
233 (
234  const polyMesh& mesh,
235  const scalarField& vol
236 )
237 {
238  const labelList& own = mesh.faceOwner();
239  const labelList& nei = mesh.faceNeighbour();
240  const polyBoundaryMesh& pbm = mesh.boundaryMesh();
241 
242  tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
243  scalarField& ratio = tratio.ref();
244 
245  // Internal faces
246  forAll(nei, facei)
247  {
248  scalar volOwn = vol[own[facei]];
249  scalar volNei = vol[nei[facei]];
250 
251  ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
252  }
253 
254 
255  // Coupled faces
256 
257  scalarField neiVol;
258  syncTools::swapBoundaryCellList(mesh, vol, neiVol);
259 
260  forAll(pbm, patchi)
261  {
262  const polyPatch& pp = pbm[patchi];
263  if (pp.coupled())
264  {
265  forAll(pp, i)
266  {
267  label facei = pp.start() + i;
268  label bFacei = facei - mesh.nInternalFaces();
269 
270  scalar volOwn = vol[own[facei]];
271  scalar volNei = neiVol[bFacei];
272 
273  ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
274  }
275  }
276  }
277 
278  return tratio;
279 }
280 
281 
282 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:453
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1243
dimensionedTensor skew(const dimensionedTensor &dt)
label nFaces() const
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
static tmp< scalarField > faceWeights(const polyMesh &mesh, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate interpolation factors field.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:316
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
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:90
static tmp< scalarField > volRatio(const polyMesh &mesh, const scalarField &vol)
Generate volume ratio field.
Foam::polyBoundaryMesh.
volScalarField scalarField(fieldObject, mesh)
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:35
label patchi
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66