directAMI.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-2014 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 "directAMI.H"
27 
28 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
29 
30 template<class SourcePatch, class TargetPatch>
32 (
33  labelList& mapFlag,
34  labelList& srcTgtSeed,
35  DynamicList<label>& srcSeeds,
36  DynamicList<label>& nonOverlapFaces,
37  label& srcFaceI,
38  label& tgtFaceI
39 ) const
40 {
41  const labelList& srcNbr = this->srcPatch_.faceFaces()[srcFaceI];
42  const labelList& tgtNbr = this->tgtPatch_.faceFaces()[tgtFaceI];
43 
44  const pointField& srcPoints = this->srcPatch_.points();
45  const pointField& tgtPoints = this->tgtPatch_.points();
46 
47  const vectorField& srcCf = this->srcPatch_.faceCentres();
48 
49  forAll(srcNbr, i)
50  {
51  label srcI = srcNbr[i];
52 
53  if ((mapFlag[srcI] == 0) && (srcTgtSeed[srcI] == -1))
54  {
55  // first attempt: match by comparing face centres
56  const face& srcF = this->srcPatch_[srcI];
57  const point& srcC = srcCf[srcI];
58 
59  scalar tol = GREAT;
60  forAll(srcF, fpI)
61  {
62  const point& p = srcPoints[srcF[fpI]];
63  scalar d2 = magSqr(p - srcC);
64  if (d2 < tol)
65  {
66  tol = d2;
67  }
68  }
69  tol = max(SMALL, 0.0001*sqrt(tol));
70 
71  bool found = false;
72  forAll(tgtNbr, j)
73  {
74  label tgtI = tgtNbr[j];
75  const face& tgtF = this->tgtPatch_[tgtI];
76  const point tgtC = tgtF.centre(tgtPoints);
77 
78  if (mag(srcC - tgtC) < tol)
79  {
80  // new match - append to lists
81  found = true;
82 
83  srcTgtSeed[srcI] = tgtI;
84  srcSeeds.append(srcI);
85 
86  break;
87  }
88  }
89 
90  // second attempt: match by shooting a ray into the tgt face
91  if (!found)
92  {
93  const vector srcN = srcF.normal(srcPoints);
94 
95  forAll(tgtNbr, j)
96  {
97  label tgtI = tgtNbr[j];
98  const face& tgtF = this->tgtPatch_[tgtI];
99  pointHit ray = tgtF.ray(srcCf[srcI], srcN, tgtPoints);
100 
101  if (ray.hit())
102  {
103  // new match - append to lists
104  found = true;
105 
106  srcTgtSeed[srcI] = tgtI;
107  srcSeeds.append(srcI);
108 
109  break;
110  }
111  }
112  }
113 
114  // no match available for source face srcI
115  if (!found)
116  {
117  mapFlag[srcI] = -1;
118  nonOverlapFaces.append(srcI);
119 
120  if (debug)
121  {
122  Pout<< "source face not found: id=" << srcI
123  << " centre=" << srcCf[srcI]
124  << " normal=" << srcF.normal(srcPoints)
125  << " points=" << srcF.points(srcPoints)
126  << endl;
127 
128  Pout<< "target neighbours:" << nl;
129  forAll(tgtNbr, j)
130  {
131  label tgtI = tgtNbr[j];
132  const face& tgtF = this->tgtPatch_[tgtI];
133 
134  Pout<< "face id: " << tgtI
135  << " centre=" << tgtF.centre(tgtPoints)
136  << " normal=" << tgtF.normal(tgtPoints)
137  << " points=" << tgtF.points(tgtPoints)
138  << endl;
139  }
140  }
141  }
142  }
143  }
144 
145  if (srcSeeds.size())
146  {
147  srcFaceI = srcSeeds.remove();
148  tgtFaceI = srcTgtSeed[srcFaceI];
149  }
150  else
151  {
152  srcFaceI = -1;
153  tgtFaceI = -1;
154  }
155 }
156 
157 
158 template<class SourcePatch, class TargetPatch>
160 (
161  labelList& mapFlag,
162  DynamicList<label>& nonOverlapFaces,
163  label& srcFaceI,
164  label& tgtFaceI
165 ) const
166 {
167  forAll(mapFlag, faceI)
168  {
169  if (mapFlag[faceI] == 0)
170  {
171  tgtFaceI = this->findTargetFace(faceI);
172 
173  if (tgtFaceI < 0)
174  {
175  mapFlag[faceI] = -1;
176  nonOverlapFaces.append(faceI);
177  }
178  else
179  {
180  srcFaceI = faceI;
181  break;
182  }
183  }
184  }
185 }
186 
187 
188 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
189 
190 template<class SourcePatch, class TargetPatch>
192 (
193  const SourcePatch& srcPatch,
194  const TargetPatch& tgtPatch,
195  const scalarField& srcMagSf,
196  const scalarField& tgtMagSf,
198  const bool reverseTarget,
199  const bool requireMatch
200 )
201 :
203  (
204  srcPatch,
205  tgtPatch,
206  srcMagSf,
207  tgtMagSf,
208  triMode,
209  reverseTarget,
210  requireMatch
211  )
212 {}
213 
214 
215 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
216 
217 template<class SourcePatch, class TargetPatch>
219 {}
220 
221 
222 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
223 
224 template<class SourcePatch, class TargetPatch>
226 (
227  labelListList& srcAddress,
228  scalarListList& srcWeights,
229  labelListList& tgtAddress,
230  scalarListList& tgtWeights,
231  label srcFaceI,
232  label tgtFaceI
233 )
234 {
235  bool ok =
236  this->initialise
237  (
238  srcAddress,
239  srcWeights,
240  tgtAddress,
241  tgtWeights,
242  srcFaceI,
243  tgtFaceI
244  );
245 
246  if (!ok)
247  {
248  return;
249  }
250 
251 
252  // temporary storage for addressing and weights
253  List<DynamicList<label> > srcAddr(this->srcPatch_.size());
254  List<DynamicList<label> > tgtAddr(this->tgtPatch_.size());
255 
256 
257  // construct weights and addressing
258  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259 
260  // list of faces currently visited for srcFaceI to avoid multiple hits
261  DynamicList<label> srcSeeds(10);
262 
263  // list to keep track of tgt faces used to seed src faces
264  labelList srcTgtSeed(srcAddr.size(), -1);
265  srcTgtSeed[srcFaceI] = tgtFaceI;
266 
267  // list to keep track of whether src face can be mapped
268  // 1 = mapped, 0 = untested, -1 = cannot map
269  labelList mapFlag(srcAddr.size(), 0);
270 
271  label nTested = 0;
272  DynamicList<label> nonOverlapFaces;
273  do
274  {
275  srcAddr[srcFaceI].append(tgtFaceI);
276  tgtAddr[tgtFaceI].append(srcFaceI);
277 
278  mapFlag[srcFaceI] = 1;
279 
280  nTested++;
281 
282  // Do advancing front starting from srcFaceI, tgtFaceI
283  appendToDirectSeeds
284  (
285  mapFlag,
286  srcTgtSeed,
287  srcSeeds,
288  nonOverlapFaces,
289  srcFaceI,
290  tgtFaceI
291  );
292 
293  if (srcFaceI < 0 && nTested < this->srcPatch_.size())
294  {
295  restartAdvancingFront(mapFlag, nonOverlapFaces, srcFaceI, tgtFaceI);
296  }
297 
298  } while (srcFaceI >= 0);
299 
300  if (nonOverlapFaces.size() != 0)
301  {
302  Pout<< " AMI: " << nonOverlapFaces.size()
303  << " non-overlap faces identified"
304  << endl;
305 
306  this->srcNonOverlap_.transfer(nonOverlapFaces);
307  }
308 
309  // transfer data to persistent storage
310  forAll(srcAddr, i)
311  {
312  scalar magSf = this->srcMagSf_[i];
313 // srcWeights[i] = scalarList(srcAddr[i].size(), magSf);
314  srcWeights[i] = scalarList(1, magSf);
315  srcAddress[i].transfer(srcAddr[i]);
316  }
317  forAll(tgtAddr, i)
318  {
319  scalar magSf = this->tgtMagSf_[i];
320 // tgtWeights[i] = scalarList(tgtAddr[i].size(), magSf);
321  tgtWeights[i] = scalarList(1, magSf);
322  tgtAddress[i].transfer(tgtAddr[i]);
323  }
324 }
325 
326 
327 // ************************************************************************* //
dimensionedScalar sqrt(const dimensionedScalar &ds)
vector point
Point is a vector.
Definition: point.H:41
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
Base class for Arbitrary Mesh Interface (AMI) methods.
Definition: AMIMethod.H:55
virtual ~directAMI()
Destructor.
Definition: directAMI.C:218
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 size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
PointHit< point > pointHit
Definition: pointHit.H:41
virtual void calculate(labelListList &srcAddress, scalarListList &srcWeights, labelListList &tgtAddress, scalarListList &tgtWeights, label srcFaceI=-1, label tgtFaceI=-1)
Update addressing and weights.
Definition: directAMI.C:226
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
volVectorField vectorField(fieldObject, mesh)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAll(list, i)
Definition: UList.H:421
Direct mapped Arbitrary Mesh Interface (AMI) method.
Definition: directAMI.H:49
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt > > &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:106
List< label > labelList
A List of labels.
Definition: labelList.H:56
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53