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-2010 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
  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 "PatchToPatchInterpolation.H"
  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 = point::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 // ************************************************************************* //