detachInterface.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-2015 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 
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 faceZoneMesh& zoneMesh = 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 = zoneMesh[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 " << zoneMesh[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 = zoneMesh[faceZoneID_.index()]();
104  const pointField& points = mesh.points();
105  const labelListList& meshEdgeFaces = mesh.edgeFaces();
106 
107  const labelList& mp = masterFaceLayer.meshPoints();
108  const edgeList& zoneLocalEdges = masterFaceLayer.edges();
109 
110  const labelList& meshEdges = zoneMesh[faceZoneID_.index()].meshEdges();
111 
112  // Create the points
113 
114  labelList addedPoints(mp.size(), -1);
115 
116  // Go through boundary edges of the master patch. If all the faces from
117  // this patch are internal, mark the points in the addedPoints lookup
118  // with their original labels to stop duplication
119  label nIntEdges = masterFaceLayer.nInternalEdges();
120 
121  for (label curEdgeID = nIntEdges; curEdgeID < meshEdges.size(); curEdgeID++)
122  {
123  const labelList& curFaces = meshEdgeFaces[meshEdges[curEdgeID]];
124 
125  bool edgeIsInternal = true;
126 
127  forAll(curFaces, faceI)
128  {
129  if (!mesh.isInternalFace(curFaces[faceI]))
130  {
131  // The edge belongs to a boundary face
132  edgeIsInternal = false;
133  break;
134  }
135  }
136 
137  if (edgeIsInternal)
138  {
139  const edge& e = zoneLocalEdges[curEdgeID];
140 
141  // Reset the point creation
142  addedPoints[e.start()] = mp[e.start()];
143  addedPoints[e.end()] = mp[e.end()];
144  }
145  }
146 // Pout<< "addedPoints before point creation: " << addedPoints << endl;
147 
148  // Create new points for face zone
149  forAll(addedPoints, pointI)
150  {
151  if (addedPoints[pointI] < 0)
152  {
153  addedPoints[pointI] =
154  ref.setAction
155  (
157  (
158  points[mp[pointI]], // point
159  mp[pointI], // master point
160  -1, // zone ID
161  true // supports a cell
162  )
163  );
164  //Pout<< "Adding point " << addedPoints[pointI]
165  // << " coord1:" << points[mp[pointI]]
166  // << " coord2:" << masterFaceLayer.localPoints()[pointI]
167  // << " for original point " << mp[pointI] << endl;
168  }
169  }
170 
171  // Modify faces in the master zone and duplicate for the slave zone
172 
173  const labelList& mf = zoneMesh[faceZoneID_.index()];
174  const boolList& mfFlip = zoneMesh[faceZoneID_.index()].flipMap();
175  const faceList& zoneFaces = masterFaceLayer.localFaces();
176 
177  const faceList& faces = mesh.faces();
178  const labelList& own = mesh.faceOwner();
179  const labelList& nei = mesh.faceNeighbour();
180 
181  forAll(mf, faceI)
182  {
183  const label curFaceID = mf[faceI];
184 
185  // Build the face for the slave patch by renumbering
186  const face oldFace = zoneFaces[faceI].reverseFace();
187 
188  face newFace(oldFace.size());
189 
190  forAll(oldFace, pointI)
191  {
192  newFace[pointI] = addedPoints[oldFace[pointI]];
193  }
194 
195  if (mfFlip[faceI])
196  {
197  // Face needs to be flipped for the master patch
198  ref.setAction
199  (
201  (
202  faces[curFaceID].reverseFace(), // modified face
203  curFaceID, // label of face being modified
204  nei[curFaceID], // owner
205  -1, // neighbour
206  true, // face flip
207  masterPatchID_.index(), // patch for face
208  false, // remove from zone
209  faceZoneID_.index(), // zone for face
210  !mfFlip[faceI] // face flip in zone
211  )
212  );
213 
214  // Add renumbered face into the slave patch
215  //label addedFaceI =
216  ref.setAction
217  (
219  (
220  newFace, // face
221  own[curFaceID], // owner
222  -1, // neighbour
223  -1, // master point
224  -1, // master edge
225  curFaceID, // master face
226  false, // flip flux
227  slavePatchID_.index(), // patch to add the face to
228  -1, // zone for face
229  false // zone flip
230  )
231  );
232  //{
233  // pointField newPts(ref.points());
234  //Pout<< "Flip. Modifying face: " << ref.faces()[curFaceID]
235  // << " fc:" << ref.faces()[curFaceID].centre(newPts)
236  // << " next to cell: " << nei[curFaceID]
237  // << " and adding face: " << newFace
238  // << " fc:" << ref.faces()[addedFaceI].centre(newPts)
239  // << " next to cell: " << own[curFaceID] << endl;
240  //}
241  }
242  else
243  {
244  // No flip
245  ref.setAction
246  (
248  (
249  faces[curFaceID], // modified face
250  curFaceID, // label of face being modified
251  own[curFaceID], // owner
252  -1, // neighbour
253  false, // face flip
254  masterPatchID_.index(), // patch for face
255  false, // remove from zone
256  faceZoneID_.index(), // zone for face
257  mfFlip[faceI] // face flip in zone
258  )
259  );
260 
261  // Add renumbered face into the slave patch
262  //label addedFaceI =
263  ref.setAction
264  (
266  (
267  newFace, // face
268  nei[curFaceID], // owner
269  -1, // neighbour
270  -1, // master point
271  -1, // master edge
272  curFaceID, // master face
273  true, // flip flux
274  slavePatchID_.index(), // patch to add the face to
275  -1, // zone for face
276  false // face flip in zone
277  )
278  );
279  //{
280  // pointField newPts(ref.points());
281  //Pout<< "No flip. Modifying face: " << ref.faces()[curFaceID]
282  // << " fc:" << ref.faces()[curFaceID].centre(newPts)
283  // << " next to cell: " << own[curFaceID]
284  // << " and adding face: " << newFace
285  // << " fc:" << ref.faces()[addedFaceI].centre(newPts)
286  // << " next to cell: " << nei[curFaceID] << endl;
287  //}
288  }
289  }
290 
291  // Modify the remaining faces of the master cells to reconnect to the new
292  // layer of faces.
293 
294  // Algorithm: Go through all the cells of the master zone and make
295  // a map of faces to avoid duplicates. Do not insert the faces in
296  // the master patch (as they have already been dealt with). Make
297  // a master layer point renumbering map, which for every point in
298  // the master layer gives its new label. Loop through all faces in
299  // the map and attempt to renumber them using the master layer
300  // point renumbering map. Once the face is renumbered, compare it
301  // with the original face; if they are the same, the face has not
302  // changed; if not, modify the face but keep all of its old
303  // attributes (apart from the vertex numbers).
304 
305  // Create the map of faces in the master cell layer
306  const labelList& mc =
307  mesh.faceZones()[faceZoneID_.index()].masterCells();
308 
309  labelHashSet masterCellFaceMap(6*mc.size());
310 
311  const cellList& cells = mesh.cells();
312 
313  forAll(mc, cellI)
314  {
315  const labelList& curFaces = cells[mc[cellI]];
316 
317  forAll(curFaces, faceI)
318  {
319  // Check if the face belongs to the master patch; if not add it
320  if (zoneMesh.whichZone(curFaces[faceI]) != faceZoneID_.index())
321  {
322  masterCellFaceMap.insert(curFaces[faceI]);
323  }
324  }
325  }
326 
327  // Extend the map to include first neighbours of the master cells to
328  // deal with multiple corners.
329  { // Protection and memory management
330  // Make a map of master cells for quick reject
331  labelHashSet mcMap(2*mc.size());
332 
333  forAll(mc, mcI)
334  {
335  mcMap.insert(mc[mcI]);
336  }
337 
338  // Go through all the faces in the masterCellFaceMap. If the
339  // cells around them are not already used, add all of their
340  // faces to the map
341  const labelList mcf = masterCellFaceMap.toc();
342 
343  forAll(mcf, mcfI)
344  {
345  // Do the owner side
346  const label ownCell = own[mcf[mcfI]];
347 
348  if (!mcMap.found(ownCell))
349  {
350  // Cell not found. Add its faces to the map
351  const cell& curFaces = cells[ownCell];
352 
353  forAll(curFaces, faceI)
354  {
355  masterCellFaceMap.insert(curFaces[faceI]);
356  }
357  }
358 
359  // Do the neighbour side if face is internal
360  if (mesh.isInternalFace(mcf[mcfI]))
361  {
362  const label neiCell = nei[mcf[mcfI]];
363 
364  if (!mcMap.found(neiCell))
365  {
366  // Cell not found. Add its faces to the map
367  const cell& curFaces = cells[neiCell];
368 
369  forAll(curFaces, faceI)
370  {
371  masterCellFaceMap.insert(curFaces[faceI]);
372  }
373  }
374  }
375  }
376  }
377 
378  // Create the master layer point map
379  Map<label> masterLayerPointMap(2*mp.size());
380 
381  forAll(mp, pointI)
382  {
383  masterLayerPointMap.insert
384  (
385  mp[pointI],
386  addedPoints[pointI]
387  );
388  }
389 
390  // Grab the list of faces of the master layer
391  const labelList masterCellFaces = masterCellFaceMap.toc();
392 
393  forAll(masterCellFaces, faceI)
394  {
395  // Attempt to renumber the face using the masterLayerPointMap.
396  // Missing point remain the same
397 
398  const label curFaceID = masterCellFaces[faceI];
399 
400  const face& oldFace = faces[curFaceID];
401 
402  face newFace(oldFace.size());
403 
404  bool changed = false;
405 
406  forAll(oldFace, pointI)
407  {
408  if (masterLayerPointMap.found(oldFace[pointI]))
409  {
410  changed = true;
411 
412  newFace[pointI] = masterLayerPointMap.find(oldFace[pointI])();
413  }
414  else
415  {
416  newFace[pointI] = oldFace[pointI];
417  }
418  }
419 
420  // If the face has changed, create a modification entry
421  if (changed)
422  {
423  if (mesh.isInternalFace(curFaceID))
424  {
425  ref.setAction
426  (
428  (
429  newFace, // face
430  curFaceID, // master face
431  own[curFaceID], // owner
432  nei[curFaceID], // neighbour
433  false, // flip flux
434  -1, // patch for face
435  false, // remove from zone
436  -1, // zone for face
437  false // face zone flip
438  )
439  );
440 
441  // Pout<< "modifying stick-out face. Internal Old face: "
442  // << oldFace
443  // << " new face: " << newFace
444  // << " own: " << own[curFaceID]
445  // << " nei: " << nei[curFaceID]
446  // << endl;
447  }
448  else
449  {
450  ref.setAction
451  (
453  (
454  newFace, // face
455  curFaceID, // master face
456  own[curFaceID], // owner
457  -1, // neighbour
458  false, // flip flux
459  mesh.boundaryMesh().whichPatch(curFaceID), // patch
460  false, // remove from zone
461  -1, // zone for face
462  false // face zone flip
463  )
464  );
465 
466  // Pout<< "modifying stick-out face. Boundary Old face: "
467  // << oldFace
468  // << " new face: " << newFace
469  // << " own: " << own[curFaceID]
470  // << " patch: "
471  // << mesh.boundaryMesh().whichPatch(curFaceID)
472  // << endl;
473  }
474  }
475  }
476 
477  if (debug)
478  {
479  Pout<< "void attachDetach::detachInterface("
480  << "polyTopoChange& ref) const "
481  << " for object " << name() << " : "
482  << "Finished detaching interface" << endl;
483  }
484 }
485 
486 
487 // ************************************************************************* //
Foam::constant::atomic::mp
const dimensionedScalar mp
Proton mass.
Foam::PrimitivePatch::edges
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
Definition: PrimitivePatchTemplate.C:212
Foam::HashTable::toc
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:201
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
polyAddFace.H
polyTopoChanger.H
polyTopoChange.H
Foam::polyTopoChange
Direct mesh changes based on v1.3 polyTopoChange syntax.
Definition: polyTopoChange.H:97
Foam::edge
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
Foam::Map< label >
Foam::polyModifyFace
Class describing modification of a face.
Definition: polyModifyFace.H:48
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
polyMesh.H
Foam::HashSet< label, Hash< label > >
Foam::polyMesh
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
Foam::label
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
Foam::Field
Pre-declare SubField and related Field type.
Definition: Field.H:57
Foam::ZoneMesh< faceZone, polyMesh >
polyAddPoint.H
Foam::ZoneMesh::whichZone
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:226
Foam::attachDetach::detachInterface
void detachInterface(polyTopoChange &) const
Detach interface.
Definition: detachInterface.C:38
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:18
Foam::HashTable::found
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:109
polyModifyFace.H
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Foam::PrimitivePatch::nInternalEdges
label nInternalEdges() const
Number of internal edges.
Definition: PrimitivePatchTemplate.C:232
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:318
Foam::polyTopoChange::setAction
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
Definition: polyTopoChange.C:2530
Foam::Pout
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:94
attachDetach.H
Foam::HashSet::insert
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116
Foam::polyAddPoint
Class containing data for point addition.
Definition: polyAddPoint.H:47
Foam::polyAddFace
A face addition data class. A face can be inflated either from a point or from another face and can e...
Definition: polyAddFace.H:51
Foam::PrimitivePatch::localFaces
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
Definition: PrimitivePatchTemplate.C:372
Foam::face
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::List::size
void size(const label)
Override size to be inconsistent with allocated storage.
Foam::PrimitivePatch::meshPoints
const labelList & meshPoints() const
Return labelList of mesh points in patch.
Definition: PrimitivePatchTemplate.C:392
Foam::cell
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
Foam::help::reverseFace
faceType reverseFace(const faceType &f)
reverse the face
Definition: helperFunctionsTopologyManipulationI.H:113
Foam::name
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatchTemplate.H:88