FreeFOAM The Cross-Platform CFD Toolkit
twoDPointCorrector.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  Class applies a two-dimensional correction to mesh motion point field.
26  The correction guarantees that the mesh does not get twisted during motion
27  and thus introduce a third dimension into a 2-D problem.
28 
29 \*---------------------------------------------------------------------------*/
30 
31 #include "twoDPointCorrector.H"
32 #include <OpenFOAM/polyMesh.H>
35 
36 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 
41 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
42 
43 const scalar twoDPointCorrector::edgeOrthogonalityTol = 1.0 - 1e-4;
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 void twoDPointCorrector::calcAddressing() const
48 {
49  // Find geometry normal
50  planeNormalPtr_ = new vector(0, 0, 0);
51  vector& pn = *planeNormalPtr_;
52 
53  bool isWedge = false;
54 
55  // Algorithm:
56  // Attempt to find wedge patch and work out the normal from it.
57  // If not found, find an empty patch with faces in it and use the
58  // normal of the first face. If neither is found, declare an
59  // error.
60 
61  // Try and find a wedge patch
62  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
63 
64  forAll (patches, patchI)
65  {
66  if (isA<wedgePolyPatch>(patches[patchI]))
67  {
68  isWedge = true;
69 
70  pn = refCast<const wedgePolyPatch>(patches[patchI]).centreNormal();
71 
72  if (polyMesh::debug)
73  {
74  Pout << "Found normal from wedge patch " << patchI;
75  }
76 
77  break;
78  }
79  }
80 
81  // Try to find an empty patch with faces
82  if (!isWedge)
83  {
84  forAll (patches, patchI)
85  {
86  if (isA<emptyPolyPatch>(patches[patchI]) && patches[patchI].size())
87  {
88  pn = patches[patchI].faceAreas()[0];
89 
90  if (polyMesh::debug)
91  {
92  Pout << "Found normal from empty patch " << patchI;
93  }
94 
95  break;
96  }
97  }
98  }
99 
100 
101  if (mag(pn) < VSMALL)
102  {
104  (
105  "twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh, "
106  "const vector& n)"
107  ) << "Cannot determine normal vector from patches."
108  << abort(FatalError);
109  }
110  else
111  {
112  pn /= mag(pn);
113  }
114 
115  if (polyMesh::debug)
116  {
117  Pout << " twoDPointCorrector normal: " << pn << endl;
118  }
119 
120  // Select edges to be included in check.
121  normalEdgeIndicesPtr_ = new labelList(mesh_.nEdges());
122  labelList& neIndices = *normalEdgeIndicesPtr_;
123 
124  const edgeList& meshEdges = mesh_.edges();
125 
126  const pointField& meshPoints = mesh_.points();
127 
128  label nNormalEdges = 0;
129 
130  forAll (meshEdges, edgeI)
131  {
132  vector edgeVector =
133  meshEdges[edgeI].vec(meshPoints)/
134  (meshEdges[edgeI].mag(meshPoints) + VSMALL);
135 
136  if (mag(edgeVector & pn) > edgeOrthogonalityTol)
137  {
138  // this edge is normal to the plane. Add it to the list
139  neIndices[nNormalEdges++] = edgeI;
140  }
141  }
142 
143  neIndices.setSize(nNormalEdges);
144 
145  // Construction check: number of points in a read 2-D or wedge geometry
146  // should be odd and the number of edges normal to the plane should be
147  // exactly half the number of points
148  if (!isWedge)
149  {
150  if (meshPoints.size() % 2 != 0)
151  {
152  WarningIn
153  (
154  "twoDPointCorrector::twoDPointCorrector("
155  "const polyMesh& mesh, const vector& n)"
156  ) << "the number of vertices in the geometry "
157  << "is odd - this should not be the case for a 2-D case. "
158  << "Please check the geometry."
159  << endl;
160  }
161 
162  if (2*nNormalEdges != meshPoints.size())
163  {
164  WarningIn
165  (
166  "twoDPointCorrector::twoDPointCorrector("
167  "const polyMesh& mesh, const vector& n)"
168  ) << "The number of points in the mesh is "
169  << "not equal to twice the number of edges normal to the plane "
170  << "- this may be OK only for wedge geometries.\n"
171  << " Please check the geometry or adjust "
172  << "the orthogonality tolerance.\n" << endl
173  << "Number of normal edges: " << nNormalEdges
174  << " number of points: " << meshPoints.size()
175  << endl;
176  }
177  }
178 }
179 
180 
181 void twoDPointCorrector::clearAddressing() const
182 {
183  deleteDemandDrivenData(planeNormalPtr_);
184  deleteDemandDrivenData(normalEdgeIndicesPtr_);
185 }
186 
187 
188 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
189 
190 twoDPointCorrector::twoDPointCorrector(const polyMesh& mesh)
191 :
192  mesh_(mesh),
193  required_(mesh_.nGeometricD() == 2),
194  planeNormalPtr_(NULL),
195  normalEdgeIndicesPtr_(NULL)
196 {}
197 
198 
199 
200 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
201 
203 {
204  clearAddressing();
205 }
206 
207 
208 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
209 
211 {
212  const vector& pn = planeNormal();
213 
214  if (mag(pn.x()) >= edgeOrthogonalityTol)
215  {
216  return vector::X;
217  }
218  else if (mag(pn.y()) >= edgeOrthogonalityTol)
219  {
220  return vector::Y;
221  }
222  else if (mag(pn.z()) >= edgeOrthogonalityTol)
223  {
224  return vector::Z;
225  }
226  else
227  {
228  FatalErrorIn("direction twoDPointCorrector::normalDir() const")
229  << "Plane normal not aligned with the coordinate system" << nl
230  << " pn = " << pn
231  << abort(FatalError);
232 
233  return vector::Z;
234  }
235 }
236 
237 
238 // Return plane normal
240 {
241  if (!planeNormalPtr_)
242  {
243  calcAddressing();
244  }
245 
246  return *planeNormalPtr_;
247 }
248 
249 
250 // Return indices of normal edges.
252 {
253  if (!normalEdgeIndicesPtr_)
254  {
255  calcAddressing();
256  }
257 
258  return *normalEdgeIndicesPtr_;
259 }
260 
261 
263 {
264  if (!required_) return;
265 
266  // Algorithm:
267  // Loop through all edges. Calculate the average point position A for
268  // the front and the back. Correct the position of point P (in two planes)
269  // such that vectors AP and planeNormal are parallel
270 
271  // Get reference to edges
272  const edgeList& meshEdges = mesh_.edges();
273 
274  const labelList& neIndices = normalEdgeIndices();
275  const vector& pn = planeNormal();
276 
277  forAll (neIndices, edgeI)
278  {
279  point& pStart = p[meshEdges[neIndices[edgeI]].start()];
280 
281  point& pEnd = p[meshEdges[neIndices[edgeI]].end()];
282 
283  // calculate average point position
284  const point A = 0.5*(pStart + pEnd);
285 
286  // correct point locations
287  pStart = A + pn*(pn & (pStart - A));
288  pEnd = A + pn*(pn & (pEnd - A));
289  }
290 }
291 
292 
294 {
295  clearAddressing();
296 }
297 
298 
299 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
300 
301 } // End namespace Foam
302 
303 // ************************ vim: set sw=4 sts=4 et: ************************ //