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