FreeFOAM The Cross-Platform CFD Toolkit
enrichedPatchFaces.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 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "enrichedPatch.H"
29 #include <OpenFOAM/DynamicList.H>
30 
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 const Foam::label Foam::enrichedPatch::enrichedFaceRatio_ = 3;
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
39 (
40  const labelListList& pointsIntoMasterEdges,
41  const labelListList& pointsIntoSlaveEdges,
42  const pointField& projectedSlavePoints
43 )
44 {
45  if (enrichedFacesPtr_)
46  {
48  (
49  "void enrichedPatch::calcEnrichedFaces\n"
50  "(\n"
51  " const labelListList& pointsIntoMasterEdges,\n"
52  " const labelListList& pointsIntoSlaveEdges,\n"
53  " const pointField& projectedSlavePoints\n"
54  ")"
55  ) << "Enriched faces already calculated."
56  << abort(FatalError);
57  }
58 
59  // Create a list of enriched faces
60  // Algorithm:
61  // 1) Grab the original face and start from point zero.
62  // 2) If the point has been merged away, grab the merge label;
63  // otherwise, keep the original label.
64  // 3) Go to the next edge. Collect all the points to be added along
65  // the edge; order them in the necessary direction and insert onto the
66  // face.
67  // 4) Grab the next point and return on step 2.
68  enrichedFacesPtr_ = new faceList(masterPatch_.size() + slavePatch_.size());
69  faceList& enrichedFaces = *enrichedFacesPtr_;
70 
71  label nEnrichedFaces = 0;
72 
73  const pointField& masterLocalPoints = masterPatch_.localPoints();
74  const faceList& masterLocalFaces = masterPatch_.localFaces();
75  const labelListList& masterFaceEdges = masterPatch_.faceEdges();
76 
77  const faceList& slaveLocalFaces = slavePatch_.localFaces();
78  const labelListList& slaveFaceEdges = slavePatch_.faceEdges();
79 
80  // For correct functioning of the enrichedPatch class, the slave
81  // faces need to be inserted first. See comments in
82  // enrichedPatch.H
83 
84  // Get reference to the point merge map
85  const Map<label>& pmm = pointMergeMap();
86 
87  // Add slave faces into the enriched faces list
88 
89  forAll (slavePatch_, faceI)
90  {
91  const face oldFace = slavePatch_[faceI];
92  const face oldLocalFace = slaveLocalFaces[faceI];
93 // Info << "old slave face " << faceI << ": " << oldFace << endl;
94  const labelList& curEdges = slaveFaceEdges[faceI];
95 
96  DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
97 
98  // Note: The number of points and edges in a face is always identical
99  // so both can be done is the same loop
100  forAll (oldFace, i)
101  {
102  // Add the point
104  pmm.find(oldFace[i]);
105 
106  if (mpIter == pmm.end())
107  {
108  // Point not mapped
109  newFace.append(oldFace[i]);
110 
111  // Add the projected point into the patch support
112  pointMap().insert
113  (
114  oldFace[i], // Global label of point
115  projectedSlavePoints[oldLocalFace[i]] // Projected position
116  );
117  }
118  else
119  {
120  // Point mapped
121  newFace.append(mpIter());
122 
123  // Add the projected point into the patch support
124  pointMap().insert
125  (
126  mpIter(), // Merged global label of point
127  projectedSlavePoints[oldLocalFace[i]] // Projected position
128  );
129  }
130 
131  // Grab the edge points
132 
133  const labelList& slavePointsOnEdge =
134  pointsIntoSlaveEdges[curEdges[i]];
135 // Info << "slavePointsOnEdge for " << curEdges[i] << ": " << slavePointsOnEdge << endl;
136  // If there are no points on the edge, skip everything
137  // If there is only one point, no need for sorting
138  if (slavePointsOnEdge.size())
139  {
140  // Sort edge points in order
141  scalarField edgePointWeights(slavePointsOnEdge.size());
142  const point& startPoint = projectedSlavePoints[oldLocalFace[i]];
143 
144  vector e =
145  projectedSlavePoints[oldLocalFace.nextLabel(i)]
146  - startPoint;
147 
148  scalar magSqrE = magSqr(e);
149 
150  if (magSqrE > SMALL)
151  {
152  e /= magSqrE;
153  }
154  else
155  {
157  (
158  "void enrichedPatch::calcEnrichedFaces\n"
159  "(\n"
160  " const labelListList& pointsIntoMasterEdges,\n"
161  " const labelListList& pointsIntoSlaveEdges,\n"
162  " const pointField& projectedSlavePoints\n"
163  ")"
164  ) << "Zero length edge in slave patch for face " << i
165  << ". This is not allowed."
166  << abort(FatalError);
167  }
168 
169  pointField slavePosOnEdge(slavePointsOnEdge.size());
170 
171  forAll (slavePointsOnEdge, edgePointI)
172  {
173  slavePosOnEdge[edgePointI] =
174  pointMap().find(slavePointsOnEdge[edgePointI])();
175 
176  edgePointWeights[edgePointI] =
177  (e & (slavePosOnEdge[edgePointI] - startPoint));
178  }
179 
180  if (debug)
181  {
182  // Check weights: all new points should be on the edge
183  if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
184  {
186  (
187  "void enrichedPatch::calcEnrichedFaces\n"
188  "(\n"
189  " const labelListList& pointsIntoMasterEdges,\n"
190  " const labelListList& pointsIntoSlaveEdges,\n"
191  " const pointField& projectedSlavePoints\n"
192  ")"
193  ) << "Invalid point edge weights. Some of points are"
194  << " not on the edge for edge " << curEdges[i]
195  << " of face " << faceI << " in slave patch." << nl
196  << "Min weight: " << min(edgePointWeights)
197  << " Max weight: " << max(edgePointWeights)
198  << abort(FatalError);
199  }
200  }
201 
202  // Go through the points and collect them based on
203  // weights from lower to higher. This gives the
204  // correct order of points along the edge.
205  for (label passI = 0; passI < edgePointWeights.size(); passI++)
206  {
207  // Max weight can only be one, so the sorting is
208  // done by elimination.
209  label nextPoint = -1;
210  scalar dist = 2;
211 
212  forAll (edgePointWeights, wI)
213  {
214  if (edgePointWeights[wI] < dist)
215  {
216  dist = edgePointWeights[wI];
217  nextPoint = wI;
218  }
219  }
220 
221  // Insert the next point and reset its weight to exclude it
222  // from future picks
223  newFace.append(slavePointsOnEdge[nextPoint]);
224  edgePointWeights[nextPoint] = GREAT;
225 
226  // Add the point into patch support
227  pointMap().insert
228  (
229  slavePointsOnEdge[nextPoint],
230  slavePosOnEdge[nextPoint]
231  );
232  }
233  }
234  }
235  // Info<< "New slave face " << faceI << ": " << newFace << endl;
236 
237  // Add the new face to the list
238  enrichedFaces[nEnrichedFaces].transfer(newFace);
239  nEnrichedFaces++;
240  }
241 
242  // Add master faces into the enriched faces list
243 
244  forAll (masterPatch_, faceI)
245  {
246  const face& oldFace = masterPatch_[faceI];
247  const face& oldLocalFace = masterLocalFaces[faceI];
248 // Info << "old master face: " << oldFace << endl;
249  const labelList& curEdges = masterFaceEdges[faceI];
250 
251  DynamicList<label> newFace(oldFace.size()*enrichedFaceRatio_);
252 
253  // Note: The number of points and edges in a face is always identical
254  // so both can be done is the same loop
255  forAll (oldFace, i)
256  {
257  // Add the point
259  pmm.find(oldFace[i]);
260 
261  if (mpIter == pmm.end())
262  {
263  // Point not mapped
264  newFace.append(oldFace[i]);
265 
266  // Add the point into patch support
267  pointMap().insert
268  (
269  oldFace[i],
270  masterLocalPoints[oldLocalFace[i]]
271  );
272  }
273  else
274  {
275  // Point mapped
276  newFace.append(mpIter());
277 
278  // Add the point into support
279  pointMap().insert(mpIter(), masterLocalPoints[oldLocalFace[i]]);
280  }
281 
282  // Grab the edge points
283 
284  const labelList& masterPointsOnEdge =
285  pointsIntoMasterEdges[curEdges[i]];
286 
287  // If there are no points on the edge, skip everything
288  // If there is only one point, no need for sorting
289  if (masterPointsOnEdge.size())
290  {
291  // Sort edge points in order
292  scalarField edgePointWeights(masterPointsOnEdge.size());
293  const point& startPoint = masterLocalPoints[oldLocalFace[i]];
294 
295  vector e =
296  masterLocalPoints[oldLocalFace.nextLabel(i)]
297  - startPoint;
298 
299  scalar magSqrE = magSqr(e);
300 
301  if (magSqrE > SMALL)
302  {
303  e /= magSqrE;
304  }
305  else
306  {
308  (
309  "void enrichedPatch::calcEnrichedFaces\n"
310  "(\n"
311  " const labelListList& pointsIntoMasterEdges,\n"
312  " const labelListList& pointsIntoSlaveEdges,\n"
313  " const pointField& projectedSlavePoints\n"
314  ")"
315  ) << "Zero length edge in master patch for face " << i
316  << ". This is not allowed."
317  << abort(FatalError);
318  }
319 
320  pointField masterPosOnEdge(masterPointsOnEdge.size());
321 
322  forAll (masterPointsOnEdge, edgePointI)
323  {
324  masterPosOnEdge[edgePointI] =
325  pointMap().find(masterPointsOnEdge[edgePointI])();
326 
327  edgePointWeights[edgePointI] =
328  (e & (masterPosOnEdge[edgePointI] - startPoint));
329  }
330 
331  if (debug)
332  {
333  // Check weights: all new points should be on the edge
334  if (min(edgePointWeights) < 0 || max(edgePointWeights) > 1)
335  {
337  (
338  "void enrichedPatch::calcEnrichedFaces\n"
339  "(\n"
340  " const labelListList& pointsIntoMasterEdges,\n"
341  " const labelListList& pointsIntoSlaveEdges,\n"
342  " const pointField& projectedSlavePoints\n"
343  ")"
344  ) << "Invalid point edge weights. Some of points are"
345  << " not on the edge for edge " << curEdges[i]
346  << " of face " << faceI << " in master patch." << nl
347  << "Min weight: " << min(edgePointWeights)
348  << " Max weight: " << max(edgePointWeights)
349  << abort(FatalError);
350  }
351  }
352 
353  // Go through the points and collect them based on
354  // weights from lower to higher. This gives the
355  // correct order of points along the edge.
356  for (label pass = 0; pass < edgePointWeights.size(); pass++)
357  {
358  // Max weight can only be one, so the sorting is
359  // done by elimination.
360  label nextPoint = -1;
361  scalar dist = 2;
362 
363  forAll (edgePointWeights, wI)
364  {
365  if (edgePointWeights[wI] < dist)
366  {
367  dist = edgePointWeights[wI];
368  nextPoint = wI;
369  }
370  }
371 
372  // Insert the next point and reset its weight to exclude it
373  // from future picks
374  newFace.append(masterPointsOnEdge[nextPoint]);
375  edgePointWeights[nextPoint] = GREAT;
376 
377  // Add the point into patch support
378  pointMap().insert
379  (
380  masterPointsOnEdge[nextPoint],
381  masterPosOnEdge[nextPoint]
382  );
383  }
384  }
385  }
386  // Info<< "New master face: " << newFace << endl;
387 
388  // Add the new face to the list
389  enrichedFaces[nEnrichedFaces].transfer(newFace);
390  nEnrichedFaces++;
391  }
392 
393  // Check the support for the enriched patch
394  if (debug)
395  {
396  if (!checkSupport())
397  {
398  Info<< "Enriched patch support OK. Slave faces: "
399  << slavePatch_.size() << " Master faces: "
400  << masterPatch_.size() << endl;
401  }
402  else
403  {
405  (
406  "void enrichedPatch::calcEnrichedFaces\n"
407  "(\n"
408  " const labelListList& pointsIntoMasterEdges,\n"
409  " const labelListList& pointsIntoSlaveEdges,\n"
410  " const pointField& projectedSlavePoints\n"
411  ")"
412  ) << "Error in enriched patch support"
413  << abort(FatalError);
414  }
415  }
416 }
417 
418 
419 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
420 
422 {
423  if (!enrichedFacesPtr_)
424  {
425  FatalErrorIn("const faceList& enrichedPatch::enrichedFaces() const")
426  << "Enriched faces not available yet. Please use "
427  << "void enrichedPatch::calcEnrichedFaces\n"
428  << "(\n"
429  << " const labelListList& pointsIntoMasterEdges,\n"
430  << " const labelListList& pointsIntoSlaveEdges,\n"
431  << " const pointField& projectedSlavePoints\n"
432  << ")"
433  << " before trying to access faces."
434  << abort(FatalError);
435  }
436 
437  return *enrichedFacesPtr_;
438 }
439 
440 
441 // ************************ vim: set sw=4 sts=4 et: ************************ //