src/OpenFOAM/interpolations/primitivePatchInterpolation/PrimitivePatchInterpolation.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 "PrimitivePatchInterpolation.H"
  28 #include "faceList.H"
  29 #include "demandDrivenData.H"
  30
  31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  32
  33 namespace Foam
  34 {
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 template<class Patch>
  39 const scalarListList&
  40 PrimitivePatchInterpolation<Patch>::faceToPointWeights() const
  41 {
  42     if (!faceToPointWeightsPtr_)
  43     {
  44         makeFaceToPointWeights();
  45     }
  46
  47     return *faceToPointWeightsPtr_;
  48 }
  49
  50
  51 template<class Patch>
  52 void PrimitivePatchInterpolation<Patch>::makeFaceToPointWeights() const
  53 {
  54     if (faceToPointWeightsPtr_)
  55     {
  56         FatalErrorIn
  57         (
  58             "PrimitivePatchInterpolation<Patch>::makeFaceToPointWeights() const"
  59         )   << "Face-to-edge weights already calculated"
  60             << abort(FatalError);
  61     }
  62
  63     const pointField& points = patch_.localPoints();
  64     const faceList& faces = patch_.localFaces();
  65
  66     faceToPointWeightsPtr_ = new scalarListList(points.size());
  67     scalarListList& weights = *faceToPointWeightsPtr_;
  68
  69     // get reference to addressing
  70     const labelListList& pointFaces = patch_.pointFaces();
  71
  72     forAll(pointFaces, pointi)
  73     {
  74         const labelList& curFaces = pointFaces[pointi];
  75
  76         scalarList& pw = weights[pointi];
  77         pw.setSize(curFaces.size());
  78
  79         scalar sumw = 0.0;
  80
  81         forAll(curFaces, facei)
  82         {
  83             pw[facei] =
  84                 1.0/mag(faces[curFaces[facei]].centre(points) - points[pointi]);
  85             sumw += pw[facei];
  86         }
  87
  88         forAll(curFaces, facei)
  89         {
  90             pw[facei] /= sumw;
  91         }
  92     }
  93 }
  94
  95
  96 template<class Patch>
  97 const scalarList&
  98 PrimitivePatchInterpolation<Patch>::faceToEdgeWeights() const
  99 {
 100     if (!faceToEdgeWeightsPtr_)
 101     {
 102         makeFaceToEdgeWeights();
 103     }
 104
 105     return *faceToEdgeWeightsPtr_;
 106 }
 107
 108
 109 template<class Patch>
 110 void PrimitivePatchInterpolation<Patch>::makeFaceToEdgeWeights() const
 111 {
 112     if (faceToEdgeWeightsPtr_)
 113     {
 114         FatalErrorIn
 115         (
 116             "PrimitivePatchInterpolation<Patch>::makeFaceToEdgeWeights() const"
 117         )   << "Face-to-edge weights already calculated"
 118             << abort(FatalError);
 119     }
 120
 121     const pointField& points = patch_.localPoints();
 122     const faceList& faces = patch_.localFaces();
 123     const edgeList& edges = patch_.edges();
 124     const labelListList& edgeFaces = patch_.edgeFaces();
 125
 126     faceToEdgeWeightsPtr_ = new scalarList(patch_.nInternalEdges());
 127     scalarList& weights = *faceToEdgeWeightsPtr_;
 128
 129     for (label edgei = 0; edgei < weights.size(); edgei++)
 130     {
 131         vector P = faces[edgeFaces[edgei][0]].centre(points);
 132         vector N = faces[edgeFaces[edgei][1]].centre(points);
 133         vector S = points[edges[edgei].start()];
 134         vector e = edges[edgei].vec(points);
 135
 136         scalar alpha =
 137             -(((N - P)^(S - P))&((N - P)^e))/(((N - P)^e )&((N - P)^e));
 138
 139         vector E = S + alpha*e;
 140
 141         weights[edgei] = mag(N - E)/(mag(N - E) + mag(E - P));
 142     }
 143 }
 144
 145
 146 template<class Patch>
 147 void PrimitivePatchInterpolation<Patch>::clearWeights()
 148 {
 149     deleteDemandDrivenData(faceToPointWeightsPtr_);
 150     deleteDemandDrivenData(faceToEdgeWeightsPtr_);
 151 }
 152
 153
 154 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 155
 156 template<class Patch>
 157 PrimitivePatchInterpolation<Patch>::PrimitivePatchInterpolation(const Patch& p)
 158 :
 159     patch_(p),
 160     faceToPointWeightsPtr_(NULL),
 161     faceToEdgeWeightsPtr_(NULL)
 162 {}
 163
 164
 165 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
 166
 167 template<class Patch>
 168 PrimitivePatchInterpolation<Patch>::~PrimitivePatchInterpolation()
 169 {
 170     clearWeights();
 171 }
 172
 173
 174 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 175
 176 template<class Patch>
 177 template<class Type>
 178 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::faceToPointInterpolate
 179 (
 180     const Field<Type>& ff
 181 ) const
 182 {
 183     // Check size of the given field
 184     if (ff.size() != patch_.size())
 185     {
 186         FatalErrorIn
 187         (
 188             "tmp<Field<Type> > PrimitivePatchInterpolation::"
 189             "faceToPointInterpolate(const Field<Type> ff)"
 190         )   << "given field does not correspond to patch. Patch size: "
 191             << patch_.size() << " field size: " << ff.size()
 192             << abort(FatalError);
 193     }
 194
 195     tmp<Field<Type> > tresult
 196     (
 197         new Field<Type>
 198         (
 199             patch_.nPoints(), pTraits<Type>::zero
 200         )
 201     );
 202
 203     Field<Type>& result = tresult();
 204
 205     const labelListList& pointFaces = patch_.pointFaces();
 206     const scalarListList& weights = faceToPointWeights();
 207
 208     forAll(pointFaces, pointi)
 209     {
 210         const labelList& curFaces = pointFaces[pointi];
 211         const scalarList& w = weights[pointi];
 212
 213         forAll(curFaces, facei)
 214         {
 215             result[pointi] += w[facei]*ff[curFaces[facei]];
 216         }
 217     }
 218
 219     return tresult;
 220 }
 221
 222
 223 template<class Patch>
 224 template<class Type>
 225 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::faceToPointInterpolate
 226 (
 227     const tmp<Field<Type> >& tff
 228 ) const
 229 {
 230     tmp<Field<Type> > tint = faceToPointInterpolate(tff());
 231     tff.clear();
 232     return tint;
 233 }
 234
 235
 236 template<class Patch>
 237 template<class Type>
 238 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::pointToFaceInterpolate
 239 (
 240     const Field<Type>& pf
 241 ) const
 242 {
 243     if (pf.size() != patch_.nPoints())
 244     {
 245         FatalErrorIn
 246         (
 247             "tmp<Field<Type> > PrimitivePatchInterpolation::"
 248             "pointToFaceInterpolate(const Field<Type> pf)"
 249         )   << "given field does not correspond to patch. Patch size: "
 250             << patch_.nPoints() << " field size: " << pf.size()
 251             << abort(FatalError);
 252     }
 253
 254     tmp<Field<Type> > tresult
 255     (
 256         new Field<Type>
 257         (
 258             patch_.size(),
 259             pTraits<Type>::zero
 260         )
 261     );
 262
 263     Field<Type>& result = tresult();
 264
 265     const faceList& localFaces = patch_.localFaces();
 266
 267     forAll(result, facei)
 268     {
 269         const labelList& curPoints = localFaces[facei];
 270
 271         forAll(curPoints, pointi)
 272         {
 273             result[facei] += pf[curPoints[pointi]];
 274         }
 275
 276         result[facei] /= curPoints.size();
 277     }
 278
 279     return tresult;
 280 }
 281
 282
 283 template<class Patch>
 284 template<class Type>
 285 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::pointToFaceInterpolate
 286 (
 287     const tmp<Field<Type> >& tpf
 288 ) const
 289 {
 290     tmp<Field<Type> > tint = pointToFaceInterpolate(tpf());
 291     tpf.clear();
 292     return tint;
 293 }
 294
 295
 296 template<class Patch>
 297 template<class Type>
 298 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::faceToEdgeInterpolate
 299 (
 300     const Field<Type>& pf
 301 ) const
 302 {
 303     // Check size of the given field
 304     if (pf.size() != patch_.size())
 305     {
 306         FatalErrorIn
 307         (
 308             "tmp<Field<Type> > PrimitivePatchInterpolation::"
 309             "faceToEdgeInterpolate(const Field<Type> ff)"
 310         )   << "given field does not correspond to patch. Patch size: "
 311             << patch_.size() << " field size: " << pf.size()
 312             << abort(FatalError);
 313     }
 314
 315     tmp<Field<Type> > tresult
 316     (
 317         new Field<Type>(patch_.nEdges(), pTraits<Type>::zero)
 318     );
 319
 320     Field<Type>& result = tresult();
 321
 322     const pointField& points = patch_.localPoints();
 323     const faceList& faces = patch_.localFaces();
 324     const edgeList& edges = patch_.edges();
 325     const labelListList& edgeFaces = patch_.edgeFaces();
 326
 327     const scalarList& weights = faceToEdgeWeights();
 328
 329     for (label edgei = 0; edgei < patch_.nInternalEdges(); edgei++)
 330     {
 331         result[edgei] =
 332             weights[edgei]*pf[edgeFaces[edgei][0]]
 333           + (1.0 - weights[edgei])*pf[edgeFaces[edgei][1]];
 334     }
 335
 336     for (label edgei = patch_.nInternalEdges(); edgei < edges.size(); edgei++)
 337     {
 338         result[edgei] = pf[edgeFaces[edgei][0]];
 339     }
 340
 341     return tresult;
 342 }
 343
 344
 345 template<class Patch>
 346 template<class Type>
 347 tmp<Field<Type> > PrimitivePatchInterpolation<Patch>::faceToEdgeInterpolate
 348 (
 349     const tmp<Field<Type> >& tpf
 350 ) const
 351 {
 352     tmp<Field<Type> > tint = faceToEdgeInterpolate(tpf());
 353     tpf.clear();
 354     return tint;
 355 }
 356
 357
 358 template<class Patch>
 359 bool PrimitivePatchInterpolation<Patch>::movePoints()
 360 {
 361     clearWeights();
 362
 363     return true;
 364 }
 365
 366
 367 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 368
 369 } // End namespace Foam
 370
 371 // ************************************************************************* //