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) 2011-2016 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 "directions.H"
27 #include "polyMesh.H"
28 #include "twoDPointCorrector.H"
29 #include "directionInfo.H"
30 #include "MeshWave.H"
31 #include "OFstream.H"
32 #include "meshTools.H"
33 #include "hexMatcher.H"
34 #include "Switch.H"
35 #include "globalMeshData.H"
36 
37 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 
39 namespace Foam
40 {
41  template<>
42  const char* Foam::NamedEnum
43  <
45  3
46  >::names[] =
47  {
48  "tan1",
49  "tan2",
50  "normal"
51  };
52 }
53 
55  Foam::directions::directionTypeNames_;
56 
57 
58 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
59 
60 void Foam::directions::writeOBJ(Ostream& os, const point& pt)
61 {
62  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
63 }
64 
65 
66 void Foam::directions::writeOBJ
67 (
68  Ostream& os,
69  const point& pt0,
70  const point& pt1,
71  label& vertI
72 )
73 {
74  writeOBJ(os, pt0);
75  writeOBJ(os, pt1);
76 
77  os << "l " << vertI + 1 << ' ' << vertI + 2 << endl;
78 
79  vertI += 2;
80 }
81 
82 
83 void Foam::directions::writeOBJ
84 (
85  const fileName& fName,
86  const primitiveMesh& mesh,
87  const vectorField& dirs
88 )
89 {
90  Pout<< "Writing cell info to " << fName << " as vectors at the cellCentres"
91  << endl << endl;
92 
93  OFstream xDirStream(fName);
94 
95  label vertI = 0;
96 
97  forAll(dirs, celli)
98  {
99  const point& ctr = mesh.cellCentres()[celli];
100 
101  // Calculate local length scale
102  scalar minDist = GREAT;
103 
104  const labelList& nbrs = mesh.cellCells()[celli];
105 
106  forAll(nbrs, nbrI)
107  {
108  minDist = min(minDist, mag(mesh.cellCentres()[nbrs[nbrI]] - ctr));
109  }
110 
111  scalar scale = 0.5*minDist;
112 
113  writeOBJ(xDirStream, ctr, ctr + scale*dirs[celli], vertI);
114  }
115 }
116 
117 
118 void Foam::directions::check2D
119 (
120  const twoDPointCorrector* correct2DPtr,
121  const vector& vec
122 )
123 {
124  if (correct2DPtr)
125  {
126  if (mag(correct2DPtr->planeNormal() & vec) > 1e-6)
127  {
129  << "is not normal to plane defined in dynamicMeshDict."
130  << endl
131  << "Either make case 3D or adjust vector."
132  << exit(FatalError);
133  }
134  }
135 }
136 
137 
138 Foam::vectorField Foam::directions::propagateDirection
139 (
140  const polyMesh& mesh,
141  const bool useTopo,
142  const polyPatch& pp,
143  const vectorField& ppField,
144  const vector& defaultDir
145 )
146 {
147  // Seed all faces on patch
148  labelList changedFaces(pp.size());
149  List<directionInfo> changedFacesInfo(pp.size());
150 
151  if (useTopo)
152  {
153  forAll(pp, patchFacei)
154  {
155  label meshFacei = pp.start() + patchFacei;
156 
157  label celli = mesh.faceOwner()[meshFacei];
158 
159  if (!hexMatcher().isA(mesh, celli))
160  {
162  << "useHexTopology specified but cell " << celli
163  << " on face " << patchFacei << " of patch " << pp.name()
164  << " is not a hex" << exit(FatalError);
165  }
166 
167  const vector& cutDir = ppField[patchFacei];
168 
169  // Get edge(bundle) on cell most in direction of cutdir
170  label edgeI = meshTools::cutDirToEdge(mesh, celli, cutDir);
171 
172  // Convert edge into index on face
173  label faceIndex =
175  (
176  mesh,
177  celli,
178  meshFacei,
179  edgeI
180  );
181 
182  // Set initial face and direction
183  changedFaces[patchFacei] = meshFacei;
184  changedFacesInfo[patchFacei] =
186  (
187  faceIndex,
188  cutDir
189  );
190  }
191  }
192  else
193  {
194  forAll(pp, patchFacei)
195  {
196  changedFaces[patchFacei] = pp.start() + patchFacei;
197  changedFacesInfo[patchFacei] =
199  (
200  -2, // Geometric information only
201  ppField[patchFacei]
202  );
203  }
204  }
205 
206  MeshWave<directionInfo> directionCalc
207  (
208  mesh,
209  changedFaces,
210  changedFacesInfo,
211  mesh.globalData().nTotalCells()+1
212  );
213 
214  const List<directionInfo>& cellInfo = directionCalc.allCellInfo();
215 
216  vectorField dirField(cellInfo.size());
217 
218  label nUnset = 0;
219  label nGeom = 0;
220  label nTopo = 0;
221 
222  forAll(cellInfo, celli)
223  {
224  label index = cellInfo[celli].index();
225 
226  if (index == -3)
227  {
228  // Never visited
230  << "Cell " << celli << " never visited to determine "
231  << "local coordinate system" << endl
232  << "Using direction " << defaultDir << " instead" << endl;
233 
234  dirField[celli] = defaultDir;
235 
236  nUnset++;
237  }
238  else if (index == -2)
239  {
240  // Geometric direction
241  dirField[celli] = cellInfo[celli].n();
242 
243  nGeom++;
244  }
245  else if (index == -1)
246  {
248  << "Illegal index " << index << endl
249  << "Value is only allowed on faces" << abort(FatalError);
250  }
251  else
252  {
253  // Topological edge cut. Convert into average cut direction.
254  dirField[celli] = meshTools::edgeToCutDir(mesh, celli, index);
255 
256  nTopo++;
257  }
258  }
259 
260  Pout<< "Calculated local coords for " << defaultDir
261  << endl
262  << " Geometric cut cells : " << nGeom << endl
263  << " Topological cut cells : " << nTopo << endl
264  << " Unset cells : " << nUnset << endl
265  << endl;
266 
267  return dirField;
268 }
269 
270 
271 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
272 
273 Foam::directions::directions
274 (
275  const polyMesh& mesh,
276  const dictionary& dict,
277  const twoDPointCorrector* correct2DPtr
278 )
279 :
280  List<vectorField>(wordList(dict.lookup("directions")).size())
281 {
282  const wordList wantedDirs(dict.lookup("directions"));
283  const word coordSystem(dict.lookup("coordinateSystem"));
284 
285  bool wantNormal = false;
286  bool wantTan1 = false;
287  bool wantTan2 = false;
288  label nDirs = 0;
289 
290  if (coordSystem != "fieldBased")
291  {
292  forAll(wantedDirs, i)
293  {
294  directionType wantedDir = directionTypeNames_[wantedDirs[i]];
295 
296  if (wantedDir == NORMAL)
297  {
298  wantNormal = true;
299  }
300  else if (wantedDir == TAN1)
301  {
302  wantTan1 = true;
303  }
304  else if (wantedDir == TAN2)
305  {
306  wantTan2 = true;
307  }
308  }
309  }
310 
311 
312  if (coordSystem == "global")
313  {
314  const dictionary& globalDict = dict.subDict("globalCoeffs");
315 
316  vector tan1(globalDict.lookup("tan1"));
317  check2D(correct2DPtr, tan1);
318 
319  vector tan2(globalDict.lookup("tan2"));
320  check2D(correct2DPtr, tan2);
321 
322  vector normal = tan1 ^ tan2;
323  normal /= mag(normal);
324 
325  Pout<< "Global Coordinate system:" << endl
326  << " normal : " << normal << endl
327  << " tan1 : " << tan1 << endl
328  << " tan2 : " << tan2
329  << endl << endl;
330 
331  if (wantNormal)
332  {
333  operator[](nDirs++) = vectorField(1, normal);
334  }
335  if (wantTan1)
336  {
337  operator[](nDirs++) = vectorField(1, tan1);
338  }
339  if (wantTan2)
340  {
341  operator[](nDirs++) = vectorField(1, tan2);
342  }
343  }
344  else if (coordSystem == "patchLocal")
345  {
346  const dictionary& patchDict = dict.subDict("patchLocalCoeffs");
347 
348  const word patchName(patchDict.lookup("patch"));
349 
350  const label patchi = mesh.boundaryMesh().findPatchID(patchName);
351 
352  if (patchi == -1)
353  {
355  << "Cannot find patch "
356  << patchName
357  << exit(FatalError);
358  }
359 
360  // Take zeroth face on patch
361  const polyPatch& pp = mesh.boundaryMesh()[patchi];
362 
363  vector tan1(patchDict.lookup("tan1"));
364 
365  const vector& n0 = pp.faceNormals()[0];
366 
367  if (correct2DPtr)
368  {
369  tan1 = correct2DPtr->planeNormal() ^ n0;
370 
372  << "Discarding user specified tan1 since 2D case." << endl
373  << "Recalculated tan1 from face normal and planeNormal as "
374  << tan1 << endl << endl;
375  }
376 
377  Switch useTopo(dict.lookup("useHexTopology"));
378 
379  vectorField normalDirs;
380  vectorField tan1Dirs;
381 
382  if (wantNormal || wantTan2)
383  {
384  normalDirs =
385  propagateDirection
386  (
387  mesh,
388  useTopo,
389  pp,
390  pp.faceNormals(),
391  n0
392  );
393 
394  if (wantNormal)
395  {
396  this->operator[](nDirs++) = normalDirs;
397  }
398  }
399 
400  if (wantTan1 || wantTan2)
401  {
402  tan1Dirs =
403  propagateDirection
404  (
405  mesh,
406  useTopo,
407  pp,
408  vectorField(pp.size(), tan1),
409  tan1
410  );
411 
412 
413  if (wantTan1)
414  {
415  this->operator[](nDirs++) = tan1Dirs;
416  }
417  }
418  if (wantTan2)
419  {
420  tmp<vectorField> tan2Dirs = normalDirs ^ tan1Dirs;
421 
422  this->operator[](nDirs++) = tan2Dirs;
423  }
424  }
425  else if (coordSystem == "fieldBased")
426  {
427  forAll(wantedDirs, i)
428  {
429  operator[](nDirs++) =
431  (
432  IOobject
433  (
434  mesh.instance()/wantedDirs[i],
435  mesh,
438  )
439  );
440  }
441  }
442  else
443  {
445  << "Unknown coordinate system "
446  << coordSystem << endl
447  << "Known types are global, patchLocal and fieldBased"
448  << exit(FatalError);
449  }
450 }
451 
452 
453 // ************************************************************************* //
bool isA(const Type &t)
Check if a dynamic_cast to typeid is possible.
Definition: typeInfo.H:134
vector edgeToCutDir(const primitiveMesh &, const label celli, const label edgeI)
Given edge on hex find all &#39;parallel&#39; (i.e. non-connected)
Definition: meshTools.C:747
const Cmpt & z() const
Definition: VectorI.H:87
const labelListList & cellCells() const
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A class for handling file names.
Definition: fileName.H:69
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const Cmpt & x() const
Definition: VectorI.H:75
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
const vector & planeNormal() const
Return plane normal.
Output to file stream.
Definition: OFstream.H:81
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:633
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
A cellMatcher for hex cells.
Definition: hexMatcher.H:51
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
Holds information regarding type of cell. Used in inside/outside determination in cellClassification...
Definition: cellInfo.H:63
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:52
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
label cutDirToEdge(const primitiveMesh &, const label celli, const vector &cutDir)
Reverse of edgeToCutDir: given direction find edge bundle and.
Definition: meshTools.C:795
dynamicFvMesh & mesh
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
Class applies a two-dimensional correction to mesh motion point field.
A class for handling words, derived from string.
Definition: word.H:59
FaceCellWave plus data.
Definition: MeshWave.H:56
static label edgeToFaceIndex(const primitiveMesh &mesh, const label celli, const label facei, const label edgeI)
Given edge on hex cell find corresponding edge on face. Is either.
Definition: directionInfo.C:93
const Cmpt & y() const
Definition: VectorI.H:81
const vectorField & cellCentres() const
const word & name() const
Return name.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1140
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
Holds direction in which to split cell (in fact a local coordinate axes). Information is a label and ...
Definition: directionInfo.H:84
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
IOField< vector > vectorIOField
vectorField with IO.
Definition: vectorIOField.H:42
A normal distribution model.
List< word > wordList
A List of words.
Definition: fileName.H:54
label patchi
const Field< PointType > & faceNormals() const
Return face normals for patch.
#define WarningInFunction
Report a warning using Foam::Warning.
directionType
Enumeration listing the possible coordinate directions.
Definition: directions.H:74
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A class for managing temporary objects.
Definition: PtrList.H:54
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label findPatchID(const word &patchName) const
Find patch index given a name.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
const fileName & instance() const
Definition: IOobject.H:337
label nTotalCells() const
Return total number of cells in decomposed mesh.
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451