globalIndexAndTransform.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-2020 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 #include "DynamicField.H"
29 #include "globalMeshData.H"
30 
31 // * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(globalIndexAndTransform, 0);
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 Foam::label Foam::globalIndexAndTransform::matchTransform
42 (
43  const List<transformer>& refTransforms,
44  label& matchedRefTransformI,
45  const transformer& testTransform,
46  scalar tolerance,
47  bool checkBothSigns
48 ) const
49 {
50  matchedRefTransformI = -1;
51 
52  forAll(refTransforms, i)
53  {
54  const transformer& refTransform = refTransforms[i];
55 
56  const scalar maxVectorMag = sqrt
57  (
58  max(magSqr(testTransform.t()), magSqr(refTransform.t()))
59  );
60 
61  // Test the difference between vector parts to see if it is
62  // less than tolerance times the larger vector part magnitude.
63 
64  scalar vectorDiff =
65  mag(refTransform.t() - testTransform.t())
66  /(maxVectorMag + vSmall)
67  /tolerance;
68 
69  // Test the difference between tensor parts to see if it is
70  // less than the tolerance. sqrt(3.0) factor used to scale
71  // differences as this is magnitude of a rotation tensor. If
72  // neither transform has a rotation, then the test is not
73  // necessary.
74 
75  scalar tensorDiff = 0;
76 
77  if (refTransform.transforms() || testTransform.transforms())
78  {
79  tensorDiff =
80  mag(refTransform.T() - testTransform.T())
81  /sqrt(3.0)
82  /tolerance;
83  }
84 
85  // ...Diff result is < 1 if the test part matches the ref part
86  // within tolerance
87 
88  if (vectorDiff < 1 && tensorDiff < 1)
89  {
90  matchedRefTransformI = i;
91 
92  return +1;
93  }
94 
95  if (checkBothSigns)
96  {
97  // Test the inverse transform differences too
98 
99  const transformer testInvTransform = inv(testTransform);
100 
101  vectorDiff =
102  mag(refTransform.t() - testInvTransform.t())
103  /(maxVectorMag + vSmall)
104  /tolerance;
105 
106  tensorDiff = 0;
107 
108  if (refTransform.transforms() || testTransform.transforms())
109  {
110  tensorDiff =
111  mag(refTransform.T() - testInvTransform.T())
112  /sqrt(3.0)
113  /tolerance;
114  }
115 
116  if (vectorDiff < 1 && tensorDiff < 1)
117  {
118  matchedRefTransformI = i;
119 
120  return -1;
121  }
122  }
123  }
124 
125  return 0;
126 }
127 
128 
129 void Foam::globalIndexAndTransform::determineTransforms()
130 {
131  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
132 
133  DynamicList<transformer> localTransforms;
134  DynamicField<scalar> localTols;
135 
136  label dummyMatch = -1;
137 
138  forAll(patches, patchi)
139  {
140  const polyPatch& pp = patches[patchi];
141 
142  if (isA<coupledPolyPatch>(pp))
143  {
144  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
145 
146  const transformer transform(inv(cpp.transform()));
147 
148  if (transform.transformsPosition())
149  {
150  if
151  (
152  matchTransform
153  (
154  localTransforms,
155  dummyMatch,
156  transform,
157  cpp.matchTolerance(),
158  false
159  ) == 0
160  )
161  {
162  localTransforms.append(transform);
163  localTols.append(cpp.matchTolerance());
164  }
165  }
166  }
167  }
168 
169 
170  // Collect transforms on master
171  List<List<transformer>> allTransforms(Pstream::nProcs());
172  allTransforms[Pstream::myProcNo()] = localTransforms;
173  Pstream::gatherList(allTransforms);
174 
175  // Collect matching tolerance on master
176  List<scalarField> allTols(Pstream::nProcs());
177  allTols[Pstream::myProcNo()] = localTols;
178  Pstream::gatherList(allTols);
179 
180  if (Pstream::master())
181  {
182  localTransforms.clear();
183 
184  forAll(allTransforms, proci)
185  {
186  const List<transformer>& procTransVecs =
187  allTransforms[proci];
188 
189  forAll(procTransVecs, pSVI)
190  {
191  const transformer& transform = procTransVecs[pSVI];
192 
193  if (transform.transformsPosition())
194  {
195  if
196  (
197  matchTransform
198  (
199  localTransforms,
200  dummyMatch,
201  transform,
202  allTols[proci][pSVI],
203  true
204  ) == 0
205  )
206  {
207  localTransforms.append(transform);
208  }
209  }
210  }
211  }
212  }
213 
214  transforms_.transfer(localTransforms);
215  Pstream::scatter(transforms_);
216 }
217 
218 
219 void Foam::globalIndexAndTransform::determineTransformPermutations()
220 {
221  label nTransformPermutations = pow(label(3), transforms_.size());
222 
223  transformPermutations_.setSize(nTransformPermutations);
224 
225  forAll(transformPermutations_, tPI)
226  {
227  transformer transform;
228 
229  label transformIndex = tPI;
230 
231  // Invert the ternary index encoding using repeated division by
232  // three
233 
234  forAll(transforms_, b)
235  {
236  const label w = (transformIndex % 3) - 1;
237 
238  transformIndex /= 3;
239 
240  if (w > 0)
241  {
242  transform = transforms_[b] & transform;
243  }
244  else if (w < 0)
245  {
246  transform = inv(transforms_[b]) & transform;
247  }
248  }
249 
250  transformPermutations_[tPI] = transform;
251  }
252 
253 
254  // Encode index for 0 sign
255  labelList permutationIndices(nIndependentTransforms(), 0);
256  nullTransformIndex_ = encodeTransformIndex(permutationIndices);
257 }
258 
259 
260 void Foam::globalIndexAndTransform::determinePatchTransformSign()
261 {
262  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
263 
264  patchTransformSign_.setSize(patches.size(), labelPair(-1, 0));
265 
266  forAll(patches, patchi)
267  {
268  const polyPatch& pp = patches[patchi];
269 
270  if (isA<coupledPolyPatch>(pp))
271  {
272  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
273 
274  const transformer transform(inv(cpp.transform()));
275 
276  if (transform.transformsPosition())
277  {
278  label matchTransI;
279  label sign = matchTransform
280  (
281  transforms_,
282  matchTransI,
283  transform,
284  cpp.matchTolerance(),
285  true
286  );
287  patchTransformSign_[patchi] = labelPair(matchTransI, sign);
288  }
289  }
290  }
291 }
292 
293 
294 bool Foam::globalIndexAndTransform::uniqueTransform
295 (
296  const point& pt,
297  labelPairList& trafos,
298  const label patchi,
299  const labelPair& patchTrafo
300 ) const
301 {
302  if (findIndex(trafos, patchTrafo) == -1)
303  {
304  // New transform. Check if already have 3
305  if (trafos.size() == 3)
306  {
307  if (patchi > -1)
308  {
310  << "Point " << pt
311  << " is on patch " << mesh_.boundaryMesh()[patchi].name();
312  }
313  else
314  {
316  << "Point " << pt << " is on a coupled patch";
317  }
318  Warning
319  << " with transformation " << patchTrafo
320  << " but also on 3 other patches with transforms "
321  << trafos << nl
322  << "This is not a space filling tiling and might"
323  << " indicate a setup problem and give problems"
324  << " for e.g. lagrangian tracking or interpolation" << endl;
325 
326  // Already warned so no need to extend more
327  trafos.clear();
328  return false;
329  }
330 
331  return true;
332  }
333  else
334  {
335  return false;
336  }
337 }
338 
339 
340 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
341 
343 :
344  mesh_(mesh),
345  transforms_(),
346  transformPermutations_(),
347  patchTransformSign_()
348 {
349  determineTransforms();
350 
351  determineTransformPermutations();
352 
353  determinePatchTransformSign();
354 
355  if (debug && transforms_.size() > 0)
356  {
357  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
358 
359  Info<< "Determined global transforms :" << endl;
360  Info<< "\t\ttranslation\trotation" << endl;
361  forAll(transforms_, i)
362  {
363  Info<< '\t' << i << '\t';
364  const transformer& trafo = transforms_[i];
365  Info<< trafo.t() << '\t' << trafo.T() << endl;
366  }
367  Info<< endl;
368 
369 
370  Info<< "\tpatch\ttransform\tsign" << endl;
371  forAll(patchTransformSign_, patchi)
372  {
373  if (patchTransformSign_[patchi].first() != -1)
374  {
375  Info<< '\t' << patches[patchi].name()
376  << '\t' << patchTransformSign_[patchi].first()
377  << '\t' << patchTransformSign_[patchi].second()
378  << endl;
379  }
380  }
381  Info<< endl;
382 
383 
384  Info<< "Permutations of transformations:" << endl
385  << "\t\ttranslation\trotation" << endl;
386  forAll(transformPermutations_, i)
387  {
388  Info<< '\t' << i << '\t';
389  const transformer& trafo = transformPermutations_[i];
390  Info<< trafo.t() << '\t' << trafo.T() << endl;
391  }
392  Info<< "nullTransformIndex:" << nullTransformIndex() << endl
393  << endl;
394  }
395 
396 
397  if (transforms_.size() > 0)
398  {
399  // Check that the transforms are space filling : any point
400  // can only use up to three transforms
401 
402  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
403 
404 
405  // 1. Collect transform&sign per point and do local check
406 
407  List<labelPairList> pointToTrafos(mesh_.nPoints());
408 
409  forAll(patches, patchi)
410  {
411  const polyPatch& pp = patches[patchi];
412 
413  const labelPair& transSign = patchTransformSign_[patchi];
414 
415  if (transSign.first() > -1)
416  {
417  const labelList& mp = pp.meshPoints();
418  forAll(mp, i)
419  {
420  labelPairList& trafos = pointToTrafos[mp[i]];
421 
422  bool newTransform = uniqueTransform
423  (
424  mesh_.points()[mp[i]],
425  trafos,
426  patchi,
427  transSign
428  );
429 
430  if (newTransform)
431  {
432  trafos.append(transSign);
433  }
434  }
435  }
436  }
437 
438 
439  // Synchronise across collocated (= untransformed) points
440  // TBD: there is a big problem in that globalMeshData uses
441  // globalIndexAndTransform. Triggers recursion.
442  if (false)
443  {
444  const globalMeshData& gmd = mesh_.globalData();
445  const indirectPrimitivePatch& cpp = gmd.coupledPatch();
446  const labelList& meshPoints = cpp.meshPoints();
447  const mapDistribute& slavesMap = gmd.globalCoPointSlavesMap();
448  const labelListList& slaves = gmd.globalCoPointSlaves();
449 
450  List<labelPairList> elems(slavesMap.constructSize());
451  forAll(meshPoints, i)
452  {
453  elems[i] = pointToTrafos[meshPoints[i]];
454  }
455 
456  // Pull slave data onto master. No need to update transformed slots.
457  slavesMap.distribute(elems, false);
458 
459  // Combine master data with slave data
460  forAll(slaves, i)
461  {
462  labelPairList& trafos = elems[i];
463 
464  const labelList& slavePoints = slaves[i];
465 
466  // Combine master with untransformed slave data
467  forAll(slavePoints, j)
468  {
469  const labelPairList& slaveTrafos = elems[slavePoints[j]];
470 
471  forAll(slaveTrafos, slaveI)
472  {
473  bool newTransform = uniqueTransform
474  (
475  mesh_.points()[meshPoints[i]],
476  trafos,
477  -1,
478  slaveTrafos[slaveI]
479  );
480 
481  if (newTransform)
482  {
483  trafos.append(slaveTrafos[slaveI]);
484  }
485  }
486  }
487  }
488  }
489  }
490 }
491 
492 
493 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
dimensionedScalar sign(const dimensionedScalar &ds)
const mapDistribute & globalCoPointSlavesMap() const
#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.
Definition: IOobject.H:303
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space...
Definition: transformer.H:83
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute data using default commsType.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
globalIndexAndTransform(const polyMesh &mesh)
Construct from components.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
const labelList & meshPoints() const
Return labelList of mesh points in patch. They are constructed.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
const vector & t() const
Return the translation vector.
Definition: transformerI.H:82
A list of faces which address into the list of points.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
dynamicFvMesh & mesh
label nullTransformIndex() const
Return the transformIndex (index in transformPermutations)
const dimensionedScalar mp
Proton mass.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1394
List< label > labelList
A List of labels.
Definition: labelList.H:56
List< labelPair > labelPairList
List of labelPairs.
Definition: labelPair.H:49
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 > &)
const tensor & T() const
Return the transformation tensor.
Definition: transformerI.H:94
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
messageStream Warning
const indirectPrimitivePatch & coupledPatch() const
Return patch of all coupled faces.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
label constructSize() const
Constructed data size.
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
label patchi
Class containing processor-to-processor mapping information.
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
const Type & first() const
Return first.
Definition: Pair.H:98
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477
const labelListList & globalCoPointSlaves() const