matchingCellsToCells.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) 2013-2024 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 "matchingCellsToCells.H"
27 #include "indexedOctree.H"
28 #include "treeDataCell.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace cellsToCellss
36 {
39 }
40 }
41 
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 bool Foam::cellsToCellss::matching::intersect
46 (
47  const polyMesh& srcMesh,
48  const polyMesh& tgtMesh,
49  const label srcCelli,
50  const label tgtCelli
51 ) const
52 {
53  return tgtMesh.pointInCell
54  (
55  srcMesh.cellCentres()[srcCelli],
56  tgtCelli,
58  );
59 }
60 
61 
62 bool Foam::cellsToCellss::matching::findInitialSeeds
63 (
64  const polyMesh& srcMesh,
65  const polyMesh& tgtMesh,
66  const labelList& srcCellIDs,
67  const boolList& mapFlag,
68  const label startSeedI,
69  label& srcSeedI,
70  label& tgtSeedI
71 ) const
72 {
73  const cellList& srcCells = srcMesh.cells();
74  const faceList& srcFaces = srcMesh.faces();
75  const pointField& srcPts = srcMesh.points();
76 
77  for (label i = startSeedI; i < srcCellIDs.size(); i++)
78  {
79  label srcI = srcCellIDs[i];
80 
81  if (mapFlag[srcI])
82  {
83  const point srcCtr(srcCells[srcI].centre(srcPts, srcFaces));
84  label tgtI = tgtMesh.cellTree().findInside(srcCtr);
85 
86  if (tgtI != -1 && intersect(srcMesh, tgtMesh, srcI, tgtI))
87  {
88  srcSeedI = srcI;
89  tgtSeedI = tgtI;
90 
91  return true;
92  }
93  }
94  }
95 
96  if (debug)
97  {
98  Pout<< "could not find starting seed" << endl;
99  }
100 
101  return false;
102 }
103 
104 
105 Foam::scalar Foam::cellsToCellss::matching::calculateAddressing
106 (
107  const polyMesh& srcMesh,
108  const polyMesh& tgtMesh,
109  labelListList& srcToTgtCellAddr,
110  scalarListList& srcToTgtCellWght,
111  labelListList& tgtToSrcCellAddr,
112  scalarListList& tgtToSrcCellWght,
113  const label srcSeedI,
114  const label tgtSeedI,
115  const labelList& srcCellIDs, // not used
116  boolList& mapFlag,
117  label& startSeedI
118 )
119 {
120  scalar V = 0;
121 
122  // store a list of src cells already mapped
123  labelList srcTgtSeed(srcMesh.nCells(), -1);
124 
125  List<DynamicList<label>> srcToTgt(srcMesh.nCells());
126  List<DynamicList<label>> tgtToSrc(tgtMesh.nCells());
127 
128  DynamicList<label> srcSeeds(10);
129 
130  const scalarField& srcVc = srcMesh.cellVolumes();
131  const scalarField& tgtVc = tgtMesh.cellVolumes();
132 
133  label srcCelli = srcSeedI;
134  label tgtCelli = tgtSeedI;
135 
136  do
137  {
138  // store src/tgt cell pair
139  srcToTgt[srcCelli].append(tgtCelli);
140  tgtToSrc[tgtCelli].append(srcCelli);
141 
142  // mark source cell srcSeedI as matched
143  mapFlag[srcCelli] = false;
144 
145  // accumulate intersection volume
146  V += srcVc[srcCelli];
147 
148  // find new source seed cell
149  appendToDirectSeeds
150  (
151  srcMesh,
152  tgtMesh,
153  mapFlag,
154  srcTgtSeed,
155  srcSeeds,
156  srcCelli,
157  tgtCelli
158  );
159  }
160  while (srcCelli >= 0);
161 
162  // transfer addressing into persistent storage
163  forAll(srcToTgtCellAddr, i)
164  {
165  srcToTgtCellWght[i] = scalarList(srcToTgt[i].size(), srcVc[i]);
166  srcToTgtCellAddr[i].transfer(srcToTgt[i]);
167  }
168 
169  forAll(tgtToSrcCellAddr, i)
170  {
171  tgtToSrcCellWght[i] = scalarList(tgtToSrc[i].size(), tgtVc[i]);
172  tgtToSrcCellAddr[i].transfer(tgtToSrc[i]);
173  }
174 
175  return V;
176 }
177 
178 
179 void Foam::cellsToCellss::matching::appendToDirectSeeds
180 (
181  const polyMesh& srcMesh,
182  const polyMesh& tgtMesh,
183  boolList& mapFlag,
184  labelList& srcTgtSeed,
185  DynamicList<label>& srcSeeds,
186  label& srcSeedI,
187  label& tgtSeedI
188 ) const
189 {
190  const labelList& srcNbr = srcMesh.cellCells()[srcSeedI];
191  const labelList& tgtNbr = tgtMesh.cellCells()[tgtSeedI];
192 
193  forAll(srcNbr, i)
194  {
195  label srcI = srcNbr[i];
196 
197  if (mapFlag[srcI] && (srcTgtSeed[srcI] == -1))
198  {
199  // source cell srcI not yet mapped
200 
201  // identify if target cell exists for source cell srcI
202  bool found = false;
203  forAll(tgtNbr, j)
204  {
205  label tgtI = tgtNbr[j];
206 
207  if (intersect(srcMesh, tgtMesh, srcI, tgtI))
208  {
209  // new match - append to lists
210  found = true;
211 
212  srcTgtSeed[srcI] = tgtI;
213  srcSeeds.append(srcI);
214 
215  break;
216  }
217  }
218 
219  if (!found)
220  {
221  // no match available for source cell srcI
222  mapFlag[srcI] = false;
223  }
224  }
225  }
226 
227  if (srcSeeds.size())
228  {
229  srcSeedI = srcSeeds.remove();
230  tgtSeedI = srcTgtSeed[srcSeedI];
231  }
232  else
233  {
234  srcSeedI = -1;
235  tgtSeedI = -1;
236  }
237 }
238 
239 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
240 
242 (
243  const polyMesh& srcMesh,
244  const polyMesh& tgtMesh
245 )
246 {
247  initialise(srcMesh, tgtMesh);
248 
249  // Determine (potentially) participating source mesh cells
250  const labelList srcCellIDs(maskCells(srcMesh, tgtMesh));
251 
252  // Initialise list to keep track of whether src cell can be mapped
253  boolList mapFlag(srcMesh.nCells(), false);
254  UIndirectList<bool>(mapFlag, srcCellIDs) = true;
255 
256  // Find initial point in tgt mesh
257  label srcSeedI = -1;
258  label tgtSeedI = -1;
259  label startSeedI = 0;
260  bool startWalk =
261  findInitialSeeds
262  (
263  srcMesh,
264  tgtMesh,
265  srcCellIDs,
266  mapFlag,
267  startSeedI,
268  srcSeedI,
269  tgtSeedI
270  );
271 
272  if (startWalk)
273  {
274  return
275  calculateAddressing
276  (
277  srcMesh,
278  tgtMesh,
279  srcLocalTgtCells_,
280  srcWeights_,
281  tgtLocalSrcCells_,
282  tgtWeights_,
283  srcSeedI,
284  tgtSeedI,
285  srcCellIDs,
286  mapFlag,
287  startSeedI
288  );
289  }
290  else
291  {
292  return 0;
293  }
294 }
295 
296 
298 (
299  const polyMesh& srcMesh,
300  labelListList& srcToTgtAddr,
301  scalarListList& srcToTgtWght
302 ) const
303 {
304  forAll(srcToTgtWght, srcCelli)
305  {
306  if (srcToTgtWght[srcCelli].size() > 1)
307  {
308  srcToTgtAddr[srcCelli].resize(1);
309  srcToTgtWght[srcCelli].resize(1);
310  srcToTgtWght[srcCelli][0] = 1;
311  }
312  }
313 }
314 
315 
316 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
317 
319 :
320  cellsToCells()
321 {}
322 
323 
324 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
325 
327 {}
328 
329 
330 // ************************************************************************* //
bool found
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
void resize(const label)
Alias for setSize(const label)
Definition: ListI.H:138
A List with indirect addressing.
Definition: UIndirectList.H:60
Class to calculate interpolative addressing and weights between the cells of two overlapping meshes.
Definition: cellsToCells.H:56
Matching, one-to-one, cells-to-cells interpolation class.
virtual scalar calculate(const polyMesh &srcMesh, const polyMesh &tgtMesh)
Calculate the addressing and weights.
virtual void normalise(const polyMesh &mesh, labelListList &localOtherCells, scalarListList &weights) const
Normalise the weights for a given mesh.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
label nCells() const
A class for handling words, derived from string.
Definition: word.H:62
defineTypeNameAndDebug(intersection, 0)
addToRunTimeSelectionTable(cellsToCells, intersection, word)
Namespace for OpenFOAM.
List< scalarList > scalarListList
Definition: scalarList.H:51
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
List< cell > cellList
list of cells
Definition: cellList.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
vector point
Point is a vector.
Definition: point.H:41
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
List< face > faceList
Definition: faceListFwd.H:41