FreeFOAM The Cross-Platform CFD Toolkit
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 | Copyright (C) 1991-2011 OpenCFD Ltd.
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 <OpenFOAM/ListOps.H>
28 #include <OpenFOAM/transform.H>
29 #include <OpenFOAM/OFstream.H>
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
34 
35 template<>
37 {
38  "unknown",
39  "rotational",
40  "translational"
41 };
42 
45 
46 Foam::scalar Foam::coupledPolyPatch::matchTol = 1E-3;
47 
48 
49 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
50 
52 {
53  os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
54 }
55 
56 
58 (
59  Ostream& os,
60  const pointField& points,
61  const labelList& pointLabels
62 )
63 {
64  forAll(pointLabels, i)
65  {
66  writeOBJ(os, points[pointLabels[i]]);
67  }
68 }
69 
70 
72 (
73  Ostream& os,
74  const point& p0,
75  const point& p1,
76  label& vertI
77 )
78 {
79  writeOBJ(os, p0);
80  vertI++;
81 
82  writeOBJ(os, p1);
83  vertI++;
84 
85  os<< "l " << vertI-1 << ' ' << vertI << nl;
86 }
87 
88 
90 (
91  const fileName& fName,
92  const UList<face>& faces,
93  const pointField& points
94 )
95 {
96  OFstream os(fName);
97 
98  Map<label> foamToObj(4*faces.size());
99 
100  label vertI = 0;
101 
102  forAll(faces, i)
103  {
104  const face& f = faces[i];
105 
106  forAll(f, fp)
107  {
108  if (foamToObj.insert(f[fp], vertI))
109  {
110  writeOBJ(os, points[f[fp]]);
111  vertI++;
112  }
113  }
114 
115  os << 'l';
116  forAll(f, fp)
117  {
118  os << ' ' << foamToObj[f[fp]]+1;
119  }
120  os << ' ' << foamToObj[f[0]]+1 << nl;
121  }
122 }
123 
124 
125 Foam::pointField Foam::coupledPolyPatch::calcFaceCentres
126 (
127  const UList<face>& faces,
128  const pointField& points
129 )
130 {
131  pointField ctrs(faces.size());
132 
133  forAll(faces, faceI)
134  {
135  ctrs[faceI] = faces[faceI].centre(points);
136  }
137 
138  return ctrs;
139 }
140 
141 
143 (
144  const UList<face>& faces,
145  const pointField& points
146 )
147 {
148  pointField anchors(faces.size());
149 
150  forAll(faces, faceI)
151  {
152  anchors[faceI] = points[faces[faceI][0]];
153  }
154 
155  return anchors;
156 }
157 
158 
160 (
161  const labelList& oldToNew,
162  const label oldFaceI
163 ) const
164 {
165  label faceI = oldToNew[oldFaceI];
166 
167  return faceI >= start() && faceI < start()+size();
168 }
169 
170 
172 (
173  const labelList& patchStarts,
174  const label faceI
175 )
176 {
177  forAll(patchStarts, patchI)
178  {
179  if (patchStarts[patchI] <= faceI)
180  {
181  if (patchI == patchStarts.size()-1)
182  {
183  return patchI;
184  }
185  else if (patchStarts[patchI+1] > faceI)
186  {
187  return patchI;
188  }
189  }
190  }
191 
192  return -1;
193 }
194 
195 
197 (
198  const UList<face>& faces,
199  const pointField& points,
200  const pointField& faceCentres
201 )
202 {
203  // Calculate typical distance per face
204  scalarField tols(faces.size());
205 
206  forAll(faces, faceI)
207  {
208  const point& cc = faceCentres[faceI];
209 
210  const face& f = faces[faceI];
211 
212  // 1. calculate a typical size of the face. Use maximum distance
213  // to face centre
214  scalar maxLenSqr = -GREAT;
215  // 2. as measure of truncation error when comparing two coordinates
216  // use SMALL * maximum component
217  scalar maxCmpt = -GREAT;
218 
219  forAll(f, fp)
220  {
221  const point& pt = points[f[fp]];
222  maxLenSqr = max(maxLenSqr, magSqr(pt - cc));
223  maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt)));
224  }
225  tols[faceI] = max(SMALL*maxCmpt, matchTol*Foam::sqrt(maxLenSqr));
226  }
227  return tols;
228 }
229 
230 
232 (
233  const pointField& points,
234  const face& f,
235  const point& anchor,
236  const scalar tol
237 )
238 {
239  label anchorFp = -1;
240  scalar minDistSqr = GREAT;
241 
242  forAll(f, fp)
243  {
244  scalar distSqr = magSqr(anchor - points[f[fp]]);
245 
246  if (distSqr < minDistSqr)
247  {
248  minDistSqr = distSqr;
249  anchorFp = fp;
250  }
251  }
252 
253  if (anchorFp == -1 || mag(minDistSqr) > tol)
254  {
255  return -1;
256  }
257  else
258  {
259  // Positive rotation
260  return (f.size() - anchorFp) % f.size();
261  }
262 }
263 
264 
266 (
267  const vectorField& Cf,
268  const vectorField& Cr,
269  const vectorField& nf,
270  const vectorField& nr,
271  const scalarField& smallDist,
272  const scalar absTol,
274 ) const
275 {
276  if (debug)
277  {
278  Pout<< "coupledPolyPatch::calcTransformTensors : " << name() << endl
279  << " transform:" << transformTypeNames[transform] << nl
280  << " (half)size:" << Cf.size() << nl
281  << " absTol:" << absTol << nl
282  //<< " smallDist:" << smallDist << nl
283  << " sum(mag(nf & nr)):" << sum(mag(nf & nr)) << endl;
284  }
285 
286  // Tolerance calculation.
287  // - normal calculation: assume absTol is the absolute error in a
288  // single normal/transformation calculation. Consists both of numerical
289  // precision (on the order of SMALL and of writing precision
290  // (from e.g. decomposition)
291  // Then the overall error of summing the normals is sqrt(size())*absTol
292  // - separation calculation: pass in from the outside an allowable error.
293 
294  if (size() == 0)
295  {
296  // Dummy geometry.
297  separation_.setSize(0);
298  forwardT_ = I;
299  reverseT_ = I;
300  }
301  else
302  {
303  scalar error = absTol*Foam::sqrt(1.0*Cf.size());
304 
305  if (debug)
306  {
307  Pout<< " error:" << error << endl;
308  }
309 
310  if
311  (
312  transform == ROTATIONAL
313  || (
314  transform != TRANSLATIONAL
315  && (sum(mag(nf & nr)) < Cf.size()-error)
316  )
317  )
318  {
319  // Type is rotation or unknown and normals not aligned
320 
321  separation_.setSize(0);
322 
323  forwardT_.setSize(Cf.size());
324  reverseT_.setSize(Cf.size());
325 
326  forAll (forwardT_, facei)
327  {
328  forwardT_[facei] = rotationTensor(-nr[facei], nf[facei]);
329  reverseT_[facei] = rotationTensor(nf[facei], -nr[facei]);
330  }
331 
332  if (debug)
333  {
334  Pout<< " sum(mag(forwardT_ - forwardT_[0])):"
335  << sum(mag(forwardT_ - forwardT_[0]))
336  << endl;
337  }
338 
339  if (sum(mag(forwardT_ - forwardT_[0])) < error)
340  {
341  forwardT_.setSize(1);
342  reverseT_.setSize(1);
343 
344  if (debug)
345  {
346  Pout<< " difference in rotation less than"
347  << " local tolerance "
348  << error << ". Assuming uniform rotation." << endl;
349  }
350  }
351  }
352  else
353  {
354  // Translational or unknown and normals aligned.
355 
356  forwardT_.setSize(0);
357  reverseT_.setSize(0);
358 
359  separation_ = Cr - Cf;
360 
361  // Three situations:
362  // - separation is zero. No separation.
363  // - separation is same. Single separation vector.
364  // - separation differs per face. Separation vectorField.
365 
366  // Check for different separation per face
367  bool sameSeparation = true;
368 
369  forAll(separation_, facei)
370  {
371  scalar smallSqr = sqr(smallDist[facei]);
372 
373  if (magSqr(separation_[facei] - separation_[0]) > smallSqr)
374  {
375  if (debug)
376  {
377  Pout<< " separation " << separation_[facei]
378  << " at " << facei
379  << " differs from separation[0] " << separation_[0]
380  << " by more than local tolerance "
381  << smallDist[facei]
382  << ". Assuming non-uniform separation." << endl;
383  }
384  sameSeparation = false;
385  break;
386  }
387  }
388 
389  if (sameSeparation)
390  {
391  // Check for zero separation (at 0 so everywhere)
392  if (magSqr(separation_[0]) < sqr(smallDist[0]))
393  {
394  if (debug)
395  {
396  Pout<< " separation " << mag(separation_[0])
397  << " less than local tolerance " << smallDist[0]
398  << ". Assuming zero separation." << endl;
399  }
400 
401  separation_.setSize(0);
402  }
403  else
404  {
405  if (debug)
406  {
407  Pout<< " separation " << mag(separation_[0])
408  << " more than local tolerance " << smallDist[0]
409  << ". Assuming uniform separation." << endl;
410  }
411 
412  separation_.setSize(1);
413  }
414  }
415  }
416  }
417 
418  if (debug)
419  {
420  Pout<< " separation_:" << separation_.size() << nl
421  << " forwardT size:" << forwardT_.size() << endl;
422  }
423 }
424 
425 
426 // * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * //
427 
429 (
430  const word& name,
431  const label size,
432  const label start,
433  const label index,
434  const polyBoundaryMesh& bm
435 )
436 :
437  polyPatch(name, size, start, index, bm)
438 {}
439 
440 
442 (
443  const word& name,
444  const dictionary& dict,
445  const label index,
446  const polyBoundaryMesh& bm
447 )
448 :
449  polyPatch(name, dict, index, bm)
450 {}
451 
452 
454 (
455  const coupledPolyPatch& pp,
456  const polyBoundaryMesh& bm
457 )
458 :
459  polyPatch(pp, bm)
460 {}
461 
462 
464 (
465  const coupledPolyPatch& pp,
466  const polyBoundaryMesh& bm,
467  const label index,
468  const label newSize,
469  const label newStart
470 )
471 :
472  polyPatch(pp, bm, index, newSize, newStart)
473 {}
474 
475 
476 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
477 
479 {}
480 
481 
482 // ************************ vim: set sw=4 sts=4 et: ************************ //