setLayerPairing.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-2018 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  << "Problem with layer pairing data"
81  << abort(FatalError);
82  }
83 
84  pointsPairingPtr_ = new labelList(meshPoints.size(), -1);
85  labelList& ptc = *pointsPairingPtr_;
86 
87  facesPairingPtr_ = new labelList(mf.size(), -1);
88  labelList& ftc = *facesPairingPtr_;
89 
90  if (debug > 1)
91  {
92  Pout<< "meshPoints: " << meshPoints << nl
93  << "localPoints: "
94  << mesh.faceZones()[faceZoneID_.index()]().localPoints()
95  << endl;
96  }
97 
98  // For all faces, create the mapping
99  label nPointErrors = 0;
100  label nFaceErrors = 0;
101 
102  forAll(mf, facei)
103  {
104  // Get the local master face
105  face curLocalFace = mlf[facei];
106 
107  // Flip face based on flip index to recover original orientation
108  if (mfFlip[facei])
109  {
110  curLocalFace.flip();
111  }
112 
113  // Get the opposing face from the master cell
114  oppositeFace lidFace =
115  cells[mc[facei]].opposingFace(mf[facei], faces);
116 
117  if (!lidFace.found())
118  {
119  // This is not a valid layer; cannot continue
120  nFaceErrors++;
121  continue;
122  }
123 
124  if (debug > 1)
125  {
126  Pout<< "curMasterFace: " << faces[mf[facei]] << nl
127  << "cell shape: " << mesh.cellShapes()[mc[facei]] << nl
128  << "curLocalFace: " << curLocalFace << nl
129  << "lidFace: " << lidFace
130  << " master index: " << lidFace.masterIndex()
131  << " oppositeIndex: " << lidFace.oppositeIndex() << endl;
132  }
133 
134  // Grab the opposite face for face collapse addressing
135  ftc[facei] = lidFace.oppositeIndex();
136 
137  // Using the local face insert the points into the lid list
138  forAll(curLocalFace, pointi)
139  {
140  const label clp = curLocalFace[pointi];
141 
142  if (ptc[clp] == -1)
143  {
144  // Point not mapped yet. Insert the label
145  ptc[clp] = lidFace[pointi];
146  }
147  else
148  {
149  // Point mapped from some other face. Check the label
150  if (ptc[clp] != lidFace[pointi])
151  {
152  nPointErrors++;
153 
154  if (debug > 1)
155  {
156  Pout<< "Topological error in cell layer pairing. "
157  << "This mesh is either topologically incorrect "
158  << "or the master face layer is not defined "
159  << "consistently. Please check the "
160  << "face zone flip map." << nl
161  << "First index: " << ptc[clp]
162  << " new index: " << lidFace[pointi] << endl;
163  }
164  }
165  }
166  }
167  }
168 
169  reduce(nPointErrors, sumOp<label>());
170  reduce(nFaceErrors, sumOp<label>());
171 
172  if (nPointErrors > 0 || nFaceErrors > 0)
173  {
174  clearAddressing();
175 
176  return false;
177  }
178  else
179  {
180  // Valid layer
181  return true;
182  }
183 }
184 
185 
186 const Foam::labelList& Foam::layerAdditionRemoval::pointsPairing() const
187 {
188  if (!pointsPairingPtr_)
189  {
191  << "Problem with layer pairing data for object " << name()
192  << abort(FatalError);
193  }
194 
195  return *pointsPairingPtr_;
196 }
197 
198 const Foam::labelList& Foam::layerAdditionRemoval::facesPairing() const
199 {
200  if (!facesPairingPtr_)
201  {
203  << "Problem with layer pairing data for object " << name()
204  << abort(FatalError);
205  }
206 
207  return *facesPairingPtr_;
208 }
209 
210 
211 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
212 
214 (
215  pointField& motionPoints
216 ) const
217 {
218  if (debug)
219  {
220  Pout<< "void layerAdditionRemoval::modifyMotionPoints("
221  << "pointField& motionPoints) const for object "
222  << name() << " : ";
223  }
224 
225  if (debug)
226  {
227  Pout<< "No motion point adjustment" << endl;
228  }
229 }
230 
231 
232 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
virtual void modifyMotionPoints(pointField &motionPoints) const
Modify motion points to comply with the topological change.
List< label > labelList
A List of labels.
Definition: labelList.H:56
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label index() const
Return index of first matching zone.
Definition: DynamicID.H:107
static const char nl
Definition: Ostream.H:265
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const polyTopoChanger & topoChanger() const
Return reference to morph engine.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
const word & name() const
Return name of this modifier.
const polyMesh & mesh() const
Return the mesh reference.
List< cell > cellList
list of cells
Definition: cellList.H:42