calculateMeshToMesh0Addressing.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 Description
25  private member of meshToMesh0.
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 "meshToMesh0.H"
32 #include "SubField.H"
33 
34 #include "indexedOctree.H"
35 #include "treeDataCell.H"
36 #include "treeDataFace.H"
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 void Foam::meshToMesh0::calcAddressing()
41 {
42  if (debug)
43  {
45  << "Calculating mesh-to-mesh cell addressing" << endl;
46  }
47 
48  // set reference to cells
49  const cellList& fromCells = fromMesh_.cells();
50  const pointField& fromPoints = fromMesh_.points();
51 
52  // In an attempt to preserve the efficiency of linear search and prevent
53  // failure, a RESCUE mechanism will be set up. Here, we shall mark all
54  // cells next to the solid boundaries. If such a cell is found as the
55  // closest, the relationship between the origin and cell will be examined.
56  // If the origin is outside the cell, a global n-squared search is
57  // triggered.
58 
59  // SETTING UP RESCUE
60 
61  // visit all boundaries and mark the cell next to the boundary.
62 
63  if (debug)
64  {
65  InfoInFunction << "Setting up rescue" << endl;
66  }
67 
68  List<bool> boundaryCell(fromCells.size(), false);
69 
70  // set reference to boundary
71  const polyPatchList& patchesFrom = fromMesh_.boundaryMesh();
72 
73  forAll(patchesFrom, patchi)
74  {
75  // get reference to cells next to the boundary
76  const labelUList& bCells = patchesFrom[patchi].faceCells();
77 
78  forAll(bCells, facei)
79  {
80  boundaryCell[bCells[facei]] = true;
81  }
82  }
83 
84  treeBoundBox meshBb(fromPoints);
85 
86  scalar typDim = meshBb.avgDim()/(2.0*cbrt(scalar(fromCells.size())));
87 
88  treeBoundBox shiftedBb
89  (
90  meshBb.min(),
91  meshBb.max() + vector(typDim, typDim, typDim)
92  );
93 
94  if (debug)
95  {
96  Info<< "\nMesh" << endl;
97  Info<< " bounding box : " << meshBb << endl;
98  Info<< " bounding box (shifted) : " << shiftedBb << endl;
99  Info<< " typical dimension :" << shiftedBb.typDim() << endl;
100  }
101 
102  indexedOctree<treeDataCell> oc
103  (
104  treeDataCell(false, fromMesh_, polyMesh::CELL_TETS),
105  shiftedBb, // overall bounding box
106  8, // maxLevel
107  10, // leafsize
108  6.0 // duplicity
109  );
110 
111  if (debug)
112  {
113  oc.print(Pout, false, 0);
114  }
115 
116  cellAddresses
117  (
118  cellAddressing_,
119  toMesh_.cellCentres(),
120  fromMesh_,
121  boundaryCell,
122  oc
123  );
124 
125  forAll(toMesh_.boundaryMesh(), patchi)
126  {
127  const polyPatch& toPatch = toMesh_.boundaryMesh()[patchi];
128 
129  if (cuttingPatches_.found(toPatch.name()))
130  {
131  boundaryAddressing_[patchi].setSize(toPatch.size());
132 
133  cellAddresses
134  (
135  boundaryAddressing_[patchi],
136  toPatch.faceCentres(),
137  fromMesh_,
138  boundaryCell,
139  oc
140  );
141  }
142  else if
143  (
144  patchMap_.found(toPatch.name())
145  && fromMeshPatches_.found(patchMap_.find(toPatch.name())())
146  )
147  {
148  const polyPatch& fromPatch = fromMesh_.boundaryMesh()
149  [
150  fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
151  ];
152 
153  if (fromPatch.empty())
154  {
156  << "Source patch " << fromPatch.name()
157  << " has no faces. Not performing mapping for it."
158  << endl;
159  boundaryAddressing_[patchi].setSize(toPatch.size());
160  boundaryAddressing_[patchi] = -1;
161  }
162  else
163  {
164  treeBoundBox wallBb(fromPatch.localPoints());
165  scalar typDim =
166  wallBb.avgDim()/(2.0*sqrt(scalar(fromPatch.size())));
167 
168  treeBoundBox shiftedBb
169  (
170  wallBb.min(),
171  wallBb.max() + vector(typDim, typDim, typDim)
172  );
173 
174  // Note: allow more levels than in meshSearch. Assume patch
175  // is not as big as all boundary faces
176  indexedOctree<treeDataFace> oc
177  (
178  treeDataFace(false, fromPatch),
179  shiftedBb, // overall search domain
180  12, // maxLevel
181  10, // leafsize
182  6.0 // duplicity
183  );
184 
185  const vectorField::subField centresToBoundary =
186  toPatch.faceCentres();
187 
188  boundaryAddressing_[patchi].setSize(toPatch.size());
189 
190  scalar distSqr = sqr(wallBb.mag());
191 
192  forAll(toPatch, toi)
193  {
194  boundaryAddressing_[patchi][toi] = oc.findNearest
195  (
196  centresToBoundary[toi],
197  distSqr
198  ).index();
199  }
200  }
201  }
202  }
203 
204  if (debug)
205  {
207  << "Finished calculating mesh-to-mesh cell addressing" << endl;
208  }
209 }
210 
211 
212 void Foam::meshToMesh0::cellAddresses
213 (
214  labelList& cellAddressing_,
215  const pointField& points,
216  const fvMesh& fromMesh,
217  const List<bool>& boundaryCell,
218  const indexedOctree<treeDataCell>& oc
219 ) const
220 {
221  // the implemented search method is a simple neighbour array search.
222  // It starts from a cell zero, searches its neighbours and finds one
223  // which is nearer to the target point than the current position.
224  // The location of the "current position" is reset to that cell and
225  // search through the neighbours continues. The search is finished
226  // when all the neighbours of the cell are farther from the target
227  // point than the current cell
228 
229  // set curCell label to zero (start)
230  label curCell = 0;
231 
232  // set reference to cell to cell addressing
233  const vectorField& centresFrom = fromMesh.cellCentres();
234  const labelListList& cc = fromMesh.cellCells();
235 
236  forAll(points, toI)
237  {
238  // pick up target position
239  const vector& p = points[toI];
240 
241  // set the sqr-distance
242  scalar distSqr = magSqr(p - centresFrom[curCell]);
243 
244  bool closer;
245 
246  do
247  {
248  closer = false;
249 
250  // set the current list of neighbouring cells
251  const labelList& neighbours = cc[curCell];
252 
253  forAll(neighbours, nI)
254  {
255  scalar curDistSqr =
256  magSqr(p - centresFrom[neighbours[nI]]);
257 
258  // search through all the neighbours.
259  // If the cell is closer, reset current cell and distance
260  if (curDistSqr < (1 - small)*distSqr)
261  {
262  curCell = neighbours[nI];
263  distSqr = curDistSqr;
264  closer = true; // a closer neighbour has been found
265  }
266  }
267  } while (closer);
268 
269  cellAddressing_[toI] = -1;
270 
271  // Check point is actually in the nearest cell
272  if (fromMesh.pointInCell(p, curCell))
273  {
274  cellAddressing_[toI] = curCell;
275  }
276  else
277  {
278  // If curCell is a boundary cell then the point maybe either outside
279  // the domain or in an other region of the doamin, either way use
280  // the octree search to find it.
281  if (boundaryCell[curCell])
282  {
283  cellAddressing_[toI] = oc.findInside(p);
284  }
285  else
286  {
287  // If not on the boundary search the neighbours
288  bool found = false;
289 
290  // set the current list of neighbouring cells
291  const labelList& neighbours = cc[curCell];
292 
293  forAll(neighbours, nI)
294  {
295  // search through all the neighbours.
296  // If point is in neighbour reset current cell
297  if (fromMesh.pointInCell(p, neighbours[nI]))
298  {
299  cellAddressing_[toI] = neighbours[nI];
300  found = true;
301  break;
302  }
303  }
304 
305  if (!found)
306  {
307  // If still not found search the neighbour-neighbours
308 
309  // set the current list of neighbouring cells
310  const labelList& neighbours = cc[curCell];
311 
312  forAll(neighbours, nI)
313  {
314  // set the current list of neighbour-neighbouring cells
315  const labelList& nn = cc[neighbours[nI]];
316 
317  forAll(nn, nI)
318  {
319  // search through all the neighbours.
320  // If point is in neighbour reset current cell
321  if (fromMesh.pointInCell(p, nn[nI]))
322  {
323  cellAddressing_[toI] = nn[nI];
324  found = true;
325  break;
326  }
327  }
328  if (found) break;
329  }
330  }
331 
332  if (!found)
333  {
334  // Still not found so us the octree
335  cellAddressing_[toI] = oc.findInside(p);
336  }
337  }
338  }
339  }
340 }
341 
342 
343 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Definition: polyPatchList.H:45
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
const cellList & cells() const
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
UList< label > labelUList
Definition: UList.H:64
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1003
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
dimensionedScalar cbrt(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: labelList.H:56
const vectorField & cellCentres() const
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(1e-3))
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
#define WarningInFunction
Report a warning using Foam::Warning.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
messageStream Info
SubField< vector > subField
Declare type of subField.
Definition: Field.H:89
Field< vector > vectorField
Specialisation of Field<T> for vector.
List< cell > cellList
list of cells
Definition: cellList.H:42
#define InfoInFunction
Report an information message using Foam::Info.