FreeFOAM The Cross-Platform CFD Toolkit
pointPatchInterpolation.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 <finiteVolume/fvMesh.H>
28 #include <finiteVolume/volFields.H>
29 #include <OpenFOAM/pointFields.H>
34 
35 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
36 
37 namespace Foam
38 {
39 
40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
41 
42 defineTypeNameAndDebug(pointPatchInterpolation, 0);
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 void pointPatchInterpolation::makePatchPatchAddressing()
48 {
49  if (debug)
50  {
51  Info<< "pointPatchInterpolation::makePatchPatchAddressing() : "
52  << "constructing boundary addressing"
53  << endl;
54  }
55 
56  const fvBoundaryMesh& bm = fvMesh_.boundary();
57  const pointBoundaryMesh& pbm = pointMesh::New(fvMesh_).boundary();
58 
59  // first count the total number of patch-patch points
60 
61  label nPatchPatchPoints = 0;
62 
63  forAll(bm, patchi)
64  {
65  if(!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
66  {
67  nPatchPatchPoints += bm[patchi].patch().boundaryPoints().size();
68  }
69  }
70 
71 
72  // Go through all patches and mark up the external edge points
73  Map<label> patchPatchPointSet(2*nPatchPatchPoints);
74 
75  patchPatchPoints_.setSize(nPatchPatchPoints);
76 
77  List<pointConstraint> patchPatchPointConstraints(nPatchPatchPoints);
78 
79  label pppi = 0;
80 
81  forAll(bm, patchi)
82  {
83  if(!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
84  {
85  const labelList& bp = bm[patchi].patch().boundaryPoints();
86  const labelList& meshPoints = bm[patchi].patch().meshPoints();
87 
88  forAll(bp, pointi)
89  {
90  label ppp = meshPoints[bp[pointi]];
91 
92  Map<label>::iterator iter = patchPatchPointSet.find(ppp);
93 
94  if (iter == patchPatchPointSet.end())
95  {
96  patchPatchPointSet.insert(ppp, pppi);
97  patchPatchPoints_[pppi] = ppp;
98 
99  pbm[patchi].applyConstraint
100  (
101  bp[pointi],
102  patchPatchPointConstraints[pppi]
103  );
104  pppi++;
105  }
106  else
107  {
108  pbm[patchi].applyConstraint
109  (
110  bp[pointi],
111  patchPatchPointConstraints[iter()]
112  );
113  }
114  }
115  }
116  }
117 
118  nPatchPatchPoints = pppi;
119  patchPatchPoints_.setSize(nPatchPatchPoints);
120  patchPatchPointConstraints.setSize(nPatchPatchPoints);
121 
122  patchPatchPointConstraintPoints_.setSize(nPatchPatchPoints);
123  patchPatchPointConstraintTensors_.setSize(nPatchPatchPoints);
124 
125  label nConstraints = 0;
126 
127  forAll(patchPatchPointConstraints, i)
128  {
129  if (patchPatchPointConstraints[i].first() != 0)
130  {
131  patchPatchPointConstraintPoints_[nConstraints] =
132  patchPatchPoints_[i];
133 
134  patchPatchPointConstraintTensors_[nConstraints] =
135  patchPatchPointConstraints[i].constraintTransformation();
136 
137  nConstraints++;
138  }
139  }
140 
141  patchPatchPointConstraintPoints_.setSize(nConstraints);
142  patchPatchPointConstraintTensors_.setSize(nConstraints);
143 
144 
145  patchInterpolators_.clear();
146  patchInterpolators_.setSize(bm.size());
147 
148  forAll(bm, patchi)
149  {
150  patchInterpolators_.set
151  (
152  patchi,
153  new primitivePatchInterpolation(bm[patchi].patch())
154  );
155  }
156 
157  if (debug)
158  {
159  Info<< "pointPatchInterpolation::makePatchPatchAddressing() : "
160  << "finished constructing boundary addressing"
161  << endl;
162  }
163 }
164 
165 
166 void pointPatchInterpolation::makePatchPatchWeights()
167 {
168  if (debug)
169  {
170  Info<< "pointPatchInterpolation::makePatchPatchWeights() : "
171  << "constructing boundary weighting factors"
172  << endl;
173  }
174 
175  patchPatchPointWeights_.clear();
176  patchPatchPointWeights_.setSize(patchPatchPoints_.size());
177 
178  const labelListList& pf = fvMesh_.pointFaces();
179  const volVectorField& centres = fvMesh_.C();
180  const fvBoundaryMesh& bm = fvMesh_.boundary();
181 
182  pointScalarField sumWeights
183  (
184  IOobject
185  (
186  "sumWeights",
187  fvMesh_.polyMesh::instance(),
188  fvMesh_
189  ),
190  pointMesh::New(fvMesh_),
191  dimensionedScalar("zero", dimless, 0)
192  );
193 
194  forAll(patchPatchPoints_, pointi)
195  {
196  const label curPoint = patchPatchPoints_[pointi];
197  const labelList& curFaces = pf[curPoint];
198 
199  patchPatchPointWeights_[pointi].setSize(curFaces.size());
200  scalarList& pw = patchPatchPointWeights_[pointi];
201 
202  label nFacesAroundPoint = 0;
203 
204  const vector& pointLoc = fvMesh_.points()[curPoint];
205 
206  forAll(curFaces, facei)
207  {
208  if (!fvMesh_.isInternalFace(curFaces[facei]))
209  {
210  label patchi =
211  fvMesh_.boundaryMesh().whichPatch(curFaces[facei]);
212 
213  if (!isA<emptyFvPatch>(bm[patchi]) && !bm[patchi].coupled())
214  {
215  vector d =
216  pointLoc
217  - centres.boundaryField()[patchi]
218  [bm[patchi].patch().whichFace(curFaces[facei])];
219 
220  pw[nFacesAroundPoint] = 1.0/(mag(d)+VSMALL);
221 
222  nFacesAroundPoint++;
223  }
224  }
225  }
226 
227  // Reset the sizes of the local weights
228  pw.setSize(nFacesAroundPoint);
229 
230  // Collect the sum of weights for parallel correction
231  sumWeights[curPoint] += sum(pw);
232  }
233 
234  // Do parallel correction of weights
235 
236  // Update coupled boundaries
237  forAll(sumWeights.boundaryField(), patchi)
238  {
239  if (sumWeights.boundaryField()[patchi].coupled())
240  {
241  refCast<coupledPointPatchScalarField>
242  (sumWeights.boundaryField()[patchi]).initSwapAdd
243  (
244  sumWeights.internalField()
245  );
246  }
247  }
248 
249  forAll(sumWeights.boundaryField(), patchi)
250  {
251  if (sumWeights.boundaryField()[patchi].coupled())
252  {
253  refCast<coupledPointPatchScalarField>
254  (sumWeights.boundaryField()[patchi]).swapAdd
255  (
256  sumWeights.internalField()
257  );
258  }
259  }
260 
261 
262  // Re-scale the weights for the current point
263  forAll(patchPatchPoints_, pointi)
264  {
265  scalarList& pw = patchPatchPointWeights_[pointi];
266  scalar sumw = sumWeights[patchPatchPoints_[pointi]];
267 
268  forAll(pw, facei)
269  {
270  pw[facei] /= sumw;
271  }
272  }
273 
274 
275  if (debug)
276  {
277  Info<< "pointPatchInterpolation::makePatchPatchWeights() : "
278  << "finished constructing boundary weighting factors"
279  << endl;
280  }
281 }
282 
283 
284 // * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * //
285 
286 pointPatchInterpolation::pointPatchInterpolation(const fvMesh& vm)
287 :
288  fvMesh_(vm)
289 {
290  updateMesh();
291 }
292 
293 
294 // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * //
295 
297 {}
298 
299 
300 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
301 
303 {
304  makePatchPatchAddressing();
305  makePatchPatchWeights();
306 }
307 
308 
310 {
311  forAll(patchInterpolators_, patchi)
312  {
313  patchInterpolators_[patchi].movePoints();
314  }
315 
316  makePatchPatchWeights();
317 
318  return true;
319 }
320 
321 
322 // Specialisaion of applyCornerConstraints for scalars because
323 // no constraint need be applied
324 template<>
325 void pointPatchInterpolation::applyCornerConstraints<scalar>
326 (
328 ) const
329 {}
330 
331 
332 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
333 
334 } // End namespace Foam
335 
336 // ************************ vim: set sw=4 sts=4 et: ************************ //