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 | Copyright (C) 2012-2013 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();
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();
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();
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();
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 // ************************************************************************* //
label nFaces() const
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
dimensionedTensor skew(const dimensionedTensor &dt)
dimensioned< scalar > mag(const dimensioned< Type > &)
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
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
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
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Foam::polyBoundaryMesh.
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
static tmp< scalarField > faceWeights(const polyMesh &mesh, const vectorField &fCtrs, const vectorField &fAreas, const vectorField &cellCtrs)
Generate interpolation factors field.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
label nInternalFaces() const
static tmp< scalarField > volRatio(const polyMesh &mesh, const scalarField &vol)
Generate volume ratio field.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
volScalarField scalarField(fieldObject, mesh)
A class for managing temporary objects.
Definition: PtrList.H:118