FreeFOAM The Cross-Platform CFD Toolkit
calculateMeshToMeshAddressing.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 Description
25  private member of meshToMesh.
26  Calculates mesh to mesh addressing pattern (for each cell from one mesh
27  find the closest cell centre in the other mesh).
28 
29 \*---------------------------------------------------------------------------*/
30 
31 #include "meshToMesh.H"
32 #include <OpenFOAM/SubField.H>
33 
34 #include <meshTools/octree.H>
37 
38 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 void meshToMesh::calcAddressing()
46 {
47  if (debug)
48  {
49  Info<< "meshToMesh::calculateAddressing() : "
50  << "calculating mesh-to-mesh cell addressing" << endl;
51  }
52 
53  // set reference to cells
54  const cellList& fromCells = fromMesh_.cells();
55  const pointField& fromPoints = fromMesh_.points();
56 
57  // In an attempt to preserve the efficiency of linear search and prevent
58  // failure, a RESCUE mechanism will be set up. Here, we shall mark all
59  // cells next to the solid boundaries. If such a cell is found as the
60  // closest, the relationship between the origin and cell will be examined.
61  // If the origin is outside the cell, a global n-squared search is
62  // triggered.
63 
64  // SETTING UP RESCUE
65 
66  // visit all boundaries and mark the cell next to the boundary.
67 
68  if (debug)
69  {
70  Info<< "meshToMesh::calculateAddressing() : "
71  << "Setting up rescue" << endl;
72  }
73 
74  List<bool> boundaryCell(fromCells.size(), false);
75 
76  // set reference to boundary
77  const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
78 
79  forAll (patchesFrom, patchI)
80  {
81  // get reference to cells next to the boundary
82  const unallocLabelList& bCells = patchesFrom[patchI].faceCells();
83 
84  forAll (bCells, faceI)
85  {
86  boundaryCell[bCells[faceI]] = true;
87  }
88  }
89 
90  treeBoundBox meshBb(fromPoints);
91 
92  scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
93 
94  treeBoundBox shiftedBb
95  (
96  meshBb.min(),
97  meshBb.max() + vector(typDim, typDim, typDim)
98  );
99 
100  if (debug)
101  {
102  Info<< "\nMesh" << endl;
103  Info<< " bounding box : " << meshBb << endl;
104  Info<< " bounding box (shifted) : " << shiftedBb << endl;
105  Info<< " typical dimension :" << shiftedBb.typDim() << endl;
106  }
107 
108  // Wrap indices and mesh information into helper object
109  octreeDataCell shapes(fromMesh_);
110 
111  octree<octreeDataCell> oc
112  (
113  shiftedBb, // overall bounding box
114  shapes, // all information needed to do checks on cells
115  1, // min. levels
116  10.0, // max. size of leaves
117  10.0 // maximum ratio of cubes v.s. cells
118  );
119 
120  if (debug)
121  {
122  oc.printStats(Info);
123  }
124 
125  cellAddresses
126  (
127  cellAddressing_,
128  toMesh_.cellCentres(),
129  fromMesh_,
130  boundaryCell,
131  oc
132  );
133 
134  forAll (toMesh_.boundaryMesh(), patchi)
135  {
136  const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
137 
138  if (cuttingPatches_.found(toPatch.name()))
139  {
140  boundaryAddressing_[patchi].setSize(toPatch.size());
141 
142  cellAddresses
143  (
144  boundaryAddressing_[patchi],
145  toPatch.faceCentres(),
146  fromMesh_,
147  boundaryCell,
148  oc
149  );
150  }
151  else if
152  (
153  patchMap_.found(toPatch.name())
154  && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
155  )
156  {
157  const polyPatch& fromPatch = fromMesh_.boundaryMesh()
158  [
159  fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
160  ];
161 
162  if (fromPatch.empty())
163  {
164  WarningIn("meshToMesh::calcAddressing()")
165  << "Source patch " << fromPatch.name()
166  << " has no faces. Not performing mapping for it."
167  << endl;
168  boundaryAddressing_[patchi] = -1;
169  }
170  else
171  {
172  treeBoundBox wallBb(fromPatch.localPoints());
173  scalar typDim =
174  wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
175 
176  treeBoundBox shiftedBb
177  (
178  wallBb.min(),
179  wallBb.max() + vector(typDim, typDim, typDim)
180  );
181 
182  // Wrap data for octree into container
183  octreeDataFace shapes(fromPatch);
184 
185  octree<octreeDataFace> oc
186  (
187  shiftedBb, // overall search domain
188  shapes, // all information needed to do checks on cells
189  1, // min levels
190  20.0, // maximum ratio of cubes v.s. cells
191  2.0
192  );
193 
194 
195  const vectorField::subField centresToBoundary =
196  toPatch.faceCentres();
197 
198  boundaryAddressing_[patchi].setSize(toPatch.size());
199 
200  treeBoundBox tightest;
201  scalar tightestDist;
202 
203  forAll(toPatch, toi)
204  {
205  tightest = wallBb; // starting search bb
206  tightestDist = Foam::GREAT; // starting max distance
207 
208  boundaryAddressing_[patchi][toi] = oc.findNearest
209  (
210  centresToBoundary[toi],
211  tightest,
212  tightestDist
213  );
214  }
215  }
216  }
217  }
218 
219  if (debug)
220  {
221  Info<< "meshToMesh::calculateAddressing() : "
222  << "finished calculating mesh-to-mesh acell ddressing" << endl;
223  }
224 }
225 
226 
227 void meshToMesh::cellAddresses
228 (
229  labelList& cellAddressing_,
230  const pointField& points,
231  const fvMesh& fromMesh,
232  const List<bool>& boundaryCell,
233  const octree<octreeDataCell>& oc
234 ) const
235 {
236  // the implemented search method is a simple neighbour array search.
237  // It starts from a cell zero, searches its neighbours and finds one
238  // which is nearer to the target point than the current position.
239  // The location of the "current position" is reset to that cell and
240  // search through the neighbours continues. The search is finished
241  // when all the neighbours of the cell are farther from the target
242  // point than the current cell
243 
244  // set curCell label to zero (start)
245  register label curCell = 0;
246 
247  // set reference to cell to cell addressing
248  const vectorField& centresFrom = fromMesh.cellCentres();
249  const labelListList& cc = fromMesh.cellCells();
250 
251  forAll (points, toI)
252  {
253  // pick up target position
254  const vector& p = points[toI];
255 
256  // set the sqr-distance
257  scalar distSqr = magSqr(p - centresFrom[curCell]);
258 
259  bool closer;
260 
261  do
262  {
263  closer = false;
264 
265  // set the current list of neighbouring cells
266  const labelList& neighbours = cc[curCell];
267 
268  forAll (neighbours, nI)
269  {
270  scalar curDistSqr =
271  magSqr(p - centresFrom[neighbours[nI]]);
272 
273  // search through all the neighbours.
274  // If the cell is closer, reset current cell and distance
275  if (curDistSqr < (1 - SMALL)*distSqr)
276  {
277  curCell = neighbours[nI];
278  distSqr = curDistSqr;
279  closer = true; // a closer neighbour has been found
280  }
281  }
282  } while (closer);
283 
284  cellAddressing_[toI] = -1;
285 
286  // Check point is actually in the nearest cell
287  if (fromMesh.pointInCell(p, curCell))
288  {
289  cellAddressing_[toI] = curCell;
290  }
291  else
292  {
293  // If curCell is a boundary cell then the point maybe either outside
294  // the domain or in an other region of the doamin, either way use
295  // the octree search to find it.
296  if (boundaryCell[curCell])
297  {
298  cellAddressing_[toI] = oc.find(p);
299  }
300  else
301  {
302  // If not on the boundary search the neighbours
303  bool found = false;
304 
305  // set the current list of neighbouring cells
306  const labelList& neighbours = cc[curCell];
307 
308  forAll (neighbours, nI)
309  {
310  // search through all the neighbours.
311  // If point is in neighbour reset current cell
312  if (fromMesh.pointInCell(p, neighbours[nI]))
313  {
314  cellAddressing_[toI] = neighbours[nI];
315  found = true;
316  break;
317  }
318  }
319 
320  if (!found)
321  {
322  // If still not found search the neighbour-neighbours
323 
324  // set the current list of neighbouring cells
325  const labelList& neighbours = cc[curCell];
326 
327  forAll (neighbours, nI)
328  {
329  // set the current list of neighbour-neighbouring cells
330  const labelList& nn = cc[neighbours[nI]];
331 
332  forAll (nn, nI)
333  {
334  // search through all the neighbours.
335  // If point is in neighbour reset current cell
336  if (fromMesh.pointInCell(p, nn[nI]))
337  {
338  cellAddressing_[toI] = nn[nI];
339  found = true;
340  break;
341  }
342  }
343  if (found) break;
344  }
345  }
346 
347  if (!found)
348  {
349  // Still not found so us the octree
350  cellAddressing_[toI] = oc.find(p);
351  }
352  }
353  }
354  }
355 }
356 
357 
358 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
359 
360 } // End namespace Foam
361 
362 // ************************ vim: set sw=4 sts=4 et: ************************ //