coupledPolyPatch.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2016 OpenFOAM Foundation
9  Copyright (C) 2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "coupledPolyPatch.H"
30 #include "ListOps.H"
31 #include "transform.H"
32 #include "OFstream.H"
33 
34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
35 
36 namespace Foam
37 {
38  defineTypeNameAndDebug(coupledPolyPatch, 0);
39 
40  const scalar coupledPolyPatch::defaultMatchTol_ = 1e-4;
41 }
42 
43 
44 const Foam::Enum
45 <
47 >
49 ({
50  { transformType::UNKNOWN, "unknown" },
51  { transformType::ROTATIONAL, "rotational" },
52  { transformType::TRANSLATIONAL, "translational" },
53  { transformType::COINCIDENTFULLMATCH, "coincidentFullMatch" },
54  { transformType::NOORDERING, "noOrdering" },
55 });
56 
57 
58 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
59 
61 {
62  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
63 }
64 
65 
67 (
68  Ostream& os,
69  const pointField& points,
70  const labelList& pointLabels
71 )
72 {
73  for (const label pointi : pointLabels)
74  {
75  writeOBJ(os, points[pointi]);
76  }
77 }
78 
79 
81 (
82  Ostream& os,
83  const point& p0,
84  const point& p1,
85  label& vertI
86 )
87 {
88  writeOBJ(os, p0);
89  vertI++;
90 
91  writeOBJ(os, p1);
92  vertI++;
93 
94  os << "l " << vertI-1 << ' ' << vertI << nl;
95 }
96 
97 
99 (
100  const fileName& fName,
101  const UList<face>& faces,
102  const pointField& points
103 )
104 {
105  OFstream os(fName);
106 
107  Map<label> foamToObj(4*faces.size());
108 
109  label vertI = 0;
110 
111  for (const face& f : faces)
112  {
113  forAll(f, fp)
114  {
115  if (foamToObj.insert(f[fp], vertI))
116  {
117  writeOBJ(os, points[f[fp]]);
118  vertI++;
119  }
120  }
121 
122  os << 'l';
123  forAll(f, fp)
124  {
125  os << ' ' << foamToObj[f[fp]]+1;
126  }
127  os << ' ' << foamToObj[f[0]]+1 << nl;
128  }
129 }
130 
131 
133 (
134  const UList<face>& faces,
135  const pointField& points,
137 )
138 {
139  pointField anchors(faces.size());
140 
141  if (transform != COINCIDENTFULLMATCH)
142  {
143  // Return the first point
144  forAll(faces, facei)
145  {
146  anchors[facei] = points[faces[facei][0]];
147  }
148  }
149  else
150  {
151  // Make anchor point unique
152  forAll(faces, facei)
153  {
154  const face& f = faces[facei];
155 
156  bool unique = true;
157 
158  forAll(f, fp1)
159  {
160  const point& p1 = points[f[fp1]];
161 
162  unique = true;
163 
164  for (label fp2 = 0; fp2 < f.size(); ++fp2)
165  {
166  if (f[fp1] == f[fp2])
167  {
168  continue;
169  }
170 
171  const point& p2 = points[f[fp2]];
172 
173  // TODO: Change to a tolerance and possibly select closest
174  // point to the origin
175  if (p1 == p2)
176  {
177  unique = false;
178  break;
179  }
180  }
181 
182  if (unique)
183  {
184  anchors[facei] = p1;
185  break;
186  }
187  }
188 
189  if (!unique)
190  {
191  anchors[facei] = points[faces[facei][0]];
192  }
193  }
194  }
195 
196  return anchors;
197 }
198 
199 
201 (
202  const UList<face>& faces,
203  const pointField& points,
204  const pointField& faceCentres
205 )
206 {
207  // Calculate typical distance per face
208  scalarField tols(faces.size());
209 
210  forAll(faces, facei)
211  {
212  const point& cc = faceCentres[facei];
213 
214  const face& f = faces[facei];
215 
216  // 1. calculate a typical size of the face. Use maximum distance
217  // to face centre
218  scalar maxLenSqr = -GREAT;
219  // 2. as measure of truncation error when comparing two coordinates
220  // use SMALL * maximum component
221  scalar maxCmpt = -GREAT;
222 
223  forAll(f, fp)
224  {
225  const point& pt = points[f[fp]];
226  maxLenSqr = max(maxLenSqr, magSqr(pt - cc));
227  maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt)));
228  }
229 
230  tols[facei] = max
231  (
232  SMALL,
233  max(SMALL*maxCmpt, Foam::sqrt(maxLenSqr))
234  );
235  }
236  return tols;
237 }
238 
239 
241 (
242  const pointField& points,
243  const face& f,
244  const point& anchor,
245  const scalar tol
246 )
247 {
248  label anchorFp = -1;
249  scalar minDistSqr = GREAT;
250 
251  forAll(f, fp)
252  {
253  scalar distSqr = magSqr(anchor - points[f[fp]]);
254 
255  if (distSqr < minDistSqr)
256  {
257  minDistSqr = distSqr;
258  anchorFp = fp;
259  }
260  }
261 
262  if (anchorFp == -1 || Foam::sqrt(minDistSqr) > tol)
263  {
264  return -1;
265  }
266  else
267  {
268  // Check that anchor is unique.
269  forAll(f, fp)
270  {
271  scalar distSqr = magSqr(anchor - points[f[fp]]);
272 
273  if (distSqr == minDistSqr && fp != anchorFp)
274  {
276  << "Cannot determine unique anchor point on face "
278  << endl
279  << "Both at index " << anchorFp << " and " << fp
280  << " the vertices have the same distance "
281  << Foam::sqrt(minDistSqr)
282  << " to the anchor " << anchor
283  << ". Continuing but results might be wrong."
284  << nl << endl;
285  }
286  }
287 
288  // Positive rotation
289  return (f.size() - anchorFp) % f.size();
290  }
291 }
292 
293 
295 (
296  const vectorField& Cf,
297  const vectorField& Cr,
298  const vectorField& nf,
299  const vectorField& nr,
300  const scalarField& smallDist,
301  const scalar absTol,
303 ) const
304 {
305  if (debug)
306  {
307  Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
308  << " transform:" << transformTypeNames[transform] << nl
309  << " (half)size:" << Cf.size() << nl
310  << " absTol:" << absTol << nl
311  << " smallDist min:" << min(smallDist) << nl
312  << " smallDist max:" << max(smallDist) << nl
313  << " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
314  }
315 
316  // Tolerance calculation.
317  // - normal calculation: assume absTol is the absolute error in a
318  // single normal/transformation calculation. Consists both of numerical
319  // precision (on the order of SMALL and of writing precision
320  // (from e.g. decomposition)
321  // Then the overall error of summing the normals is sqrt(size())*absTol
322  // - separation calculation: pass in from the outside an allowable error.
323 
324  if (Cf.size() == 0)
325  {
326  // Dummy geometry. Assume non-separated, parallel.
327  separation_.setSize(0);
328  forwardT_.clear();
329  reverseT_.clear();
330  collocated_.setSize(0);
331  }
332  else
333  {
334  scalar error = absTol*Foam::sqrt(1.0*Cf.size());
335 
336  if (debug)
337  {
338  Pout<< " error:" << error << endl;
339  }
340 
341  if
342  (
343  transform == ROTATIONAL
344  || (
345  transform != TRANSLATIONAL
346  && transform != COINCIDENTFULLMATCH
347  && (sum(mag(nf & nr)) < Cf.size() - error)
348  )
349  )
350  {
351  // Type is rotation or unknown and normals not aligned
352 
353  // Assume per-face differing transformation, correct later
354 
355  separation_.setSize(0);
356 
357  forwardT_.setSize(Cf.size());
358  reverseT_.setSize(Cf.size());
359  collocated_.setSize(Cf.size());
360  collocated_ = false;
361 
362  forAll(forwardT_, facei)
363  {
364  forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]);
365  reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]);
366  }
367 
368  if (debug)
369  {
370  Pout<< " sum(mag(forwardT_ - forwardT_[0])):"
371  << sum(mag(forwardT_ - forwardT_[0]))
372  << endl;
373  }
374 
375  if (sum(mag(forwardT_ - forwardT_[0])) < error)
376  {
377  forwardT_.setSize(1);
378  reverseT_.setSize(1);
379  collocated_.setSize(1);
380 
381  if (debug)
382  {
383  Pout<< " difference in rotation less than"
384  << " local tolerance "
385  << error << ". Assuming uniform rotation." << endl;
386  }
387  }
388  }
389  else
390  {
391  // Translational or (unknown and normals aligned)
392 
393  forwardT_.setSize(0);
394  reverseT_.setSize(0);
395 
396  separation_ = Cr - Cf;
397 
398  collocated_.setSize(separation_.size());
399 
400  // Three situations:
401  // - separation is zero. No separation.
402  // - separation is same. Single separation vector.
403  // - separation differs per face. Separation vectorField.
404 
405  // Check for different separation per face
406  bool sameSeparation = true;
407  bool doneWarning = false;
408 
409  forAll(separation_, facei)
410  {
411  scalar smallSqr = sqr(smallDist[facei]);
412 
413  collocated_[facei] = (magSqr(separation_[facei]) < smallSqr);
414 
415  // Check if separation differing w.r.t. face 0.
416  if (magSqr(separation_[facei] - separation_[0]) > smallSqr)
417  {
418  sameSeparation = false;
419 
420  if (!doneWarning && debug)
421  {
422  doneWarning = true;
423 
424  Pout<< " separation " << separation_[facei]
425  << " at " << facei
426  << " differs from separation[0] " << separation_[0]
427  << " by more than local tolerance "
428  << smallDist[facei]
429  << ". Assuming non-uniform separation." << endl;
430  }
431  }
432  }
433 
434  if (sameSeparation)
435  {
436  // Check for zero separation (at 0 so everywhere)
437  if (collocated_[0])
438  {
439  if (debug)
440  {
441  Pout<< " separation " << mag(separation_[0])
442  << " less than local tolerance " << smallDist[0]
443  << ". Assuming zero separation." << endl;
444  }
445 
446  separation_.setSize(0);
447  collocated_ = boolList(1, true);
448  }
449  else
450  {
451  if (debug)
452  {
453  Pout<< " separation " << mag(separation_[0])
454  << " more than local tolerance " << smallDist[0]
455  << ". Assuming uniform separation." << endl;
456  }
457 
458  separation_.setSize(1);
459  collocated_ = boolList(1, false);
460  }
461  }
462  }
463  }
464 
465  if (debug)
466  {
467  Pout<< " separation_:" << separation_.size() << nl
468  << " forwardT size:" << forwardT_.size() << endl;
469  }
470 }
471 
472 
473 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
474 
476 (
477  const word& name,
478  const label size,
479  const label start,
480  const label index,
481  const polyBoundaryMesh& bm,
482  const word& patchType,
484 )
485 :
486  polyPatch(name, size, start, index, bm, patchType),
487  matchTolerance_(defaultMatchTol_),
488  transform_(transform)
489 {}
490 
491 
493 (
494  const word& name,
495  const dictionary& dict,
496  const label index,
497  const polyBoundaryMesh& bm,
498  const word& patchType
499 )
500 :
501  polyPatch(name, dict, index, bm, patchType),
502  matchTolerance_(dict.getOrDefault("matchTolerance", defaultMatchTol_)),
503  transform_
504  (
505  transformTypeNames.getOrDefault
506  (
507  "transform",
508  dict,
509  transformType::UNKNOWN
510  )
511  )
512 {}
513 
514 
516 (
517  const coupledPolyPatch& pp,
518  const polyBoundaryMesh& bm
519 )
520 :
521  polyPatch(pp, bm),
522  matchTolerance_(pp.matchTolerance_),
523  transform_(pp.transform_)
524 {}
525 
526 
528 (
529  const coupledPolyPatch& pp,
530  const labelList& faceCells
531 )
532 :
533  polyPatch(pp, faceCells),
534  matchTolerance_(pp.matchTolerance_),
535  transform_(pp.transform_)
536 {}
537 
538 
540 (
541  const coupledPolyPatch& pp,
542  const polyBoundaryMesh& bm,
543  const label index,
544  const label newSize,
545  const label newStart
546 )
547 :
548  polyPatch(pp, bm, index, newSize, newStart),
549  matchTolerance_(pp.matchTolerance_),
550  transform_(pp.transform_)
551 {}
552 
553 
555 (
556  const coupledPolyPatch& pp,
557  const polyBoundaryMesh& bm,
558  const label index,
559  const labelUList& mapAddressing,
560  const label newStart
561 )
562 :
563  polyPatch(pp, bm, index, mapAddressing, newStart),
564  matchTolerance_(pp.matchTolerance_),
565  transform_(pp.transform_)
566 {}
567 
568 
569 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
570 
572 {}
573 
574 
575 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
576 
578 {
580  //if (matchTolerance_ != defaultMatchTol_)
581  {
582  os.writeEntry("matchTolerance", matchTolerance_);
583  os.writeEntry("transform", transformTypeNames[transform_]);
584  }
585 }
586 
587 
588 // ************************************************************************* //
Foam::coupledPolyPatch::getAnchorPoints
static pointField getAnchorPoints(const UList< face > &, const pointField &, const transformType)
Get a unique anchor point for all faces.
Definition: coupledPolyPatch.C:133
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
Foam::Vector::x
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
Foam::Enum
Enum is a wrapper around a list of names/values that represent particular enumeration (or int) values...
Definition: IOstreamOption.H:57
Foam::word
A class for handling words, derived from Foam::string.
Definition: word.H:65
Foam::fileName
A class for handling file names.
Definition: fileName.H:73
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO.
Definition: polyBoundaryMesh.H:63
Foam::meshTools::writeOBJ
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Foam::coupledPolyPatch::transformTypeNames
static const Enum< transformType > transformTypeNames
Definition: coupledPolyPatch.H:70
Foam::Map< label >
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::coupledPolyPatch
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
Definition: coupledPolyPatch.H:54
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:521
Foam::Vector::z
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
OFstream.H
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::coupledPolyPatch::write
virtual void write(Ostream &os) const
Write the polyPatch data as a dictionary.
Definition: coupledPolyPatch.C:577
Foam::cmptMag
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:400
Foam::cmptMax
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
Definition: FieldFieldFunctions.C:253
coupledPolyPatch.H
Foam::Field< vector >
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:68
Foam::coupledPolyPatch::getRotation
static label getRotation(const pointField &points, const face &f, const point &anchor, const scalar tol)
Get the number of vertices face f needs to be rotated such that.
Definition: coupledPolyPatch.C:241
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::polyPatch::write
virtual void write(Ostream &os) const
Write the polyPatch data as a dictionary.
Definition: polyPatch.C:413
os
OBJstream os(runTime.globalPath()/outputName)
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::coupledPolyPatch::calcFaceTol
static scalarField calcFaceTol(const UList< face > &faces, const pointField &points, const pointField &faceCentres)
Calculate typical tolerance per face. Is currently max distance.
Definition: coupledPolyPatch.C:201
Foam::OFstream
Output to file stream, using an OSstream.
Definition: OFstream.H:53
Foam::Vector::y
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::nl
constexpr char nl
Definition: Ostream.H:404
f
labelList f(nPoints)
Foam::Vector< scalar >
Foam::coupledPolyPatch::~coupledPolyPatch
virtual ~coupledPolyPatch()
Destructor.
Definition: coupledPolyPatch.C:571
Foam::List< label >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::UList< face >
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::sum
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:327
Foam::coupledPolyPatch::calcTransformTensors
void calcTransformTensors(const vectorField &Cf, const vectorField &Cr, const vectorField &nf, const vectorField &nr, const scalarField &smallDist, const scalar absTol, const transformType=UNKNOWN) const
Calculate the transformation tensors.
Definition: coupledPolyPatch.C:295
Foam::Ostream::writeEntry
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
Definition: Ostream.H:236
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::UIndirectList
A List with indirect addressing.
Definition: faMatrix.H:60
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:72
ListOps.H
Various functions to operate on Lists.
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Foam::Ostream
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:56
p0
const volScalarField & p0
Definition: EEqn.H:36
transform.H
3D tensor transformation operations.
Foam::rotationTensor
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:51
Foam::coupledPolyPatch::writeOBJ
static void writeOBJ(Ostream &os, const point &pt)
Write point in OBJ format.
Definition: coupledPolyPatch.C:60
Foam::coupledPolyPatch::transformType
transformType
Definition: coupledPolyPatch.H:60
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::coupledPolyPatch::coupledPolyPatch
coupledPolyPatch(const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm, const word &patchType, const transformType transform)
Construct from components.
Definition: coupledPolyPatch.C:476
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328
Foam::error
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:73
Foam::faceCells
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:56
pointLabels
labelList pointLabels(nPoints, -1)