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) 2011 OpenFOAM Foundation
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 "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  {
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 List<typename Patch::FaceType>& 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 List<typename Patch::FaceType>& 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  forAll(weights, 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 List<typename Patch::FaceType>& 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 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>
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>
357 {
358  clearWeights();
359 
360  return true;
361 }
362 
363 
364 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
365 
366 } // End namespace Foam
367 
368 // ************************************************************************* //
List< scalarList > scalarListList
Definition: scalarList.H:51
dimensioned< scalar > mag(const dimensioned< Type > &)
void deleteDemandDrivenData(DataPtr &dataPtr)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Interpolation class within a primitive patch. Allows interpolation from points to faces and vice vers...
tmp< Field< Type > > faceToEdgeInterpolate(const Field< Type > &ff) const
Interpolate from faces to edges.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
volScalarField & p
Definition: createFields.H:51
#define forAll(list, i)
Definition: UList.H:421
Template functions to aid in the implementation of demand driven data.
tmp< Field< Type > > faceToPointInterpolate(const Field< Type > &ff) const
Interpolate from faces to points.
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
Pre-declare SubField and related Field type.
Definition: Field.H:57
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt > > &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:106
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
List< edge > edgeList
Definition: edgeList.H:38
Traits class for primitives.
Definition: pTraits.H:50
error FatalError
tmp< Field< Type > > pointToFaceInterpolate(const Field< Type > &pf) const
Interpolate from points to faces.
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool movePoints()
Do what is neccessary if the mesh has moved.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
A class for managing temporary objects.
Definition: PtrList.H:118