FreeFOAM The Cross-Platform CFD Toolkit
localPointRegion.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-2010 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 "localPointRegion.H"
27 #include <OpenFOAM/syncTools.H>
28 #include <OpenFOAM/polyMesh.H>
29 #include <OpenFOAM/mapPolyMesh.H>
30 #include <OpenFOAM/globalIndex.H>
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 
38 defineTypeNameAndDebug(localPointRegion, 0);
39 
40 // Reduction class to get minimum value over face.
42 {
43 public:
44 
45  void operator()(face& x, const face& y) const
46  {
47  if (x.size())
48  {
49  label j = 0;
50  forAll(x, i)
51  {
52  x[i] = min(x[i], y[j]);
53 
54  j = y.rcIndex(j);
55  }
56  }
57  };
58 };
59 
60 
61 // Dummy transform for faces. Used in synchronisation
62 void transformList
63 (
64  const tensorField& rotTensor,
65  UList<face>& field
66 )
67 {};
68 
69 }
70 
71 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
72 
73 // Are two lists identical either in forward or in reverse order.
74 bool Foam::localPointRegion::isDuplicate
75 (
76  const face& f0,
77  const face& f1,
78  const bool forward
79 )
80 {
81  label fp1 = findIndex(f1, f0[0]);
82 
83  if (fp1 == -1)
84  {
85  return false;
86  }
87 
88  forAll(f0, fp0)
89  {
90  if (f0[fp0] != f1[fp1])
91  {
92  return false;
93  }
94 
95  if (forward)
96  {
97  fp1 = f1.fcIndex(fp1);
98  }
99  else
100  {
101  fp1 = f1.rcIndex(fp1);
102  }
103  }
104  return true;
105 }
106 
107 
108 // Count regions per point
109 void Foam::localPointRegion::countPointRegions
110 (
111  const polyMesh& mesh,
112  const boolList& candidatePoint,
113  const Map<label>& candidateFace,
114  faceList& minRegion
115 )
116 {
117  // Almost all will have only one so only
118  // populate Map if more than one.
119  labelList minPointRegion(mesh.nPoints(), -1);
120  // From global point to local (multi-region) point numbering
121  meshPointMap_.resize(candidateFace.size()/100);
122  // From local (multi-region) point to regions
123  DynamicList<labelList> pointRegions(meshPointMap_.size());
124 
125  // From faces with any duplicated point on it to local face
126  meshFaceMap_.resize(meshPointMap_.size());
127 
128  forAllConstIter(Map<label>, candidateFace, iter)
129  {
130  label faceI = iter.key();
131 
132  if (!mesh.isInternalFace(faceI))
133  {
134  const face& f = mesh.faces()[faceI];
135 
136  if (minRegion[faceI].empty())
137  {
138  FatalErrorIn("localPointRegion::countPointRegions(..)")
139  << "Face from candidateFace without minRegion set." << endl
140  << "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
141  << " verts:" << f << abort(FatalError);
142  }
143 
144  forAll(f, fp)
145  {
146  label pointI = f[fp];
147 
148  // Even points which were not candidates for splitting might
149  // be on multiple baffles that are being split so check.
150 
151  if (candidatePoint[pointI])
152  {
153  label region = minRegion[faceI][fp];
154 
155  if (minPointRegion[pointI] == -1)
156  {
157  minPointRegion[pointI] = region;
158  }
159  else if (minPointRegion[pointI] != region)
160  {
161  // Multiple regions for this point. Add.
162  Map<label>::iterator iter = meshPointMap_.find(pointI);
163  if (iter != meshPointMap_.end())
164  {
165  labelList& regions = pointRegions[iter()];
166  if (findIndex(regions, region) == -1)
167  {
168  label sz = regions.size();
169  regions.setSize(sz+1);
170  regions[sz] = region;
171  }
172  }
173  else
174  {
175  label localPointI = meshPointMap_.size();
176  meshPointMap_.insert(pointI, localPointI);
177  labelList regions(2);
178  regions[0] = minPointRegion[pointI];
179  regions[1] = region;
180  pointRegions.append(regions);
181  }
182 
183  label meshFaceMapI = meshFaceMap_.size();
184  meshFaceMap_.insert(faceI, meshFaceMapI);
185  }
186  }
187  }
188  }
189  }
190  minPointRegion.clear();
191 
192  // Add internal faces that use any duplicated point. Can only have one
193  // region!
194  forAllConstIter(Map<label>, candidateFace, iter)
195  {
196  label faceI = iter.key();
197 
198  if (mesh.isInternalFace(faceI))
199  {
200  const face& f = mesh.faces()[faceI];
201 
202  forAll(f, fp)
203  {
204  // Note: candidatePoint test not really necessary but
205  // speeds up rejection.
206  if (candidatePoint[f[fp]] && meshPointMap_.found(f[fp]))
207  {
208  label meshFaceMapI = meshFaceMap_.size();
209  meshFaceMap_.insert(faceI, meshFaceMapI);
210  }
211  }
212  }
213  }
214 
215 
216  // Transfer to member data
217  pointRegions.shrink();
218  pointRegions_.setSize(pointRegions.size());
219  forAll(pointRegions, i)
220  {
221  pointRegions_[i].transfer(pointRegions[i]);
222  }
223 
224  // Compact minRegion
225  faceRegions_.setSize(meshFaceMap_.size());
226  forAllConstIter(Map<label>, meshFaceMap_, iter)
227  {
228  faceRegions_[iter()].labelList::transfer(minRegion[iter.key()]);
229 
231  //{
232  // label faceI = iter.key();
233  // const face& f = mesh.faces()[faceI];
234  // Pout<< "Face:" << faceI << " fc:" << mesh.faceCentres()[faceI]
235  // << " verts:" << f << endl;
236  // forAll(f, fp)
237  // {
238  // Pout<< " " << f[fp] << " min:" << faceRegions_[iter()][fp]
239  // << endl;
240  // }
241  // Pout<< endl;
242  //}
243  }
244 
245  // Compact region numbering
246  // ? TBD.
247 }
248 
249 
250 void Foam::localPointRegion::calcPointRegions
251 (
252  const polyMesh& mesh,
253  boolList& candidatePoint
254 )
255 {
256  label nBnd = mesh.nFaces()-mesh.nInternalFaces();
257  const labelList& faceOwner = mesh.faceOwner();
258  const labelList& faceNeighbour = mesh.faceNeighbour();
259 
260 
262  (
263  mesh,
264  candidatePoint,
265  orEqOp<bool>(),
266  false, // nullValue
267  false // applySeparation
268  );
269 
270 
271  // Mark any face/boundaryFace/cell with a point on a candidate point.
272  // - candidateFace does not necessary have to be a baffle!
273  // - candidateFace is synchronised (since candidatePoint is)
274  Map<label> candidateFace(2*nBnd);
275  label candidateFaceI = 0;
276 
277  Map<label> candidateCell(nBnd);
278  label candidateCellI = 0;
279 
280  forAll(mesh.faces(), faceI)
281  {
282  const face& f = mesh.faces()[faceI];
283 
284  forAll(f, fp)
285  {
286  if (candidatePoint[f[fp]])
287  {
288  // Mark face
289  if (candidateFace.insert(faceI, candidateFaceI))
290  {
291  candidateFaceI++;
292  }
293 
294  // Mark cells
295  if (candidateCell.insert(faceOwner[faceI], candidateCellI))
296  {
297  candidateCellI++;
298  }
299 
300  if (mesh.isInternalFace(faceI))
301  {
302  label nei = faceNeighbour[faceI];
303  if (candidateCell.insert(nei, candidateCellI))
304  {
305  candidateCellI++;
306  }
307  }
308 
309  break;
310  }
311  }
312  }
313 
314 
315 
316  // Get global indices for cells
317  globalIndex globalCells(mesh.nCells());
318 
319 
320  // Determine for every candidate face per point the minimum region
321  // (global cell) it is connected to. (candidateFaces are the
322  // only ones using a
323  // candidate point so the only ones that can be affected)
324  faceList minRegion(mesh.nFaces());
325  forAllConstIter(Map<label>, candidateFace, iter)
326  {
327  label faceI = iter.key();
328  const face& f = mesh.faces()[faceI];
329 
330  if (mesh.isInternalFace(faceI))
331  {
332  label globOwn = globalCells.toGlobal(faceOwner[faceI]);
333  label globNei = globalCells.toGlobal(faceNeighbour[faceI]);
334  minRegion[faceI].setSize(f.size(), min(globOwn, globNei));
335  }
336  else
337  {
338  label globOwn = globalCells.toGlobal(faceOwner[faceI]);
339  minRegion[faceI].setSize(f.size(), globOwn);
340  }
341  }
342 
343  // Now minimize over all faces that are connected through internal
344  // faces or through cells. This loop iterates over the max number of
345  // cells connected to a point (=8 for hex mesh)
346 
347  while (true)
348  {
349  // Transport minimum from face across cell
350  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
351 
352  Map<label> minPointValue(100);
353  label nChanged = 0;
354  forAllConstIter(Map<label>, candidateCell, iter)
355  {
356  minPointValue.clear();
357 
358  label cellI = iter.key();
359  const cell& cFaces = mesh.cells()[cellI];
360 
361  // Determine minimum per point
362  forAll(cFaces, cFaceI)
363  {
364  label faceI = cFaces[cFaceI];
365 
366  if (minRegion[faceI].size())
367  {
368  const face& f = mesh.faces()[faceI];
369 
370  forAll(f, fp)
371  {
372  label pointI = f[fp];
373  Map<label>::iterator iter = minPointValue.find(pointI);
374 
375  if (iter == minPointValue.end())
376  {
377  minPointValue.insert(pointI, minRegion[faceI][fp]);
378  }
379  else
380  {
381  label currentMin = iter();
382  iter() = min(currentMin, minRegion[faceI][fp]);
383  }
384  }
385  }
386  }
387 
388  // Set face minimum from point minimum
389  forAll(cFaces, cFaceI)
390  {
391  label faceI = cFaces[cFaceI];
392 
393  if (minRegion[faceI].size())
394  {
395  const face& f = mesh.faces()[faceI];
396 
397  forAll(f, fp)
398  {
399  label minVal = minPointValue[f[fp]];
400 
401  if (minVal != minRegion[faceI][fp])
402  {
403  minRegion[faceI][fp] = minVal;
404  nChanged++;
405  }
406  }
407  }
408  }
409  }
410 
411  //Pout<< "nChanged:" << nChanged << endl;
412 
413  if (returnReduce(nChanged, sumOp<label>()) == 0)
414  {
415  break;
416  }
417 
418 
419  // Transport minimum across coupled faces
420  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
421 
423  (
424  mesh,
425  minRegion,
426  minEqOpFace(),
427  false // applySeparation
428  );
429  }
430 
431 
432  // Count regions per point
433  countPointRegions(mesh, candidatePoint, candidateFace, minRegion);
434  minRegion.clear();
435 
436 
438  //forAllConstIter(Map<label>, meshPointMap_, iter)
439  //{
440  // Pout<< "point:" << iter.key()
441  // << " coord:" << mesh.points()[iter.key()]
442  // << " regions:" << pointRegions_[iter()] << endl;
443  //}
444 }
445 
446 
447 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
448 
450 :
451  //nRegions_(0),
452  meshPointMap_(0),
453  pointRegions_(0),
454  meshFaceMap_(0),
455  faceRegions_(0)
456 {
457  const polyBoundaryMesh& patches = mesh.boundaryMesh();
458 
459  // Get any point on the outside which is on a non-coupled boundary
460  boolList candidatePoint(mesh.nPoints(), false);
461 
462  forAll(patches, patchI)
463  {
464  if (!patches[patchI].coupled())
465  {
466  const polyPatch& pp = patches[patchI];
467 
468  forAll(pp.meshPoints(), i)
469  {
470  candidatePoint[pp.meshPoints()[i]] = true;
471  }
472  }
473  }
474 
475  calcPointRegions(mesh, candidatePoint);
476 }
477 
478 
480 (
481  const polyMesh& mesh,
482  const labelList& candidatePoints
483 )
484 :
485  //nRegions_(0),
486  meshPointMap_(0),
487  pointRegions_(0),
488  meshFaceMap_(0),
489  faceRegions_(0)
490 {
491  boolList candidatePoint(mesh.nPoints(), false);
492 
493  forAll(candidatePoints, i)
494  {
495  candidatePoint[candidatePoints[i]] = true;
496  }
497 
498  calcPointRegions(mesh, candidatePoint);
499 }
500 
501 
502 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
503 
504 // Return a list (in allPatch indices) with either -1 or the face label
505 // of the face that uses the same vertices.
507 (
508  const primitiveMesh& mesh,
509  const labelList& boundaryFaces
510 )
511 {
512  // Addressing engine for all boundary faces.
513  indirectPrimitivePatch allPatch
514  (
515  IndirectList<face>(mesh.faces(), boundaryFaces),
516  mesh.points()
517  );
518 
519  labelList duplicateFace(allPatch.size(), -1);
520  label nDuplicateFaces = 0;
521 
522  // Find all duplicate faces.
523  forAll(allPatch, bFaceI)
524  {
525  const face& f = allPatch.localFaces()[bFaceI];
526 
527  // Get faces connected to f[0].
528  // Check whether share all points with f.
529  const labelList& pFaces = allPatch.pointFaces()[f[0]];
530 
531  forAll(pFaces, i)
532  {
533  label otherFaceI = pFaces[i];
534 
535  if (otherFaceI > bFaceI)
536  {
537  const face& otherF = allPatch.localFaces()[otherFaceI];
538 
539  if (isDuplicate(f, otherF, true))
540  {
542  (
543  "findDuplicateFaces(const primitiveMesh&"
544  ", const labelList&)"
545  ) << "Face:" << bFaceI + mesh.nInternalFaces()
546  << " has local points:" << f
547  << " which are in same order as face:"
548  << otherFaceI + mesh.nInternalFaces()
549  << " with local points:" << otherF
550  << abort(FatalError);
551  }
552  else if (isDuplicate(f, otherF, false))
553  {
554  label meshFace0 = bFaceI + mesh.nInternalFaces();
555  label meshFace1 = otherFaceI + mesh.nInternalFaces();
556 
557  if
558  (
559  duplicateFace[bFaceI] != -1
560  || duplicateFace[otherFaceI] != -1
561  )
562  {
564  (
565  "findDuplicateFaces(const primitiveMesh&"
566  ", const labelList&)"
567  ) << "One of two duplicate faces already marked"
568  << " as duplicate." << nl
569  << "This means that three or more faces share"
570  << " the same points and this is illegal." << nl
571  << "Face:" << meshFace0
572  << " with local points:" << f
573  << " which are in same order as face:"
574  << meshFace1
575  << " with local points:" << otherF
576  << abort(FatalError);
577  }
578 
579  duplicateFace[bFaceI] = otherFaceI;
580  duplicateFace[otherFaceI] = bFaceI;
581  nDuplicateFaces++;
582  }
583  }
584  }
585  }
586 
587  return duplicateFace;
588 }
589 
590 
592 {
593  {
594  Map<label> newMap(meshFaceMap_.size());
595 
596  forAllConstIter(Map<label>, meshFaceMap_, iter)
597  {
598  label newFaceI = map.reverseFaceMap()[iter.key()];
599 
600  if (newFaceI >= 0)
601  {
602  newMap.insert(newFaceI, iter());
603  }
604  }
605  meshFaceMap_.transfer(newMap);
606  }
607  {
608  Map<label> newMap(meshPointMap_.size());
609 
610  forAllConstIter(Map<label>, meshPointMap_, iter)
611  {
612  label newPointI = map.reversePointMap()[iter.key()];
613 
614  if (newPointI >= 0)
615  {
616  newMap.insert(newPointI, iter());
617  }
618  }
619 
620  meshPointMap_.transfer(newMap);
621  }
622 }
623 
624 
625 // ************************ vim: set sw=4 sts=4 et: ************************ //