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