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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2017 OpenFOAM Foundation
9 Copyright (C) 2017-2018 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "mirrorFvMesh.H"
30#include "Time.H"
31#include "plane.H"
32
33// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
34
36:
37 mirrorFvMesh
38 (
39 io,
40 IOdictionary
41 (
42 IOobject
43 (
44 "mirrorMeshDict",
45 io.time().system(),
46 io.time(),
47 IOobject::MUST_READ_IF_MODIFIED,
48 IOobject::NO_WRITE
49 )
50 )
51 )
52{}
53
54
56(
57 const IOobject& io,
58 const IOdictionary& mirrorDict
59)
60:
61 fvMesh(io)
62{
63 plane mirrorPlane(mirrorDict);
64
65 const scalar planeTolerance
66 (
67 mirrorDict.get<scalar>("planeTolerance")
68 );
69
70 const pointField& oldPoints = points();
71 const faceList& oldFaces = faces();
72 const cellList& oldCells = cells();
73 const label nOldInternalFaces = nInternalFaces();
74 const polyPatchList& oldPatches = boundaryMesh();
75
76 Info<< "Mirroring mesh at origin:" << mirrorPlane.origin()
77 << " normal:" << mirrorPlane.normal() << nl;
78
79 // Mirror the points
80 Info<< "Mirroring points. Old points: " << oldPoints.size();
81
82 pointField newPoints(2*oldPoints.size());
83 label nNewPoints = 0;
84
85 labelList mirrorPointLookup(oldPoints.size(), -1);
86
87 // Grab the old points
88 for (const point& pt : oldPoints)
89 {
90 newPoints[nNewPoints] = pt;
91 ++nNewPoints;
92 }
93
94 forAll(oldPoints, pointi)
95 {
96 const scalar alpha =
97 mirrorPlane.normalIntersect
98 (
99 oldPoints[pointi],
100 mirrorPlane.normal()
101 );
102
103 // Check plane on tolerance
104 if (mag(alpha) > planeTolerance)
105 {
106 // The point gets mirrored
107 newPoints[nNewPoints] =
108 oldPoints[pointi] + 2.0*alpha*mirrorPlane.normal();
109
110 // remember the point correspondence
111 mirrorPointLookup[pointi] = nNewPoints;
112 nNewPoints++;
113 }
114 else
115 {
116 // The point is on the plane and does not get mirrored
117 // Adjust plane location
118 newPoints[nNewPoints] =
119 oldPoints[pointi] + alpha*mirrorPlane.normal();
120
121 mirrorPointLookup[pointi] = pointi;
122 }
123 }
124
125 // Reset the size of the point list
126 Info<< " New points: " << nNewPoints << endl;
127 newPoints.setSize(nNewPoints);
128
129 // Construct new to old map
130 pointMapPtr_.reset(new labelList(newPoints.size()));
131 labelList& pointMap = pointMapPtr_();
132 // Insert old points
133 forAll(oldPoints, oldPointi)
134 {
135 pointMap[oldPointi] = oldPointi;
136 }
137 forAll(mirrorPointLookup, oldPointi)
138 {
139 pointMap[mirrorPointLookup[oldPointi]] = oldPointi;
140 }
141
142
143 Info<< "Mirroring faces. Old faces: " << oldFaces.size();
144
145 // Algorithm:
146 // During mirroring, the faces that were previously boundary faces
147 // in the mirror plane may become internal faces. In order to
148 // deal with the ordering of the faces, the algorithm is split
149 // into two parts. For original faces, the internal faces are
150 // distributed to their owner cells. Once all internal faces are
151 // distributed, the boundary faces are visited and if they are in
152 // the mirror plane they are added to the master cells (the future
153 // boundary faces are not touched). After the first phase, the
154 // internal faces are collected in the cell order and numbering
155 // information is added. Then, the internal faces are mirrored
156 // and the face numbering data is stored for the mirrored section.
157 // Once all the internal faces are mirrored, the boundary faces
158 // are added by mirroring the faces patch by patch.
159
160 // Distribute internal faces
161 labelListList newCellFaces(oldCells.size());
162
163 const labelUList& oldOwnerStart = lduAddr().ownerStartAddr();
164
165 forAll(newCellFaces, celli)
166 {
167 labelList& curFaces = newCellFaces[celli];
168
169 const label s = oldOwnerStart[celli];
170 const label e = oldOwnerStart[celli + 1];
171
172 curFaces.setSize(e - s);
173
174 forAll(curFaces, i)
175 {
176 curFaces[i] = s + i;
177 }
178 }
179
180 // Distribute boundary faces. Remember the faces that have been inserted
181 // as internal
182 boolListList insertedBouFace(oldPatches.size());
183
184 forAll(oldPatches, patchi)
185 {
186 const polyPatch& curPatch = oldPatches[patchi];
187
188 if (curPatch.coupled())
189 {
191 << "Found coupled patch " << curPatch.name() << endl
192 << " Mirroring faces on coupled patches destroys"
193 << " the ordering. This might be fixed by running a dummy"
194 << " createPatch afterwards." << endl;
195 }
196
197 boolList& curInsBouFace = insertedBouFace[patchi];
198
199 curInsBouFace.setSize(curPatch.size());
200 curInsBouFace = false;
201
202 // Get faceCells for face insertion
203 const labelUList& curFaceCells = curPatch.faceCells();
204 const label curStart = curPatch.start();
205
206 forAll(curPatch, facei)
207 {
208 // Find out if the mirrored face is identical to the
209 // original. If so, the face needs to become internal and
210 // added to its owner cell
211 const face& origFace = curPatch[facei];
212
213 face mirrorFace(origFace.size());
214 forAll(mirrorFace, pointi)
215 {
216 mirrorFace[pointi] = mirrorPointLookup[origFace[pointi]];
217 }
218
219 if (origFace == mirrorFace)
220 {
221 // The mirror is identical to current face. This will
222 // become an internal face
223 const label oldSize = newCellFaces[curFaceCells[facei]].size();
224
225 newCellFaces[curFaceCells[facei]].setSize(oldSize + 1);
226 newCellFaces[curFaceCells[facei]][oldSize] = curStart + facei;
227
228 curInsBouFace[facei] = true;
229 }
230 }
231 }
232
233 // Construct the new list of faces. Boundary faces are added
234 // last, cush that each patch is mirrored separately. The
235 // addressing is stored in two separate arrays: first for the
236 // original cells (face order has changed) and then for the
237 // mirrored cells.
238 labelList masterFaceLookup(oldFaces.size(), -1);
239 labelList mirrorFaceLookup(oldFaces.size(), -1);
240
241 faceList newFaces(2*oldFaces.size());
242 label nNewFaces = 0;
243
244 // Insert original (internal) faces
245 forAll(newCellFaces, celli)
246 {
247 const labelList& curCellFaces = newCellFaces[celli];
248
249 forAll(curCellFaces, cfI)
250 {
251 newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
252 masterFaceLookup[curCellFaces[cfI]] = nNewFaces;
253
254 nNewFaces++;
255 }
256 }
257
258 // Mirror internal faces
259 for (label facei = 0; facei < nOldInternalFaces; facei++)
260 {
261 const face& oldFace = oldFaces[facei];
262 face& nf = newFaces[nNewFaces];
263 nf.setSize(oldFace.size());
264
265 nf[0] = mirrorPointLookup[oldFace[0]];
266
267 for (label i = 1; i < oldFace.size(); i++)
268 {
269 nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
270 }
271
272 mirrorFaceLookup[facei] = nNewFaces;
273 nNewFaces++;
274 }
275
276 // Mirror boundary faces patch by patch
277 labelList newPatchSizes(boundary().size(), Zero);
278 labelList newPatchStarts(boundary().size(), -1);
279
280 forAll(boundaryMesh(), patchi)
281 {
282 const label curPatchSize = boundaryMesh()[patchi].size();
283 const label curPatchStart = boundaryMesh()[patchi].start();
284 const boolList& curInserted = insertedBouFace[patchi];
285
286 newPatchStarts[patchi] = nNewFaces;
287
288 // Master side
289 for (label facei = 0; facei < curPatchSize; facei++)
290 {
291 // Check if the face has already been added. If not, add it and
292 // insert the numbering details.
293 if (!curInserted[facei])
294 {
295 newFaces[nNewFaces] = oldFaces[curPatchStart + facei];
296
297 masterFaceLookup[curPatchStart + facei] = nNewFaces;
298 nNewFaces++;
299 }
300 }
301
302 // Mirror side
303 for (label facei = 0; facei < curPatchSize; facei++)
304 {
305 // Check if the face has already been added. If not, add it and
306 // insert the numbering details.
307 if (!curInserted[facei])
308 {
309 const face& oldFace = oldFaces[curPatchStart + facei];
310 face& nf = newFaces[nNewFaces];
311 nf.setSize(oldFace.size());
312
313 nf[0] = mirrorPointLookup[oldFace[0]];
314
315 for (label i = 1; i < oldFace.size(); i++)
316 {
317 nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
318 }
319
320 mirrorFaceLookup[curPatchStart + facei] = nNewFaces;
321 nNewFaces++;
322 }
323 else
324 {
325 // Grab the index of the master face for the mirror side
326 mirrorFaceLookup[curPatchStart + facei] =
327 masterFaceLookup[curPatchStart + facei];
328 }
329 }
330
331 // If patch exists, grab the name and type of the original patch
332 if (nNewFaces > newPatchStarts[patchi])
333 {
334 newPatchSizes[patchi] =
335 nNewFaces - newPatchStarts[patchi];
336 }
337 }
338
339 // Tidy up the lists
340 newFaces.setSize(nNewFaces);
341 Info<< " New faces: " << nNewFaces << endl;
342
343 Info<< "Mirroring patches. Old patches: " << boundary().size()
344 << " New patches: " << boundary().size() << endl;
345
346 Info<< "Mirroring cells. Old cells: " << oldCells.size()
347 << " New cells: " << 2*oldCells.size() << endl;
348
349 cellList newCells(2*oldCells.size());
350 label nNewCells = 0;
351
352 // Construct new to old cell map
353 cellMapPtr_.reset(new labelList(newCells.size()));
354 labelList& cellMap = cellMapPtr_();
355
356 // Grab the original cells. Take care of face renumbering.
357 forAll(oldCells, celli)
358 {
359 const cell& oc = oldCells[celli];
360
361 cell& nc = newCells[nNewCells];
362 nc.setSize(oc.size());
363
364 forAll(oc, i)
365 {
366 nc[i] = masterFaceLookup[oc[i]];
367 }
368
369 cellMap[nNewCells] = celli;
370
371 nNewCells++;
372 }
373
374 // Mirror the cells
375 forAll(oldCells, celli)
376 {
377 const cell& oc = oldCells[celli];
378
379 cell& nc = newCells[nNewCells];
380 nc.setSize(oc.size());
381
382 forAll(oc, i)
383 {
384 nc[i] = mirrorFaceLookup[oc[i]];
385 }
386
387 cellMap[nNewCells] = celli;
388
389 nNewCells++;
390 }
391
392 // Mirror the cell shapes
393 Info<< "Mirroring cell shapes." << endl;
394
395 Info<< nl << "Creating new mesh" << endl;
396 mirrorMeshPtr_ = autoPtr<fvMesh>::New
397 (
398 io,
399 std::move(newPoints),
400 std::move(newFaces),
401 std::move(newCells)
402 );
403
404 fvMesh& pMesh = mirrorMeshPtr_();
405
406 // Add the boundary patches
407 List<polyPatch*> p(newPatchSizes.size());
408
409 forAll(p, patchi)
410 {
411 p[patchi] = boundaryMesh()[patchi].clone
412 (
413 pMesh.boundaryMesh(),
414 patchi,
415 newPatchSizes[patchi],
416 newPatchStarts[patchi]
417 ).ptr();
418 }
419
420 pMesh.addPatches(p);
421}
422
423
424// ************************************************************************* //
void setSize(const label n)
Alias for resize()
Definition: List.H:218
volScalarField & p
faceListList boundary
const pointField & points
const cellShapeList & cells
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
IOobject io("surfaceFilmProperties", mesh.time().constant(), mesh, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE, false)
#define WarningInFunction
Report a warning using Foam::Warning.
PtrList< polyPatch > polyPatchList
Store lists of polyPatch as a PtrList.
Definition: polyPatch.H:63
int system(const std::string &command, const bool bg=false)
Execute the specified command via the shell.
Definition: MSwindows.C:1158
List< label > labelList
A List of labels.
Definition: List.H:66
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
List< List< bool > > boolListList
A List of boolList.
Definition: boolList.H:48
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
List< bool > boolList
A List of bools.
Definition: List.H:64
List< face > faceList
A List of faces.
Definition: faceListFwd.H:47
UList< label > labelUList
A UList of labels.
Definition: UList.H:85
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & alpha
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333