setLayerPairing.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-2014 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 Description
25  Remove a layer of cells and prepare addressing data
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "layerAdditionRemoval.H"
30 #include "polyMesh.H"
31 #include "primitiveMesh.H"
32 #include "polyTopoChange.H"
33 #include "oppositeFace.H"
34 #include "polyTopoChanger.H"
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 bool Foam::layerAdditionRemoval::setLayerPairing() const
39 {
40  // Note:
41  // This is also the most complex part of the topological change.
42  // Therefore it will be calculated here and stored as temporary
43  // data until the actual topological change, after which it will
44  // be cleared.
45 
46  // Algorithm for point collapse
47  // 1) Go through the master cell layer and for every face of
48  // the face zone find the opposite face in the master cell.
49  // Check the direction of the opposite face and adjust as
50  // necessary. Check other faces to find an edge defining
51  // relative orientation of the two faces and adjust the face
52  // as necessary. Once the face is adjusted, record the
53  // addressing between the master and slave vertex layer.
54 
55  const polyMesh& mesh = topoChanger().mesh();
56 
57  const labelList& mc =
58  mesh.faceZones()[faceZoneID_.index()].masterCells();
59 
60  const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
61 
62  const boolList& mfFlip =
63  mesh.faceZones()[faceZoneID_.index()].flipMap();
64 
65  const faceList& faces = mesh.faces();
66  const cellList& cells = mesh.cells();
67 
68  // Grab the local faces from the master zone
69  const faceList& mlf =
70  mesh.faceZones()[faceZoneID_.index()]().localFaces();
71 
72  const labelList& meshPoints =
73  mesh.faceZones()[faceZoneID_.index()]().meshPoints();
74 
75  // Create a list of points to collapse for every point of
76  // the master patch
77  if (pointsPairingPtr_ || facesPairingPtr_)
78  {
80  (
81  "void Foam::layerAdditionRemoval::setLayerPairing() const"
82  ) << "Problem with layer pairing data"
83  << abort(FatalError);
84  }
85 
86  pointsPairingPtr_ = new labelList(meshPoints.size(), -1);
87  labelList& ptc = *pointsPairingPtr_;
88 
89  facesPairingPtr_ = new labelList(mf.size(), -1);
90  labelList& ftc = *facesPairingPtr_;
91 
92  if (debug > 1)
93  {
94  Pout<< "meshPoints: " << meshPoints << nl
95  << "localPoints: "
96  << mesh.faceZones()[faceZoneID_.index()]().localPoints()
97  << endl;
98  }
99 
100  // For all faces, create the mapping
101  label nPointErrors = 0;
102  label nFaceErrors = 0;
103 
104  forAll(mf, faceI)
105  {
106  // Get the local master face
107  face curLocalFace = mlf[faceI];
108 
109  // Flip face based on flip index to recover original orientation
110  if (mfFlip[faceI])
111  {
112  curLocalFace.flip();
113  }
114 
115  // Get the opposing face from the master cell
116  oppositeFace lidFace =
117  cells[mc[faceI]].opposingFace(mf[faceI], faces);
118 
119  if (!lidFace.found())
120  {
121  // This is not a valid layer; cannot continue
122  nFaceErrors++;
123  continue;
124  }
125 
126  if (debug > 1)
127  {
128  Pout<< "curMasterFace: " << faces[mf[faceI]] << nl
129  << "cell shape: " << mesh.cellShapes()[mc[faceI]] << nl
130  << "curLocalFace: " << curLocalFace << nl
131  << "lidFace: " << lidFace
132  << " master index: " << lidFace.masterIndex()
133  << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
134  }
135 
136  // Grab the opposite face for face collapse addressing
137  ftc[faceI] = lidFace.oppositeIndex();
138 
139  // Using the local face insert the points into the lid list
140  forAll(curLocalFace, pointI)
141  {
142  const label clp = curLocalFace[pointI];
143 
144  if (ptc[clp] == -1)
145  {
146  // Point not mapped yet. Insert the label
147  ptc[clp] = lidFace[pointI];
148  }
149  else
150  {
151  // Point mapped from some other face. Check the label
152  if (ptc[clp] != lidFace[pointI])
153  {
154  nPointErrors++;
155 
156  if (debug > 1)
157  {
158  Pout<< "Topological error in cell layer pairing. "
159  << "This mesh is either topologically incorrect "
160  << "or the master face layer is not defined "
161  << "consistently. Please check the "
162  << "face zone flip map." << nl
163  << "First index: " << ptc[clp]
164  << " new index: " << lidFace[pointI] << endl;
165  }
166  }
167  }
168  }
169  }
170 
171  reduce(nPointErrors, sumOp<label>());
172  reduce(nFaceErrors, sumOp<label>());
173 
174  if (nPointErrors > 0 || nFaceErrors > 0)
175  {
176  clearAddressing();
177 
178  return false;
179  }
180  else
181  {
182  // Valid layer
183  return true;
184  }
185 }
186 
187 
188 const Foam::labelList& Foam::layerAdditionRemoval::pointsPairing() const
189 {
190  if (!pointsPairingPtr_)
191  {
193  (
194  "const labelList& layerAdditionRemoval::pointsPairing() const"
195  ) << "Problem with layer pairing data for object " << name()
196  << abort(FatalError);
197  }
198 
199  return *pointsPairingPtr_;
200 }
201 
202 const Foam::labelList& Foam::layerAdditionRemoval::facesPairing() const
203 {
204  if (!facesPairingPtr_)
205  {
207  (
208  "const labelList& layerAdditionRemoval::facesPairing() const"
209  ) << "Problem with layer pairing data for object " << name()
210  << abort(FatalError);
211  }
212 
213  return *facesPairingPtr_;
214 }
215 
216 
217 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
218 
220 (
221  pointField& motionPoints
222 ) const
223 {
224  if (debug)
225  {
226  Pout<< "void layerAdditionRemoval::modifyMotionPoints("
227  << "pointField& motionPoints) const for object "
228  << name() << " : ";
229  }
230 
231  if (debug)
232  {
233  Pout<< "No motion point adjustment" << endl;
234  }
235 }
236 
237 
238 // ************************************************************************* //
virtual void modifyMotionPoints(pointField &motionPoints) const
Modify motion points to comply with the topological change.
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
errorManip< error > abort(error &err)
Definition: errorManip.H:132
word name(const complex &)
Return a string representation of a complex.
List< label > labelList
A List of labels.
Definition: labelList.H:56
#define forAll(list, i)
Definition: UList.H:421
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:249
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
List< cell > cellList
list of cells
Definition: cellList.H:42
error FatalError
Definition: error.H:303
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
static const char nl
Definition: Ostream.H:258
#define FatalErrorIn(fn)
Definition: error.H:318
List< face > faceList
Definition: faceListFwd.H:43
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
const polyMesh & mesh() const
Return the mesh reference.