cyclicTransform.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) 2020-2024 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 "cyclicTransform.H"
27 #include "unitConversion.H"
28 #include "IOmanip.H"
29 #include "stringOps.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 
36 template<class Type>
37 Type sum(const Type& x, const bool global)
38 {
39  return global ? returnReduce(x, sumOp<Type>()) : x;
40 }
41 
42 template<class Type>
43 Type sum(const Field<Type>& x, const bool global)
44 {
45  return global ? gSum(x) : sum(x);
46 }
47 
48 template<class Type>
49 Type sum(const tmp<Field<Type>>& x, const bool global)
50 {
51  const Type s = sum(x(), global);
52  x.clear();
53  return s;
54 }
55 
56 }
57 
58 
59 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
60 
61 namespace Foam
62 {
63  template<>
65  {
66  "unspecified",
67  "none",
68  "rotational",
69  "translational"
70  };
71 
74 
76 }
77 
78 
79 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
80 
81 void Foam::cyclicTransform::update()
82 {
83  if (!transformComplete_)
84  {
85  return;
86  }
87 
88  switch (transformType_)
89  {
90  case UNSPECIFIED:
91  break;
92 
93  case NONE:
94  transform_ = transformer();
95  break;
96 
97  case ROTATIONAL:
98  if (rotationAngle_ == 0)
99  {
100  transform_ = transformer();
101  }
102  else
103  {
104  const tensor R =
105  quaternion(rotationAxis_, rotationAngle_).R();
106 
107  if (mag(rotationCentre_) == 0)
108  {
109  transform_ = transformer::rotation(R);
110  }
111  else
112  {
113  transform_ =
114  transformer::translation(rotationCentre_)
116  & transformer::translation(- rotationCentre_);
117  }
118  }
119  break;
120 
121  case TRANSLATIONAL:
122  if (mag(separation_) == 0)
123  {
124  transform_ = transformer();
125  }
126  else
127  {
128  transform_ = transformer::translation(separation_);
129  }
130  break;
131  }
132 }
133 
134 
135 bool Foam::cyclicTransform::set
136 (
137  const cyclicTransform& t,
138  const scalar lengthScale,
139  const scalar matchTolerance
140 )
141 {
142  // If the supplied transform is unspecified then there is nothing to do
143  if (t.transformType_ == UNSPECIFIED)
144  {
145  return true;
146  }
147 
148  // If this transform is specified then we need to check that it is
149  // compatible with the supplied transform
150  if (transformType_ != UNSPECIFIED)
151  {
152  // If the transforms are of different types then they are incompatible
153  if (transformType_ != t.transformType_)
154  {
155  return false;
156  }
157 
158  // Length-tolerance
159  const scalar lengthTolerance = lengthScale*matchTolerance;
160 
161  // If the transforms are both rotational then the axes must be in the
162  // same direction, the centre points must lie on the same line, and the
163  // angles (if available) must be opposing.
164  if (transformType_ == ROTATIONAL)
165  {
166  const scalar dot = rotationAxis_ & t.rotationAxis_;
167 
168  if (mag(dot) < 1 - matchTolerance)
169  {
170  return false;
171  }
172 
173  if
174  (
175  (rotationAxis_ & (rotationCentre_ - t.rotationCentre_))
176  > lengthTolerance
177  )
178  {
179  return false;
180  }
181 
182  if (transformComplete_ && t.transformComplete_)
183  {
184  if
185  (
186  mag(rotationAngle_ - sign(dot)*t.rotationAngle_)
187  > matchTolerance
188  )
189  {
190  return false;
191  }
192  }
193  }
194 
195  // If the transforms are both translational then the separations must
196  // be opposing
197  if (transformType_ == TRANSLATIONAL)
198  {
199  if (transformComplete_ && t.transformComplete_)
200  {
201  if (mag(separation_ - t.separation_) > lengthTolerance)
202  {
203  return false;
204  }
205  }
206  }
207  }
208 
209  // If the supplied transform is more complete then overwrite this with it
210  if (!transformComplete_ && t.transformComplete_)
211  {
212  *this = t;
213  }
214 
215  return true;
216 }
217 
218 
219 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
220 
222 :
223  cyclicTransform(true)
224 {}
225 
226 
228 (
229  const bool defaultIsNone
230 )
231 :
232  transformType_(defaultIsNone ? NONE : UNSPECIFIED),
233  rotationAxis_(vector::uniform(NaN)),
234  rotationCentre_(vector::uniform(NaN)),
235  rotationAngle_(NaN),
236  separation_(vector::uniform(NaN)),
237  transformComplete_(transformType_ == NONE),
238  transform_()
239 {}
240 
241 
243 (
244  const dictionary& dict,
245  const bool defaultIsNone
246 )
247 :
248  transformType_
249  (
250  transformTypeNames
251  [
252  dict.lookupOrDefaultBackwardsCompatible<word>
253  (
254  {"transformType", "transform"},
255  transformTypeNames[defaultIsNone ? NONE : UNSPECIFIED]
256  )
257  ]
258  ),
259  rotationAxis_
260  (
261  transformType_ == ROTATIONAL
262  ? normalised(dict.lookup<vector>("rotationAxis"))
263  : vector::uniform(NaN)
264  ),
265  rotationCentre_
266  (
267  transformType_ == ROTATIONAL
268  ? dict.lookup<point>("rotationCentre")
269  : point::uniform(NaN)
270  ),
271  rotationAngle_
272  (
273  dict.lookupOrDefault<scalar>("rotationAngle", unitDegrees, NaN)
274  ),
275  separation_
276  (
277  transformType_ == TRANSLATIONAL
278  ? (
280  (
281  {"separation", "separationVector"}
282  )
283  )
284  : vector::uniform(NaN)
285  ),
286  transformComplete_
287  (
288  transformType_ == NONE
289  || (
290  transformType_ == ROTATIONAL
291  && dict.found("rotationAngle")
292  )
293  || (
294  transformType_ == TRANSLATIONAL
295  && (dict.found("separation") || dict.found("separationVector"))
296  )
297  ),
298  transform_()
299 {
300  if (transformComplete_)
301  {
302  update();
303  }
304 }
305 
306 
308 (
309  const word& name,
310  const vectorField& areas,
311  const cyclicTransform& transform,
312  const word& nbrName,
313  const cyclicTransform& nbrTransform,
314  const scalar matchTolerance,
315  const bool global
316 )
317 :
319 {
320  // Calculate the total (vector) areas for the supplied patch data
321  const vector area = sum(areas, global);
322 
323  // Calculate patch length scales
324  const scalar lengthScale = sqrt(mag(area));
325 
326  // Copy as much data from the neighbour as possible
327  if (!transformComplete_ && nbrTransform.transformType_ != UNSPECIFIED)
328  {
329  if (!set(inv(nbrTransform), lengthScale, matchTolerance))
330  {
332  << "Patch " << name
333  << " and neighbour patch " << nbrName
334  << " have incompatible transforms:" << nl << nl << incrIndent;
335 
337  << indent << name << nl << indent << token::BEGIN_BLOCK << nl
338  << incrIndent;
339 
341 
343  << decrIndent << indent << token::END_BLOCK << nl << nl;
344 
346  << indent << nbrName << nl << indent << token::BEGIN_BLOCK << nl
347  << incrIndent;
348 
349  nbrTransform.write(FatalError);
350 
352  << decrIndent << indent << token::END_BLOCK << nl << nl;
353 
355  << decrIndent << exit(FatalError);
356  }
357  }
358 }
359 
360 
362 (
363  const word& name,
364  const pointField& ctrs,
365  const vectorField& areas,
366  const cyclicTransform& transform,
367  const word& nbrName,
368  const pointField& nbrCtrs,
369  const vectorField& nbrAreas,
370  const cyclicTransform& nbrTransform,
371  const scalar matchTolerance,
372  const bool global
373 )
374 :
376  (
377  name,
378  areas,
379  transform,
380  nbrName,
381  nbrTransform,
382  matchTolerance
383  )
384 {
385  // If there is no geometry from which to calculate the transform then
386  // nothing can be calculated
387  if (sum(areas.size(), global) == 0 || sum(nbrAreas.size(), global) == 0)
388  {
389  return;
390  }
391 
392  // Calculate the total (vector) areas for the supplied patch data
393  const vector area = sum(areas, global);
394  const vector nbrArea = sum(nbrAreas, global);
395 
396  // Calculate the centroids for the supplied patch data
397  const scalarField magAreas(mag(areas));
398  const scalarField magNbrAreas(mag(nbrAreas));
399  const scalar sumMagAreas = sum(magAreas, global);
400  const scalar sumMagNbrAreas = sum(magNbrAreas, global);
401  const point ctr = sum(ctrs*magAreas, global)/sumMagAreas;
402  const point nbrCtr = sum(nbrCtrs*magNbrAreas, global)/sumMagNbrAreas;
403 
404  // Calculate patch length scales
405  const scalar lengthScale = sqrt(sumMagAreas);
406 
407  // Calculate the transformation from the patch geometry
408  if (!transformComplete_)
409  {
410  // Store the old transformation type
411  const transformTypes oldTransformType = transformType_;
412 
413  // Calculate the average patch normals
414  const vector normal = normalised(area);
415  const vector negNbrNormal = - normalised(nbrArea);
416 
417  // Calculate the angle and distance separations
418  const scalar dot = normal & negNbrNormal;
419  const vector delta = ctr - nbrCtr;
420 
421  // Determine the type of transformation if it has not been specified
422  if (transformType_ == UNSPECIFIED)
423  {
424  transformType_ =
425  dot < 1 - rootSmall
426  ? ROTATIONAL
427  : mag(delta) > lengthScale*rootSmall
428  ? TRANSLATIONAL
429  : NONE;
430  }
431 
432  // If the transformation is known to be rotational, then we need to
433  // calculate the angle. If the transformation was previously
434  // unspecified then we also need to calculate the axis and the centre
435  // of rotation.
436  if (transformType_ == ROTATIONAL)
437  {
438  // Calculate the axis, if necessary
439  if (transformType_ != oldTransformType)
440  {
441  const vector midNormal = normalised(normal + negNbrNormal);
442  const vector axis =
443  (ctr - nbrCtr)
444  ^ (
445  normal*(negNbrNormal & midNormal)
446  - negNbrNormal*(normal & midNormal)
447  );
448  const vector axis180 =
449  (ctr - nbrCtr)
450  ^ (normal - negNbrNormal);
451 
452  rotationAxis_ =
453  normalised
454  (
455  mag(axis) > lengthScale*rootSmall
456  ? axis
457  : axis180
458  );
459  }
460 
461  const tensor PerpA = tensor::I - sqr(rotationAxis_);
462  const vector normalPerpA = normalised(PerpA & normal);
463  const vector negNbrNormalPerpA = normalised(PerpA & negNbrNormal);
464  const scalar theta =
465  acos(min(max(normalPerpA & negNbrNormalPerpA, -1), 1));
466  rotationAngle_ =
467  - sign((normalPerpA ^ negNbrNormalPerpA) & rotationAxis_)*theta;
468 
469  // Calculate the angle
470  // Calculate the centre of rotation, if necessary
471  if (transformType_ != oldTransformType)
472  {
473  const tensor R = quaternion(rotationAxis_, theta).R();
474  tensor A = tensor::I - R;
475  vector b = ctr - (R & nbrCtr);
476  const label i = findMax(cmptMag(rotationAxis_));
477  forAll(b, j)
478  {
479  A(i, j) = j == i;
480  }
481  b[i] = 0;
482  rotationCentre_ = inv(A) & b;
483  }
484  }
485 
486  // If the transformation is known to be translational then we just need
487  // to set the separation.
488  if (transformType_ == TRANSLATIONAL)
489  {
490  separation_ = delta;
491  }
492 
493  // Update the transform object
494  transformComplete_ = true;
495  update();
496 
497  // Print results of calculation
498  if (debug)
499  {
500  Info<< "Transformation calculated between patches " << name
501  << " and " << nbrName << ":" << nl << token::BEGIN_BLOCK << nl
502  << incrIndent;
503 
505 
507  }
508  }
509 
510  // Check the transformation is correct to within the matching tolerance
511  const point nbrCtrT =
512  transform_.transformPosition(nbrCtr);
513 
514  const scalar ctrNbrCtrTDistance = mag(ctr - nbrCtrT);
515 
516  if (ctrNbrCtrTDistance > lengthScale*matchTolerance)
517  {
518  OStringStream str;
519  str << "Patches " << name << " and " << nbrName << " are potentially "
520  << "not geometrically similar enough to be coupled." << nl << nl
521  << "The distance between the transformed centres of these patches "
522  << "is " << ctrNbrCtrTDistance << ", which is greater than the "
523  << "patch length scale (" << lengthScale << ") multiplied by the "
524  << "match tolerance (" << matchTolerance << ")." << nl << nl
525  << "Check that the patches are geometrically similar and that any "
526  << "transformations defined between them are correct." << nl << nl
527  << "If the patches and their transformations are defined correctly "
528  << "but small irregularities in the mesh mean this geometric test "
529  << "is failing, then it might be appropriate to relax the failure "
530  << "criteria by increasing the \"matchTolerance\" setting for "
531  << "these patches in the \"polyMesh/boundary\" file.";
533  << nl << stringOps::breakIntoIndentedLines(str.str()).c_str()
534  << exit(FatalError);
535  }
536 }
537 
538 
539 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
540 
542 {}
543 
544 
545 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
546 
548 {
549  const label oldPrecision = os.precision();
550 
551  os.precision(16);
552 
553  if (transformType_ != UNSPECIFIED)
554  {
555  writeEntry(os, "transformType", transformTypeNames[transformType_]);
556  }
557 
558  if (transformType_ == ROTATIONAL)
559  {
560  writeEntry(os, "rotationAxis", rotationAxis_);
561  writeEntry(os, "rotationCentre", rotationCentre_);
562 
563  if (transformComplete_)
564  {
565  writeEntry(os, "rotationAngle", unitDegrees, rotationAngle_);
566  }
567  }
568 
569  if (transformType_ == TRANSLATIONAL)
570  {
571  if (transformComplete_)
572  {
573  writeEntry(os, "separation", separation_);
574  }
575  }
576 
577  os.precision(oldPrecision);
578 }
579 
580 
581 // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
582 
583 Foam::cyclicTransform Foam::operator&
584 (
585  const transformer& t,
586  const cyclicTransform& c0
587 )
588 {
589  cyclicTransform c1(c0);
590 
591  if (c1.transformType_ == cyclicTransform::ROTATIONAL)
592  {
593  c1.rotationAxis_ = normalised(t.transform(c1.rotationAxis_));
594  c1.rotationCentre_ = t.transformPosition(c1.rotationCentre_);
595  }
596 
597  if (c1.transformType_ == cyclicTransform::TRANSLATIONAL)
598  {
599  if (c1.transformComplete_)
600  {
601  c1.separation_ = t.transform(c1.separation_);
602  }
603  }
604 
605  if (c1.transformComplete_)
606  {
607  c1.update();
608  }
609 
610  return c1;
611 }
612 
613 
615 {
616  cyclicTransform c1(c0);
617 
618  if (c1.transformType_ == cyclicTransform::ROTATIONAL)
619  {
620  if (c1.transformComplete_)
621  {
622  c1.rotationAngle_ = - c1.rotationAngle_;
623  }
624  }
625 
626  if (c1.transformType_ == cyclicTransform::TRANSLATIONAL)
627  {
628  if (c1.transformComplete_)
629  {
630  c1.separation_ = - c1.separation_;
631  }
632  }
633 
634  if (c1.transformComplete_)
635  {
636  c1.update();
637  }
638 
639  return c1;
640 }
641 
642 
643 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
Istream and Ostream manipulators taking arguments.
scalar delta
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Pre-declare SubField and related Field type.
Definition: Field.H:83
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:54
Output to memory buffer stream.
Definition: OStringStream.H:52
string str() const
Return the string.
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:57
virtual int precision() const =0
Get precision of output field.
static const Tensor I
Definition: Tensor.H:83
static Form uniform(const Cmpt &s)
Return a VectorSpace with all elements = s.
Definition: VectorSpaceI.H:168
Cyclic plane transformation.
static const NamedEnum< transformTypes, 4 > transformTypeNames
friend cyclicTransform inv(const cyclicTransform &c)
void write(Ostream &os) const
Write the data to a dictionary.
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T, if not found return the given default.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:710
ITstream & lookupBackwardsCompatible(const wordList &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream, trying a list of keywords.
Definition: dictionary.C:721
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:509
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:61
tensor R() const
The rotation tensor corresponding the quaternion.
Definition: quaternionI.H:323
A class for managing temporary objects.
Definition: tmp.H:55
@ BEGIN_BLOCK
Definition: token.H:113
@ END_BLOCK
Definition: token.H:114
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space.
Definition: transformer.H:84
static transformer rotation(const tensor &T)
Construct a pure rotation transformer.
Definition: transformerI.H:43
static transformer translation(const vector &t)
Construct a pure translation transformer.
Definition: transformerI.H:31
vector transformPosition(const vector &v) const
Transform the given position.
Definition: transformerI.H:153
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
volScalarField & b
Definition: createFields.H:25
const dimensionedScalar c1
First radiation constant: default SI units: [W/m^2].
string breakIntoIndentedLines(const string &str, const string::size_type nLength=80, const string::size_type nIndent=0)
Break a string up into indented lines.
Definition: stringOps.C:963
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Type gSum(const FieldField< Field, Type > &f)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:241
dimensionedScalar sign(const dimensionedScalar &ds)
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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
word name(const bool)
Return a word representation of a bool.
Definition: boolIO.C:39
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
messageStream Info
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:234
void dot(FieldField< Field1, typename innerProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void writeEntry(Ostream &os, const HashTable< T, Key, Hash > &ht)
Definition: HashTableIO.C:96
dimensionedScalar sqrt(const dimensionedScalar &ds)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:504
dimensionSet normalised(const dimensionSet &)
Definition: dimensionSet.C:510
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:296
error FatalError
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:227
static const char nl
Definition: Ostream.H:266
const unitConversion unitDegrees
dimensionedScalar acos(const dimensionedScalar &ds)
dictionary dict