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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2017-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "slidingInterface.H"
30#include "polyMesh.H"
31#include "mapPolyMesh.H"
32#include "polyTopoChanger.H"
33
34// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
35
36void Foam::slidingInterface::calcAttachedAddressing() const
37{
38 if (debug)
39 {
41 << " for object " << name() << " : "
42 << "Calculating zone face-cell addressing."
43 << endl;
44 }
45
46 if (!attached_)
47 {
48 // Clear existing addressing
49 clearAttachedAddressing();
50
51 const polyMesh& mesh = topoChanger().mesh();
52 const labelList& own = mesh.faceOwner();
53 const labelList& nei = mesh.faceNeighbour();
54 const faceZoneMesh& faceZones = mesh.faceZones();
55
56 // Master side
57
58 const primitiveFacePatch& masterPatch =
59 faceZones[masterFaceZoneID_.index()]();
60
61 const labelList& masterPatchFaces =
62 faceZones[masterFaceZoneID_.index()];
63
64 const boolList& masterFlip =
65 faceZones[masterFaceZoneID_.index()].flipMap();
66
67 masterFaceCellsPtr_.reset(new labelList(masterPatchFaces.size()));
68 auto& mfc = *masterFaceCellsPtr_;
69
70 forAll(masterPatchFaces, facei)
71 {
72 if (masterFlip[facei])
73 {
74 mfc[facei] = nei[masterPatchFaces[facei]];
75 }
76 else
77 {
78 mfc[facei] = own[masterPatchFaces[facei]];
79 }
80 }
81
82 // Slave side
83
84 const primitiveFacePatch& slavePatch =
85 faceZones[slaveFaceZoneID_.index()]();
86
87 const labelList& slavePatchFaces =
88 faceZones[slaveFaceZoneID_.index()];
89
90 const boolList& slaveFlip =
91 faceZones[slaveFaceZoneID_.index()].flipMap();
92
93 slaveFaceCellsPtr_.reset(new labelList(slavePatchFaces.size()));
94 auto& sfc = *slaveFaceCellsPtr_;
95
96 forAll(slavePatchFaces, facei)
97 {
98 if (slaveFlip[facei])
99 {
100 sfc[facei] = nei[slavePatchFaces[facei]];
101 }
102 else
103 {
104 sfc[facei] = own[slavePatchFaces[facei]];
105 }
106 }
107
108 // Check that the addressing is valid
109 if (min(mfc) < 0 || min(sfc) < 0)
110 {
111 if (debug)
112 {
113 forAll(mfc, facei)
114 {
115 if (mfc[facei] < 0)
116 {
117 Pout<< "No cell next to master patch face " << facei
118 << ". Global face no: " << mfc[facei]
119 << " own: " << own[masterPatchFaces[facei]]
120 << " nei: " << nei[masterPatchFaces[facei]]
121 << " flip: " << masterFlip[facei] << endl;
122 }
123 }
124
125 forAll(sfc, facei)
126 {
127 if (sfc[facei] < 0)
128 {
129 Pout<< "No cell next to slave patch face " << facei
130 << ". Global face no: " << sfc[facei]
131 << " own: " << own[slavePatchFaces[facei]]
132 << " nei: " << nei[slavePatchFaces[facei]]
133 << " flip: " << slaveFlip[facei] << endl;
134 }
135 }
136 }
137
139 << "decoupled mesh or sliding interface definition."
140 << abort(FatalError);
141 }
142
143 // Calculate stick-out faces
144 const labelListList& pointFaces = mesh.pointFaces();
145
146 labelHashSet stickOutFaceMap
147 (
149 * max(masterPatch.size(), slavePatch.size())
150 );
151
152 // Master side
153 const labelList& masterMeshPoints = masterPatch.meshPoints();
154
155 stickOutFaceMap.clear();
156
157 for (const label pointi : masterMeshPoints)
158 {
159 for (const label facei : pointFaces[pointi])
160 {
161 const label zoneIdx = faceZones.whichZone(facei);
162
163 // Add if face not already part of master or slave face zone
164 // This handles partially attached faces.
165 if
166 (
167 zoneIdx != masterFaceZoneID_.index()
168 && zoneIdx != slaveFaceZoneID_.index()
169 )
170 {
171 stickOutFaceMap.insert(facei);
172 }
173 }
174 }
175
176 masterStickOutFacesPtr_.reset(new labelList(stickOutFaceMap.toc()));
177
178 // Sort in debug mode for easier diagnostics
179 if (debug)
180 {
181 Foam::sort(*masterStickOutFacesPtr_);
182 }
183
184 // Slave side
185 const labelList& slaveMeshPoints = slavePatch.meshPoints();
186
187 stickOutFaceMap.clear();
188
189 for (const label pointi : slaveMeshPoints)
190 {
191 for (const label facei : pointFaces[pointi])
192 {
193 const label zoneIdx = faceZones.whichZone(facei);
194
195 // Add if face not already part of master or slave face zone
196 // This handles partially attached faces.
197 if
198 (
199 zoneIdx != masterFaceZoneID_.index()
200 && zoneIdx != slaveFaceZoneID_.index()
201 )
202 {
203 stickOutFaceMap.insert(facei);
204 }
205 }
206 }
207
208 slaveStickOutFacesPtr_.reset(new labelList(stickOutFaceMap.toc()));
209
210 // Sort in debug mode for easier diagnostics
211 if (debug)
212 {
213 Foam::sort(*slaveStickOutFacesPtr_);
214 }
215
216 stickOutFaceMap.clear();
217
218 // Retired point addressing does not exist at this stage.
219 // It will be filled when the interface is coupled.
220 retiredPointMapPtr_.reset
221 (
222 new Map<label>
223 (
224 2*faceZones[slaveFaceZoneID_.index()]().nPoints()
225 )
226 );
227
228 // Ditto for cut point edge map. This is a rough guess of its size
229 cutPointEdgePairMapPtr_.reset
230 (
231 new Map<Pair<edge>>
232 (
233 faceZones[slaveFaceZoneID_.index()]().nEdges()
234 )
235 );
236 }
237 else
238 {
240 << "cannot be assembled for object " << name()
241 << abort(FatalError);
242 }
243
244 if (debug)
245 {
247 << " for object " << name() << " : "
248 << "Finished calculating zone face-cell addressing."
249 << endl;
250 }
251}
252
253
254void Foam::slidingInterface::clearAttachedAddressing() const
255{
256 masterFaceCellsPtr_.reset(nullptr);
257 slaveFaceCellsPtr_.reset(nullptr);
258
259 masterStickOutFacesPtr_.reset(nullptr);
260 slaveStickOutFacesPtr_.reset(nullptr);
261
262 retiredPointMapPtr_.reset(nullptr);
263 cutPointEdgePairMapPtr_.reset(nullptr);
264}
265
266
267void Foam::slidingInterface::renumberAttachedAddressing
268(
269 const mapPolyMesh& m
270) const
271{
272 // Get reference to reverse cell renumbering
273 // The renumbering map is needed the other way around, i.e. giving
274 // the new cell number for every old cell next to the interface.
275 const labelList& reverseCellMap = m.reverseCellMap();
276
277 const labelList& mfc = masterFaceCells();
278 const labelList& sfc = slaveFaceCells();
279
280 // Master side
281 unique_ptr<labelList> newMfcPtr(new labelList(mfc.size(), -1));
282 auto& newMfc = *newMfcPtr;
283
284 const labelList& mfzRenumber =
285 m.faceZoneFaceMap()[masterFaceZoneID_.index()];
286
287 forAll(mfc, facei)
288 {
289 label newCelli = reverseCellMap[mfc[mfzRenumber[facei]]];
290
291 if (newCelli >= 0)
292 {
293 newMfc[facei] = newCelli;
294 }
295 }
296
297 // Slave side
298 unique_ptr<labelList> newSfcPtr(new labelList(sfc.size(), -1));
299 auto& newSfc = *newSfcPtr;
300
301 const labelList& sfzRenumber =
302 m.faceZoneFaceMap()[slaveFaceZoneID_.index()];
303
304 forAll(sfc, facei)
305 {
306 label newCelli = reverseCellMap[sfc[sfzRenumber[facei]]];
307
308 if (newCelli >= 0)
309 {
310 newSfc[facei] = newCelli;
311 }
312 }
313
314 if (debug)
315 {
316 // Check if all the mapped cells are live
317 if (min(newMfc) < 0 || min(newSfc) < 0)
318 {
320 << "Error in cell renumbering for object " << name()
321 << ". Some of master cells next "
322 << "to the interface have been removed."
323 << abort(FatalError);
324 }
325 }
326
327 // Renumber stick-out faces
328
329 const labelList& reverseFaceMap = m.reverseFaceMap();
330
331 // Master side
332 const labelList& msof = masterStickOutFaces();
333
334 unique_ptr<labelList> newMsofPtr(new labelList(msof.size(), -1));
335 auto& newMsof = *newMsofPtr;
336
337 forAll(msof, facei)
338 {
339 label newFacei = reverseFaceMap[msof[facei]];
340
341 if (newFacei >= 0)
342 {
343 newMsof[facei] = newFacei;
344 }
345 }
346// Pout<< "newMsof: " << newMsof << endl;
347 // Slave side
348 const labelList& ssof = slaveStickOutFaces();
349
350 unique_ptr<labelList> newSsofPtr(new labelList(ssof.size(), -1));
351 auto& newSsof = *newSsofPtr;
352
353 forAll(ssof, facei)
354 {
355 label newFacei = reverseFaceMap[ssof[facei]];
356
357 if (newFacei >= 0)
358 {
359 newSsof[facei] = newFacei;
360 }
361 }
362// Pout<< "newSsof: " << newSsof << endl;
363 if (debug)
364 {
365 // Check if all the mapped cells are live
366 if (min(newMsof) < 0 || min(newSsof) < 0)
367 {
369 << "Error in face renumbering for object " << name()
370 << ". Some of stick-out next "
371 << "to the interface have been removed."
372 << abort(FatalError);
373 }
374 }
375
376 // Renumber the retired point map. Need to take a copy!
377 const Map<label> rpm = retiredPointMap();
378
379 unique_ptr<Map<label>> newRpmPtr(new Map<label>(rpm.size()));
380 auto& newRpm = *newRpmPtr;
381
382 // Get reference to point renumbering
383 const labelList& reversePointMap = m.reversePointMap();
384
385 forAllConstIters(rpm, iter)
386 {
387 const label key = reversePointMap[iter.key()];
388 const label val = reversePointMap[iter.val()];
389
390 if (debug)
391 {
392 // Check if all the mapped cells are live
393 if (key < 0 || val < 0)
394 {
396 << "Error in retired point numbering for object "
397 << name() << ". Some of master "
398 << "points have been removed."
399 << abort(FatalError);
400 }
401 }
402
403 newRpm.insert(key, val);
404 }
405
406 // Renumber the cut point edge pair map. Need to take a copy!
407 const Map<Pair<edge>> cpepm = cutPointEdgePairMap();
408
409 unique_ptr<Map<Pair<edge>>> newCpepmPtr(new Map<Pair<edge>>(cpepm.size()));
410 auto& newCpepm = *newCpepmPtr;
411
412 forAllConstIters(cpepm, iter)
413 {
414 const label key = reversePointMap[iter.key()];
415
416 const Pair<edge>& oldPe = iter.val();
417
418 // Re-do the edges in global addressing
419 const label ms = reversePointMap[oldPe.first().start()];
420 const label me = reversePointMap[oldPe.first().end()];
421
422 const label ss = reversePointMap[oldPe.second().start()];
423 const label se = reversePointMap[oldPe.second().end()];
424
425 if (debug)
426 {
427 // Check if all the mapped cells are live
428 if (key < 0 || ms < 0 || me < 0 || ss < 0 || se < 0)
429 {
431 << "Error in cut point edge pair map numbering for object "
432 << name() << ". Some of master points have been removed."
433 << abort(FatalError);
434 }
435 }
436
437 newCpepm.insert(key, Pair<edge>(edge(ms, me), edge(ss, se)));
438 }
439
440 if (!projectedSlavePointsPtr_)
441 {
443 << "Error in projected point numbering for object " << name()
444 << abort(FatalError);
445 }
446
447 // Renumber the projected slave zone points
448 const pointField& projectedSlavePoints = *projectedSlavePointsPtr_;
449
450 unique_ptr<pointField> newProjectedSlavePointsPtr
451 (
452 new pointField(projectedSlavePoints.size())
453 );
454 auto& newProjectedSlavePoints = *newProjectedSlavePointsPtr;
455
456 const labelList& sfzPointRenumber =
457 m.faceZonePointMap()[slaveFaceZoneID_.index()];
458
459 forAll(newProjectedSlavePoints, pointi)
460 {
461 if (sfzPointRenumber[pointi] > -1)
462 {
463 newProjectedSlavePoints[pointi] =
464 projectedSlavePoints[sfzPointRenumber[pointi]];
465 }
466 }
467
468 // Re-set the lists
469 clearAttachedAddressing();
470
471 projectedSlavePointsPtr_.reset(nullptr);
472
473 masterFaceCellsPtr_ = std::move(newMfcPtr);
474 slaveFaceCellsPtr_ = std::move(newSfcPtr);
475
476 masterStickOutFacesPtr_ = std::move(newMsofPtr);
477 slaveStickOutFacesPtr_ = std::move(newSsofPtr);
478
479 retiredPointMapPtr_ = std::move(newRpmPtr);
480 cutPointEdgePairMapPtr_ = std::move(newCpepmPtr);
481 projectedSlavePointsPtr_ = std::move(newProjectedSlavePointsPtr);
482}
483
484
485const Foam::labelList& Foam::slidingInterface::masterFaceCells() const
486{
487 if (!masterFaceCellsPtr_)
488 {
490 << "Master zone face-cell addressing not available for object "
491 << name()
492 << abort(FatalError);
493 }
494
495 return *masterFaceCellsPtr_;
496}
497
498
499const Foam::labelList& Foam::slidingInterface::slaveFaceCells() const
500{
501 if (!slaveFaceCellsPtr_)
502 {
504 << "Slave zone face-cell addressing not available for object "
505 << name()
506 << abort(FatalError);
507 }
508
509 return *slaveFaceCellsPtr_;
510}
511
512
513const Foam::labelList& Foam::slidingInterface::masterStickOutFaces() const
514{
515 if (!masterStickOutFacesPtr_)
516 {
518 << "Master zone stick-out face addressing not available for object "
519 << name()
520 << abort(FatalError);
521 }
522
523 return *masterStickOutFacesPtr_;
524}
525
526
527const Foam::labelList& Foam::slidingInterface::slaveStickOutFaces() const
528{
529 if (!slaveStickOutFacesPtr_)
530 {
532 << "Slave zone stick-out face addressing not available for object "
533 << name()
534 << abort(FatalError);
535 }
536
537 return *slaveStickOutFacesPtr_;
538}
539
540
541const Foam::Map<Foam::label>& Foam::slidingInterface::retiredPointMap() const
542{
543 if (!retiredPointMapPtr_)
544 {
546 << "Retired point map not available for object " << name()
547 << abort(FatalError);
548 }
549
550 return *retiredPointMapPtr_;
551}
552
553
555Foam::slidingInterface::cutPointEdgePairMap() const
556{
557 if (!cutPointEdgePairMapPtr_)
558 {
560 << "Retired point map not available for object " << name()
561 << abort(FatalError);
562 }
563
564 return *cutPointEdgePairMapPtr_;
565}
566
567
568// ************************************************************************* //
label index() const
The index of the first matching items, -1 if no matches.
Definition: DynamicID.H:123
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
A HashTable to objects of type <T> with a label key.
Definition: Map.H:60
const word & name() const
Return name of this modifier.
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
const faceZoneMesh & faceZones() const noexcept
Return face zone mesh.
Definition: polyMesh.H:498
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1127
const polyMesh & mesh() const
Return the mesh reference.
const labelListList & pointFaces() const
static const unsigned facesPerCell_
Estimated number of faces per cell.
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define FUNCTION_NAME
const dimensionedScalar me
Electron mass.
auto key(const Type &t) -> typename std::enable_if< std::is_enum< Type >::value, typename std::underlying_type< Type >::type >::type
Definition: foamGltfBase.H:108
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
List< label > labelList
A List of labels.
Definition: List.H:66
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void sort(UList< T > &list)
Sort the list.
Definition: UList.C:342
errorManip< error > abort(error &err)
Definition: errorManip.H:144
List< bool > boolList
A List of bools.
Definition: List.H:64
HashSet< label, Hash< label > > labelHashSet
A HashSet of labels, uses label hasher.
Definition: HashSet.H:85
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
A PrimitivePatch with List storage for the faces, const reference for the point field.
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278