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