cellsToCells.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) 2012-2023 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 "cellsToCells.H"
27 #include "globalIndex.H"
28 #include "PatchTools.H"
29 #include "patchToPatchTools.H"
30 #include "emptyPolyPatch.H"
31 #include "wedgePolyPatch.H"
32 #include "processorPolyPatch.H"
33 #include "Time.H"
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
41 }
42 
43 
44 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
45 
47 (
48  const polyMesh& srcMesh,
49  const polyMesh& tgtMesh
50 )
51 {
52  srcLocalTgtCells_.setSize(srcMesh.nCells());
53  srcWeights_.setSize(srcMesh.nCells());
54  forAll(srcLocalTgtCells_, srcCelli)
55  {
56  srcLocalTgtCells_[srcCelli].clear();
57  srcWeights_[srcCelli].clear();
58  }
59 
60  tgtLocalSrcCells_.setSize(tgtMesh.nCells());
61  tgtWeights_.setSize(tgtMesh.nCells());
62  forAll(tgtLocalSrcCells_, tgtCelli)
63  {
64  tgtLocalSrcCells_[tgtCelli].clear();
65  tgtWeights_[tgtCelli].clear();
66  }
67 }
68 
69 
71 (
72  const polyMesh& srcMesh,
73  const polyMesh& tgtMesh
74 ) const
75 {
77  (
78  max(srcMesh.bounds().min(), tgtMesh.bounds().min()),
79  min(srcMesh.bounds().max(), tgtMesh.bounds().max())
80  );
81 
82  meshBb.inflate(0.01);
83 
84  const cellList& srcCells = srcMesh.cells();
85  const faceList& srcFaces = srcMesh.faces();
86  const pointField& srcPts = srcMesh.points();
87 
88  DynamicList<label> resultDyn(srcMesh.nCells());
89  forAll(srcCells, srcCelli)
90  {
91  const boundBox cellBb
92  (
93  srcCells[srcCelli].points(srcFaces, srcPts),
94  false
95  );
96 
97  if (meshBb.overlaps(cellBb))
98  {
99  resultDyn.append(srcCelli);
100  }
101  }
102 
103  labelList result;
104  result.transfer(resultDyn);
105  return result;
106 }
107 
108 
110 (
111  const label celli,
112  const polyMesh& mesh,
113  const DynamicList<label>& visitedCells,
114  DynamicList<label>& nbrCells
115 ) const
116 {
117  // Get all cell-cells
118  const labelList& allNbrCells = mesh.cellCells()[celli];
119 
120  // Filter out cells already visited
121  forAll(allNbrCells, i)
122  {
123  const label nbrCelli = allNbrCells[i];
124 
125  if
126  (
127  findIndex(visitedCells, nbrCelli) == -1
128  && findIndex(nbrCells, nbrCelli) == -1
129  )
130  {
131  nbrCells.append(nbrCelli);
132  }
133  }
134 }
135 
136 
137 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
138 
140 :
141  singleProcess_(-1),
142  srcLocalTgtCells_(),
143  tgtLocalSrcCells_(),
144  srcWeights_(),
145  tgtWeights_(),
146  srcMapPtr_(nullptr),
147  tgtMapPtr_(nullptr),
148  localSrcProcCellsPtr_(nullptr),
149  localTgtProcCellsPtr_(nullptr),
150  localTgtMeshPtr_(nullptr)
151 {}
152 
153 
154 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
155 
157 {}
158 
159 
160 // * * * * * * * * * * * * * * * * Selector * * * * * * * * * * * * * * * * //
161 
163 (
164  const word& cellsToCellsType
165 )
166 {
167  wordConstructorTable::iterator cstrIter =
168  wordConstructorTablePtr_->find(cellsToCellsType);
169 
170  if (cstrIter == wordConstructorTablePtr_->end())
171  {
173  << "Unknown " << typeName << " type "
174  << cellsToCellsType << endl << endl
175  << "Valid " << typeName << " types are : " << endl
176  << wordConstructorTablePtr_->sortedToc()
177  << exit(FatalError);
178  }
179 
180  return cstrIter()();
181 }
182 
183 
184 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
185 
187 {
188  PackedBoolList result(srcLocalTgtCells_.size());
189  forAll(srcLocalTgtCells_, srcCelli)
190  {
191  result[srcCelli] = !srcLocalTgtCells_[srcCelli].empty();
192  }
193  return result;
194 }
195 
196 
198 {
199  PackedBoolList result(tgtLocalSrcCells_.size());
200  forAll(tgtLocalSrcCells_, tgtCelli)
201  {
202  result[tgtCelli] = !tgtLocalSrcCells_[tgtCelli].empty();
203  }
204  return result;
205 }
206 
207 
209 (
210  const polyMesh& tgtMesh,
211  const label srcCelli,
212  const point& p
213 ) const
214 {
215  forAll(srcLocalTgtCells_[srcCelli], i)
216  {
217  const label tgtCelli = srcLocalTgtCells_[srcCelli][i];
218 
219  const polyMesh& localTgtMesh =
220  singleProcess_ == -1 ? localTgtMeshPtr_() : tgtMesh;
221 
222  if (localTgtMesh.pointInCell(p, tgtCelli))
223  {
224  return
225  singleProcess_ == -1
226  ? localTgtProcCellsPtr_()[tgtCelli]
227  : remote(Pstream::myProcNo(), tgtCelli);
228  }
229  }
230 
231  return remote();
232 }
233 
234 
236 (
237  const polyMesh& srcMesh,
238  const polyMesh& tgtMesh
239 )
240 {
241  cpuTime time;
242 
243  // Determine numbers of faces on both sides, report, and quit if either
244  // side is empty
245  const label srcTotalSize = returnReduce(srcMesh.nCells(), sumOp<label>());
246  const label tgtTotalSize = returnReduce(tgtMesh.nCells(), sumOp<label>());
247  if (srcTotalSize == 0 || tgtTotalSize == 0)
248  {
249  return 0;
250  }
251 
252  Info<< indent << typeName << ": Calculating couplings between "
253  << srcTotalSize << " source cells and " << tgtTotalSize
254  << " target cells" << incrIndent << endl;
255 
256  singleProcess_ =
258  (
259  srcMesh.nCells(),
260  tgtMesh.nCells()
261  );
262 
263  scalar V = 0;
264 
265  if (isSingleProcess())
266  {
267  // Do the intersection
268  V = calculate(srcMesh, tgtMesh);
269 
270  // Normalise the weights
271  normalise(srcMesh, srcLocalTgtCells_, srcWeights_);
272  normalise(tgtMesh, tgtLocalSrcCells_, tgtWeights_);
273  }
274  else
275  {
276  // Create the target map of overlapping cells. This map gets remote
277  // parts of the target mesh so that everything needed to compute an
278  // intersection is available locally to the source. Use it to create a
279  // source-local target mesh.
280  tgtMapPtr_ =
282  (
283  tgtMeshSendCells(srcMesh, tgtMesh)
284  );
285  localTgtProcCellsPtr_.reset
286  (
287  new List<remote>
288  (
289  distributeMesh
290  (
291  tgtMapPtr_(),
292  tgtMesh,
293  localTgtMeshPtr_
294  )
295  )
296  );
297  const polyMesh& localTgtMesh = localTgtMeshPtr_();
298 
299  if (debug > 1)
300  {
301  Pout<< "Writing local target mesh: "
302  << localTgtMesh.name() << endl;
303  localTgtMesh.write();
304  }
305 
306  // Do the intersection
307  V = calculate(srcMesh, localTgtMesh);
308 
309  // Trim the local target mesh
310  trimLocalTgt();
311 
312  if (debug > 1)
313  {
314  Pout<< "Writing trimmed local target mesh: "
315  << localTgtMesh.name() << endl;
316  localTgtMesh.write();
317  }
318 
319  // Construct the source map
320  srcMapPtr_ =
322  (
324  (
325  tgtLocalSrcCells_,
326  localTgtProcCellsPtr_()
327  )
328  );
329  localSrcProcCellsPtr_.reset
330  (
331  new List<remote>
332  (
334  )
335  );
336 
337  // Collect the addressing on the target
339  (
340  tgtMesh.nCells(),
341  tgtMapPtr_(),
342  localSrcProcCellsPtr_(),
343  tgtLocalSrcCells_
344  );
345 
346  // Collect the weights on the target
348  (
349  tgtMesh.nCells(),
350  tgtMapPtr_(),
351  tgtWeights_
352  );
353 
354  // Normalise the weights
355  normalise(srcMesh, srcLocalTgtCells_, srcWeights_);
356  normalise(tgtMesh, tgtLocalSrcCells_, tgtWeights_);
357 
358  // Collect volume intersection contributions
359  reduce(V, sumOp<scalar>());
360  }
361 
362  label nCouples = 0;
363  forAll(srcLocalTgtCells_, srcCelli)
364  {
365  nCouples += srcLocalTgtCells_[srcCelli].size();
366  }
367  forAll(tgtLocalSrcCells_, tgtCelli)
368  {
369  nCouples += tgtLocalSrcCells_[tgtCelli].size();
370  }
371  reduce(nCouples, sumOp<label>());
372 
373  if (nCouples != 0)
374  {
375  Info<< indent << "Overlapping volume = " << V << endl
376  << indent << nCouples << " couplings calculated in "
377  << time.cpuTimeIncrement() << 's' << endl;
378  }
379  else
380  {
381  Info<< indent << "No couplings found" << endl;
382  }
383 
384  Info<< decrIndent;
385 
386  return V;
387 }
388 
389 
390 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
const word & name() const
Return name.
Definition: IOobject.H:310
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
void setSize(const label)
Reset size of List.
Definition: List.C:281
A bit-packed bool list.
bool empty() const
Return true if the list is empty (ie, size() is zero).
Definition: PackedListI.H:718
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:59
const point & min() const
Minimum point defining the bounding box.
Definition: boundBoxI.H:60
const point & max() const
Maximum point defining the bounding box.
Definition: boundBoxI.H:66
Class to calculate interpolative addressing and weights between the cells of two overlapping meshes.
Definition: cellsToCells.H:56
cellsToCells()
Construct null.
Definition: cellsToCells.C:139
labelListList tgtLocalSrcCells_
For each target cell, the coupled local source cells.
Definition: cellsToCells.H:68
virtual ~cellsToCells()
Destructor.
Definition: cellsToCells.C:156
labelList maskCells(const polyMesh &srcMesh, const polyMesh &tgtMesh) const
Return src cell IDs for the overlap region.
Definition: cellsToCells.C:71
scalar update(const polyMesh &srcMesh, const polyMesh &tgtMesh)
Update addressing and weights for the given meshes. Returns the.
Definition: cellsToCells.C:236
void appendNbrCells(const label tgtCelli, const polyMesh &mesh, const DynamicList< label > &visitedTgtCells, DynamicList< label > &nbrTgtCellIDs) const
Append target cell neighbour cells to cellIDs list.
Definition: cellsToCells.C:110
remote srcToTgtPoint(const polyMesh &tgtMesh, const label srcCelli, const point &p) const
Find the target processor and cell associated with a point in a.
Definition: cellsToCells.C:209
PackedBoolList srcCoupled() const
Return a list indicating which source cells are coupled.
Definition: cellsToCells.C:186
static autoPtr< cellsToCells > New(const word &cellsToCellsType)
Select from name.
Definition: cellsToCells.C:163
scalarListList tgtWeights_
For each target cell, the coupled source weights.
Definition: cellsToCells.H:74
labelListList srcLocalTgtCells_
For each source cell, the coupled local target cells.
Definition: cellsToCells.H:65
scalarListList srcWeights_
For each source cell, the coupled target weights.
Definition: cellsToCells.H:71
PackedBoolList tgtCoupled() const
Return a list indicating which target cells are coupled.
Definition: cellsToCells.C:197
void initialise(const polyMesh &srcMesh, const polyMesh &tgtMesh)
Initialise the addressing and weights.
Definition: cellsToCells.C:47
Starts timing CPU usage and return elapsed time from start.
Definition: cpuTime.H:55
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:60
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1326
bool pointInCell(const point &p, label celli, const cellDecomposition=CELL_TETS) const
Test if point p is in the celli.
Definition: polyMesh.C:1617
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1313
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:410
label nCells() const
const labelListList & cellCells() const
const cellList & cells() const
virtual bool write(const bool write=true) const
Write using setting from DB.
Struct for keeping processor, element (cell, face, point) index.
Definition: remote.H:57
A class for handling words, derived from string.
Definition: word.H:62
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
const pointField & points
label calculate(const fvMesh &mesh, const labelHashSet &patchIDs, const scalar minFaceFraction, GeometricField< scalar, PatchField, GeoMesh > &distance)
Calculate distance data from patches.
label singleProcess(const label srcSize, const label tgtSize)
Determine whether this intersection is confined to a single processor or.
autoPtr< distributionMap > constructDistributionMap(const labelListList &procSendIndices)
Turn a list of send-to-processor indices into a distribution map.
List< remote > distributeAddressing(const distributionMap &map)
Construct local addressing from the given distribution map. The result is a.
void rDistributeTgtAddressing(const label tgtSize, const distributionMap &tgtMap, const List< remote > &localSrcProcFaces, labelListList &tgtLocalSrcFaces)
Reverse distribute a set of target addressing.
static void rDistributeListList(const label size, const distributionMap &map, List< List< Type >> &data)
Reverse distribute a list-list given the map.
labelListList procSendIndices(const labelListList &tgtLocalSrcFaces, const List< remote > &localTgtProcFaces)
Given a local intersection addressing, determine what elements need sending.
Namespace for OpenFOAM.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:241
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
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
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
messageStream Info
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:234
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
error FatalError
quaternion normalise(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:603
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:227
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(1e-3))
volScalarField & p