directions.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 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 "FaceCellWave.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  "e1",
49  "e2",
50  "e3"
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  List<directionInfo> faceInfo(mesh.nFaces()), cellInfo(mesh.nCells());
207  FaceCellWave<directionInfo> directionCalc
208  (
209  mesh,
210  changedFaces,
211  changedFacesInfo,
212  faceInfo,
213  cellInfo,
214  mesh.globalData().nTotalCells() + 1
215  );
216 
217  vectorField dirField(cellInfo.size());
218 
219  label nUnset = 0;
220  label nGeom = 0;
221  label nTopo = 0;
222 
223  forAll(cellInfo, celli)
224  {
225  label index = cellInfo[celli].index();
226 
227  if (index == -3)
228  {
229  // Never visited
231  << "Cell " << celli << " never visited to determine "
232  << "local coordinate system" << endl
233  << "Using direction " << defaultDir << " instead" << endl;
234 
235  dirField[celli] = defaultDir;
236 
237  nUnset++;
238  }
239  else if (index == -2)
240  {
241  // Geometric direction
242  dirField[celli] = cellInfo[celli].n();
243 
244  nGeom++;
245  }
246  else if (index == -1)
247  {
249  << "Illegal index " << index << endl
250  << "Value is only allowed on faces" << abort(FatalError);
251  }
252  else
253  {
254  // Topological edge cut. Convert into average cut direction.
255  dirField[celli] = meshTools::edgeToCutDir(mesh, celli, index);
256 
257  nTopo++;
258  }
259  }
260 
261  reduce(nGeom, sumOp<label>());
262  reduce(nTopo, sumOp<label>());
263  reduce(nUnset, sumOp<label>());
264 
265  Info<< "Calculated local coords for " << defaultDir
266  << endl
267  << " Geometric cut cells : " << nGeom << endl
268  << " Topological cut cells : " << nTopo << endl
269  << " Unset cells : " << nUnset << endl
270  << endl;
271 
272  return dirField;
273 }
274 
275 
276 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
277 
279 (
280  const polyMesh& mesh,
281  const dictionary& dict,
282  const twoDPointCorrector* correct2DPtr
283 )
284 :
285  List<vectorField>(wordList(dict.lookup("directions")).size())
286 {
287  const wordList wantedDirs(dict.lookup("directions"));
288  const word coordSystem(dict.lookup("coordinateSystem"));
289 
290  bool wantE3 = false;
291  bool wantE1 = false;
292  bool wantE2 = false;
293  label nDirs = 0;
294 
295  if (coordSystem != "fieldBased")
296  {
297  forAll(wantedDirs, i)
298  {
299  directionType wantedDir = directionTypeNames_[wantedDirs[i]];
300 
301  if (wantedDir == e3)
302  {
303  wantE3 = true;
304  }
305  else if (wantedDir == e1)
306  {
307  wantE1 = true;
308  }
309  else if (wantedDir == e2)
310  {
311  wantE2 = true;
312  }
313  }
314  }
315 
316 
317  if (coordSystem == "global")
318  {
319  const dictionary& globalDict = dict.subDict("globalCoeffs");
320 
321  vector e1(globalDict.lookup("e1"));
322  check2D(correct2DPtr, e1);
323 
324  vector e2(globalDict.lookup("e2"));
325  check2D(correct2DPtr, e2);
326 
327  vector e3 = e1 ^ e2;
328  e3 /= mag(e3);
329 
330  Info<< "Global Coordinate system:" << endl
331  << " e3 : " << e3 << endl
332  << " e1 : " << e1 << endl
333  << " e2 : " << e2
334  << endl << endl;
335 
336  if (wantE3)
337  {
338  operator[](nDirs++) = vectorField(1, e3);
339  }
340  if (wantE1)
341  {
342  operator[](nDirs++) = vectorField(1, e1);
343  }
344  if (wantE2)
345  {
346  operator[](nDirs++) = vectorField(1, e2);
347  }
348  }
349  else if (coordSystem == "patchLocal")
350  {
351  const dictionary& patchDict = dict.subDict("patchLocalCoeffs");
352 
353  const word patchName(patchDict.lookup("patch"));
354 
355  const label patchi = mesh.boundaryMesh().findPatchID(patchName);
356 
357  if (patchi == -1)
358  {
360  << "Cannot find patch "
361  << patchName
362  << exit(FatalError);
363  }
364 
365  // Take zeroth face on patch
366  const polyPatch& pp = mesh.boundaryMesh()[patchi];
367 
368  vector e1(patchDict.lookup("e1"));
369 
370  const vector& n0 = pp.faceNormals()[0];
371 
372  if (correct2DPtr)
373  {
374  e1 = correct2DPtr->planeNormal() ^ n0;
375 
377  << "Discarding user specified e1 since 2D case." << endl
378  << "Recalculated e1 from face normal and planeNormal as "
379  << e1 << endl << endl;
380  }
381 
382  Switch useTopo(dict.lookup("useHexTopology"));
383 
384  vectorField e3Dirs;
385  vectorField e1Dirs;
386 
387  if (wantE3 || wantE2)
388  {
389  e3Dirs =
390  propagateDirection
391  (
392  mesh,
393  useTopo,
394  pp,
395  pp.faceNormals(),
396  n0
397  );
398 
399  if (wantE3)
400  {
401  this->operator[](nDirs++) = e3Dirs;
402  }
403  }
404 
405  if (wantE1 || wantE2)
406  {
407  e1Dirs =
408  propagateDirection
409  (
410  mesh,
411  useTopo,
412  pp,
413  vectorField(pp.size(), e1),
414  e1
415  );
416 
417 
418  if (wantE1)
419  {
420  this->operator[](nDirs++) = e1Dirs;
421  }
422  }
423  if (wantE2)
424  {
425  tmp<vectorField> e2Dirs = e3Dirs ^ e1Dirs;
426 
427  this->operator[](nDirs++) = e2Dirs;
428  }
429  }
430  else if (coordSystem == "fieldBased")
431  {
432  forAll(wantedDirs, i)
433  {
434  operator[](nDirs++) =
436  (
437  IOobject
438  (
439  mesh.instance()/wantedDirs[i],
440  mesh,
443  )
444  );
445  }
446  }
447  else
448  {
450  << "Unknown coordinate system "
451  << coordSystem << endl
452  << "Known types are global, patchLocal and fieldBased"
453  << exit(FatalError);
454  }
455 }
456 
457 
458 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H: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
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
const word & name() const
Return name.
A class for handling file names.
Definition: fileName.H:79
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
label nFaces() const
Output to file stream.
Definition: OFstream.H:82
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition: FaceCellWave.H:76
label nTotalCells() const
Return total number of cells in decomposed mesh.
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/any.
Definition: Switch.H:60
label findPatchID(const word &patchName) const
Find patch index given a name.
const Cmpt & z() const
Definition: VectorI.H:87
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:66
fvMesh & mesh
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
scalar minDist(const List< pointIndexHit > &hitList)
const Cmpt & y() const
Definition: VectorI.H:81
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:1002
label cutDirToEdge(const primitiveMesh &, const label celli, const vector &cutDir)
Reverse of edgeToCutDir: given direction find edge bundle and.
Definition: meshTools.C:795
Class applies a two-dimensional correction to mesh motion point field.
A class for handling words, derived from string.
Definition: word.H:59
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
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1237
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1435
const Field< PointType > & faceNormals() const
Return face normals for patch.
const vectorField & cellCentres() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:54
const vector & planeNormal() const
Return plane normal.
Holds direction in which to split cell (in fact a local coordinate axes). Information is a label and ...
Definition: directionInfo.H:85
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
IOField< vector > vectorIOField
vectorField with IO.
Definition: vectorIOField.H:42
List< word > wordList
A List of words.
Definition: fileName.H:54
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
directionType
Enumeration listing the possible coordinate directions.
Definition: directions.H:74
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
fileName & instance() const
Return the instance directory, constant, system, <time> etc.
Definition: IOobject.C:355
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
directions(const polyMesh &mesh, const dictionary &dict, const twoDPointCorrector *correct2DPtr=nullptr)
Construct from mesh and dictionary and optional 2D corrector.
Definition: directions.C:279
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:98
const labelListList & cellCells() const
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864