FreeFOAM The Cross-Platform CFD Toolkit
directions.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 \*---------------------------------------------------------------------------*/
25 
26 #include "directions.H"
27 #include <OpenFOAM/polyMesh.H>
30 #include <OpenFOAM/MeshWave.H>
31 #include <OpenFOAM/OFstream.H>
32 #include <meshTools/meshTools.H>
33 #include <OpenFOAM/hexMatcher.H>
34 #include <OpenFOAM/Switch.H>
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 template<>
41 {
42  "tan1",
43  "tan2",
44  "normal"
45 };
46 
48  Foam::directions::directionTypeNames_;
49 
50 
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 // For debugging
55 void Foam::directions::writeOBJ(Ostream& os, const point& pt)
56 {
57  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
58 }
59 
60 // For debugging
62 (
63  Ostream& os,
64  const point& pt0,
65  const point& pt1,
66  label& vertI
67 )
68 {
69  writeOBJ(os, pt0);
70  writeOBJ(os, pt1);
71 
72  os << "l " << vertI + 1 << ' ' << vertI + 2 << endl;
73 
74  vertI += 2;
75 }
76 
77 
78 // Dump to file.
80 (
81  const fileName& fName,
82  const primitiveMesh& mesh,
83  const vectorField& dirs
84 )
85 {
86  Pout<< "Writing cell info to " << fName << " as vectors at the cellCentres"
87  << endl << endl;
88 
89  OFstream xDirStream(fName);
90 
91  label vertI = 0;
92 
93  forAll(dirs, cellI)
94  {
95  const point& ctr = mesh.cellCentres()[cellI];
96 
97  // Calculate local length scale
98  scalar minDist = GREAT;
99 
100  const labelList& nbrs = mesh.cellCells()[cellI];
101 
102  forAll(nbrs, nbrI)
103  {
104  minDist = min(minDist, mag(mesh.cellCentres()[nbrs[nbrI]] - ctr));
105  }
106 
107  scalar scale = 0.5*minDist;
108 
109  writeOBJ(xDirStream, ctr, ctr + scale*dirs[cellI], vertI);
110  }
111 }
112 
113 
114 void Foam::directions::check2D
115 (
116  const twoDPointCorrector* correct2DPtr,
117  const vector& vec
118 )
119 {
120  if (correct2DPtr)
121  {
122  if (mag(correct2DPtr->planeNormal() & vec) > 1E-6)
123  {
124  FatalErrorIn("check2D") << "Specified vector " << vec
125  << "is not normal to plane defined in dynamicMeshDict."
126  << endl
127  << "Either make case 3D or adjust vector."
128  << exit(FatalError);
129  }
130  }
131 }
132 
133 
134 // Get direction on all cells
135 Foam::vectorField Foam::directions::propagateDirection
136 (
137  const polyMesh& mesh,
138  const bool useTopo,
139  const polyPatch& pp,
140  const vectorField& ppField,
141  const vector& defaultDir
142 )
143 {
144  // Seed all faces on patch
145  labelList changedFaces(pp.size());
146  List<directionInfo> changedFacesInfo(pp.size());
147 
148  if (useTopo)
149  {
150  forAll(pp, patchFaceI)
151  {
152  label meshFaceI = pp.start() + patchFaceI;
153 
154  label cellI = mesh.faceOwner()[meshFaceI];
155 
156  if (!hexMatcher().isA(mesh, cellI))
157  {
158  FatalErrorIn("propagateDirection")
159  << "useHexTopology specified but cell " << cellI
160  << " on face " << patchFaceI << " of patch " << pp.name()
161  << " is not a hex" << exit(FatalError);
162  }
163 
164  const vector& cutDir = ppField[patchFaceI];
165 
166  // Get edge(bundle) on cell most in direction of cutdir
167  label edgeI = meshTools::cutDirToEdge(mesh, cellI, cutDir);
168 
169  // Convert edge into index on face
170  label faceIndex =
172  (
173  mesh,
174  cellI,
175  meshFaceI,
176  edgeI
177  );
178 
179  // Set initial face and direction
180  changedFaces[patchFaceI] = meshFaceI;
181  changedFacesInfo[patchFaceI] =
182  directionInfo
183  (
184  faceIndex,
185  cutDir
186  );
187  }
188  }
189  else
190  {
191  forAll(pp, patchFaceI)
192  {
193  changedFaces[patchFaceI] = pp.start() + patchFaceI;
194  changedFacesInfo[patchFaceI] =
195  directionInfo
196  (
197  -2, // Geometric information only
198  ppField[patchFaceI]
199  );
200  }
201  }
202 
203  MeshWave<directionInfo> directionCalc
204  (
205  mesh,
206  changedFaces,
207  changedFacesInfo,
208  mesh.globalData().nTotalCells()+1
209  );
210 
211  const List<directionInfo>& cellInfo = directionCalc.allCellInfo();
212 
213  vectorField dirField(cellInfo.size());
214 
215  label nUnset = 0;
216  label nGeom = 0;
217  label nTopo = 0;
218 
219  forAll(cellInfo, cellI)
220  {
221  label index = cellInfo[cellI].index();
222 
223  if (index == -3)
224  {
225  // Never visited
226  WarningIn("propagateDirection")
227  << "Cell " << cellI << " never visited to determine "
228  << "local coordinate system" << endl
229  << "Using direction " << defaultDir << " instead" << endl;
230 
231  dirField[cellI] = defaultDir;
232 
233  nUnset++;
234  }
235  else if (index == -2)
236  {
237  // Geometric direction
238  dirField[cellI] = cellInfo[cellI].n();
239 
240  nGeom++;
241  }
242  else if (index == -1)
243  {
244  FatalErrorIn("propagateDirection")
245  << "Illegal index " << index << endl
246  << "Value is only allowed on faces" << abort(FatalError);
247  }
248  else
249  {
250  // Topological edge cut. Convert into average cut direction.
251  dirField[cellI] = meshTools::edgeToCutDir(mesh, cellI, index);
252 
253  nTopo++;
254  }
255  }
256 
257  Pout<< "Calculated local coords for " << defaultDir
258  << endl
259  << " Geometric cut cells : " << nGeom << endl
260  << " Topological cut cells : " << nTopo << endl
261  << " Unset cells : " << nUnset << endl
262  << endl;
263 
264  return dirField;
265 }
266 
267 
268 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
269 
270 // Construct from dictionary
271 Foam::directions::directions
272 (
273  const polyMesh& mesh,
274  const dictionary& dict,
275  const twoDPointCorrector* correct2DPtr
276 )
277 :
278  List<vectorField>(wordList(dict.lookup("directions")).size())
279 {
280  const wordList wantedDirs(dict.lookup("directions"));
281 
282  bool wantNormal = false;
283  bool wantTan1 = false;
284  bool wantTan2 = false;
285 
286  forAll(wantedDirs, i)
287  {
288  directionType wantedDir = directionTypeNames_[wantedDirs[i]];
289 
290  if (wantedDir == NORMAL)
291  {
292  wantNormal = true;
293  }
294  else if (wantedDir == TAN1)
295  {
296  wantTan1 = true;
297  }
298  else if (wantedDir == TAN2)
299  {
300  wantTan2 = true;
301  }
302  }
303 
304 
305  label nDirs = 0;
306 
307  word coordSystem(dict.lookup("coordinateSystem"));
308 
309  if (coordSystem == "global")
310  {
311  const dictionary& globalDict = dict.subDict("globalCoeffs");
312 
313  vector tan1(globalDict.lookup("tan1"));
314  check2D(correct2DPtr, tan1);
315 
316  vector tan2(globalDict.lookup("tan2"));
317  check2D(correct2DPtr, tan2);
318 
319  vector normal = tan1 ^ tan2;
320  normal /= mag(normal);
321 
322  Pout<< "Global Coordinate system:" << endl
323  << " normal : " << normal << endl
324  << " tan1 : " << tan1 << endl
325  << " tan2 : " << tan2
326  << endl << endl;
327 
328  if (wantNormal)
329  {
330  operator[](nDirs++) = vectorField(1, normal);
331  }
332  if (wantTan1)
333  {
334  operator[](nDirs++) = vectorField(1, tan1);
335  }
336  if (wantTan2)
337  {
338  operator[](nDirs++) = vectorField(1, tan2);
339  }
340  }
341  else if (coordSystem == "patchLocal")
342  {
343  const dictionary& patchDict = dict.subDict("patchLocalCoeffs");
344 
345  word patchName(patchDict.lookup("patch"));
346 
347  label patchI = mesh.boundaryMesh().findPatchID(patchName);
348 
349  if (patchI == -1)
350  {
352  (
353  "directions::directions(const polyMesh&, const dictionary&,"
354  "const twoDPointCorrector*"
355  ) << "Cannot find patch "
356  << patchName << exit(FatalError);
357  }
358 
359  // Take zeroth face on patch
360  const polyPatch& pp = mesh.boundaryMesh()[patchI];
361 
362  vector tan1(patchDict.lookup("tan1"));
363 
364  const vector& n0 = pp.faceNormals()[0];
365 
366  if (correct2DPtr)
367  {
368  tan1 = correct2DPtr->planeNormal() ^ n0;
369 
370  WarningIn
371  (
372  "directions::directions(const polyMesh&, const dictionary&,"
373  "const twoDPointCorrector*"
374  ) << "Discarding user specified tan1 since 2D case." << endl
375  << "Recalculated tan1 from face normal and planeNormal as "
376  << tan1 << endl << endl;
377  }
378 
379  Switch useTopo(dict.lookup("useHexTopology"));
380 
381  vectorField normalDirs;
382  vectorField tan1Dirs;
383 
384  if (wantNormal || wantTan2)
385  {
386  normalDirs =
387  propagateDirection
388  (
389  mesh,
390  useTopo,
391  pp,
392  pp.faceNormals(),
393  n0
394  );
395 
396  if (wantNormal)
397  {
398  operator[](nDirs++) = normalDirs;
399 
401  //writeOBJ("normal.obj", mesh, normalDirs);
402  }
403  }
404 
405  if (wantTan1 || wantTan2)
406  {
407  tan1Dirs =
408  propagateDirection
409  (
410  mesh,
411  useTopo,
412  pp,
413  vectorField(pp.size(), tan1),
414  tan1
415  );
416 
417 
418  if (wantTan1)
419  {
420  operator[](nDirs++) = tan1Dirs;
421 
423  //writeOBJ("tan1.obj", mesh, tan1Dirs);
424  }
425  }
426  if (wantTan2)
427  {
428  vectorField tan2Dirs = normalDirs ^ tan1Dirs;
429 
430  operator[](nDirs++) = tan2Dirs;
431 
433  //writeOBJ("tan2.obj", mesh, tan2Dirs);
434  }
435  }
436  else
437  {
439  (
440  "directions::directions(const polyMesh&, const dictionary&,"
441  "const twoDPointCorrector*"
442  ) << "Unknown coordinate system "
443  << coordSystem << endl
444  << "Known types are global and patchLocal" << exit(FatalError);
445  }
446 }
447 
448 
449 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
450 
451 
452 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
453 
454 
455 // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
456 
457 
458 // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
459 
460 
461 // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
462 
463 
464 // ************************ vim: set sw=4 sts=4 et: ************************ //