nearestFaceAMI.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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2020-2022 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "nearestFaceAMI.H"
30
31// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32
33namespace Foam
34{
38}
39
40// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41
42Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcFaceMap
43(
44 const List<nearestAndDist>& localInfo,
45 const primitivePatch& srcPatch,
46 const primitivePatch& tgtPatch
47) const
48{
49 // Generate the list of processor bounding boxes for tgtPatch
50 List<boundBox> procBbs(Pstream::nProcs());
51 procBbs[Pstream::myProcNo()] =
52 boundBox(tgtPatch.points(), tgtPatch.meshPoints(), true);
54
55 // Identify which of my local src faces intersect with each processor
56 // tgtPatch bb within the current match's search distance
57 const pointField& srcCcs = srcPatch.faceCentres();
58 List<DynamicList<label>> dynSendMap(Pstream::nProcs());
59
60 forAll(localInfo, srcFacei)
61 {
62 // Test local srcPatch face centres against remote processor tgtPatch bb
63 // using distance from initial pass
64
65 const scalar r2 = localInfo[srcFacei].second();
66
67 forAll(procBbs, proci)
68 {
69 if (proci != Pstream::myProcNo())
70 {
71 if (procBbs[proci].overlaps(srcCcs[srcFacei], r2))
72 {
73 dynSendMap[proci].append(srcFacei);
74 }
75 }
76 }
77 }
78
79 // Convert dynamicList to labelList
81 forAll(sendMap, proci)
82 {
83 dynSendMap[proci].shrink();
84 sendMap[proci].transfer(dynSendMap[proci]);
85
86 if (debug)
87 {
88 Pout<< "send map - to proc " << proci << " sending "
89 << sendMap[proci].size() << " elements" << endl;
90 }
91 }
92
93 return autoPtr<mapDistribute>::New(std::move(sendMap));
94}
95
96
97Foam::autoPtr<Foam::mapDistribute> Foam::nearestFaceAMI::calcDistributed
98(
99 const primitivePatch& src,
100 const primitivePatch& tgt,
101 labelListList& srcToTgtAddr,
102 scalarListList& srcToTgtWght
103) const
104{
105 autoPtr<indexedOctree<treeType>> tgtTreePtr;
106 if (tgt.size())
107 {
108 tgtTreePtr = this->createTree(tgt);
109 }
110
111 // Create global indexing for tgtPatch
112 globalIndex globalTgtCells(tgt.size());
113
114
115 // First pass
116 // ==========
117 // For each srcPatch face, determine local match on tgtPatch
118
119 // Identify local nearest matches
120 pointField srcCcs(src.faceCentres());
121
122 List<nearestAndDist> localInfo(src.size());
123 if (tgt.size())
124 {
125 const auto& tgtTree = tgtTreePtr();
126
127 forAll(srcCcs, srcCelli)
128 {
129 const point& srcCc = srcCcs[srcCelli];
130
131 pointIndexHit& test = localInfo[srcCelli].first();
132 test = tgtTree.findNearest(srcCc, GREAT);
133
134 if (test.hit())
135 {
136 // With a search radius2 of GREAT all cells should receive a hit
137 localInfo[srcCelli].second() = magSqr(srcCc - test.hitPoint());
138 test.setIndex(globalTgtCells.toGlobal(test.index()));
139 }
140 }
141 }
142 else
143 {
144 // No local tgtPatch faces - initialise nearest distance for all
145 // srcPatch faces to GREAT so that they [should be] set by remote
146 // tgtPatch faces
147 for (auto& info : localInfo)
148 {
149 info.second() = GREAT;
150 }
151 }
152
153
154 // Second pass
155 // ===========
156 // Determine remote matches
157
158 // Map returns labels of src patch faces to send to each proc
159 autoPtr<mapDistribute> mapPtr = calcFaceMap(localInfo, src, tgt);
160 mapDistribute& map = mapPtr();
161
162 List<nearestAndDist> remoteInfo(localInfo);
163 map.distribute(remoteInfo);
164
165 // Note: re-using srcCcs
166 map.distribute(srcCcs);
167
168 if (tgt.size())
169 {
170 const auto& tgtTree = tgtTreePtr();
171
172 // Test remote srcPatch faces against local tgtPatch faces
173 nearestAndDist testInfo;
174 pointIndexHit& test = testInfo.first();
175 forAll(remoteInfo, i)
176 {
177 test = tgtTree.findNearest(srcCcs[i], remoteInfo[i].second());
178 if (test.hit())
179 {
180 test.setIndex(globalTgtCells.toGlobal(test.index()));
181 testInfo.second() = magSqr(test.hitPoint() - srcCcs[i]);
182 nearestEqOp()(remoteInfo[i], testInfo);
183 }
184 }
185 }
186
187 // Send back to originating processor. Choose best if sent to multiple
188 // processors. Note that afterwards all unused entries have the unique
189 // value nearestZero (distance < 0). This is used later on to see if
190 // the sample was sent to any processor.
192 mapDistributeBase::distribute
193 (
196 src.size(),
197 map.constructMap(),
198 map.constructHasFlip(),
199 map.subMap(),
200 map.subHasFlip(),
201 remoteInfo,
203 nearestEqOp(),
204 identityOp() // No flipping
205 );
206
207
208 // Third pass
209 // ==========
210 // Combine local and remote info and filter out any connections that are
211 // further away than threshold distance squared
212
213 srcToTgtAddr.setSize(src.size());
214 srcToTgtWght.setSize(src.size());
215 forAll(srcToTgtAddr, srcFacei)
216 {
217 nearestEqOp()(localInfo[srcFacei], remoteInfo[srcFacei]);
218 if (localInfo[srcFacei].second() < maxDistance2_)
219 {
220 const label tgtFacei = localInfo[srcFacei].first().index();
221 srcToTgtAddr[srcFacei] = labelList(1, tgtFacei);
222 srcToTgtWght[srcFacei] = scalarList(1, 1.0);
223 }
224 }
225
226 List<Map<label>> cMap;
227 return autoPtr<mapDistribute>::New(globalTgtCells, srcToTgtAddr, cMap);
228}
229
230
231// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
232
234(
235 const dictionary& dict,
236 const bool reverseTarget
237)
238:
239 AMIInterpolation(dict, reverseTarget),
240 maxDistance2_(dict.getOrDefault<scalar>("maxDistance2", GREAT))
241{}
242
243
245(
246 const bool requireMatch,
247 const bool reverseTarget,
248 const scalar lowWeightCorrection
249)
250:
251 AMIInterpolation(requireMatch, reverseTarget, lowWeightCorrection),
252 maxDistance2_(GREAT)
253{}
254
255
257:
258 AMIInterpolation(ami),
259 maxDistance2_(ami.maxDistance2_)
260{}
261
262
263// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
264
266(
267 const primitivePatch& srcPatch,
268 const primitivePatch& tgtPatch,
269 const autoPtr<searchableSurface>& surfPtr
270)
271{
272 if (upToDate_)
273 {
274 return false;
275 }
276
277 AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr);
278
279 const auto& src = this->srcPatch0();
280 const auto& tgt = this->tgtPatch0();
281
282 // Set face area magnitudes
283 srcMagSf_ = mag(src.faceAreas());
284 tgtMagSf_ = mag(tgt.faceAreas());
285
286 // TODO: implement symmetric calculation controls; assume yes for now
287 bool symmetric_ = true;
288
289 if (this->distributed())
290 {
291 tgtMapPtr_ =
292 calcDistributed
293 (
294 src,
295 tgt,
296 srcAddress_,
297 srcWeights_
298 );
299
300 if (symmetric_)
301 {
302 srcMapPtr_ =
303 calcDistributed
304 (
305 tgt,
306 src,
307 tgtAddress_,
308 tgtWeights_
309 );
310 }
311 }
312 else
313 {
314 srcAddress_.setSize(src.size());
315 srcWeights_.setSize(src.size());
316
317 if (symmetric_)
318 {
319 tgtAddress_.setSize(tgt.size());
320 tgtWeights_.setSize(tgt.size());
321 }
322
323 const pointField& srcCcs = src.faceCentres();
324 const pointField& tgtCcs = tgt.faceCentres();
325
326 const auto tgtTreePtr = this->createTree(tgtPatch);
327 const auto& tgtTree = tgtTreePtr();
328
329 forAll(srcCcs, srcFacei)
330 {
331 const point& srcCc = srcCcs[srcFacei];
332 const pointIndexHit hit = tgtTree.findNearest(srcCc, GREAT);
333
334 if
335 (
336 hit.hit()
337 && (magSqr(srcCc - tgtCcs[hit.index()]) < maxDistance2_)
338 )
339 {
340 label tgtFacei = hit.index();
341 srcAddress_[srcFacei] = labelList(1, tgtFacei);
342 srcWeights_[srcFacei] = scalarList(1, 1.0);
343
344 if (symmetric_)
345 {
346 tgtAddress_[tgtFacei] = labelList(1, srcFacei);
347 tgtWeights_[tgtFacei] = scalarList(1, 1.0);
348 }
349 }
350 else
351 {
352 if (debug)
353 {
355 << "Unable to find target face for source face "
356 << srcFacei << endl;
357 }
358 }
359 }
360
361 if (symmetric_)
362 {
363 const auto srcTreePtr = this->createTree(srcPatch);
364 const auto& srcTree = srcTreePtr();
365
366 // Check that all source cells have connections and populate any
367 // missing entries
368 forAll(tgtWeights_, tgtCelli)
369 {
370 if (tgtAddress_[tgtCelli].empty())
371 {
372 const point& tgtCc = tgtCcs[tgtCelli];
373 pointIndexHit hit = srcTree.findNearest(tgtCc, GREAT);
374
375 if
376 (
377 hit.hit()
378 && (magSqr(tgtCc - srcCcs[hit.index()]) < maxDistance2_)
379 )
380 {
381 tgtAddress_[tgtCelli] = labelList(1, hit.index());
382 tgtWeights_[tgtCelli] = scalarList(1, 1.0);
383 }
384 }
385 else
386 {
387 if (debug)
388 {
390 << "Unable to find source face for target face "
391 << tgtCelli << endl;
392 }
393 }
394 }
395 }
396 }
397
398 srcWeightsSum_.setSize(srcWeights_.size(), 1);
399 tgtWeightsSum_.setSize(tgtWeights_.size(), 1);
400
401 upToDate_ = true;
402
403 return upToDate_;
404}
405
406
408{
410 os.writeEntryIfDifferent<scalar>("maxDistance2", GREAT, maxDistance2_);
411}
412
413
414// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
static const List< T > & null()
Return a null List.
Definition: ListI.H:109
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Ostream & writeEntryIfDifferent(const word &key, const T &value1, const T &value2)
Write a keyword/value entry only when the two values differ.
Definition: Ostream.H:251
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:66
label index() const noexcept
Return the hit index.
bool hit() const noexcept
Is there a hit?
A list of faces which address into the list of points.
label nProcs() const noexcept
Number of ranks associated with PstreamBuffers.
static void allGatherList(const List< commsStruct > &comms, List< T > &values, const int tag, const label comm)
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
@ nonBlocking
"nonBlocking"
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
virtual bool write()
Write the output fields.
Nearest-face Arbitrary Mesh Interface (AMI) method.
virtual bool calculate(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const autoPtr< searchableSurface > &surfPtr=nullptr)
Update addressing and weights.
int myProcNo() const noexcept
Return processor number.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
OBJstream os(runTime.globalPath()/outputName)
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
PointIndexHit< point > pointIndexHit
A PointIndexHit for 3D points.
Definition: pointIndexHit.H:46
void component(FieldField< Field, typename FieldField< Field, Type >::cmptType > &sf, const FieldField< Field, Type > &f, const direction d)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Tuple2< pointIndexHit, scalar > nearestAndDist
Combine operator for nearest.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:64
vector point
Point is a vector.
Definition: point.H:43
const nearestAndDist nearestZero(nearestAndDist(pointIndexHit(), -GREAT))
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
A PrimitivePatch with a SubList addressing for the faces, const reference for the point field.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
List< scalarList > scalarListList
A List of scalarList.
Definition: scalarList.H:66
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333