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