globalIndexAndTransform.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-2012 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 
27 #include "cyclicPolyPatch.H"
28 
29 // * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 defineTypeNameAndDebug(globalIndexAndTransform, 0);
34 const label globalIndexAndTransform::base_ = 32;
35 }
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 Foam::label Foam::globalIndexAndTransform::matchTransform
41 (
42  const List<vectorTensorTransform>& refTransforms,
43  label& matchedRefTransformI,
44  const vectorTensorTransform& testTransform,
45  scalar tolerance,
46  bool checkBothSigns
47 ) const
48 {
49  matchedRefTransformI = -1;
50 
51  forAll(refTransforms, i)
52  {
53  const vectorTensorTransform& refTransform = refTransforms[i];
54 
55  scalar maxVectorMag = sqrt
56  (
57  max(magSqr(testTransform.t()), magSqr(refTransform.t()))
58  );
59 
60  // Test the difference between vector parts to see if it is
61  // less than tolerance times the larger vector part magnitude.
62 
63  scalar vectorDiff =
64  mag(refTransform.t() - testTransform.t())
65  /(maxVectorMag + VSMALL)
66  /tolerance;
67 
68  // Test the difference between tensor parts to see if it is
69  // less than the tolerance. sqrt(3.0) factor used to scale
70  // differnces as this is magnitude of a rotation tensor. If
71  // neither transform has a rotation, then the test is not
72  // necessary.
73 
74  scalar tensorDiff = 0;
75 
76  if (refTransform.hasR() || testTransform.hasR())
77  {
78  tensorDiff =
79  mag(refTransform.R() - testTransform.R())
80  /sqrt(3.0)
81  /tolerance;
82  }
83 
84  // ...Diff result is < 1 if the test part matches the ref part
85  // within tolerance
86 
87  if (vectorDiff < 1 && tensorDiff < 1)
88  {
89  matchedRefTransformI = i;
90 
91  return +1;
92  }
93 
94  if (checkBothSigns)
95  {
96  // Test the inverse transform differences too
97 
98  vectorDiff =
99  mag(refTransform.t() + testTransform.t())
100  /(maxVectorMag + VSMALL)
101  /tolerance;
102 
103  tensorDiff = 0;
104 
105  if (refTransform.hasR() || testTransform.hasR())
106  {
107  tensorDiff =
108  mag(refTransform.R() - testTransform.R().T())
109  /sqrt(3.0)
110  /tolerance;
111  }
112 
113  if (vectorDiff < 1 && tensorDiff < 1)
114  {
115  matchedRefTransformI = i;
116 
117  return -1;
118  }
119  }
120  }
121 
122  return 0;
123 }
124 
125 
126 void Foam::globalIndexAndTransform::determineTransforms()
127 {
128  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
129 
130  transforms_ = List<vectorTensorTransform>(6);
131  scalarField maxTol(6);
132 
133  label nextTrans = 0;
134 
135  label dummyMatch = -1;
136 
137  forAll(patches, patchI)
138  {
139  const polyPatch& pp = patches[patchI];
140 
141  // Note: special check for unordered cyclics. These are in fact
142  // transform bcs and should probably be split off.
143  if
144  (
145  isA<coupledPolyPatch>(pp)
146  && !(
147  isA<cyclicPolyPatch>(pp)
148  && (
149  refCast<const cyclicPolyPatch>(pp).transform()
151  )
152  )
153  )
154  {
155  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
156 
157  if (cpp.separated())
158  {
159  const vectorField& sepVecs = cpp.separation();
160 
161  forAll(sepVecs, sVI)
162  {
163  const vector& sepVec = sepVecs[sVI];
164 
165  if (mag(sepVec) > SMALL)
166  {
167  vectorTensorTransform transform(sepVec);
168 
169  if
170  (
171  matchTransform
172  (
173  transforms_,
174  dummyMatch,
175  transform,
176  cpp.matchTolerance(),
177  false
178  ) == 0
179  )
180  {
181  if (nextTrans == 6)
182  {
184  (
185  "void Foam::globalIndexAndTransform::"
186  "determineTransforms()"
187  ) << "More than six unsigned transforms"
188  << " detected:" << nl << transforms_
189  << exit(FatalError);
190  }
191  transforms_[nextTrans] = transform;
192  maxTol[nextTrans++] = cpp.matchTolerance();
193  }
194  }
195  }
196  }
197  else if (!cpp.parallel())
198  {
199  const tensorField& transTensors = cpp.reverseT();
200 
201  forAll(transTensors, tTI)
202  {
203  const tensor& transT = transTensors[tTI];
204 
205  if (mag(transT - I) > SMALL)
206  {
207  vectorTensorTransform transform(transT);
208 
209  if
210  (
211  matchTransform
212  (
213  transforms_,
214  dummyMatch,
215  transform,
216  cpp.matchTolerance(),
217  false
218  ) == 0
219  )
220  {
221  if (nextTrans == 6)
222  {
224  (
225  "void Foam::globalIndexAndTransform::"
226  "determineTransforms()"
227  ) << "More than six unsigned transforms"
228  << " detected:" << nl << transforms_
229  << exit(FatalError);
230  }
231  transforms_[nextTrans] = transform;
232  maxTol[nextTrans++] = cpp.matchTolerance();
233  }
234  }
235  }
236  }
237  }
238  }
239 
240 
241  // Collect transforms on master
242 
243  List<List<vectorTensorTransform> > allTransforms(Pstream::nProcs());
244  allTransforms[Pstream::myProcNo()] = transforms_;
245  Pstream::gatherList(allTransforms);
246 
247  // Collect matching tolerance on master
248  List<scalarField> allTols(Pstream::nProcs());
249  allTols[Pstream::myProcNo()] = maxTol;
250  Pstream::gatherList(allTols);
251 
252  if (Pstream::master())
253  {
254  transforms_ = List<vectorTensorTransform>(3);
255 
256  label nextTrans = 0;
257 
258  forAll(allTransforms, procI)
259  {
260  const List<vectorTensorTransform>& procTransVecs =
261  allTransforms[procI];
262 
263  forAll(procTransVecs, pSVI)
264  {
265  const vectorTensorTransform& transform = procTransVecs[pSVI];
266 
267  if (mag(transform.t()) > SMALL || transform.hasR())
268  {
269  if
270  (
271  matchTransform
272  (
273  transforms_,
274  dummyMatch,
275  transform,
276  allTols[procI][pSVI],
277  true
278  ) == 0
279  )
280  {
281  transforms_[nextTrans++] = transform;
282  }
283 
284  if (nextTrans > 3)
285  {
287  (
288  "void Foam::globalIndexAndTransform::"
289  "determineTransforms()"
290  )
291  << "More than three independent basic "
292  << "transforms detected:" << nl
293  << allTransforms
294  << transforms_
295  << exit(FatalError);
296  }
297  }
298  }
299  }
300 
301  transforms_.setSize(nextTrans);
302  }
303 
304  Pstream::scatter(transforms_);
305 
306  if (transforms_.size() > 3)
307  {
308  WarningIn
309  (
310  "void globalIndexAndTransform::determineTransforms()"
311  ) << "More than three independent basic "
312  << "transforms detected:" << nl
313  << transforms_ << nl
314  << "This is not a space filling tiling and will probably"
315  << " give problems for e.g. lagrangian tracking or interpolation"
316  << endl;
317  }
318 }
319 
320 
321 void Foam::globalIndexAndTransform::determineTransformPermutations()
322 {
323  label nTransformPermutations = pow(label(3), transforms_.size());
324 
325  transformPermutations_.setSize(nTransformPermutations);
326 
327  forAll(transformPermutations_, tPI)
328  {
329  vectorTensorTransform transform;
330 
331  label transformIndex = tPI;
332 
333  // Invert the ternary index encoding using repeated division by
334  // three
335 
336  forAll(transforms_, b)
337  {
338  const label w = (transformIndex % 3) - 1;
339 
340  transformIndex /= 3;
341 
342  if (w > 0)
343  {
344  transform &= transforms_[b];
345  }
346  else if (w < 0)
347  {
348  transform &= inv(transforms_[b]);
349  }
350  }
351 
352  transformPermutations_[tPI] = transform;
353  }
354 
355 
356  // Encode index for 0 sign
357  labelList permutationIndices(nIndependentTransforms(), 0);
358  nullTransformIndex_ = encodeTransformIndex(permutationIndices);
359 }
360 
361 
362 void Foam::globalIndexAndTransform::determinePatchTransformSign()
363 {
364  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
365 
366  patchTransformSign_.setSize(patches.size(), Pair<label>(-1, 0));
367 
368  label matchTransI = -1;
369 
370  forAll(patches, patchI)
371  {
372  const polyPatch& pp = patches[patchI];
373 
374  // Pout<< nl << patchI << " " << pp.name() << endl;
375 
376  // Note: special check for unordered cyclics. These are in fact
377  // transform bcs and should probably be split off.
378  if
379  (
380  isA<coupledPolyPatch>(pp)
381  && !(
382  isA<cyclicPolyPatch>(pp)
383  && (
384  refCast<const cyclicPolyPatch>(pp).transform()
386  )
387  )
388  )
389  {
390  const coupledPolyPatch& cpp =
391  refCast<const coupledPolyPatch>(pp);
392 
393  if (cpp.separated())
394  {
395  const vectorField& sepVecs = cpp.separation();
396 
397  // Pout<< "sepVecs " << sepVecs << endl;
398 
399  // This loop is implicitly expecting only a single
400  // value for separation()
401  forAll(sepVecs, sVI)
402  {
403  const vector& sepVec = sepVecs[sVI];
404 
405  if (mag(sepVec) > SMALL)
406  {
407  vectorTensorTransform t(sepVec);
408 
409  label sign = matchTransform
410  (
411  transforms_,
412  matchTransI,
413  t,
414  cpp.matchTolerance(),
415  true
416  );
417 
418  // Pout<< sign << " " << matchTransI << endl;
419 
420  // List<label> permutation(transforms_.size(), 0);
421 
422  // permutation[matchTransI] = sign;
423 
424  // Pout<< encodeTransformIndex(permutation) << nl
425  // << transformPermutations_
426  // [
427  // encodeTransformIndex(permutation)
428  // ]
429  // << endl;
430 
431  patchTransformSign_[patchI] =
432  Pair<label>(matchTransI, sign);
433  }
434  }
435 
436  }
437  else if (!cpp.parallel())
438  {
439  const tensorField& transTensors = cpp.reverseT();
440 
441  // Pout<< "transTensors " << transTensors << endl;
442 
443  // This loop is implicitly expecting only a single
444  // value for reverseT()
445  forAll(transTensors, tTI)
446  {
447  const tensor& transT = transTensors[tTI];
448 
449  if (mag(transT - I) > SMALL)
450  {
451  vectorTensorTransform t(transT);
452 
453  label sign = matchTransform
454  (
455  transforms_,
456  matchTransI,
457  t,
458  cpp.matchTolerance(),
459  true
460  );
461 
462  // Pout<< sign << " " << matchTransI << endl;
463 
464  // List<label> permutation(transforms_.size(), 0);
465 
466  // permutation[matchTransI] = sign;
467 
468  // Pout<< encodeTransformIndex(permutation) << nl
469  // << transformPermutations_
470  // [
471  // encodeTransformIndex(permutation)
472  // ]
473  // << endl;
474 
475  patchTransformSign_[patchI] =
476  Pair<label>(matchTransI, sign);
477  }
478  }
479  }
480  }
481  }
482 
483  // Pout<< patchTransformSign_ << endl;
484 }
485 
486 
487 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
488 
489 Foam::globalIndexAndTransform::globalIndexAndTransform
490 (
491  const polyMesh& mesh
492 )
493 :
494  mesh_(mesh),
495  transforms_(),
496  transformPermutations_(),
497  patchTransformSign_()
498 {
499  determineTransforms();
500 
501  determineTransformPermutations();
502 
503  determinePatchTransformSign();
504 
505  if (debug && transforms_.size() > 0)
506  {
507  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
508 
509  Info<< "Determined global transforms :" << endl;
510  Info<< "\t\ttranslation\trotation" << endl;
511  forAll(transforms_, i)
512  {
513  Info<< '\t' << i << '\t';
514  const vectorTensorTransform& trafo = transforms_[i];
515  if (trafo.hasR())
516  {
517  Info<< trafo.t() << '\t' << trafo.R();
518  }
519  else
520  {
521  Info<< trafo.t() << '\t' << "---";
522  }
523  Info<< endl;
524  }
525  Info<< endl;
526 
527 
528  Info<< "\tpatch\ttransform\tsign" << endl;
529  forAll(patchTransformSign_, patchI)
530  {
531  if (patchTransformSign_[patchI].first() != -1)
532  {
533  Info<< '\t' << patches[patchI].name()
534  << '\t' << patchTransformSign_[patchI].first()
535  << '\t' << patchTransformSign_[patchI].second()
536  << endl;
537  }
538  }
539  Info<< endl;
540 
541 
542  Info<< "Permutations of transformations:" << endl
543  << "\t\ttranslation\trotation" << endl;
544  forAll(transformPermutations_, i)
545  {
546  Info<< '\t' << i << '\t';
547  const vectorTensorTransform& trafo = transformPermutations_[i];
548  if (trafo.hasR())
549  {
550  Info<< trafo.t() << '\t' << trafo.R();
551  }
552  else
553  {
554  Info<< trafo.t() << '\t' << "---";
555  }
556  Info<< endl;
557  }
558  Info<< "nullTransformIndex:" << nullTransformIndex() << endl
559  << endl;
560  }
561 }
562 
563 
564 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
565 
567 {}
568 
569 
570 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
dimensioned< scalar > mag(const dimensioned< Type > &)
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const sphericalTensor I(1)
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
dimensionedScalar sign(const dimensionedScalar &ds)
messageStream Info
Namespace for OpenFOAM.
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:386
static const char nl
Definition: Ostream.H:260
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
Foam::polyBoundaryMesh.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define forAll(list, i)
Definition: UList.H:421
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
const word & name() const
Return name.
Definition: IOobject.H:260
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Vector-tensor class used to perform translations and rotations in 3D space.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:398
defineTypeNameAndDebug(combustionModel, 0)