PrimitivePatchInterpolation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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  << "Face-to-edge weights already calculated"
57  << abort(FatalError);
58  }
59 
60  const pointField& points = patch_.localPoints();
61  const List<typename Patch::FaceType>& faces = patch_.localFaces();
62 
63  faceToPointWeightsPtr_ = new scalarListList(points.size());
64  scalarListList& weights = *faceToPointWeightsPtr_;
65 
66  // get reference to addressing
67  const labelListList& pointFaces = patch_.pointFaces();
68 
69  forAll(pointFaces, pointi)
70  {
71  const labelList& curFaces = pointFaces[pointi];
72 
73  scalarList& pw = weights[pointi];
74  pw.setSize(curFaces.size());
75 
76  scalar sumw = 0.0;
77 
78  forAll(curFaces, facei)
79  {
80  pw[facei] =
81  1.0/mag(faces[curFaces[facei]].centre(points) - points[pointi]);
82  sumw += pw[facei];
83  }
84 
85  forAll(curFaces, facei)
86  {
87  pw[facei] /= sumw;
88  }
89  }
90 }
91 
92 
93 template<class Patch>
94 const scalarList&
95 PrimitivePatchInterpolation<Patch>::faceToEdgeWeights() const
96 {
97  if (!faceToEdgeWeightsPtr_)
98  {
99  makeFaceToEdgeWeights();
100  }
101 
102  return *faceToEdgeWeightsPtr_;
103 }
104 
105 
106 template<class Patch>
107 void PrimitivePatchInterpolation<Patch>::makeFaceToEdgeWeights() const
108 {
109  if (faceToEdgeWeightsPtr_)
110  {
112  << "Face-to-edge weights already calculated"
113  << abort(FatalError);
114  }
115 
116  const pointField& points = patch_.localPoints();
117  const List<typename Patch::FaceType>& faces = patch_.localFaces();
118  const edgeList& edges = patch_.edges();
119  const labelListList& edgeFaces = patch_.edgeFaces();
120 
121  faceToEdgeWeightsPtr_ = new scalarList(patch_.nInternalEdges());
122  scalarList& weights = *faceToEdgeWeightsPtr_;
123 
124  forAll(weights, edgei)
125  {
126  vector P = faces[edgeFaces[edgei][0]].centre(points);
127  vector N = faces[edgeFaces[edgei][1]].centre(points);
128  vector S = points[edges[edgei].start()];
129  vector e = edges[edgei].vec(points);
130 
131  scalar alpha =
132  -(((N - P)^(S - P))&((N - P)^e))/(((N - P)^e )&((N - P)^e));
133 
134  vector E = S + alpha*e;
135 
136  weights[edgei] = mag(N - E)/(mag(N - E) + mag(E - P));
137  }
138 }
139 
140 
141 template<class Patch>
142 void PrimitivePatchInterpolation<Patch>::clearWeights()
143 {
144  deleteDemandDrivenData(faceToPointWeightsPtr_);
145  deleteDemandDrivenData(faceToEdgeWeightsPtr_);
146 }
147 
148 
149 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
150 
151 template<class Patch>
153 :
154  patch_(p),
155  faceToPointWeightsPtr_(nullptr),
156  faceToEdgeWeightsPtr_(nullptr)
157 {}
158 
159 
160 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
161 
162 template<class Patch>
164 {
165  clearWeights();
166 }
167 
168 
169 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
170 
171 template<class Patch>
172 template<class Type>
174 (
175  const Field<Type>& ff
176 ) const
177 {
178  // Check size of the given field
179  if (ff.size() != patch_.size())
180  {
182  << "given field does not correspond to patch. Patch size: "
183  << patch_.size() << " field size: " << ff.size()
184  << abort(FatalError);
185  }
186 
187  tmp<Field<Type>> tresult
188  (
189  new Field<Type>
190  (
191  patch_.nPoints(), Zero
192  )
193  );
194 
195  Field<Type>& result = tresult.ref();
196 
197  const labelListList& pointFaces = patch_.pointFaces();
198  const scalarListList& weights = faceToPointWeights();
199 
200  forAll(pointFaces, pointi)
201  {
202  const labelList& curFaces = pointFaces[pointi];
203  const scalarList& w = weights[pointi];
204 
205  forAll(curFaces, facei)
206  {
207  result[pointi] += w[facei]*ff[curFaces[facei]];
208  }
209  }
210 
211  return tresult;
212 }
213 
214 
215 template<class Patch>
216 template<class Type>
218 (
219  const tmp<Field<Type>>& tff
220 ) const
221 {
222  tmp<Field<Type>> tint = faceToPointInterpolate(tff());
223  tff.clear();
224  return tint;
225 }
226 
227 
228 template<class Patch>
229 template<class Type>
231 (
232  const Field<Type>& pf
233 ) const
234 {
235  if (pf.size() != patch_.nPoints())
236  {
238  << "given field does not correspond to patch. Patch size: "
239  << patch_.nPoints() << " field size: " << pf.size()
240  << abort(FatalError);
241  }
242 
243  tmp<Field<Type>> tresult
244  (
245  new Field<Type>
246  (
247  patch_.size(),
248  Zero
249  )
250  );
251 
252  Field<Type>& result = tresult.ref();
253 
254  const List<typename Patch::FaceType>& localFaces = patch_.localFaces();
255 
256  forAll(result, facei)
257  {
258  const labelList& curPoints = localFaces[facei];
259 
260  forAll(curPoints, pointi)
261  {
262  result[facei] += pf[curPoints[pointi]];
263  }
264 
265  result[facei] /= curPoints.size();
266  }
267 
268  return tresult;
269 }
270 
271 
272 template<class Patch>
273 template<class Type>
275 (
276  const tmp<Field<Type>>& tpf
277 ) const
278 {
279  tmp<Field<Type>> tint = pointToFaceInterpolate(tpf());
280  tpf.clear();
281  return tint;
282 }
283 
284 
285 template<class Patch>
286 template<class Type>
288 (
289  const Field<Type>& pf
290 ) const
291 {
292  // Check size of the given field
293  if (pf.size() != patch_.size())
294  {
296  << "given field does not correspond to patch. Patch size: "
297  << patch_.size() << " field size: " << pf.size()
298  << abort(FatalError);
299  }
300 
301  tmp<Field<Type>> tresult
302  (
303  new Field<Type>(patch_.nEdges(), Zero)
304  );
305 
306  Field<Type>& result = tresult.ref();
307 
308  const edgeList& edges = patch_.edges();
309  const labelListList& edgeFaces = patch_.edgeFaces();
310 
311  const scalarList& weights = faceToEdgeWeights();
312 
313  for (label edgei = 0; edgei < patch_.nInternalEdges(); edgei++)
314  {
315  result[edgei] =
316  weights[edgei]*pf[edgeFaces[edgei][0]]
317  + (1.0 - weights[edgei])*pf[edgeFaces[edgei][1]];
318  }
319 
320  for (label edgei = patch_.nInternalEdges(); edgei < edges.size(); edgei++)
321  {
322  result[edgei] = pf[edgeFaces[edgei][0]];
323  }
324 
325  return tresult;
326 }
327 
328 
329 template<class Patch>
330 template<class Type>
332 (
333  const tmp<Field<Type>>& tpf
334 ) const
335 {
336  tmp<Field<Type>> tint = faceToEdgeInterpolate(tpf());
337  tpf.clear();
338  return tint;
339 }
340 
341 
342 template<class Patch>
344 {
345  clearWeights();
346 
347  return true;
348 }
349 
350 
351 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
352 
353 } // End namespace Foam
354 
355 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Interpolation class within a primitive patch. Allows interpolation from points to faces and vice vers...
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
T & ref() const
Return non-const reference or generate a fatal error.
Definition: tmpI.H:181
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
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:116
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< edge > edgeList
Definition: edgeList.H:38
PrimitivePatchInterpolation(const Patch &p)
Construct from PrimitivePatch.
Pre-declare SubField and related Field type.
Definition: Field.H:56
List< scalarList > scalarListList
Definition: scalarList.H:51
List< scalar > scalarList
A List of scalars.
Definition: scalarList.H:50
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void setSize(const label)
Reset size of List.
Definition: List.C:281
Template functions to aid in the implementation of demand driven data.
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
volScalarField & p
A class for managing temporary objects.
Definition: PtrList.H:53
void deleteDemandDrivenData(DataPtr &dataPtr)
Namespace for OpenFOAM.