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-2021 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  "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  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  reduce(nGeom, sumOp<label>());
261  reduce(nTopo, sumOp<label>());
262  reduce(nUnset, sumOp<label>());
263 
264  Info<< "Calculated local coords for " << defaultDir
265  << endl
266  << " Geometric cut cells : " << nGeom << endl
267  << " Topological cut cells : " << nTopo << endl
268  << " Unset cells : " << nUnset << endl
269  << endl;
270 
271  return dirField;
272 }
273 
274 
275 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
276 
278 (
279  const polyMesh& mesh,
280  const dictionary& dict,
281  const twoDPointCorrector* correct2DPtr
282 )
283 :
284  List<vectorField>(wordList(dict.lookup("directions")).size())
285 {
286  const wordList wantedDirs(dict.lookup("directions"));
287  const word coordSystem(dict.lookup("coordinateSystem"));
288 
289  bool wantE3 = false;
290  bool wantE1 = false;
291  bool wantE2 = false;
292  label nDirs = 0;
293 
294  if (coordSystem != "fieldBased")
295  {
296  forAll(wantedDirs, i)
297  {
298  directionType wantedDir = directionTypeNames_[wantedDirs[i]];
299 
300  if (wantedDir == e3)
301  {
302  wantE3 = true;
303  }
304  else if (wantedDir == e1)
305  {
306  wantE1 = true;
307  }
308  else if (wantedDir == e2)
309  {
310  wantE2 = true;
311  }
312  }
313  }
314 
315 
316  if (coordSystem == "global")
317  {
318  const dictionary& globalDict = dict.subDict("globalCoeffs");
319 
320  vector e1(globalDict.lookup("e1"));
321  check2D(correct2DPtr, e1);
322 
323  vector e2(globalDict.lookup("e2"));
324  check2D(correct2DPtr, e2);
325 
326  vector e3 = e1 ^ e2;
327  e3 /= mag(e3);
328 
329  Info<< "Global Coordinate system:" << endl
330  << " e3 : " << e3 << endl
331  << " e1 : " << e1 << endl
332  << " e2 : " << e2
333  << endl << endl;
334 
335  if (wantE3)
336  {
337  operator[](nDirs++) = vectorField(1, e3);
338  }
339  if (wantE1)
340  {
341  operator[](nDirs++) = vectorField(1, e1);
342  }
343  if (wantE2)
344  {
345  operator[](nDirs++) = vectorField(1, e2);
346  }
347  }
348  else if (coordSystem == "patchLocal")
349  {
350  const dictionary& patchDict = dict.subDict("patchLocalCoeffs");
351 
352  const word patchName(patchDict.lookup("patch"));
353 
354  const label patchi = mesh.boundaryMesh().findPatchID(patchName);
355 
356  if (patchi == -1)
357  {
359  << "Cannot find patch "
360  << patchName
361  << exit(FatalError);
362  }
363 
364  // Take zeroth face on patch
365  const polyPatch& pp = mesh.boundaryMesh()[patchi];
366 
367  vector e1(patchDict.lookup("e1"));
368 
369  const vector& n0 = pp.faceNormals()[0];
370 
371  if (correct2DPtr)
372  {
373  e1 = correct2DPtr->planeNormal() ^ n0;
374 
376  << "Discarding user specified e1 since 2D case." << endl
377  << "Recalculated e1 from face normal and planeNormal as "
378  << e1 << endl << endl;
379  }
380 
381  Switch useTopo(dict.lookup("useHexTopology"));
382 
383  vectorField e3Dirs;
384  vectorField e1Dirs;
385 
386  if (wantE3 || wantE2)
387  {
388  e3Dirs =
389  propagateDirection
390  (
391  mesh,
392  useTopo,
393  pp,
394  pp.faceNormals(),
395  n0
396  );
397 
398  if (wantE3)
399  {
400  this->operator[](nDirs++) = e3Dirs;
401  }
402  }
403 
404  if (wantE1 || wantE2)
405  {
406  e1Dirs =
407  propagateDirection
408  (
409  mesh,
410  useTopo,
411  pp,
412  vectorField(pp.size(), e1),
413  e1
414  );
415 
416 
417  if (wantE1)
418  {
419  this->operator[](nDirs++) = e1Dirs;
420  }
421  }
422  if (wantE2)
423  {
424  tmp<vectorField> e2Dirs = e3Dirs ^ e1Dirs;
425 
426  this->operator[](nDirs++) = e2Dirs;
427  }
428  }
429  else if (coordSystem == "fieldBased")
430  {
431  forAll(wantedDirs, i)
432  {
433  operator[](nDirs++) =
435  (
436  IOobject
437  (
438  mesh.instance()/wantedDirs[i],
439  mesh,
442  )
443  );
444  }
445  }
446  else
447  {
449  << "Unknown coordinate system "
450  << coordSystem << endl
451  << "Known types are global, patchLocal and fieldBased"
452  << exit(FatalError);
453  }
454 }
455 
456 
457 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
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
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
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:323
Cell-face mesh analysis engine.
Definition: primitiveMesh.H:74
Output to file stream.
Definition: OFstream.H:82
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
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:63
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:982
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
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
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1394
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:84
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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
const fileName & instance() const
Definition: IOobject.H:390
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:309
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
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
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:278
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
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:844