FreeFOAM The Cross-Platform CFD Toolkit
primitiveMeshCheckMotion.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) 1991-2010 OpenCFD Ltd.
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 Description
25  Given a set of points, find out if the mesh resulting from point motion will
26  be valid without actually changing the mesh.
27 
28 \*---------------------------------------------------------------------------*/
29 
30 #include <OpenFOAM/primitiveMesh.H>
32 #include <OpenFOAM/cell.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
41 (
42  const pointField& newPoints,
43  const bool report
44 ) const
45 {
46  if (debug || report)
47  {
48  Pout<< "bool primitiveMesh::checkMeshMotion("
49  << "const pointField& newPoints, const bool report) const: "
50  << "checking mesh motion" << endl;
51  }
52 
53  bool error = false;
54 
55  const faceList& f = faces();
56 
57  const labelList& own = faceOwner();
58  const labelList& nei = faceNeighbour();
59 
60  vectorField fCtrs(nFaces());
61  vectorField fAreas(nFaces());
62 
63  makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
64 
65  // Check cell volumes and calculate new cell centres
66  vectorField cellCtrs(nCells());
67  scalarField cellVols(nCells());
68 
69  makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
70 
71  scalar minVolume = GREAT;
72  label nNegVols = 0;
73 
74  forAll (cellVols, cellI)
75  {
76  if (cellVols[cellI] < VSMALL)
77  {
78  if (debug || report)
79  {
80  Pout<< "Zero or negative cell volume detected for cell "
81  << cellI << ". Volume = " << cellVols[cellI] << endl;
82  }
83 
84  nNegVols++;
85  }
86 
87  minVolume = min(minVolume, cellVols[cellI]);
88  }
89 
90  if (nNegVols > 0)
91  {
92  error = true;
93 
94  Pout<< "Zero or negative cell volume in mesh motion in " << nNegVols
95  << " cells. Min volume: " << minVolume << endl;
96  }
97  else
98  {
99  if (debug || report)
100  {
101  Pout<< "Min volume = " << minVolume
102  << ". Total volume = " << sum(cellVols)
103  << ". Cell volumes OK." << endl;
104  }
105  }
106 
107  // Check face areas
108 
109  scalar minArea = GREAT;
110  label nNegAreas = 0;
111  label nPyrErrors = 0;
112  label nDotProductErrors = 0;
113 
114  forAll (f, faceI)
115  {
116  const scalar a = Foam::mag(fAreas[faceI]);
117 
118  if (a < VSMALL)
119  {
120  if (debug || report)
121  {
122  if (isInternalFace(faceI))
123  {
124  Pout<< "Zero or negative face area detected for "
125  << "internal face "<< faceI << " between cells "
126  << own[faceI] << " and " << nei[faceI]
127  << ". Face area magnitude = " << a << endl;
128  }
129  else
130  {
131  Pout<< "Zero or negative face area detected for "
132  << "boundary face " << faceI << " next to cell "
133  << own[faceI] << ". Face area magnitude = "
134  << a << endl;
135  }
136  }
137 
138  nNegAreas++;
139  }
140 
141  minArea = min(minArea, a);
142 
143  // Create the owner pyramid - it will have negative volume
144  scalar pyrVol =
145  pyramidPointFaceRef(f[faceI], cellCtrs[own[faceI]]).mag(newPoints);
146 
147  if (pyrVol > SMALL)
148  {
149  if (debug || report)
150  {
151  Pout<< "Negative pyramid volume: " << -pyrVol
152  << " for face " << faceI << " " << f[faceI]
153  << " and owner cell: " << own[faceI] << endl
154  << "Owner cell vertex labels: "
155  << cells()[own[faceI]].labels(f)
156  << endl;
157  }
158 
159  nPyrErrors++;
160  }
161 
162  if (isInternalFace(faceI))
163  {
164  // Create the neighbour pyramid - it will have positive volume
165  scalar pyrVol =
167  (
168  f[faceI],
169  cellCtrs[nei[faceI]]
170  ).mag(newPoints);
171 
172  if (pyrVol < -SMALL)
173  {
174  if (debug || report)
175  {
176  Pout<< "Negative pyramid volume: " << pyrVol
177  << " for face " << faceI << " " << f[faceI]
178  << " and neighbour cell: " << nei[faceI] << nl
179  << "Neighbour cell vertex labels: "
180  << cells()[nei[faceI]].labels(f)
181  << endl;
182  }
183 
184  nPyrErrors++;
185  }
186 
187  const vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
188  const vector& s = fAreas[faceI];
189  scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
190 
191  // Only write full message the first time
192  if (dDotS < SMALL && nDotProductErrors == 0)
193  {
194  // Non-orthogonality greater than 90 deg
195  WarningIn
196  (
197  "primitiveMesh::checkMeshMotion"
198  "(const pointField& newPoints, const bool report) const"
199  ) << "Severe non-orthogonality in mesh motion for face "
200  << faceI
201  << " between cells " << own[faceI] << " and " << nei[faceI]
202  << ": Angle = " << ::acos(dDotS)/mathematicalConstant::pi*180.0
203  << " deg." << endl;
204 
205  nDotProductErrors++;
206  }
207  }
208  }
209 
210  if (nNegAreas > 0)
211  {
212  error = true;
213 
214  WarningIn
215  (
216  "primitiveMesh::checkMeshMotion"
217  "(const pointField& newPoints, const bool report) const"
218  ) << "Zero or negative face area in mesh motion in " << nNegAreas
219  << " faces. Min area: " << minArea << endl;
220  }
221  else
222  {
223  if (debug || report)
224  {
225  Pout<< "Min area = " << minArea
226  << ". Face areas OK." << endl;
227  }
228  }
229 
230  if (nPyrErrors > 0)
231  {
232  Pout<< "Detected " << nPyrErrors
233  << " negative pyramid volume in mesh motion" << endl;
234 
235  error = true;
236  }
237  else
238  {
239  if (debug || report)
240  {
241  Pout<< "Pyramid volumes OK." << endl;
242  }
243  }
244 
245  if (nDotProductErrors > 0)
246  {
247  Pout<< "Detected " << nDotProductErrors
248  << " in non-orthogonality in mesh motion." << endl;
249 
250  error = true;
251  }
252  else
253  {
254  if (debug || report)
255  {
256  Pout<< "Non-orthogonality check OK." << endl;
257  }
258  }
259 
260  if (!error && (debug || report))
261  {
262  Pout << "Mesh motion check OK." << endl;
263  }
264 
265  return error;
266 }
267 
268 
269 // ************************ vim: set sw=4 sts=4 et: ************************ //