detachInterface.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2021 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 "attachDetach.H"
27 #include "polyMesh.H"
28 #include "primitiveMesh.H"
29 #include "polyTopoChange.H"
30 #include "polyTopoChanger.H"
31 #include "polyAddPoint.H"
32 #include "polyModifyFace.H"
33 #include "polyAddFace.H"
34 
35 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
36 
37 void Foam::attachDetach::detachInterface
38 (
39  polyTopoChange& ref
40 ) const
41 {
42  // Algorithm:
43  // 1. Create new points for points of the master face zone
44  // 2. Modify all faces of the master zone, by putting them into the master
45  // patch (look for orientation) and their renumbered mirror images
46  // into the slave patch
47  // 3. Create a point renumbering list, giving a new point index for original
48  // points in the face patch
49  // 4. Grab all faces in cells on the master side and renumber them
50  // using the point renumbering list. Exclude the ones that belong to
51  // the master face zone
52  //
53  // Note on point creation:
54  // In order to take into account the issues related to partial
55  // blocking in an attach/detach mesh modifier, special treatment
56  // is required for the duplication of points on the edge of the
57  // face zone. Points which are shared only by internal edges need
58  // not to be duplicated, as this would propagate the discontinuity
59  // in the mesh beyond the face zone. Therefore, before creating
60  // the new points, check the external edge loop. For each edge
61  // check if the edge is internal (i.e. does not belong to a
62  // patch); if so, exclude both of its points from duplication.
63 
64  if (debug)
65  {
66  Pout<< "void attachDetach::detachInterface("
67  << "polyTopoChange& ref) const "
68  << " for object " << name() << " : "
69  << "Detaching interface" << endl;
70  }
71 
72  const polyMesh& mesh = topoChanger().mesh();
73  const meshFaceZones& meshZones = mesh.faceZones();
74 
75  // Check that zone is in increasing order (needed since adding faces
76  // in same order - otherwise polyTopoChange face ordering will mess up
77  // correspondence)
78  if (debug)
79  {
80  const labelList& faceLabels = meshZones[faceZoneID_.index()];
81  if (faceLabels.size() > 0)
82  {
83  for (label i = 1; i < faceLabels.size(); i++)
84  {
85  if (faceLabels[i] <= faceLabels[i-1])
86  {
88  << "faceZone " << meshZones[faceZoneID_.index()].name()
89  << " does not have mesh face labels in"
90  << " increasing order." << endl
91  << "Face label " << faceLabels[i]
92  << " at position " << i
93  << " is smaller than the previous value "
94  << faceLabels[i-1]
95  << exit(FatalError);
96  }
97  }
98  }
99  }
100 
101 
102 
103  const primitiveFacePatch& masterFaceLayer =
104  meshZones[faceZoneID_.index()]();
105  const pointField& points = mesh.points();
106  const labelListList& meshEdgeFaces = mesh.edgeFaces();
107 
108  const labelList& mp = masterFaceLayer.meshPoints();
109  const edgeList& zoneLocalEdges = masterFaceLayer.edges();
110 
111  const labelList& meshEdges = meshZones[faceZoneID_.index()].meshEdges();
112 
113  // Create the points
114 
115  labelList addedPoints(mp.size(), -1);
116 
117  // Go through boundary edges of the master patch. If all the faces from
118  // this patch are internal, mark the points in the addedPoints lookup
119  // with their original labels to stop duplication
120  label nIntEdges = masterFaceLayer.nInternalEdges();
121 
122  for (label curEdgeID = nIntEdges; curEdgeID < meshEdges.size(); curEdgeID++)
123  {
124  const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
125 
126  bool edgeIsInternal = true;
127 
128  forAll(curFaces, facei)
129  {
130  if (!mesh.isInternalFace(curFaces[facei]))
131  {
132  // The edge belongs to a boundary face
133  edgeIsInternal = false;
134  break;
135  }
136  }
137 
138  if (edgeIsInternal)
139  {
140  const edge& e = zoneLocalEdges[curEdgeID];
141 
142  // Reset the point creation
143  addedPoints[e.start()] = mp[e.start()];
144  addedPoints[e.end()] = mp[e.end()];
145  }
146  }
147 // Pout<< "addedPoints before point creation: " << addedPoints << endl;
148 
149  // Create new points for face zone
150  forAll(addedPoints, pointi)
151  {
152  if (addedPoints[pointi] < 0)
153  {
154  addedPoints[pointi] =
155  ref.setAction
156  (
157  polyAddPoint
158  (
159  points[mp[pointi]], // point
160  mp[pointi], // master point
161  -1, // zone ID
162  true // supports a cell
163  )
164  );
165  // Pout<< "Adding point " << addedPoints[pointi]
166  // << " coord1:" << points[mp[pointi]]
167  // << " coord2:" << masterFaceLayer.localPoints()[pointi]
168  // << " for original point " << mp[pointi] << endl;
169  }
170  }
171 
172  // Modify faces in the master zone and duplicate for the slave zone
173 
174  const labelList& mf = meshZones[faceZoneID_.index()];
175  const boolList& mfFlip = meshZones[faceZoneID_.index()].flipMap();
176  const faceList& zoneFaces = masterFaceLayer.localFaces();
177 
178  const faceList& faces = mesh.faces();
179  const labelList& own = mesh.faceOwner();
180  const labelList& nei = mesh.faceNeighbour();
181 
182  forAll(mf, facei)
183  {
184  const label curFaceID = mf[facei];
185 
186  // Build the face for the slave patch by renumbering
187  const face oldFace = zoneFaces[facei].reverseFace();
188 
189  face newFace(oldFace.size());
190 
191  forAll(oldFace, pointi)
192  {
193  newFace[pointi] = addedPoints[oldFace[pointi]];
194  }
195 
196  if (mfFlip[facei])
197  {
198  // Face needs to be flipped for the master patch
199  ref.setAction
200  (
201  polyModifyFace
202  (
203  faces[curFaceID].reverseFace(), // modified face
204  curFaceID, // label of face being modified
205  nei[curFaceID], // owner
206  -1, // neighbour
207  true, // face flip
208  masterPatchID_.index(), // patch for face
209  false, // remove from zone
210  faceZoneID_.index(), // zone for face
211  !mfFlip[facei] // face flip in zone
212  )
213  );
214 
215  // Add renumbered face into the slave patch
216  // label addedFacei =
217  ref.setAction
218  (
219  polyAddFace
220  (
221  newFace, // face
222  own[curFaceID], // owner
223  -1, // neighbour
224  -1, // master point
225  -1, // master edge
226  curFaceID, // master face
227  false, // flip flux
228  slavePatchID_.index(), // patch to add the face to
229  -1, // zone for face
230  false // zone flip
231  )
232  );
233  //{
234  // pointField newPts(ref.points());
235  // Pout<< "Flip. Modifying face: " << ref.faces()[curFaceID]
236  // << " fc:" << ref.faces()[curFaceID].centre(newPts)
237  // << " next to cell: " << nei[curFaceID]
238  // << " and adding face: " << newFace
239  // << " fc:" << ref.faces()[addedFacei].centre(newPts)
240  // << " next to cell: " << own[curFaceID] << endl;
241  //}
242  }
243  else
244  {
245  // No flip
246  ref.setAction
247  (
248  polyModifyFace
249  (
250  faces[curFaceID], // modified face
251  curFaceID, // label of face being modified
252  own[curFaceID], // owner
253  -1, // neighbour
254  false, // face flip
255  masterPatchID_.index(), // patch for face
256  false, // remove from zone
257  faceZoneID_.index(), // zone for face
258  mfFlip[facei] // face flip in zone
259  )
260  );
261 
262  // Add renumbered face into the slave patch
263  // label addedFacei =
264  ref.setAction
265  (
266  polyAddFace
267  (
268  newFace, // face
269  nei[curFaceID], // owner
270  -1, // neighbour
271  -1, // master point
272  -1, // master edge
273  curFaceID, // master face
274  true, // flip flux
275  slavePatchID_.index(), // patch to add the face to
276  -1, // zone for face
277  false // face flip in zone
278  )
279  );
280  //{
281  // pointField newPts(ref.points());
282  // Pout<< "No flip. Modifying face: " << ref.faces()[curFaceID]
283  // << " fc:" << ref.faces()[curFaceID].centre(newPts)
284  // << " next to cell: " << own[curFaceID]
285  // << " and adding face: " << newFace
286  // << " fc:" << ref.faces()[addedFacei].centre(newPts)
287  // << " next to cell: " << nei[curFaceID] << endl;
288  //}
289  }
290  }
291 
292  // Modify the remaining faces of the master cells to reconnect to the new
293  // layer of faces.
294 
295  // Algorithm: Go through all the cells of the master zone and make
296  // a map of faces to avoid duplicates. Do not insert the faces in
297  // the master patch (as they have already been dealt with). Make
298  // a master layer point renumbering map, which for every point in
299  // the master layer gives its new label. Loop through all faces in
300  // the map and attempt to renumber them using the master layer
301  // point renumbering map. Once the face is renumbered, compare it
302  // with the original face; if they are the same, the face has not
303  // changed; if not, modify the face but keep all of its old
304  // attributes (apart from the vertex numbers).
305 
306  // Create the map of faces in the master cell layer
307  const labelList& mc =
308  mesh.faceZones()[faceZoneID_.index()].masterCells();
309 
310  labelHashSet masterCellFaceMap(6*mc.size());
311 
312  const cellList& cells = mesh.cells();
313 
314  forAll(mc, celli)
315  {
316  const labelList& curFaces = cells[mc[celli]];
317 
318  forAll(curFaces, facei)
319  {
320  // Check if the face belongs to the master patch; if not add it
321  if (meshZones.whichZone(curFaces[facei]) != faceZoneID_.index())
322  {
323  masterCellFaceMap.insert(curFaces[facei]);
324  }
325  }
326  }
327 
328  // Extend the map to include first neighbours of the master cells to
329  // deal with multiple corners.
330  { // Protection and memory management
331  // Make a map of master cells for quick reject
332  labelHashSet mcMap(2*mc.size());
333 
334  forAll(mc, mcI)
335  {
336  mcMap.insert(mc[mcI]);
337  }
338 
339  // Go through all the faces in the masterCellFaceMap. If the
340  // cells around them are not already used, add all of their
341  // faces to the map
342  const labelList mcf = masterCellFaceMap.toc();
343 
344  forAll(mcf, mcfI)
345  {
346  // Do the owner side
347  const label ownCell = own[mcf[mcfI]];
348 
349  if (!mcMap.found(ownCell))
350  {
351  // Cell not found. Add its faces to the map
352  const cell& curFaces = cells[ownCell];
353 
354  forAll(curFaces, facei)
355  {
356  masterCellFaceMap.insert(curFaces[facei]);
357  }
358  }
359 
360  // Do the neighbour side if face is internal
361  if (mesh.isInternalFace(mcf[mcfI]))
362  {
363  const label neiCell = nei[mcf[mcfI]];
364 
365  if (!mcMap.found(neiCell))
366  {
367  // Cell not found. Add its faces to the map
368  const cell& curFaces = cells[neiCell];
369 
370  forAll(curFaces, facei)
371  {
372  masterCellFaceMap.insert(curFaces[facei]);
373  }
374  }
375  }
376  }
377  }
378 
379  // Create the master layer point map
380  Map<label> masterLayerPointMap(2*mp.size());
381 
382  forAll(mp, pointi)
383  {
384  masterLayerPointMap.insert
385  (
386  mp[pointi],
387  addedPoints[pointi]
388  );
389  }
390 
391  // Grab the list of faces of the master layer
392  const labelList masterCellFaces = masterCellFaceMap.toc();
393 
394  forAll(masterCellFaces, facei)
395  {
396  // Attempt to renumber the face using the masterLayerPointMap.
397  // Missing point remain the same
398 
399  const label curFaceID = masterCellFaces[facei];
400 
401  const face& oldFace = faces[curFaceID];
402 
403  face newFace(oldFace.size());
404 
405  bool changed = false;
406 
407  forAll(oldFace, pointi)
408  {
409  if (masterLayerPointMap.found(oldFace[pointi]))
410  {
411  changed = true;
412 
413  newFace[pointi] = masterLayerPointMap.find(oldFace[pointi])();
414  }
415  else
416  {
417  newFace[pointi] = oldFace[pointi];
418  }
419  }
420 
421  // If the face has changed, create a modification entry
422  if (changed)
423  {
424  if (mesh.isInternalFace(curFaceID))
425  {
426  ref.setAction
427  (
428  polyModifyFace
429  (
430  newFace, // face
431  curFaceID, // master face
432  own[curFaceID], // owner
433  nei[curFaceID], // neighbour
434  false, // flip flux
435  -1, // patch for face
436  false, // remove from zone
437  -1, // zone for face
438  false // face zone flip
439  )
440  );
441 
442  // Pout<< "modifying stick-out face. Internal Old face: "
443  // << oldFace
444  // << " new face: " << newFace
445  // << " own: " << own[curFaceID]
446  // << " nei: " << nei[curFaceID]
447  // << endl;
448  }
449  else
450  {
451  ref.setAction
452  (
453  polyModifyFace
454  (
455  newFace, // face
456  curFaceID, // master face
457  own[curFaceID], // owner
458  -1, // neighbour
459  false, // flip flux
460  mesh.boundaryMesh().whichPatch(curFaceID), // patch
461  false, // remove from zone
462  -1, // zone for face
463  false // face zone flip
464  )
465  );
466 
467  // Pout<< "modifying stick-out face. Boundary Old face: "
468  // << oldFace
469  // << " new face: " << newFace
470  // << " own: " << own[curFaceID]
471  // << " patch: "
472  // << mesh.boundaryMesh().whichPatch(curFaceID)
473  // << endl;
474  }
475  }
476  }
477 
478  if (debug)
479  {
480  Pout<< "void attachDetach::detachInterface("
481  << "polyTopoChange& ref) const "
482  << " for object " << name() << " : "
483  << "Finished detaching interface" << endl;
484  }
485 }
486 
487 
488 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
const word & name() const
Return name of this modifier.
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
const polyMesh & mesh() const
Return the mesh reference.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const pointField & points
const cellShapeList & cells
const dimensionedScalar mp
Proton mass.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const doubleScalar e
Definition: doubleScalar.H:105
List< label > labelList
A List of labels.
Definition: labelList.H:56
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
List< cell > cellList
list of cells
Definition: cellList.H:42
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:211
PrimitivePatch< List< face >, const pointField & > primitiveFacePatch
Foam::primitiveFacePatch.
error FatalError
List< edge > edgeList
Definition: edgeList.H:38
MeshZones< faceZone, polyMesh > meshFaceZones
A MeshZones with the type faceZone.
List< face > faceList
Definition: faceListFwd.H:43
trDeltaTf ref()