FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/Time.H>
28 #include <OpenFOAM/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,
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 unallocLabelList& 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  {
158  WarningIn("mirrorFvMesh::mirrorFvMesh(const IOobject&)")
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 unallocLabelList& 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 
246  wordList newPatchTypes(boundary().size());
247  wordList newPatchNames(boundary().size());
248  labelList newPatchSizes(boundary().size(), -1);
249  labelList newPatchStarts(boundary().size(), -1);
250  label nNewPatches = 0;
251 
252  forAll (boundaryMesh(), patchI)
253  {
254  const label curPatchSize = boundaryMesh()[patchI].size();
255  const label curPatchStart = boundaryMesh()[patchI].start();
256  const boolList& curInserted = insertedBouFace[patchI];
257 
258  newPatchStarts[nNewPatches] = nNewFaces;
259 
260  // Master side
261  for (label faceI = 0; faceI < curPatchSize; faceI++)
262  {
263  // Check if the face has already been added. If not, add it and
264  // insert the numbering details.
265  if (!curInserted[faceI])
266  {
267  newFaces[nNewFaces] = oldFaces[curPatchStart + faceI];
268 
269  masterFaceLookup[curPatchStart + faceI] = nNewFaces;
270  nNewFaces++;
271  }
272  }
273 
274  // Mirror side
275  for (label faceI = 0; faceI < curPatchSize; faceI++)
276  {
277  // Check if the face has already been added. If not, add it and
278  // insert the numbering details.
279  if (!curInserted[faceI])
280  {
281  const face& oldFace = oldFaces[curPatchStart + faceI];
282  face& nf = newFaces[nNewFaces];
283  nf.setSize(oldFace.size());
284 
285  nf[0] = mirrorPointLookup[oldFace[0]];
286 
287  for (label i = 1; i < oldFace.size(); i++)
288  {
289  nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
290  }
291 
292  mirrorFaceLookup[curPatchStart + faceI] = nNewFaces;
293  nNewFaces++;
294  }
295  else
296  {
297  // Grab the index of the master face for the mirror side
298  mirrorFaceLookup[curPatchStart + faceI] =
299  masterFaceLookup[curPatchStart + faceI];
300  }
301  }
302 
303  // If patch exists, grab the name and type of the original patch
304  if (nNewFaces > newPatchStarts[nNewPatches])
305  {
306  newPatchTypes[nNewPatches] = boundaryMesh()[patchI].type();
307  newPatchNames[nNewPatches] = boundaryMesh()[patchI].name();
308  newPatchSizes[nNewPatches] =
309  nNewFaces - newPatchStarts[nNewPatches];
310 
311  nNewPatches++;
312  }
313  }
314 
315  // Tidy up the lists
316  newFaces.setSize(nNewFaces);
317  Info << " New faces: " << nNewFaces << endl;
318 
319  newPatchTypes.setSize(nNewPatches);
320  newPatchNames.setSize(nNewPatches);
321  newPatchSizes.setSize(nNewPatches);
322  newPatchStarts.setSize(nNewPatches);
323 
324  Info << "Mirroring patches. Old patches: " << boundary().size()
325  << " New patches: " << nNewPatches << endl;
326 
327  Info<< "Mirroring cells. Old cells: " << oldCells.size()
328  << " New cells: " << 2*oldCells.size() << endl;
329 
330  cellList newCells(2*oldCells.size());
331  label nNewCells = 0;
332 
333  // Grab the original cells. Take care of face renumbering.
334  forAll (oldCells, cellI)
335  {
336  const cell& oc = oldCells[cellI];
337 
338  cell& nc = newCells[nNewCells];
339  nc.setSize(oc.size());
340 
341  forAll (oc, i)
342  {
343  nc[i] = masterFaceLookup[oc[i]];
344  }
345 
346  nNewCells++;
347  }
348 
349  // Mirror the cells
350  forAll (oldCells, cellI)
351  {
352  const cell& oc = oldCells[cellI];
353 
354  cell& nc = newCells[nNewCells];
355  nc.setSize(oc.size());
356 
357  forAll (oc, i)
358  {
359  nc[i] = mirrorFaceLookup[oc[i]];
360  }
361 
362  nNewCells++;
363  }
364 
365  // Mirror the cell shapes
366  Info << "Mirroring cell shapes." << endl;
367 
368  Info << nl << "Creating new mesh" << endl;
369  mirrorMeshPtr_ = new fvMesh
370  (
371  io,
372  xferMove(newPoints),
373  xferMove(newFaces),
374  xferMove(newCells)
375  );
376 
377  fvMesh& pMesh = *mirrorMeshPtr_;
378 
379  // Add the boundary patches
380  List<polyPatch*> p(newPatchTypes.size());
381 
382  forAll (p, patchI)
383  {
384  p[patchI] = polyPatch::New
385  (
386  newPatchTypes[patchI],
387  newPatchNames[patchI],
388  newPatchSizes[patchI],
389  newPatchStarts[patchI],
390  patchI,
391  pMesh.boundaryMesh()
392  ).ptr();
393  }
394 
395  pMesh.addPatches(p);
396 }
397 
398 
399 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
400 
402 {}
403 
404 
405 // ************************ vim: set sw=4 sts=4 et: ************************ //