mirrorFvMesh.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 "mirrorFvMesh.H"
27 #include "Time.H"
28 #include "plane.H"
29 
30 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
31 
32 Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io)
33 :
34  fvMesh(io),
35  mirrorMeshDict_
36  (
37  IOobject
38  (
39  "mirrorMeshDict",
40  time().system(),
41  *this,
42  IOobject::MUST_READ_IF_MODIFIED,
43  IOobject::NO_WRITE
44  )
45  ),
46  mirrorMeshPtr_(NULL)
47 {
48  plane mirrorPlane(mirrorMeshDict_);
49 
50  scalar planeTolerance
51  (
52  readScalar(mirrorMeshDict_.lookup("planeTolerance"))
53  );
54 
55  const pointField& oldPoints = points();
56  const faceList& oldFaces = faces();
57  const cellList& oldCells = cells();
58  const label nOldInternalFaces = nInternalFaces();
59  const polyPatchList& oldPatches = boundaryMesh();
60 
61  // Mirror the points
62  Info<< "Mirroring points. Old points: " << oldPoints.size();
63 
64  pointField newPoints(2*oldPoints.size());
65  label nNewPoints = 0;
66 
67  labelList mirrorPointLookup(oldPoints.size(), -1);
68 
69  // Grab the old points
70  forAll(oldPoints, pointI)
71  {
72  newPoints[nNewPoints] = oldPoints[pointI];
73  nNewPoints++;
74  }
75 
76  forAll(oldPoints, pointI)
77  {
78  scalar alpha =
79  mirrorPlane.normalIntersect
80  (
81  oldPoints[pointI],
82  mirrorPlane.normal()
83  );
84 
85  // Check plane on tolerance
86  if (mag(alpha) > planeTolerance)
87  {
88  // The point gets mirrored
89  newPoints[nNewPoints] =
90  oldPoints[pointI] + 2.0*alpha*mirrorPlane.normal();
91 
92  // remember the point correspondence
93  mirrorPointLookup[pointI] = nNewPoints;
94  nNewPoints++;
95  }
96  else
97  {
98  // The point is on the plane and does not get mirrored
99  // Adjust plane location
100  newPoints[nNewPoints] =
101  oldPoints[pointI] + alpha*mirrorPlane.normal();
102 
103  mirrorPointLookup[pointI] = pointI;
104  }
105  }
106 
107  // Reset the size of the point list
108  Info<< " New points: " << nNewPoints << endl;
109  newPoints.setSize(nNewPoints);
110 
111  Info<< "Mirroring faces. Old faces: " << oldFaces.size();
112 
113  // Algorithm:
114  // During mirroring, the faces that were previously boundary faces
115  // in the mirror plane may become ineternal faces. In order to
116  // deal with the ordering of the faces, the algorithm is split
117  // into two parts. For original faces, the internal faces are
118  // distributed to their owner cells. Once all internal faces are
119  // distributed, the boundary faces are visited and if they are in
120  // the mirror plane they are added to the master cells (the future
121  // boundary faces are not touched). After the first phase, the
122  // internal faces are collected in the cell order and numbering
123  // information is added. Then, the internal faces are mirrored
124  // and the face numbering data is stored for the mirrored section.
125  // Once all the internal faces are mirrored, the boundary faces
126  // are added by mirroring the faces patch by patch.
127 
128  // Distribute internal faces
129  labelListList newCellFaces(oldCells.size());
130 
131  const labelUList& oldOwnerStart = lduAddr().ownerStartAddr();
132 
133  forAll(newCellFaces, cellI)
134  {
135  labelList& curFaces = newCellFaces[cellI];
136 
137  const label s = oldOwnerStart[cellI];
138  const label e = oldOwnerStart[cellI + 1];
139 
140  curFaces.setSize(e - s);
141 
142  forAll(curFaces, i)
143  {
144  curFaces[i] = s + i;
145  }
146  }
147 
148  // Distribute boundary faces. Remember the faces that have been inserted
149  // as internal
150  boolListList insertedBouFace(oldPatches.size());
151 
152  forAll(oldPatches, patchI)
153  {
154  const polyPatch& curPatch = oldPatches[patchI];
155 
156  if (curPatch.coupled())
157  {
159  << "Found coupled patch " << curPatch.name() << endl
160  << " Mirroring faces on coupled patches destroys"
161  << " the ordering. This might be fixed by running a dummy"
162  << " createPatch afterwards." << endl;
163  }
164 
165  boolList& curInsBouFace = insertedBouFace[patchI];
166 
167  curInsBouFace.setSize(curPatch.size());
168  curInsBouFace = false;
169 
170  // Get faceCells for face insertion
171  const labelUList& curFaceCells = curPatch.faceCells();
172  const label curStart = curPatch.start();
173 
174  forAll(curPatch, faceI)
175  {
176  // Find out if the mirrored face is identical to the
177  // original. If so, the face needs to become internal and
178  // added to its owner cell
179  const face& origFace = curPatch[faceI];
180 
181  face mirrorFace(origFace.size());
182  forAll(mirrorFace, pointI)
183  {
184  mirrorFace[pointI] = mirrorPointLookup[origFace[pointI]];
185  }
186 
187  if (origFace == mirrorFace)
188  {
189  // The mirror is identical to current face. This will
190  // become an internal face
191  const label oldSize = newCellFaces[curFaceCells[faceI]].size();
192 
193  newCellFaces[curFaceCells[faceI]].setSize(oldSize + 1);
194  newCellFaces[curFaceCells[faceI]][oldSize] = curStart + faceI;
195 
196  curInsBouFace[faceI] = true;
197  }
198  }
199  }
200 
201  // Construct the new list of faces. Boundary faces are added
202  // last, cush that each patch is mirrored separately. The
203  // addressing is stored in two separate arrays: first for the
204  // original cells (face order has changed) and then for the
205  // mirrored cells.
206  labelList masterFaceLookup(oldFaces.size(), -1);
207  labelList mirrorFaceLookup(oldFaces.size(), -1);
208 
209  faceList newFaces(2*oldFaces.size());
210  label nNewFaces = 0;
211 
212  // Insert original (internal) faces
213  forAll(newCellFaces, cellI)
214  {
215  const labelList& curCellFaces = newCellFaces[cellI];
216 
217  forAll(curCellFaces, cfI)
218  {
219  newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
220  masterFaceLookup[curCellFaces[cfI]] = nNewFaces;
221 
222  nNewFaces++;
223  }
224  }
225 
226  // Mirror internal faces
227  for (label faceI = 0; faceI < nOldInternalFaces; faceI++)
228  {
229  const face& oldFace = oldFaces[faceI];
230  face& nf = newFaces[nNewFaces];
231  nf.setSize(oldFace.size());
232 
233  nf[0] = mirrorPointLookup[oldFace[0]];
234 
235  for (label i = 1; i < oldFace.size(); i++)
236  {
237  nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
238  }
239 
240  mirrorFaceLookup[faceI] = nNewFaces;
241  nNewFaces++;
242  }
243 
244  // Mirror boundary faces patch by patch
245  labelList newToOldPatch(boundary().size(), -1);
246  labelList newPatchSizes(boundary().size(), -1);
247  labelList newPatchStarts(boundary().size(), -1);
248  label nNewPatches = 0;
249 
250  forAll(boundaryMesh(), patchI)
251  {
252  const label curPatchSize = boundaryMesh()[patchI].size();
253  const label curPatchStart = boundaryMesh()[patchI].start();
254  const boolList& curInserted = insertedBouFace[patchI];
255 
256  newPatchStarts[nNewPatches] = nNewFaces;
257 
258  // Master side
259  for (label faceI = 0; faceI < curPatchSize; faceI++)
260  {
261  // Check if the face has already been added. If not, add it and
262  // insert the numbering details.
263  if (!curInserted[faceI])
264  {
265  newFaces[nNewFaces] = oldFaces[curPatchStart + faceI];
266 
267  masterFaceLookup[curPatchStart + faceI] = nNewFaces;
268  nNewFaces++;
269  }
270  }
271 
272  // Mirror side
273  for (label faceI = 0; faceI < curPatchSize; faceI++)
274  {
275  // Check if the face has already been added. If not, add it and
276  // insert the numbering details.
277  if (!curInserted[faceI])
278  {
279  const face& oldFace = oldFaces[curPatchStart + faceI];
280  face& nf = newFaces[nNewFaces];
281  nf.setSize(oldFace.size());
282 
283  nf[0] = mirrorPointLookup[oldFace[0]];
284 
285  for (label i = 1; i < oldFace.size(); i++)
286  {
287  nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
288  }
289 
290  mirrorFaceLookup[curPatchStart + faceI] = nNewFaces;
291  nNewFaces++;
292  }
293  else
294  {
295  // Grab the index of the master face for the mirror side
296  mirrorFaceLookup[curPatchStart + faceI] =
297  masterFaceLookup[curPatchStart + faceI];
298  }
299  }
300 
301  // If patch exists, grab the name and type of the original patch
302  if (nNewFaces > newPatchStarts[nNewPatches])
303  {
304  newToOldPatch[nNewPatches] = patchI;
305 
306  newPatchSizes[nNewPatches] =
307  nNewFaces - newPatchStarts[nNewPatches];
308 
309  nNewPatches++;
310  }
311  }
312 
313  // Tidy up the lists
314  newFaces.setSize(nNewFaces);
315  Info<< " New faces: " << nNewFaces << endl;
316 
317  newToOldPatch.setSize(nNewPatches);
318  newPatchSizes.setSize(nNewPatches);
319  newPatchStarts.setSize(nNewPatches);
320 
321  Info<< "Mirroring patches. Old patches: " << boundary().size()
322  << " New patches: " << nNewPatches << endl;
323 
324  Info<< "Mirroring cells. Old cells: " << oldCells.size()
325  << " New cells: " << 2*oldCells.size() << endl;
326 
327  cellList newCells(2*oldCells.size());
328  label nNewCells = 0;
329 
330  // Grab the original cells. Take care of face renumbering.
331  forAll(oldCells, cellI)
332  {
333  const cell& oc = oldCells[cellI];
334 
335  cell& nc = newCells[nNewCells];
336  nc.setSize(oc.size());
337 
338  forAll(oc, i)
339  {
340  nc[i] = masterFaceLookup[oc[i]];
341  }
342 
343  nNewCells++;
344  }
345 
346  // Mirror the cells
347  forAll(oldCells, cellI)
348  {
349  const cell& oc = oldCells[cellI];
350 
351  cell& nc = newCells[nNewCells];
352  nc.setSize(oc.size());
353 
354  forAll(oc, i)
355  {
356  nc[i] = mirrorFaceLookup[oc[i]];
357  }
358 
359  nNewCells++;
360  }
361 
362  // Mirror the cell shapes
363  Info<< "Mirroring cell shapes." << endl;
364 
365  Info<< nl << "Creating new mesh" << endl;
366  mirrorMeshPtr_ = new fvMesh
367  (
368  io,
369  xferMove(newPoints),
370  xferMove(newFaces),
371  xferMove(newCells)
372  );
373 
374  fvMesh& pMesh = *mirrorMeshPtr_;
375 
376  // Add the boundary patches
377  List<polyPatch*> p(newPatchSizes.size());
378 
379  forAll(p, patchI)
380  {
381  p[patchI] = boundaryMesh()[newToOldPatch[patchI]].clone
382  (
383  pMesh.boundaryMesh(),
384  patchI,
385  newPatchSizes[patchI],
386  newPatchStarts[patchI]
387  ).ptr();
388  }
389 
390  pMesh.addPatches(p);
391 }
392 
393 
394 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
395 
397 {}
398 
399 
400 // ************************************************************************* //
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
p
p
Definition: pEqn.H:62
Foam::labelList
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:406
Foam::constant::atomic::alpha
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
Definition: readThermalProperties.H:216
boundary
faceListList boundary(nPatches)
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
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::boolListList
List< List< bool > > boolListList
Definition: boolList.H:51
plane.H
Foam::nl
static const char nl
Definition: Ostream.H:260
Foam::Info
messageStream Info
Foam::polyPatchList
PtrList< polyPatch > polyPatchList
container classes for polyPatch
Definition: polyPatchList.H:45
Foam::cellList
List< cell > cellList
list of cells
Definition: cellList.H:42
Foam::boolList
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
s
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
readScalar
#define readScalar
Definition: doubleScalar.C:38
Foam::List::setSize
void setSize(const label)
Reset size of List.
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Foam::xferMove
Xfer< T > xferMove(T &)
Construct by transferring the contents of the arg.
Foam::faceList
List< face > faceList
Definition: faceListFwd.H:43
Foam::mirrorFvMesh::mirrorFvMesh
mirrorFvMesh(const mirrorFvMesh &)
Disallow default bitwise copy construct.
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:94
List
Definition: Test.C:19
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::labelUList
UList< label > labelUList
Definition: UList.H:63
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:259
Foam::system
int system(const std::string &command)
Execute the specified command.
Definition: POSIX.C:1155
Foam::mirrorFvMesh::~mirrorFvMesh
~mirrorFvMesh()
Destructor.
mirrorFvMesh.H