CalcPatchToPatchWeights.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) 2011-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 
27 #include "objectHit.H"
28 #include "pointHit.H"
29 
30 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
37 template<class FromPatch, class ToPatch>
38 scalar PatchToPatchInterpolation<FromPatch, ToPatch>::projectionTol_ = 0.05;
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 template<class FromPatch, class ToPatch>
44 void PatchToPatchInterpolation<FromPatch, ToPatch>::calcPointAddressing() const
45 {
46  // Calculate pointWeights
47 
48  pointWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.nPoints());
49  FieldField<Field, scalar>& pointWeights = *pointWeightsPtr_;
50 
51  pointDistancePtr_ = new scalarField(toPatch_.nPoints(), great);
52  scalarField& pointDistance = *pointDistancePtr_;
53 
54  const pointField& fromPatchPoints = fromPatch_.localPoints();
55  const List<typename FromPatch::FaceType>& fromPatchFaces =
56  fromPatch_.localFaces();
57 
58  const pointField& toPatchPoints = toPatch_.localPoints();
59  const vectorField& projectionDirection = toPatch_.pointNormals();
60  const edgeList& toPatchEdges = toPatch_.edges();
61  const labelListList& toPatchPointEdges = toPatch_.pointEdges();
62 
63  if (debug)
64  {
65  Info<< "projecting points" << endl;
66  }
67 
68  List<objectHit> proj =
69  toPatch_.projectPoints(fromPatch_, projectionDirection, alg_, dir_);
70 
71  pointAddressingPtr_ = new labelList(proj.size(), -1);
72  labelList& pointAddressing = *pointAddressingPtr_;
73 
74  bool doWeights = false;
75 
76  forAll(pointAddressing, pointi)
77  {
78  doWeights = false;
79 
80  const typename FromPatch::FaceType& hitFace =
81  fromPatchFaces[proj[pointi].hitObject()];
82 
83  point hitPoint = Zero;
84 
85  if (proj[pointi].hit())
86  {
87  // A hit exists
88  doWeights = true;
89 
90  pointAddressing[pointi] = proj[pointi].hitObject();
91 
92  pointHit curHit =
93  hitFace.ray
94  (
95  toPatchPoints[pointi],
96  projectionDirection[pointi],
97  fromPatchPoints,
98  alg_,
99  dir_
100  );
101 
102  // Grab distance to target
103  if (dir_ == intersection::CONTACT_SPHERE)
104  {
105  pointDistance[pointi] =
106  hitFace.contactSphereDiameter
107  (
108  toPatchPoints[pointi],
109  projectionDirection[pointi],
110  fromPatchPoints
111  );
112  }
113  else
114  {
115  pointDistance[pointi] = curHit.distance();
116  }
117 
118  // Grab hit point
119  hitPoint = curHit.hitPoint();
120  }
121  else if (projectionTol_ > small)
122  {
123  // Check for a near miss
124  pointHit ph =
125  hitFace.ray
126  (
127  toPatchPoints[pointi],
128  projectionDirection[pointi],
129  fromPatchPoints,
130  alg_,
131  dir_
132  );
133 
134  scalar dist =
135  Foam::mag
136  (
137  toPatchPoints[pointi]
138  + projectionDirection[pointi]*ph.distance()
139  - ph.missPoint()
140  );
141 
142  // Calculate the local tolerance
143  scalar minEdgeLength = great;
144 
145  // Do shortest edge of hit object
146  edgeList hitFaceEdges =
147  fromPatchFaces[proj[pointi].hitObject()].edges();
148 
149  forAll(hitFaceEdges, edgeI)
150  {
151  minEdgeLength =
152  min
153  (
154  minEdgeLength,
155  hitFaceEdges[edgeI].mag(fromPatchPoints)
156  );
157  }
158 
159  const labelList& curEdges = toPatchPointEdges[pointi];
160 
161  forAll(curEdges, edgeI)
162  {
163  minEdgeLength =
164  min
165  (
166  minEdgeLength,
167  toPatchEdges[curEdges[edgeI]].mag(toPatchPoints)
168  );
169  }
170 
171  if (dist < minEdgeLength*projectionTol_)
172  {
173  // This point is being corrected
174  doWeights = true;
175 
176  pointAddressing[pointi] = proj[pointi].hitObject();
177 
178  // Grab nearest point on face as hit point
179  hitPoint = ph.missPoint();
180 
181  // Grab distance to target
182  if (dir_ == intersection::CONTACT_SPHERE)
183  {
184  pointDistance[pointi] =
185  hitFace.contactSphereDiameter
186  (
187  toPatchPoints[pointi],
188  projectionDirection[pointi],
189  fromPatchPoints
190  );
191  }
192  else
193  {
194  pointDistance[pointi] =
195  (
196  projectionDirection[pointi]
197  /mag(projectionDirection[pointi])
198  )
199  & (hitPoint - toPatchPoints[pointi]);
200  }
201  }
202  }
203 
204  if (doWeights)
205  {
206  // Set interpolation pointWeights
207  pointWeights.set(pointi, new scalarField(hitFace.size()));
208 
209  pointField hitFacePoints = hitFace.points(fromPatchPoints);
210 
211  forAll(hitFacePoints, masterPointi)
212  {
213  pointWeights[pointi][masterPointi] =
214  1.0/
215  (
216  mag
217  (
218  hitFacePoints[masterPointi]
219  - hitPoint
220  )
221  + vSmall
222  );
223  }
224 
225  pointWeights[pointi] /= sum(pointWeights[pointi]);
226  }
227  else
228  {
229  pointWeights.set(pointi, new scalarField(0));
230  }
231  }
232 }
233 
234 
235 template<class FromPatch, class ToPatch>
236 void PatchToPatchInterpolation<FromPatch, ToPatch>::calcFaceAddressing() const
237 {
238  faceWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.size());
239  FieldField<Field, scalar>& faceWeights = *faceWeightsPtr_;
240 
241  faceDistancePtr_ = new scalarField(toPatch_.size(), great);
242  scalarField& faceDistance = *faceDistancePtr_;
243 
244  if (debug)
245  {
246  Info<< "projecting face centres" << endl;
247  }
248 
249  const pointField& fromPatchPoints = fromPatch_.points();
250  const typename FromPatch::FaceListType& fromPatchFaces = fromPatch_;
251  const labelListList& fromPatchFaceFaces = fromPatch_.faceFaces();
252 
253  vectorField fromPatchFaceCentres(fromPatchFaces.size());
254 
255  forAll(fromPatchFaceCentres, facei)
256  {
257  fromPatchFaceCentres[facei] =
258  fromPatchFaces[facei].centre(fromPatchPoints);
259  }
260 
261  const pointField& toPatchPoints = toPatch_.points();
262  const typename ToPatch::FaceListType& toPatchFaces = toPatch_;
263 
264  const vectorField& projectionDirection = toPatch_.faceNormals();
265 
266  List<objectHit> proj =
267  toPatch_.projectFaceCentres
268  (
269  fromPatch_,
270  projectionDirection,
271  alg_,
272  dir_
273  );
274 
275  faceAddressingPtr_ = new labelList(proj.size(), -1);
276  labelList& faceAddressing = *faceAddressingPtr_;
277 
278  forAll(faceAddressing, facei)
279  {
280  if (proj[facei].hit())
281  {
282  // A hit exists
283  faceAddressing[facei] = proj[facei].hitObject();
284 
285  const typename FromPatch::FaceType& hitFace =
286  fromPatchFaces[faceAddressing[facei]];
287 
288  pointHit curHit =
289  hitFace.ray
290  (
291  toPatchFaces[facei].centre(toPatchPoints),
292  projectionDirection[facei],
293  fromPatchPoints,
294  alg_,
295  dir_
296  );
297 
298  // grab distance to target
299  faceDistance[facei] = curHit.distance();
300 
301  // grab face centre of the hit face
302  const point& hitFaceCentre =
303  fromPatchFaceCentres[faceAddressing[facei]];
304 
305  // grab neighbours of hit face
306  const labelList& neighbours =
307  fromPatchFaceFaces[faceAddressing[facei]];
308 
309  scalar m = mag(curHit.hitPoint() - hitFaceCentre);
310 
311  if
312  (
313  m < directHitTol // Direct hit
314  || neighbours.empty()
315  )
316  {
317  faceWeights.set(facei, new scalarField(1));
318  faceWeights[facei][0] = 1.0;
319  }
320  else
321  {
322  // set interpolation faceWeights
323 
324  // The first coefficient corresponds to the centre face.
325  // The rest is ordered in the same way as the faceFaces list.
326  faceWeights.set(facei, new scalarField(neighbours.size() + 1));
327 
328  faceWeights[facei][0] = 1.0/m;
329 
330  forAll(neighbours, nI)
331  {
332  faceWeights[facei][nI + 1] =
333  1.0/
334  (
335  mag
336  (
337  fromPatchFaceCentres[neighbours[nI]]
338  - curHit.hitPoint()
339  )
340  + vSmall
341  );
342  }
343  }
344 
345  faceWeights[facei] /= sum(faceWeights[facei]);
346  }
347  else
348  {
349  faceWeights.set(facei, new scalarField(0));
350  }
351  }
352 }
353 
354 
355 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
356 
357 } // End namespace Foam
358 
359 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< edge > edgeList
Definition: edgeList.H:38
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
vector point
Point is a vector.
Definition: point.H:41
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
PointHit< point > pointHit
Definition: pointHit.H:41
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
Namespace for OpenFOAM.