slidingInterfaceAttachedAddressing.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 "slidingInterface.H"
27 #include "polyMesh.H"
28 #include "mapPolyMesh.H"
29 #include "polyTopoChanger.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 void Foam::slidingInterface::calcAttachedAddressing() const
34 {
35  if (debug)
36  {
37  Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
38  << " for object " << name() << " : "
39  << "Calculating zone face-cell addressing."
40  << endl;
41  }
42 
43  if (!attached_)
44  {
45  // Clear existing addressing
46  clearAttachedAddressing();
47 
48  const polyMesh& mesh = topoChanger().mesh();
49  const labelList& own = mesh.faceOwner();
50  const labelList& nei = mesh.faceNeighbour();
51  const faceZoneMesh& faceZones = mesh.faceZones();
52 
53  // Master side
54 
55  const primitiveFacePatch& masterPatch =
56  faceZones[masterFaceZoneID_.index()]();
57 
58  const labelList& masterPatchFaces =
59  faceZones[masterFaceZoneID_.index()];
60 
61  const boolList& masterFlip =
62  faceZones[masterFaceZoneID_.index()].flipMap();
63 
64  masterFaceCellsPtr_ = new labelList(masterPatchFaces.size());
65  labelList& mfc = *masterFaceCellsPtr_;
66 
67  forAll(masterPatchFaces, facei)
68  {
69  if (masterFlip[facei])
70  {
71  mfc[facei] = nei[masterPatchFaces[facei]];
72  }
73  else
74  {
75  mfc[facei] = own[masterPatchFaces[facei]];
76  }
77  }
78 
79  // Slave side
80 
81  const primitiveFacePatch& slavePatch =
82  faceZones[slaveFaceZoneID_.index()]();
83 
84  const labelList& slavePatchFaces =
85  faceZones[slaveFaceZoneID_.index()];
86 
87  const boolList& slaveFlip =
88  faceZones[slaveFaceZoneID_.index()].flipMap();
89 
90  slaveFaceCellsPtr_ = new labelList(slavePatchFaces.size());
91  labelList& sfc = *slaveFaceCellsPtr_;
92 
93  forAll(slavePatchFaces, facei)
94  {
95  if (slaveFlip[facei])
96  {
97  sfc[facei] = nei[slavePatchFaces[facei]];
98  }
99  else
100  {
101  sfc[facei] = own[slavePatchFaces[facei]];
102  }
103  }
104 
105  // Check that the addressing is valid
106  if (min(mfc) < 0 || min(sfc) < 0)
107  {
108  if (debug)
109  {
110  forAll(mfc, facei)
111  {
112  if (mfc[facei] < 0)
113  {
114  Pout<< "No cell next to master patch face " << facei
115  << ". Global face no: " << mfc[facei]
116  << " own: " << own[masterPatchFaces[facei]]
117  << " nei: " << nei[masterPatchFaces[facei]]
118  << " flip: " << masterFlip[facei] << endl;
119  }
120  }
121 
122  forAll(sfc, facei)
123  {
124  if (sfc[facei] < 0)
125  {
126  Pout<< "No cell next to slave patch face " << facei
127  << ". Global face no: " << sfc[facei]
128  << " own: " << own[slavePatchFaces[facei]]
129  << " nei: " << nei[slavePatchFaces[facei]]
130  << " flip: " << slaveFlip[facei] << endl;
131  }
132  }
133  }
134 
136  << "decoupled mesh or sliding interface definition."
137  << abort(FatalError);
138  }
139 
140  // Calculate stick-out faces
141  const labelListList& pointFaces = mesh.pointFaces();
142 
143  // Master side
144  labelHashSet masterStickOutFaceMap
145  (
146  primitiveMesh::facesPerCell_*(masterPatch.size())
147  );
148 
149  const labelList& masterMeshPoints = masterPatch.meshPoints();
150 
151  forAll(masterMeshPoints, pointi)
152  {
153  const labelList& curFaces = pointFaces[masterMeshPoints[pointi]];
154 
155  forAll(curFaces, facei)
156  {
157  // Check if the face belongs to the master face zone;
158  // if not add it
159  if
160  (
161  faceZones.whichZone(curFaces[facei])
162  != masterFaceZoneID_.index()
163  )
164  {
165  masterStickOutFaceMap.insert(curFaces[facei]);
166  }
167  }
168  }
169 
170  masterStickOutFacesPtr_ = new labelList(masterStickOutFaceMap.toc());
171 
172  // Slave side
173  labelHashSet slaveStickOutFaceMap
174  (
175  primitiveMesh::facesPerCell_*(slavePatch.size())
176  );
177 
178  const labelList& slaveMeshPoints = slavePatch.meshPoints();
179 
180  forAll(slaveMeshPoints, pointi)
181  {
182  const labelList& curFaces = pointFaces[slaveMeshPoints[pointi]];
183 
184  forAll(curFaces, facei)
185  {
186  // Check if the face belongs to the slave face zone;
187  // if not add it
188  if
189  (
190  faceZones.whichZone(curFaces[facei])
191  != slaveFaceZoneID_.index()
192  )
193  {
194  slaveStickOutFaceMap.insert(curFaces[facei]);
195  }
196  }
197  }
198 
199  slaveStickOutFacesPtr_ = new labelList(slaveStickOutFaceMap.toc());
200 
201 
202  // Retired point addressing does not exist at this stage.
203  // It will be filled when the interface is coupled.
204  retiredPointMapPtr_ =
205  new Map<label>
206  (
207  2*faceZones[slaveFaceZoneID_.index()]().nPoints()
208  );
209 
210  // Ditto for cut point edge map. This is a rough guess of its size
211  //
212  cutPointEdgePairMapPtr_ =
213  new Map<Pair<edge>>
214  (
215  faceZones[slaveFaceZoneID_.index()]().nEdges()
216  );
217  }
218  else
219  {
221  << "cannot be assembled for object " << name()
222  << abort(FatalError);
223  }
224 
225  if (debug)
226  {
227  Pout<< "void Foam::slidingInterface::calcAttachedAddressing() const "
228  << " for object " << name() << " : "
229  << "Finished calculating zone face-cell addressing."
230  << endl;
231  }
232 }
233 
234 
235 void Foam::slidingInterface::clearAttachedAddressing() const
236 {
237  deleteDemandDrivenData(masterFaceCellsPtr_);
238  deleteDemandDrivenData(slaveFaceCellsPtr_);
239 
240  deleteDemandDrivenData(masterStickOutFacesPtr_);
241  deleteDemandDrivenData(slaveStickOutFacesPtr_);
242 
243  deleteDemandDrivenData(retiredPointMapPtr_);
244  deleteDemandDrivenData(cutPointEdgePairMapPtr_);
245 }
246 
247 
248 void Foam::slidingInterface::renumberAttachedAddressing
249 (
250  const mapPolyMesh& m
251 ) const
252 {
253  // Get reference to reverse cell renumbering
254  // The renumbering map is needed the other way around, i.e. giving
255  // the new cell number for every old cell next to the interface.
256  const labelList& reverseCellMap = m.reverseCellMap();
257 
258  const labelList& mfc = masterFaceCells();
259  const labelList& sfc = slaveFaceCells();
260 
261  // Master side
262  labelList* newMfcPtr = new labelList(mfc.size(), -1);
263  labelList& newMfc = *newMfcPtr;
264 
265  const labelList& mfzRenumber =
266  m.faceZoneFaceMap()[masterFaceZoneID_.index()];
267 
268  forAll(mfc, facei)
269  {
270  label newCelli = reverseCellMap[mfc[mfzRenumber[facei]]];
271 
272  if (newCelli >= 0)
273  {
274  newMfc[facei] = newCelli;
275  }
276  }
277 
278  // Slave side
279  labelList* newSfcPtr = new labelList(sfc.size(), -1);
280  labelList& newSfc = *newSfcPtr;
281 
282  const labelList& sfzRenumber =
283  m.faceZoneFaceMap()[slaveFaceZoneID_.index()];
284 
285  forAll(sfc, facei)
286  {
287  label newCelli = reverseCellMap[sfc[sfzRenumber[facei]]];
288 
289  if (newCelli >= 0)
290  {
291  newSfc[facei] = newCelli;
292  }
293  }
294 
295  if (debug)
296  {
297  // Check if all the mapped cells are live
298  if (min(newMfc) < 0 || min(newSfc) < 0)
299  {
301  << "Error in cell renumbering for object " << name()
302  << ". Some of master cells next "
303  << "to the interface have been removed."
304  << abort(FatalError);
305  }
306  }
307 
308  // Renumber stick-out faces
309 
310  const labelList& reverseFaceMap = m.reverseFaceMap();
311 
312  // Master side
313  const labelList& msof = masterStickOutFaces();
314 
315  labelList* newMsofPtr = new labelList(msof.size(), -1);
316  labelList& newMsof = *newMsofPtr;
317 
318  forAll(msof, facei)
319  {
320  label newFacei = reverseFaceMap[msof[facei]];
321 
322  if (newFacei >= 0)
323  {
324  newMsof[facei] = newFacei;
325  }
326  }
327 // Pout<< "newMsof: " << newMsof << endl;
328  // Slave side
329  const labelList& ssof = slaveStickOutFaces();
330 
331  labelList* newSsofPtr = new labelList(ssof.size(), -1);
332  labelList& newSsof = *newSsofPtr;
333 
334  forAll(ssof, facei)
335  {
336  label newFacei = reverseFaceMap[ssof[facei]];
337 
338  if (newFacei >= 0)
339  {
340  newSsof[facei] = newFacei;
341  }
342  }
343 // Pout<< "newSsof: " << newSsof << endl;
344  if (debug)
345  {
346  // Check if all the mapped cells are live
347  if (min(newMsof) < 0 || min(newSsof) < 0)
348  {
350  << "Error in face renumbering for object " << name()
351  << ". Some of stick-out next "
352  << "to the interface have been removed."
353  << abort(FatalError);
354  }
355  }
356 
357  // Renumber the retired point map. Need to take a copy!
358  const Map<label> rpm = retiredPointMap();
359 
360  Map<label>* newRpmPtr = new Map<label>(rpm.size());
361  Map<label>& newRpm = *newRpmPtr;
362 
363  const labelList rpmToc = rpm.toc();
364 
365  // Get reference to point renumbering
366  const labelList& reversePointMap = m.reversePointMap();
367 
368  label key, value;
369 
370  forAll(rpmToc, rpmTocI)
371  {
372  key = reversePointMap[rpmToc[rpmTocI]];
373 
374  value = reversePointMap[rpm.find(rpmToc[rpmTocI])()];
375 
376  if (debug)
377  {
378  // Check if all the mapped cells are live
379  if (key < 0 || value < 0)
380  {
382  << "Error in retired point numbering for object "
383  << name() << ". Some of master "
384  << "points have been removed."
385  << abort(FatalError);
386  }
387  }
388 
389  newRpm.insert(key, value);
390  }
391 
392  // Renumber the cut point edge pair map. Need to take a copy!
393  const Map<Pair<edge>> cpepm = cutPointEdgePairMap();
394 
395  Map<Pair<edge>>* newCpepmPtr = new Map<Pair<edge>>(cpepm.size());
396  Map<Pair<edge>>& newCpepm = *newCpepmPtr;
397 
398  const labelList cpepmToc = cpepm.toc();
399 
400  forAll(cpepmToc, cpepmTocI)
401  {
402  key = reversePointMap[cpepmToc[cpepmTocI]];
403 
404  const Pair<edge>& oldPe = cpepm.find(cpepmToc[cpepmTocI])();
405 
406  // Re-do the edges in global addressing
407  const label ms = reversePointMap[oldPe.first().start()];
408  const label me = reversePointMap[oldPe.first().end()];
409 
410  const label ss = reversePointMap[oldPe.second().start()];
411  const label se = reversePointMap[oldPe.second().end()];
412 
413  if (debug)
414  {
415  // Check if all the mapped cells are live
416  if (key < 0 || ms < 0 || me < 0 || ss < 0 || se < 0)
417  {
419  << "Error in cut point edge pair map numbering for object "
420  << name() << ". Some of master points have been removed."
421  << abort(FatalError);
422  }
423  }
424 
425  newCpepm.insert(key, Pair<edge>(edge(ms, me), edge(ss, se)));
426  }
427 
428  if (!projectedSlavePointsPtr_)
429  {
431  << "Error in projected point numbering for object " << name()
432  << abort(FatalError);
433  }
434 
435  // Renumber the projected slave zone points
436  const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
437 
438  pointField* newProjectedSlavePointsPtr
439  (
440  new pointField(projectedSlavePoints.size())
441  );
442  pointField& newProjectedSlavePoints = *newProjectedSlavePointsPtr;
443 
444  const labelList& sfzPointRenumber =
445  m.faceZonePointMap()[slaveFaceZoneID_.index()];
446 
447  forAll(newProjectedSlavePoints, pointi)
448  {
449  if (sfzPointRenumber[pointi] > -1)
450  {
451  newProjectedSlavePoints[pointi] =
452  projectedSlavePoints[sfzPointRenumber[pointi]];
453  }
454  }
455 
456  // Re-set the lists
457  clearAttachedAddressing();
458 
459  deleteDemandDrivenData(projectedSlavePointsPtr_);
460 
461  masterFaceCellsPtr_ = newMfcPtr;
462  slaveFaceCellsPtr_ = newSfcPtr;
463 
464  masterStickOutFacesPtr_ = newMsofPtr;
465  slaveStickOutFacesPtr_ = newSsofPtr;
466 
467  retiredPointMapPtr_ = newRpmPtr;
468  cutPointEdgePairMapPtr_ = newCpepmPtr;
469  projectedSlavePointsPtr_ = newProjectedSlavePointsPtr;
470 }
471 
472 
473 const Foam::labelList& Foam::slidingInterface::masterFaceCells() const
474 {
475  if (!masterFaceCellsPtr_)
476  {
478  << "Master zone face-cell addressing not available for object "
479  << name()
480  << abort(FatalError);
481  }
482 
483  return *masterFaceCellsPtr_;
484 }
485 
486 
487 const Foam::labelList& Foam::slidingInterface::slaveFaceCells() const
488 {
489  if (!slaveFaceCellsPtr_)
490  {
492  << "Slave zone face-cell addressing not available for object "
493  << name()
494  << abort(FatalError);
495  }
496 
497  return *slaveFaceCellsPtr_;
498 }
499 
500 
501 const Foam::labelList& Foam::slidingInterface::masterStickOutFaces() const
502 {
503  if (!masterStickOutFacesPtr_)
504  {
506  << "Master zone stick-out face addressing not available for object "
507  << name()
508  << abort(FatalError);
509  }
510 
511  return *masterStickOutFacesPtr_;
512 }
513 
514 
515 const Foam::labelList& Foam::slidingInterface::slaveStickOutFaces() const
516 {
517  if (!slaveStickOutFacesPtr_)
518  {
520  << "Slave zone stick-out face addressing not available for object "
521  << name()
522  << abort(FatalError);
523  }
524 
525  return *slaveStickOutFacesPtr_;
526 }
527 
528 
529 const Foam::Map<Foam::label>& Foam::slidingInterface::retiredPointMap() const
530 {
531  if (!retiredPointMapPtr_)
532  {
534  << "Retired point map not available for object " << name()
535  << abort(FatalError);
536  }
537 
538  return *retiredPointMapPtr_;
539 }
540 
541 
543 Foam::slidingInterface::cutPointEdgePairMap() const
544 {
545  if (!cutPointEdgePairMapPtr_)
546  {
548  << "Retired point map not available for object " << name()
549  << abort(FatalError);
550  }
551 
552  return *cutPointEdgePairMapPtr_;
553 }
554 
555 
556 // ************************************************************************* //
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
#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
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
const word & name() const
Return name of this modifier.
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
const polyMesh & mesh() const
Return the mesh reference.
static const unsigned facesPerCell_
Estimated number of faces per cell.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
label nPoints
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
PrimitivePatch< face, List, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void deleteDemandDrivenData(DataPtr &dataPtr)
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49