src/foam/interpolations/patchToPatchInterpolation/CalcPatchToPatchWeights.C

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