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