AMIMethod.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-2022 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 "AMIMethod.H"
27 #include "meshTools.H"
28 #include "distributionMap.H"
29 #include "unitConversion.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(AMIMethod, 0);
36  defineRunTimeSelectionTable(AMIMethod, components);
37 }
38 
39 
40 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
41 
43 {
44  if (debug && (!srcPatch_.size() || !tgtPatch_.size()))
45  {
46  Pout<< "AMI: Patches not on processor: Source faces = "
47  << srcPatch_.size() << ", target faces = " << tgtPatch_.size()
48  << endl;
49  }
50 
51 
52  if (conformal())
53  {
54  const scalar maxBoundsError = 0.05;
55 
56  // check bounds of source and target
57  boundBox bbSrc(srcPatch_.points(), srcPatch_.meshPoints(), true);
58  boundBox bbTgt(tgtPatch_.points(), tgtPatch_.meshPoints(), true);
59 
60  boundBox bbTgtInf(bbTgt);
61  bbTgtInf.inflate(maxBoundsError);
62 
63  if (!bbTgtInf.contains(bbSrc))
64  {
66  << "Source and target patch bounding boxes are not similar"
67  << nl
68  << " source box span : " << bbSrc.span() << nl
69  << " target box span : " << bbTgt.span() << nl
70  << " source box : " << bbSrc << nl
71  << " target box : " << bbTgt << nl
72  << " inflated target box : " << bbTgtInf << endl;
73  }
74  }
75 }
76 
77 
79 (
80  labelListList& srcAddress,
81  scalarListList& srcWeights,
82  labelListList& tgtAddress,
83  scalarListList& tgtWeights,
84  label& srcFacei,
85  label& tgtFacei
86 )
87 {
88  checkPatches();
89 
90  // set initial sizes for weights and addressing - must be done even if
91  // returns false below
92  srcAddress.setSize(srcPatch_.size());
93  srcWeights.setSize(srcPatch_.size());
94  tgtAddress.setSize(tgtPatch_.size());
95  tgtWeights.setSize(tgtPatch_.size());
96 
97  // check that patch sizes are valid
98  if (!srcPatch_.size())
99  {
100  return false;
101  }
102  else if (!tgtPatch_.size())
103  {
105  << srcPatch_.size() << " source faces but no target faces" << endl;
106 
107  return false;
108  }
109 
110  // reset the octree
111  resetTree();
112 
113  // find initial face match using brute force/octree search
114  if ((srcFacei == -1) || (tgtFacei == -1))
115  {
116  srcFacei = 0;
117  tgtFacei = 0;
118  bool foundFace = false;
119  forAll(srcPatch_, facei)
120  {
121  tgtFacei = findTargetFace(facei);
122  if (tgtFacei >= 0)
123  {
124  srcFacei = facei;
125  foundFace = true;
126  break;
127  }
128  }
129 
130  if (!foundFace)
131  {
132  if (requireMatch_)
133  {
135  << "Unable to find initial target face"
136  << abort(FatalError);
137  }
138 
139  return false;
140  }
141  }
142 
143  if (debug)
144  {
145  Pout<< "AMI: initial target face = " << tgtFacei << endl;
146  }
147 
148  return true;
149 }
150 
151 
153 (
154  const scalar area,
155  const face& f1,
156  const face& f2,
157  const pointField& f1Points,
158  const pointField& f2Points
159 ) const
160 {
161  static label count = 1;
162 
163  const pointField f1pts = f1.points(f1Points);
164  const pointField f2pts = f2.points(f2Points);
165 
166  Pout<< "Face intersection area (" << count << "):" << nl
167  << " f1 face = " << f1 << nl
168  << " f1 pts = " << f1pts << nl
169  << " f2 face = " << f2 << nl
170  << " f2 pts = " << f2pts << nl
171  << " area = " << area
172  << endl;
173 
174  OFstream os("areas" + name(count) + ".obj");
175 
176  forAll(f1pts, i)
177  {
178  meshTools::writeOBJ(os, f1pts[i]);
179  }
180  os<< "l";
181  forAll(f1pts, i)
182  {
183  os<< " " << i + 1;
184  }
185  os<< " 1" << endl;
186 
187 
188  forAll(f2pts, i)
189  {
190  meshTools::writeOBJ(os, f2pts[i]);
191  }
192  os<< "l";
193  forAll(f2pts, i)
194  {
195  os<< " " << f1pts.size() + i + 1;
196  }
197  os<< " " << f1pts.size() + 1 << endl;
198 
199  count++;
200 }
201 
202 
204 {
205  // Clear the old octree
206  treePtr_.clear();
207 
208  treeBoundBox bb(tgtPatch_.points(), tgtPatch_.meshPoints());
209  bb.inflate(0.01);
210 
211  if (!treePtr_.valid())
212  {
213  treePtr_.reset
214  (
216  (
217  treeType
218  (
219  false,
220  tgtPatch_,
222  ),
223  bb, // overall search domain
224  8, // maxLevel
225  10, // leaf size
226  3.0 // duplicity
227  )
228  );
229  }
230 }
231 
232 
234 (
235  const label srcFacei
236 ) const
237 {
238  label targetFacei = -1;
239 
240  const pointField& srcPts = srcPatch_.points();
241  const face& srcFace = srcPatch_[srcFacei];
242  const point srcPt = srcFace.centre(srcPts);
243  const scalar srcFaceArea = srcMagSf_[srcFacei];
244 
245  pointIndexHit sample = treePtr_->findNearest(srcPt, 10.0*srcFaceArea);
246 
247  if (sample.hit())
248  {
249  targetFacei = sample.index();
250 
251  if (debug)
252  {
253  Pout<< "Source point = " << srcPt << ", Sample point = "
254  << sample.hitPoint() << ", Sample index = " << sample.index()
255  << endl;
256  }
257  }
258 
259  return targetFacei;
260 }
261 
262 
264 (
265  const label facei,
266  const primitivePatch& patch,
267  const DynamicList<label>& visitedFaces,
268  DynamicList<label>& faceIDs
269 ) const
270 {
271  const labelList& nbrFaces = patch.faceFaces()[facei];
272 
273  // filter out faces already visited from face neighbours
274  forAll(nbrFaces, i)
275  {
276  label nbrFacei = nbrFaces[i];
277  bool valid = true;
278  forAll(visitedFaces, j)
279  {
280  if (nbrFacei == visitedFaces[j])
281  {
282  valid = false;
283  break;
284  }
285  }
286 
287  if (valid)
288  {
289  forAll(faceIDs, j)
290  {
291  if (nbrFacei == faceIDs[j])
292  {
293  valid = false;
294  break;
295  }
296  }
297  }
298 
299  // prevent addition of face if it is not on the same plane-ish
300  if (valid)
301  {
302  const vector& n1 = patch.faceNormals()[facei];
303  const vector& n2 = patch.faceNormals()[nbrFacei];
304 
305  scalar cosI = n1 & n2;
306 
307  if (cosI > cos(maxWalkAngle()))
308  {
309  faceIDs.append(nbrFacei);
310  }
311  }
312  }
313 }
314 
315 
316 Foam::scalar Foam::AMIMethod::maxWalkAngle() const
317 {
318  return degToRad(89);
319 }
320 
321 
322 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
323 
325 (
326  const primitivePatch& srcPatch,
327  const primitivePatch& tgtPatch,
328  const scalarField& srcMagSf,
329  const scalarField& tgtMagSf,
331  const bool reverseTarget,
332  const bool requireMatch
333 )
334 :
335  srcPatch_(srcPatch),
336  tgtPatch_(tgtPatch),
337  reverseTarget_(reverseTarget),
338  requireMatch_(requireMatch),
339  srcMagSf_(srcMagSf),
340  tgtMagSf_(tgtMagSf),
341  srcNonOverlap_(),
342  triMode_(triMode)
343 {}
344 
345 
346 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
347 
349 {}
350 
351 
352 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
353 
355 {
356  return true;
357 }
358 
359 
360 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
virtual ~AMIMethod()
Destructor.
Definition: AMIMethod.C:348
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Unit conversion functions.
Output to file stream.
Definition: OFstream.H:82
AMIMethod(const primitivePatch &srcPatch, const primitivePatch &tgtPatch, const scalarField &srcMagSf, const scalarField &tgtMagSf, const faceAreaIntersect::triangulationMode &triMode, const bool reverseTarget, const bool requireMatch)
Construct from components.
Definition: AMIMethod.C:325
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label count(const ListType &l, typename ListType::const_reference x)
Count the number of occurrences of a value in a list.
A bounding box defined in terms of the points at its extremities.
Definition: boundBox.H:58
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
label findTargetFace(const label srcFacei) const
Find face on target patch that overlaps source face.
Definition: AMIMethod.C:234
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
A list of faces which address into the list of points.
const Point & hitPoint() const
Return hit point.
const labelListList & faceFaces() const
Return face-face addressing.
dimensionedScalar cos(const dimensionedScalar &ds)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Encapsulation of data needed to search on PrimitivePatches.
virtual bool conformal() const
Flag to indicate that interpolation patches are conformal.
Definition: AMIMethod.C:354
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
void resetTree()
Reset the octree for the target patch face search.
Definition: AMIMethod.C:203
pointField points(const pointField &) const
Return the points corresponding to this face.
Definition: faceI.H:64
const Field< PointType > & faceNormals() const
Return face normals for patch.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void checkPatches() const
Check AMI patch coupling.
Definition: AMIMethod.C:42
void inflate(const scalar s)
Inflate box by factor*mag(span) in all dimensions.
Definition: boundBox.C:210
static const char nl
Definition: Ostream.H:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
static vector centre(const PointField &ps)
Return centre point given face points.
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
void setSize(const label)
Reset size of List.
Definition: List.C:281
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
void appendNbrFaces(const label facei, const primitivePatch &patch, const DynamicList< label > &visitedFaces, DynamicList< label > &faceIDs) const
Add faces neighbouring facei to the ID list.
Definition: AMIMethod.C:264
#define WarningInFunction
Report a warning using Foam::Warning.
bool initialise(labelListList &srcAddress, scalarListList &srcWeights, labelListList &tgtAddress, scalarListList &tgtWeights, label &srcFacei, label &tgtFacei)
Initialise and return true if all ok.
Definition: AMIMethod.C:79
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
virtual scalar maxWalkAngle() const
The maximum edge angle that the walk will cross.
Definition: AMIMethod.C:316
label index() const
Return index.
Namespace for OpenFOAM.
void writeIntersectionOBJ(const scalar area, const face &f1, const face &f2, const pointField &f1Points, const pointField &f2Points) const
Write triangle intersection to OBJ file.
Definition: AMIMethod.C:153