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-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 
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  << "More than six unsigned transforms"
185  << " detected:" << nl << transforms_
186  << exit(FatalError);
187  }
188  transforms_[nextTrans] = transform;
189  maxTol[nextTrans++] = cpp.matchTolerance();
190  }
191  }
192  }
193  }
194  else if (!cpp.parallel())
195  {
196  const tensorField& transTensors = cpp.reverseT();
197 
198  forAll(transTensors, tTI)
199  {
200  const tensor& transT = transTensors[tTI];
201 
202  if (mag(transT - I) > SMALL)
203  {
204  vectorTensorTransform transform(transT);
205 
206  if
207  (
208  matchTransform
209  (
210  transforms_,
211  dummyMatch,
212  transform,
213  cpp.matchTolerance(),
214  false
215  ) == 0
216  )
217  {
218  if (nextTrans == 6)
219  {
221  << "More than six unsigned transforms"
222  << " detected:" << nl << transforms_
223  << exit(FatalError);
224  }
225  transforms_[nextTrans] = transform;
226  maxTol[nextTrans++] = cpp.matchTolerance();
227  }
228  }
229  }
230  }
231  }
232  }
233 
234 
235  // Collect transforms on master
236 
237  List<List<vectorTensorTransform>> allTransforms(Pstream::nProcs());
238  allTransforms[Pstream::myProcNo()] = transforms_;
239  Pstream::gatherList(allTransforms);
240 
241  // Collect matching tolerance on master
242  List<scalarField> allTols(Pstream::nProcs());
243  allTols[Pstream::myProcNo()] = maxTol;
244  Pstream::gatherList(allTols);
245 
246  if (Pstream::master())
247  {
248  transforms_ = List<vectorTensorTransform>(3);
249 
250  label nextTrans = 0;
251 
252  forAll(allTransforms, proci)
253  {
254  const List<vectorTensorTransform>& procTransVecs =
255  allTransforms[proci];
256 
257  forAll(procTransVecs, pSVI)
258  {
259  const vectorTensorTransform& transform = procTransVecs[pSVI];
260 
261  if (mag(transform.t()) > SMALL || transform.hasR())
262  {
263  if
264  (
265  matchTransform
266  (
267  transforms_,
268  dummyMatch,
269  transform,
270  allTols[proci][pSVI],
271  true
272  ) == 0
273  )
274  {
275  transforms_[nextTrans++] = transform;
276  }
277 
278  if (nextTrans > 3)
279  {
281  << "More than three independent basic "
282  << "transforms detected:" << nl
283  << allTransforms
284  << transforms_
285  << exit(FatalError);
286  }
287  }
288  }
289  }
290 
291  transforms_.setSize(nextTrans);
292  }
293 
294  Pstream::scatter(transforms_);
295 
296  if (transforms_.size() > 3)
297  {
299  << "More than three independent basic "
300  << "transforms detected:" << nl
301  << transforms_ << nl
302  << "This is not a space filling tiling and will probably"
303  << " give problems for e.g. lagrangian tracking or interpolation"
304  << endl;
305  }
306 }
307 
308 
309 void Foam::globalIndexAndTransform::determineTransformPermutations()
310 {
311  label nTransformPermutations = pow(label(3), transforms_.size());
312 
313  transformPermutations_.setSize(nTransformPermutations);
314 
315  forAll(transformPermutations_, tPI)
316  {
317  vectorTensorTransform transform;
318 
319  label transformIndex = tPI;
320 
321  // Invert the ternary index encoding using repeated division by
322  // three
323 
324  forAll(transforms_, b)
325  {
326  const label w = (transformIndex % 3) - 1;
327 
328  transformIndex /= 3;
329 
330  if (w > 0)
331  {
332  transform &= transforms_[b];
333  }
334  else if (w < 0)
335  {
336  transform &= inv(transforms_[b]);
337  }
338  }
339 
340  transformPermutations_[tPI] = transform;
341  }
342 
343 
344  // Encode index for 0 sign
345  labelList permutationIndices(nIndependentTransforms(), 0);
346  nullTransformIndex_ = encodeTransformIndex(permutationIndices);
347 }
348 
349 
350 void Foam::globalIndexAndTransform::determinePatchTransformSign()
351 {
352  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
353 
354  patchTransformSign_.setSize(patches.size(), Pair<label>(-1, 0));
355 
356  label matchTransI = -1;
357 
358  forAll(patches, patchi)
359  {
360  const polyPatch& pp = patches[patchi];
361 
362  // Pout<< nl << patchi << " " << pp.name() << endl;
363 
364  // Note: special check for unordered cyclics. These are in fact
365  // transform bcs and should probably be split off.
366  if
367  (
368  isA<coupledPolyPatch>(pp)
369  && !(
370  isA<cyclicPolyPatch>(pp)
371  && (
372  refCast<const cyclicPolyPatch>(pp).transform()
374  )
375  )
376  )
377  {
378  const coupledPolyPatch& cpp =
379  refCast<const coupledPolyPatch>(pp);
380 
381  if (cpp.separated())
382  {
383  const vectorField& sepVecs = cpp.separation();
384 
385  // Pout<< "sepVecs " << sepVecs << endl;
386 
387  // This loop is implicitly expecting only a single
388  // value for separation()
389  forAll(sepVecs, sVI)
390  {
391  const vector& sepVec = sepVecs[sVI];
392 
393  if (mag(sepVec) > SMALL)
394  {
395  vectorTensorTransform t(sepVec);
396 
397  label sign = matchTransform
398  (
399  transforms_,
400  matchTransI,
401  t,
402  cpp.matchTolerance(),
403  true
404  );
405 
406  // Pout<< sign << " " << matchTransI << endl;
407 
408  // List<label> permutation(transforms_.size(), 0);
409 
410  // permutation[matchTransI] = sign;
411 
412  // Pout<< encodeTransformIndex(permutation) << nl
413  // << transformPermutations_
414  // [
415  // encodeTransformIndex(permutation)
416  // ]
417  // << endl;
418 
419  patchTransformSign_[patchi] =
420  Pair<label>(matchTransI, sign);
421  }
422  }
423 
424  }
425  else if (!cpp.parallel())
426  {
427  const tensorField& transTensors = cpp.reverseT();
428 
429  // Pout<< "transTensors " << transTensors << endl;
430 
431  // This loop is implicitly expecting only a single
432  // value for reverseT()
433  forAll(transTensors, tTI)
434  {
435  const tensor& transT = transTensors[tTI];
436 
437  if (mag(transT - I) > SMALL)
438  {
439  vectorTensorTransform t(transT);
440 
441  label sign = matchTransform
442  (
443  transforms_,
444  matchTransI,
445  t,
446  cpp.matchTolerance(),
447  true
448  );
449 
450  // Pout<< sign << " " << matchTransI << endl;
451 
452  // List<label> permutation(transforms_.size(), 0);
453 
454  // permutation[matchTransI] = sign;
455 
456  // Pout<< encodeTransformIndex(permutation) << nl
457  // << transformPermutations_
458  // [
459  // encodeTransformIndex(permutation)
460  // ]
461  // << endl;
462 
463  patchTransformSign_[patchi] =
464  Pair<label>(matchTransI, sign);
465  }
466  }
467  }
468  }
469  }
470 
471  // Pout<< patchTransformSign_ << endl;
472 }
473 
474 
475 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
476 
477 Foam::globalIndexAndTransform::globalIndexAndTransform
478 (
479  const polyMesh& mesh
480 )
481 :
482  mesh_(mesh),
483  transforms_(),
484  transformPermutations_(),
485  patchTransformSign_()
486 {
487  determineTransforms();
488 
489  determineTransformPermutations();
490 
491  determinePatchTransformSign();
492 
493  if (debug && transforms_.size() > 0)
494  {
495  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
496 
497  Info<< "Determined global transforms :" << endl;
498  Info<< "\t\ttranslation\trotation" << endl;
499  forAll(transforms_, i)
500  {
501  Info<< '\t' << i << '\t';
502  const vectorTensorTransform& trafo = transforms_[i];
503  if (trafo.hasR())
504  {
505  Info<< trafo.t() << '\t' << trafo.R();
506  }
507  else
508  {
509  Info<< trafo.t() << '\t' << "---";
510  }
511  Info<< endl;
512  }
513  Info<< endl;
514 
515 
516  Info<< "\tpatch\ttransform\tsign" << endl;
517  forAll(patchTransformSign_, patchi)
518  {
519  if (patchTransformSign_[patchi].first() != -1)
520  {
521  Info<< '\t' << patches[patchi].name()
522  << '\t' << patchTransformSign_[patchi].first()
523  << '\t' << patchTransformSign_[patchi].second()
524  << endl;
525  }
526  }
527  Info<< endl;
528 
529 
530  Info<< "Permutations of transformations:" << endl
531  << "\t\ttranslation\trotation" << endl;
532  forAll(transformPermutations_, i)
533  {
534  Info<< '\t' << i << '\t';
535  const vectorTensorTransform& trafo = transformPermutations_[i];
536  if (trafo.hasR())
537  {
538  Info<< trafo.t() << '\t' << trafo.R();
539  }
540  else
541  {
542  Info<< trafo.t() << '\t' << "---";
543  }
544  Info<< endl;
545  }
546  Info<< "nullTransformIndex:" << nullTransformIndex() << endl
547  << endl;
548  }
549 }
550 
551 
552 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
553 
555 {}
556 
557 
558 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
#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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Vector-tensor class used to perform translations and rotations in 3D space.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:417
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:411
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static const Identity< scalar > I
Definition: Identity.H:93
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
Foam::polyBoundaryMesh.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
static const char nl
Definition: Ostream.H:262
defineTypeNameAndDebug(combustionModel, 0)
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:399
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
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
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
const word & name() const
Return name.
Definition: IOobject.H:260
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465