FreeFOAM The Cross-Platform CFD Toolkit
PrimitivePatchInterpolation_.C
Go to the documentation of this file.
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 
27 #include <OpenFOAM/faceList.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  {
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  {
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>
157 :
158  patch_(p),
159  faceToPointWeightsPtr_(NULL),
160  faceToEdgeWeightsPtr_(NULL)
161 {}
162 
163 
164 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
165 
166 template<class Patch>
168 {
169  clearWeights();
170 }
171 
172 
173 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
174 
175 template<class Patch>
176 template<class Type>
178 (
179  const Field<Type>& ff
180 ) const
181 {
182  // Check size of the given field
183  if (ff.size() != patch_.size())
184  {
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>
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>
238 (
239  const Field<Type>& pf
240 ) const
241 {
242  if (pf.size() != patch_.nPoints())
243  {
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(),
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>
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>
298 (
299  const Field<Type>& pf
300 ) const
301 {
302  // Check size of the given field
303  if (pf.size() != patch_.size())
304  {
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 pointField& points = patch_.localPoints();
322  const faceList& faces = patch_.localFaces();
323  const edgeList& edges = patch_.edges();
324  const labelListList& edgeFaces = patch_.edgeFaces();
325 
326  const scalarList& weights = faceToEdgeWeights();
327 
328  for (label edgei = 0; edgei < patch_.nInternalEdges(); edgei++)
329  {
330  result[edgei] =
331  weights[edgei]*pf[edgeFaces[edgei][0]]
332  + (1.0 - weights[edgei])*pf[edgeFaces[edgei][1]];
333  }
334 
335  for (label edgei = patch_.nInternalEdges(); edgei < edges.size(); edgei++)
336  {
337  result[edgei] = pf[edgeFaces[edgei][0]];
338  }
339 
340  return tresult;
341 }
342 
343 
344 template<class Patch>
345 template<class Type>
347 (
348  const tmp<Field<Type> >& tpf
349 ) const
350 {
351  tmp<Field<Type> > tint = faceToEdgeInterpolate(tpf());
352  tpf.clear();
353  return tint;
354 }
355 
356 
357 template<class Patch>
359 {
360  clearWeights();
361 
362  return true;
363 }
364 
365 
366 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
367 
368 } // End namespace Foam
369 
370 // ************************ vim: set sw=4 sts=4 et: ************************ //