1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "PatchToPatchInterpolation_.H"
28 #include <OpenFOAM/objectHit.H>
29 #include <OpenFOAM/pointHit.H>
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
38 template<class FromPatch, class ToPatch>
39 scalar PatchToPatchInterpolation<FromPatch, ToPatch>::projectionTol_ = 0.05;
42 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 template<class FromPatch, class ToPatch>
45 void PatchToPatchInterpolation<FromPatch, ToPatch>::calcPointAddressing() const
47 // Calculate pointWeights
49 pointWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.nPoints());
50 FieldField<Field, scalar>& pointWeights = *pointWeightsPtr_;
52 pointDistancePtr_ = new scalarField(toPatch_.nPoints(), GREAT);
53 scalarField& pointDistance = *pointDistancePtr_;
55 const pointField& fromPatchPoints = fromPatch_.localPoints();
56 const List<typename FromPatch::FaceType>& fromPatchFaces =
57 fromPatch_.localFaces();
59 const pointField& toPatchPoints = toPatch_.localPoints();
60 const vectorField& projectionDirection = toPatch_.pointNormals();
61 const edgeList& toPatchEdges = toPatch_.edges();
62 const labelListList& toPatchPointEdges = toPatch_.pointEdges();
66 Info << "projecting points" << endl;
69 List<objectHit> proj =
70 toPatch_.projectPoints(fromPatch_, projectionDirection, alg_, dir_);
72 pointAddressingPtr_ = new labelList(proj.size(), -1);
73 labelList& pointAddressing = *pointAddressingPtr_;
75 bool doWeights = false;
77 forAll (pointAddressing, pointI)
81 const typename FromPatch::FaceType& hitFace =
82 fromPatchFaces[proj[pointI].hitObject()];
84 point hitPoint = point::zero;
86 if (proj[pointI].hit())
91 pointAddressing[pointI] = proj[pointI].hitObject();
96 toPatchPoints[pointI],
97 projectionDirection[pointI],
103 // Grab distance to target
104 if (dir_ == intersection::CONTACT_SPHERE)
106 pointDistance[pointI] =
107 hitFace.contactSphereDiameter
109 toPatchPoints[pointI],
110 projectionDirection[pointI],
116 pointDistance[pointI] = curHit.distance();
120 hitPoint = curHit.hitPoint();
122 else if (projectionTol_ > SMALL)
124 // Check for a near miss
128 toPatchPoints[pointI],
129 projectionDirection[pointI],
138 toPatchPoints[pointI]
139 + projectionDirection[pointI]*ph.distance()
143 // Calculate the local tolerance
144 scalar minEdgeLength = GREAT;
146 // Do shortest edge of hit object
147 edgeList hitFaceEdges =
148 fromPatchFaces[proj[pointI].hitObject()].edges();
150 forAll (hitFaceEdges, edgeI)
156 hitFaceEdges[edgeI].mag(fromPatchPoints)
160 const labelList& curEdges = toPatchPointEdges[pointI];
162 forAll (curEdges, edgeI)
168 toPatchEdges[curEdges[edgeI]].mag(toPatchPoints)
172 if (dist < minEdgeLength*projectionTol_)
174 // This point is being corrected
177 pointAddressing[pointI] = proj[pointI].hitObject();
179 // Grab nearest point on face as hit point
180 hitPoint = ph.missPoint();
182 // Grab distance to target
183 if (dir_ == intersection::CONTACT_SPHERE)
185 pointDistance[pointI] =
186 hitFace.contactSphereDiameter
188 toPatchPoints[pointI],
189 projectionDirection[pointI],
195 pointDistance[pointI] =
197 projectionDirection[pointI]
198 /mag(projectionDirection[pointI])
200 & (hitPoint - toPatchPoints[pointI]);
207 // Set interpolation pointWeights
208 pointWeights.set(pointI, new scalarField(hitFace.size()));
210 pointField hitFacePoints = hitFace.points(fromPatchPoints);
212 forAll (hitFacePoints, masterPointI)
214 pointWeights[pointI][masterPointI] =
219 hitFacePoints[masterPointI]
226 pointWeights[pointI] /= sum(pointWeights[pointI]);
230 pointWeights.set(pointI, new scalarField(0));
236 template<class FromPatch, class ToPatch>
237 void PatchToPatchInterpolation<FromPatch, ToPatch>::calcFaceAddressing() const
239 faceWeightsPtr_ = new FieldField<Field, scalar>(toPatch_.size());
240 FieldField<Field, scalar>& faceWeights = *faceWeightsPtr_;
242 faceDistancePtr_ = new scalarField(toPatch_.size(), GREAT);
243 scalarField& faceDistance = *faceDistancePtr_;
247 Info << "projecting face centres" << endl;
250 const pointField& fromPatchPoints = fromPatch_.points();
251 const typename FromPatch::FaceListType& fromPatchFaces = fromPatch_;
252 const labelListList& fromPatchFaceFaces = fromPatch_.faceFaces();
254 vectorField fromPatchFaceCentres(fromPatchFaces.size());
256 forAll (fromPatchFaceCentres, faceI)
258 fromPatchFaceCentres[faceI] =
259 fromPatchFaces[faceI].centre(fromPatchPoints);
262 const pointField& toPatchPoints = toPatch_.points();
263 const typename ToPatch::FaceListType& toPatchFaces = toPatch_;
265 const vectorField& projectionDirection = toPatch_.faceNormals();
267 List<objectHit> proj =
268 toPatch_.projectFaceCentres
276 faceAddressingPtr_ = new labelList(proj.size(), -1);
277 labelList& faceAddressing = *faceAddressingPtr_;
279 forAll (faceAddressing, faceI)
281 if (proj[faceI].hit())
284 faceAddressing[faceI] = proj[faceI].hitObject();
286 const typename FromPatch::FaceType& hitFace =
287 fromPatchFaces[faceAddressing[faceI]];
292 toPatchFaces[faceI].centre(toPatchPoints),
293 projectionDirection[faceI],
299 // grab distance to target
300 faceDistance[faceI] = curHit.distance();
302 // grab face centre of the hit face
303 const point& hitFaceCentre =
304 fromPatchFaceCentres[faceAddressing[faceI]];
306 // grab neighbours of hit face
307 const labelList& neighbours =
308 fromPatchFaceFaces[faceAddressing[faceI]];
310 scalar m = mag(curHit.hitPoint() - hitFaceCentre);
314 m < directHitTol // Direct hit
315 || neighbours.empty()
318 faceWeights.set(faceI, new scalarField(1));
319 faceWeights[faceI][0] = 1.0;
323 // set interpolation faceWeights
325 // The first coefficient corresponds to the centre face.
326 // The rest is ordered in the same way as the faceFaces list.
327 faceWeights.set(faceI, new scalarField(neighbours.size() + 1));
329 faceWeights[faceI][0] = 1.0/m;
331 forAll (neighbours, nI)
333 faceWeights[faceI][nI + 1] =
338 fromPatchFaceCentres[neighbours[nI]]
346 faceWeights[faceI] /= sum(faceWeights[faceI]);
350 faceWeights.set(faceI, new scalarField(0));
356 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
358 } // End namespace Foam
360 // ************************ vim: set sw=4 sts=4 et: ************************ //