advancingFrontAMI.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) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2020,2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "advancingFrontAMI.H"
30#include "meshTools.H"
31#include "mapDistribute.H"
32#include "unitConversion.H"
33
34#include "findNearestMaskedOp.H"
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
41}
42
43// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
44
46{
47 const auto& src = srcPatch();
48 const auto& tgt = tgtPatch();
49
50 if (debug && (!src.size() || !tgt.size()))
51 {
52 Pout<< "AMI: Patches not on processor: Source faces = "
53 << src.size() << ", target faces = " << tgt.size()
54 << endl;
55 }
56
57
58 if (requireMatch_)
59 {
60 const scalar maxBoundsError = 0.05;
61
62 // Check bounds of source and target
63 boundBox bbSrc(src.points(), src.meshPoints(), true);
64 boundBox bbTgt(tgt.points(), tgt.meshPoints(), true);
65
66 boundBox bbTgtInf(bbTgt);
67 bbTgtInf.inflate(maxBoundsError);
68
69 if (!bbTgtInf.contains(bbSrc))
70 {
72 << "Source and target patch bounding boxes are not similar"
73 << nl
74 << " source box span : " << bbSrc.span() << nl
75 << " target box span : " << bbTgt.span() << nl
76 << " source box : " << bbSrc << nl
77 << " target box : " << bbTgt << nl
78 << " inflated target box : " << bbTgtInf << endl;
79 }
80 }
81}
82
83
85(
86 const label srcFacei,
87 const label tgtFacei
88) const
89{
90 const auto& srcPatch = this->srcPatch();
91 const auto& tgtPatch = this->tgtPatch();
92
93 if
94 (
95 (srcMagSf_[srcFacei] < ROOTVSMALL)
96 || (tgtMagSf_[tgtFacei] < ROOTVSMALL)
97 )
98 {
99 return false;
100 }
101
102 if (maxDistance2_ > 0)
103 {
104 const point& srcFc = srcPatch.faceCentres()[srcFacei];
105 const point& tgtFc = tgtPatch.faceCentres()[tgtFacei];
106 const vector& srcN = srcPatch.faceNormals()[srcFacei];
107
108 const scalar normalDist((tgtFc-srcFc)&srcN);
109 //if (magSqr(srcFc-tgtFc) >= maxDistance2_)
110 if (sqr(normalDist) >= maxDistance2_)
111 {
112 return false;
113 }
114 }
115
116 if (minCosAngle_ > -1)
117 {
118 const vector& srcN = srcPatch.faceNormals()[srcFacei];
119 vector tgtN = tgtPatch.faceNormals()[tgtFacei];
120 if (!reverseTarget_)
121 {
122 tgtN = -tgtN;
123 }
124
125 if ((srcN & tgtN) <= minCosAngle_)
126 {
127 return false;
128 }
129 }
130
131 return true;
132}
133
134
136{
137 // Create processor map of extended cells. This map gets (possibly
138 // remote) cells from the src mesh such that they (together) cover
139 // all of tgt
140 extendedTgtMapPtr_.reset(calcProcMap(srcPatch0(), tgtPatch0()));
141 const mapDistribute& map = extendedTgtMapPtr_();
142
143 // Original faces from tgtPatch
144 // Note: in globalIndexing since might be remote
145 globalIndex globalTgtFaces(tgtPatch0().size());
146 distributeAndMergePatches
147 (
148 map,
149 tgtPatch0(),
150 globalTgtFaces,
151 extendedTgtFaces_,
152 extendedTgtPoints_,
153 extendedTgtFaceIDs_
154 );
155
156 // Create a representation of the tgt patch that is extended to overlap
157 // the src patch
158 extendedTgtPatchPtr_.reset
159 (
161 (
162 SubList<face>(extendedTgtFaces_),
163 extendedTgtPoints_
164 )
165 );
166}
167
168
169bool Foam::advancingFrontAMI::initialiseWalk(label& srcFacei, label& tgtFacei)
170{
171 const auto& src = this->srcPatch();
172 const auto& tgt = this->tgtPatch();
173
174 bool foundFace = false;
175
176 // Check that patch sizes are valid
177 if (!src.size())
178 {
179 return foundFace;
180 }
181 else if (!tgt.size())
182 {
184 << src.size() << " source faces but no target faces" << endl;
185
186 return foundFace;
187 }
188
189 // Reset the octree
190 treePtr_.reset(createTree(tgt));
191
192 // Find initial face match using brute force/octree search
193 if ((srcFacei == -1) || (tgtFacei == -1))
194 {
195 srcFacei = 0;
196 tgtFacei = 0;
197 forAll(src, facei)
198 {
199 tgtFacei = findTargetFace(facei);
200 if (tgtFacei >= 0)
201 {
202 srcFacei = facei;
203 foundFace = true;
204 break;
205 }
206 }
207
208 if (!foundFace)
209 {
210 if (requireMatch_)
211 {
213 << "Unable to find initial target face"
214 << abort(FatalError);
215 }
216
217 return foundFace;
218 }
219 }
220
221 if (debug)
222 {
223 Pout<< "AMI: initial target face = " << tgtFacei << endl;
224 }
225
226 return true;
227}
228
229
231(
232 const scalar area,
233 const face& f1,
234 const face& f2,
235 const pointField& f1Points,
236 const pointField& f2Points
237) const
238{
239 static label count = 1;
240
241 const pointField f1pts = f1.points(f1Points);
242 const pointField f2pts = f2.points(f2Points);
243
244 Pout<< "Face intersection area (" << count << "):" << nl
245 << " f1 face = " << f1 << nl
246 << " f1 pts = " << f1pts << nl
247 << " f2 face = " << f2 << nl
248 << " f2 pts = " << f2pts << nl
249 << " area = " << area
250 << endl;
251
252 OFstream os("areas" + name(count) + ".obj");
253
254 for (const point& pt : f1pts)
255 {
257 }
258 os<< "l";
259 forAll(f1pts, i)
260 {
261 os<< " " << i + 1;
262 }
263 os<< " 1" << endl;
264
265 for (const point& pt : f2pts)
266 {
268 }
269 os<< "l";
270 const label n = f1pts.size();
271 forAll(f2pts, i)
272 {
273 os<< " " << n + i + 1;
274 }
275 os<< " " << n + 1 << endl;
276
277 ++count;
278}
279
280
282(
283 const label srcFacei,
284 const UList<label>& excludeFaces,
285 const label srcFacePti
286) const
287{
288 const auto& src = srcPatch();
289
290 label targetFacei = -1;
291
292 const pointField& srcPts = src.points();
293 const face& srcFace = src[srcFacei];
294
295 findNearestMaskedOp<primitivePatch> fnOp(*treePtr_, excludeFaces);
296
297 const boundBox bb(srcPts, srcFace, false);
298
299 const point srcPt =
300 srcFacePti == -1 ? bb.centre() : srcPts[srcFace[srcFacePti]];
301
302 const pointIndexHit sample =
303 treePtr_->findNearest(srcPt, magSqr(bb.max() - bb.centre()), fnOp);
304
305 if (sample.hit() && isCandidate(srcFacei, sample.index()))
306 {
307 targetFacei = sample.index();
308
309 if (debug)
310 {
311 Pout<< "Source point = " << srcPt << ", Sample point = "
312 << sample.hitPoint() << ", Sample index = " << sample.index()
313 << endl;
314 }
315 }
316
317 return targetFacei;
318}
319
320
322(
323 const label facei,
324 const primitivePatch& patch,
325 const DynamicList<label>& visitedFaces,
326 DynamicList<label>& faceIDs
327) const
328{
329 static const scalar thetaCos = Foam::cos(degToRad(89.0));
330
331 const labelList& nbrFaces = patch.faceFaces()[facei];
332
333 // Filter out faces already visited from face neighbours
334 for (const label nbrFacei : nbrFaces)
335 {
336 // Prevent addition of face if it is not on the same plane-ish
337 if (!visitedFaces.found(nbrFacei) && !faceIDs.found(nbrFacei))
338 {
339 const vector& n1 = patch.faceNormals()[facei];
340 const vector& n2 = patch.faceNormals()[nbrFacei];
341
342 const scalar cosI = n1 & n2;
343
344 if (cosI > thetaCos)
345 {
346 faceIDs.append(nbrFacei);
347 }
348 }
349 }
350}
351
352
354(
355 const primitivePatch& patch,
356 List<DynamicList<face>>& tris,
357 List<scalar>& magSf
358) const
359{
360 const pointField& points = patch.points();
361 tris.setSize(patch.size());
362 magSf.setSize(patch.size());
363
364 // Using methods that index into existing points
365 forAll(patch, facei)
366 {
367 tris[facei].clear();
368
369 switch (triMode_)
370 {
372 {
373 faceAreaIntersect::triangleFan(patch[facei], tris[facei]);
374 break;
375 }
377 {
378 patch[facei].triangles(points, tris[facei]);
379 break;
380 }
381 }
382
383 const DynamicList<face>& triFaces = tris[facei];
384 magSf[facei] = 0;
385 for (const face& f : triFaces)
386 {
387 magSf[facei] +=
389 (
390 points[f[0]],
391 points[f[1]],
392 points[f[2]]
393 ).mag();
394 }
395 }
396}
397
398
400{
401 if (!requireMatch_ && distributed())
402 {
403 scalarList newTgtMagSf(std::move(tgtMagSf_));
404
405 // Assign default sizes. Override selected values with calculated
406 // values. This is to support ACMI where some of the target faces
407 // are never used (so never get sent over and hence never assigned
408 // to)
409 tgtMagSf_ = tgtPatch0().magFaceAreas();
410
411 for (const labelList& smap : this->extendedTgtMapPtr_->subMap())
412 {
413 UIndirectList<scalar>(tgtMagSf_, smap) =
414 UIndirectList<scalar>(newTgtMagSf, smap);
415 }
416 }
417}
418
419
420// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
421
423(
424 const dictionary& dict,
425 const bool reverseTarget
426)
427:
428 AMIInterpolation(dict, reverseTarget),
429 maxDistance2_(dict.getOrDefault<scalar>("maxDistance2", -1)),
430 minCosAngle_(dict.getOrDefault<scalar>("minCosAngle", -1)),
431 srcTris_(),
432 tgtTris_(),
433 extendedTgtPatchPtr_(nullptr),
434 extendedTgtFaces_(),
435 extendedTgtPoints_(),
436 extendedTgtFaceIDs_(),
437 extendedTgtMapPtr_(nullptr),
438 srcNonOverlap_(),
439 triMode_
440 (
441 faceAreaIntersect::triangulationModeNames_.getOrDefault
442 (
443 "triMode",
444 dict,
445 faceAreaIntersect::tmMesh
446 )
447 )
448{}
449
450
452(
453 const bool requireMatch,
454 const bool reverseTarget,
455 const scalar lowWeightCorrection,
457)
458:
459 AMIInterpolation(requireMatch, reverseTarget, lowWeightCorrection),
460 maxDistance2_(-1),
461 minCosAngle_(-1),
462 srcTris_(),
463 tgtTris_(),
464 extendedTgtPatchPtr_(nullptr),
465 extendedTgtFaces_(),
466 extendedTgtPoints_(),
467 extendedTgtFaceIDs_(),
468 extendedTgtMapPtr_(nullptr),
469 srcNonOverlap_(),
470 triMode_(triMode)
471{}
472
473
475:
476 AMIInterpolation(ami),
477 maxDistance2_(ami.maxDistance2_),
478 minCosAngle_(ami.minCosAngle_),
479 srcTris_(),
480 tgtTris_(),
481 extendedTgtPatchPtr_(nullptr),
482 extendedTgtFaces_(),
483 extendedTgtPoints_(),
484 extendedTgtFaceIDs_(),
485 extendedTgtMapPtr_(nullptr),
486 srcNonOverlap_(),
487 triMode_(ami.triMode_)
488{}
489
490
491// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
492
494(
495 const primitivePatch& srcPatch,
496 const primitivePatch& tgtPatch,
497 const autoPtr<searchableSurface>& surfPtr
498)
499{
500 if (AMIInterpolation::calculate(srcPatch, tgtPatch, surfPtr))
501 {
502 // Create a representation of the target patch that covers the source
503 // patch
504 if (distributed())
505 {
506 createExtendedTgtPatch();
507 }
508
509 const auto& src = this->srcPatch();
510 const auto& tgt = this->tgtPatch();
511
512
513 if (maxDistance2_ > 0)
514 {
515 // Early trigger face centre calculation
516 (void)src.faceCentres();
517 (void)tgt.faceCentres();
518 // Early trigger face normals calculation
519 (void)src.faceNormals();
520 (void)tgt.faceNormals();
521 }
522 if (minCosAngle_ > -1)
523 {
524 // Early trigger face normals calculation
525 (void)src.faceNormals();
526 (void)tgt.faceNormals();
527 }
528
529
530 // Initialise area magnitudes
531 srcMagSf_.setSize(src.size(), 1.0);
532 tgtMagSf_.setSize(tgt.size(), 1.0);
533
534 // Source and target patch triangulations
535 triangulatePatch(src, srcTris_, srcMagSf_);
536 triangulatePatch(tgt, tgtTris_, tgtMagSf_);
537
538 checkPatches();
539
540 // Set initial sizes for weights and addressing - must be done even if
541 // returns false below
542 srcAddress_.setSize(src.size());
543 srcWeights_.setSize(src.size());
544 tgtAddress_.setSize(tgt.size());
545 tgtWeights_.setSize(tgt.size());
546
547 return true;
548 }
549
550 return false;
551}
552
553
555{
557 os.writeEntryIfDifferent<scalar>("maxDistance2", -1, maxDistance2_);
558 os.writeEntryIfDifferent<scalar>("minCosAngle", -1, minCosAngle_);
560 (
561 "triMode",
564 );
565}
566
567
568// ************************************************************************* //
label n
Interpolation class dealing with transfer of data between two primitive patches with an arbitrary mes...
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
Minimal example by using system/controlDict.functions:
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Output to file stream, using an OSstream.
Definition: OFstream.H:57
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
A list of faces which address into the list of points.
A List obtained as a section of another List.
Definition: SubList.H:70
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition: IndirectList.H:79
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Base class for Arbitrary Mesh Interface (AMI) methods.
const primitivePatch & tgtPatch() const
Return const access to the target patch.
bool initialiseWalk(label &srcFacei, label &tgtFacei)
Initialise walk and return true if all ok.
void checkPatches() const
Check AMI patch coupling.
void appendNbrFaces(const label facei, const primitivePatch &patch, const DynamicList< label > &visitedFaces, DynamicList< label > &faceIDs) const
Add faces neighbouring facei to the ID list.
virtual bool calculate(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const autoPtr< searchableSurface > &surfPtr=nullptr)
Update addressing, weights and (optional) centroids.
bool isCandidate(const label srcFacei, const label tgtFacei) const
Is source/target a valid pair (i.e. not too far/different.
virtual void nonConformalCorrection()
Correction for non-conformal interpolations, e.g. for ACMI.
void createExtendedTgtPatch()
Create a map that extends tgtPatch so that it covers srcPatch.
void writeIntersectionOBJ(const scalar area, const face &f1, const face &f2, const pointField &f1Points, const pointField &f2Points) const
Write triangle intersection to OBJ file.
const primitivePatch & srcPatch() const
Return const access to the source patch.
void triangulatePatch(const primitivePatch &patch, List< DynamicList< face > > &tris, List< scalar > &magSf) const
Helper function to decompose a patch.
label findTargetFace(const label srcFacei, const UList< label > &excludeFaces=UList< label >::null(), const label srcFacePti=-1) const
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A bounding box defined in terms of min/max extrema points.
Definition: boundBox.H:64
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:97
bool contains(const point &pt) const
Contains point? (inside or on edge)
Definition: boundBoxI.H:271
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
Definition: boundBox.C:175
vector span() const
The bounding box span (from minimum to maximum)
Definition: boundBoxI.H:127
point centre() const
The centre (midpoint) of the bounding box.
Definition: boundBoxI.H:115
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Face intersection class.
static void triangleFan(const face &f, DynamicList< face > &faces)
Decompose face into triangle fan.
static const Enum< triangulationMode > triangulationModeNames_
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
pointField points(const UList< point > &pts) const
Return the points corresponding to this face.
Definition: faceI.H:87
virtual bool write()
Write the output fields.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
Class containing processor-to-processor mapping information.
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
OBJstream os(runTime.globalPath()/outputName)
const pointField & points
#define WarningInFunction
Report a warning using Foam::Warning.
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Namespace for OpenFOAM.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
errorManip< error > abort(error &err)
Definition: errorManip.H:144
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Unit conversion functions.