src/OpenFOAM/interpolations/patchToPatchInterpolation/CalcPatchToPatchWeights.C

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