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