coupledPolyPatch.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-2018 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 "coupledPolyPatch.H"
27 #include "ListOps.H"
28 #include "transform.H"
29 #include "OFstream.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(coupledPolyPatch, 0);
36 
37  const scalar coupledPolyPatch::defaultMatchTol_ = 1e-4;
38 
39  template<>
41  {
42  "unknown",
43  "rotational",
44  "translational",
45  "coincidentFullMatch",
46  "noOrdering"
47  };
48 
51 }
52 
53 
54 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
55 
57 {
58  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
59 }
60 
61 
63 (
64  Ostream& os,
65  const pointField& points,
66  const labelList& pointLabels
67 )
68 {
69  forAll(pointLabels, i)
70  {
71  writeOBJ(os, points[pointLabels[i]]);
72  }
73 }
74 
75 
77 (
78  Ostream& os,
79  const point& p0,
80  const point& p1,
81  label& vertI
82 )
83 {
84  writeOBJ(os, p0);
85  vertI++;
86 
87  writeOBJ(os, p1);
88  vertI++;
89 
90  os << "l " << vertI-1 << ' ' << vertI << nl;
91 }
92 
93 
95 (
96  const fileName& fName,
97  const UList<face>& faces,
98  const pointField& points
99 )
100 {
101  OFstream os(fName);
102 
103  Map<label> foamToObj(4*faces.size());
104 
105  label vertI = 0;
106 
107  forAll(faces, i)
108  {
109  const face& f = faces[i];
110 
111  forAll(f, fp)
112  {
113  if (foamToObj.insert(f[fp], vertI))
114  {
115  writeOBJ(os, points[f[fp]]);
116  vertI++;
117  }
118  }
119 
120  os << 'l';
121  forAll(f, fp)
122  {
123  os << ' ' << foamToObj[f[fp]]+1;
124  }
125  os << ' ' << foamToObj[f[0]]+1 << nl;
126  }
127 }
128 
129 
131 (
132  const UList<face>& faces,
133  const pointField& points,
134  const transformType transform
135 )
136 {
137  pointField anchors(faces.size());
138 
139  if (transform != COINCIDENTFULLMATCH)
140  {
141  // Return the first point
142  forAll(faces, facei)
143  {
144  anchors[facei] = points[faces[facei][0]];
145  }
146  }
147  else
148  {
149  // Make anchor point unique
150  forAll(faces, facei)
151  {
152  const face& f = faces[facei];
153 
154  bool unique = true;
155 
156  forAll(f, fp1)
157  {
158  const point& p1 = points[f[fp1]];
159 
160  unique = true;
161 
162  for (label fp2 = 0; fp2 < f.size(); ++fp2)
163  {
164  if (f[fp1] == f[fp2])
165  {
166  continue;
167  }
168 
169  const point& p2 = points[f[fp2]];
170 
171  // TODO: Change to a tolerance and possibly select closest
172  // point to the origin
173  if (p1 == p2)
174  {
175  unique = false;
176  break;
177  }
178  }
179 
180  if (unique)
181  {
182  anchors[facei] = p1;
183  break;
184  }
185  }
186 
187  if (!unique)
188  {
189  anchors[facei] = points[faces[facei][0]];
190  }
191  }
192  }
193 
194  return anchors;
195 }
196 
197 
199 (
200  const UList<face>& faces,
201  const pointField& points,
202  const pointField& faceCentres
203 )
204 {
205  // Calculate typical distance per face
206  scalarField tols(faces.size());
207 
208  forAll(faces, facei)
209  {
210  const point& cc = faceCentres[facei];
211 
212  const face& f = faces[facei];
213 
214  // 1. calculate a typical size of the face. Use maximum distance
215  // to face centre
216  scalar maxLenSqr = -great;
217  // 2. as measure of truncation error when comparing two coordinates
218  // use small * maximum component
219  scalar maxCmpt = -great;
220 
221  forAll(f, fp)
222  {
223  const point& pt = points[f[fp]];
224  maxLenSqr = max(maxLenSqr, magSqr(pt - cc));
225  maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt)));
226  }
227 
228  tols[facei] = max
229  (
230  small,
231  max(small*maxCmpt, Foam::sqrt(maxLenSqr))
232  );
233  }
234  return tols;
235 }
236 
237 
239 (
240  const pointField& points,
241  const face& f,
242  const point& anchor,
243  const scalar tol
244 )
245 {
246  label anchorFp = -1;
247  scalar minDistSqr = great;
248 
249  forAll(f, fp)
250  {
251  scalar distSqr = magSqr(anchor - points[f[fp]]);
252 
253  if (distSqr < minDistSqr)
254  {
255  minDistSqr = distSqr;
256  anchorFp = fp;
257  }
258  }
259 
260  if (anchorFp == -1 || Foam::sqrt(minDistSqr) > tol)
261  {
262  return -1;
263  }
264  else
265  {
266  // Check that anchor is unique.
267  forAll(f, fp)
268  {
269  scalar distSqr = magSqr(anchor - points[f[fp]]);
270 
271  if (distSqr == minDistSqr && fp != anchorFp)
272  {
274  << "Cannot determine unique anchor point on face "
276  << endl
277  << "Both at index " << anchorFp << " and " << fp
278  << " the vertices have the same distance "
279  << Foam::sqrt(minDistSqr)
280  << " to the anchor " << anchor
281  << ". Continuing but results might be wrong."
282  << nl << endl;
283  }
284  }
285 
286  // Positive rotation
287  return (f.size() - anchorFp) % f.size();
288  }
289 }
290 
291 
293 (
294  const vectorField& Cf,
295  const vectorField& Cr,
296  const vectorField& nf,
297  const vectorField& nr,
298  const scalarField& smallDist,
299  const scalar absTol,
300  const transformType transform
301 ) const
302 {
303  if (debug)
304  {
305  Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
306  << " transform:" << transformTypeNames[transform] << nl
307  << " (half)size:" << Cf.size() << nl
308  << " absTol:" << absTol << nl
309  << " smallDist min:" << min(smallDist) << nl
310  << " smallDist max:" << max(smallDist) << nl
311  << " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
312  }
313 
314  // Tolerance calculation.
315  // - normal calculation: assume absTol is the absolute error in a
316  // single normal/transformation calculation. Consists both of numerical
317  // precision (on the order of small and of writing precision
318  // (from e.g. decomposition)
319  // Then the overall error of summing the normals is sqrt(size())*absTol
320  // - separation calculation: pass in from the outside an allowable error.
321 
322  if (Cf.size() == 0)
323  {
324  // Dummy geometry. Assume non-separated, parallel.
325  separation_.setSize(0);
326  forwardT_.clear();
327  reverseT_.clear();
328  collocated_.setSize(0);
329  }
330  else
331  {
332  scalar error = absTol*Foam::sqrt(1.0*Cf.size());
333 
334  if (debug)
335  {
336  Pout<< " error:" << error << endl;
337  }
338 
339  if
340  (
341  transform == ROTATIONAL
342  || (
343  transform != TRANSLATIONAL
344  && transform != COINCIDENTFULLMATCH
345  && (sum(mag(nf & nr)) < Cf.size() - error)
346  )
347  )
348  {
349  // Type is rotation or unknown and normals not aligned
350 
351  // Assume per-face differing transformation, correct later
352 
353  separation_.setSize(0);
354 
355  forwardT_.setSize(Cf.size());
356  reverseT_.setSize(Cf.size());
357  collocated_.setSize(Cf.size());
358  collocated_ = false;
359 
360  forAll(forwardT_, facei)
361  {
362  forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]);
363  reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]);
364  }
365 
366  if (debug)
367  {
368  Pout<< " sum(mag(forwardT_ - forwardT_[0])):"
369  << sum(mag(forwardT_ - forwardT_[0]))
370  << endl;
371  }
372 
373  if (sum(mag(forwardT_ - forwardT_[0])) < error)
374  {
375  forwardT_.setSize(1);
376  reverseT_.setSize(1);
377  collocated_.setSize(1);
378 
379  if (debug)
380  {
381  Pout<< " difference in rotation less than"
382  << " local tolerance "
383  << error << ". Assuming uniform rotation." << endl;
384  }
385  }
386  }
387  else
388  {
389  // Translational or (unknown and normals aligned)
390 
391  forwardT_.setSize(0);
392  reverseT_.setSize(0);
393 
394  separation_ = Cr - Cf;
395 
396  collocated_.setSize(separation_.size());
397 
398  // Three situations:
399  // - separation is zero. No separation.
400  // - separation is same. Single separation vector.
401  // - separation differs per face. Separation vectorField.
402 
403  // Check for different separation per face
404  bool sameSeparation = true;
405  bool doneWarning = false;
406 
407  forAll(separation_, facei)
408  {
409  scalar smallSqr = sqr(smallDist[facei]);
410 
411  collocated_[facei] = (magSqr(separation_[facei]) < smallSqr);
412 
413  // Check if separation differing w.r.t. face 0.
414  if (magSqr(separation_[facei] - separation_[0]) > smallSqr)
415  {
416  sameSeparation = false;
417 
418  if (!doneWarning && debug)
419  {
420  doneWarning = true;
421 
422  Pout<< " separation " << separation_[facei]
423  << " at " << facei
424  << " differs from separation[0] " << separation_[0]
425  << " by more than local tolerance "
426  << smallDist[facei]
427  << ". Assuming non-uniform separation." << endl;
428  }
429  }
430  }
431 
432  if (sameSeparation)
433  {
434  // Check for zero separation (at 0 so everywhere)
435  if (collocated_[0])
436  {
437  if (debug)
438  {
439  Pout<< " separation " << mag(separation_[0])
440  << " less than local tolerance " << smallDist[0]
441  << ". Assuming zero separation." << endl;
442  }
443 
444  separation_.setSize(0);
445  collocated_ = boolList(1, true);
446  }
447  else
448  {
449  if (debug)
450  {
451  Pout<< " separation " << mag(separation_[0])
452  << " more than local tolerance " << smallDist[0]
453  << ". Assuming uniform separation." << endl;
454  }
455 
456  separation_.setSize(1);
457  collocated_ = boolList(1, false);
458  }
459  }
460  }
461  }
462 
463  if (debug)
464  {
465  Pout<< " separation_:" << separation_.size() << nl
466  << " forwardT size:" << forwardT_.size() << endl;
467  }
468 }
469 
470 
471 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
472 
474 (
475  const word& name,
476  const label size,
477  const label start,
478  const label index,
479  const polyBoundaryMesh& bm,
480  const word& patchType,
481  const transformType transform
482 )
483 :
484  polyPatch(name, size, start, index, bm, patchType),
485  matchTolerance_(defaultMatchTol_),
486  transform_(transform)
487 {}
488 
489 
491 (
492  const word& name,
493  const dictionary& dict,
494  const label index,
495  const polyBoundaryMesh& bm,
496  const word& patchType
497 )
498 :
499  polyPatch(name, dict, index, bm, patchType),
500  matchTolerance_(dict.lookupOrDefault("matchTolerance", defaultMatchTol_)),
501  transform_
502  (
503  dict.found("transform")
504  ? transformTypeNames.read(dict.lookup("transform"))
505  : UNKNOWN
506  )
507 {}
508 
509 
511 (
512  const coupledPolyPatch& pp,
513  const polyBoundaryMesh& bm
514 )
515 :
516  polyPatch(pp, bm),
517  matchTolerance_(pp.matchTolerance_),
518  transform_(pp.transform_)
519 {}
520 
521 
523 (
524  const coupledPolyPatch& pp,
525  const polyBoundaryMesh& bm,
526  const label index,
527  const label newSize,
528  const label newStart
529 )
530 :
531  polyPatch(pp, bm, index, newSize, newStart),
532  matchTolerance_(pp.matchTolerance_),
533  transform_(pp.transform_)
534 {}
535 
536 
538 (
539  const coupledPolyPatch& pp,
540  const polyBoundaryMesh& bm,
541  const label index,
542  const labelUList& mapAddressing,
543  const label newStart
544 )
545 :
546  polyPatch(pp, bm, index, mapAddressing, newStart),
547  matchTolerance_(pp.matchTolerance_),
548  transform_(pp.transform_)
549 {}
550 
551 
552 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
553 
555 {}
556 
557 
558 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
559 
561 {
562  polyPatch::write(os);
563  // if (matchTolerance_ != defaultMatchTol_)
564  {
565  os.writeKeyword("matchTolerance") << matchTolerance_
566  << token::END_STATEMENT << nl;
567  os.writeKeyword("transform") << transformTypeNames[transform_]
568  << token::END_STATEMENT << nl;
569  }
570 }
571 
572 
573 // ************************************************************************* //
void cmptMax(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:431
#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
A class for handling file names.
Definition: fileName.H:69
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
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.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Output to file stream.
Definition: OFstream.H:82
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
const Cmpt & z() const
Definition: VectorI.H:87
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:47
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
const Cmpt & y() const
Definition: VectorI.H:81
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
Definition: polyPatch.C:357
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
Class to handle errors and exceptions in a simple, consistent stream-based manner.
Definition: error.H:66
static pointField getAnchorPoints(const UList< face > &, const pointField &, const transformType)
Get a unique anchor point for all faces.
static const NamedEnum< transformType, 5 > transformTypeNames
const pointField & points
A class for handling words, derived from string.
Definition: word.H:59
3D tensor transformation operations.
static void writeOBJ(Ostream &os, const point &pt)
Write point in OBJ format.
const Cmpt & x() const
Definition: VectorI.H:75
Foam::polyBoundaryMesh.
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const char nl
Definition: Ostream.H:265
defineTypeNameAndDebug(combustionModel, 0)
virtual void write(Ostream &) const
Write the polyPatch data as a dictionary.
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:54
labelList f(nPoints)
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.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
static scalarField calcFaceTol(const UList< face > &faces, const pointField &points, const pointField &faceCentres)
Calculate typical tolerance per face. Is currently max distance.
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.
void setSize(const label)
Reset size of List.
Definition: List.C:281
virtual ~coupledPolyPatch()
Destructor.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
Namespace for OpenFOAM.
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:576