mirrorFvMesh.C 12 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
6
7
8
9
10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11
12
13
14
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
15
16
17
18
19
20
21

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

\*---------------------------------------------------------------------------*/

#include "mirrorFvMesh.H"
#include "Time.H"
#include "plane.H"

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io)
:
    fvMesh(io),
    mirrorMeshDict_
    (
        IOobject
        (
            "mirrorMeshDict",
            time().system(),
            *this,
42
            IOobject::MUST_READ_IF_MODIFIED,
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
            IOobject::NO_WRITE
        )
    ),
    mirrorMeshPtr_(NULL)
{
    plane mirrorPlane(mirrorMeshDict_);

    scalar planeTolerance
    (
        readScalar(mirrorMeshDict_.lookup("planeTolerance"))
    );

    const pointField& oldPoints = points();
    const faceList& oldFaces = faces();
    const cellList& oldCells = cells();
    const label nOldInternalFaces = nInternalFaces();
    const polyPatchList& oldPatches = boundaryMesh();

    // Mirror the points
Mark Olesen's avatar
Mark Olesen committed
62
    Info<< "Mirroring points. Old points: " << oldPoints.size();
63
64
65
66
67
68
69

    pointField newPoints(2*oldPoints.size());
    label nNewPoints = 0;

    labelList mirrorPointLookup(oldPoints.size(), -1);

    // Grab the old points
70
    forAll(oldPoints, pointi)
71
    {
72
        newPoints[nNewPoints] = oldPoints[pointi];
73
74
75
        nNewPoints++;
    }

76
    forAll(oldPoints, pointi)
77
78
79
80
    {
        scalar alpha =
            mirrorPlane.normalIntersect
            (
81
                oldPoints[pointi],
82
83
84
85
86
87
88
89
                mirrorPlane.normal()
            );

        // Check plane on tolerance
        if (mag(alpha) > planeTolerance)
        {
            // The point gets mirrored
            newPoints[nNewPoints] =
90
                oldPoints[pointi] + 2.0*alpha*mirrorPlane.normal();
91
92

            // remember the point correspondence
93
            mirrorPointLookup[pointi] = nNewPoints;
94
95
96
97
98
99
100
            nNewPoints++;
        }
        else
        {
            // The point is on the plane and does not get mirrored
            // Adjust plane location
            newPoints[nNewPoints] =
101
                oldPoints[pointi] + alpha*mirrorPlane.normal();
102

103
            mirrorPointLookup[pointi] = pointi;
104
105
106
107
        }
    }

    // Reset the size of the point list
Mark Olesen's avatar
Mark Olesen committed
108
    Info<< " New points: " << nNewPoints << endl;
109
110
    newPoints.setSize(nNewPoints);

Mark Olesen's avatar
Mark Olesen committed
111
    Info<< "Mirroring faces. Old faces: " << oldFaces.size();
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

    // Algorithm:
    // During mirroring, the faces that were previously boundary faces
    // in the mirror plane may become ineternal faces. In order to
    // deal with the ordering of the faces, the algorithm is split
    // into two parts.  For original faces, the internal faces are
    // distributed to their owner cells.  Once all internal faces are
    // distributed, the boundary faces are visited and if they are in
    // the mirror plane they are added to the master cells (the future
    // boundary faces are not touched).  After the first phase, the
    // internal faces are collected in the cell order and numbering
    // information is added.  Then, the internal faces are mirrored
    // and the face numbering data is stored for the mirrored section.
    // Once all the internal faces are mirrored, the boundary faces
    // are added by mirroring the faces patch by patch.

    // Distribute internal faces
    labelListList newCellFaces(oldCells.size());

131
    const labelUList& oldOwnerStart = lduAddr().ownerStartAddr();
132

133
    forAll(newCellFaces, celli)
134
    {
135
        labelList& curFaces = newCellFaces[celli];
136

137
138
        const label s = oldOwnerStart[celli];
        const label e = oldOwnerStart[celli + 1];
139
140
141

        curFaces.setSize(e - s);

142
        forAll(curFaces, i)
143
144
145
146
147
148
149
150
151
        {
            curFaces[i] = s + i;
        }
    }

    // Distribute boundary faces.  Remember the faces that have been inserted
    // as internal
    boolListList insertedBouFace(oldPatches.size());

152
    forAll(oldPatches, patchi)
153
    {
154
        const polyPatch& curPatch = oldPatches[patchi];
155
156
157

        if (curPatch.coupled())
        {
158
            WarningInFunction
159
160
161
162
163
164
                << "Found coupled patch " << curPatch.name() << endl
                << "    Mirroring faces on coupled patches destroys"
                << " the ordering. This might be fixed by running a dummy"
                << " createPatch afterwards." << endl;
        }

165
        boolList& curInsBouFace = insertedBouFace[patchi];
166
167
168
169
170

        curInsBouFace.setSize(curPatch.size());
        curInsBouFace = false;

        // Get faceCells for face insertion
171
        const labelUList& curFaceCells = curPatch.faceCells();
172
173
        const label curStart = curPatch.start();

174
        forAll(curPatch, facei)
175
176
177
178
        {
            // Find out if the mirrored face is identical to the
            // original.  If so, the face needs to become internal and
            // 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
188
189
190
            }

            if (origFace == mirrorFace)
            {
                // The mirror is identical to current face.  This will
                // become an internal face
191
                const label oldSize = newCellFaces[curFaceCells[facei]].size();
192

193
194
                newCellFaces[curFaceCells[facei]].setSize(oldSize + 1);
                newCellFaces[curFaceCells[facei]][oldSize] = curStart + facei;
195

196
                curInsBouFace[facei] = true;
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
            }
        }
    }

    // Construct the new list of faces.  Boundary faces are added
    // last, cush that each patch is mirrored separately.  The
    // addressing is stored in two separate arrays: first for the
    // original cells (face order has changed) and then for the
    // mirrored cells.
    labelList masterFaceLookup(oldFaces.size(), -1);
    labelList mirrorFaceLookup(oldFaces.size(), -1);

    faceList newFaces(2*oldFaces.size());
    label nNewFaces = 0;

    // Insert original (internal) faces
213
    forAll(newCellFaces, celli)
214
    {
215
        const labelList& curCellFaces = newCellFaces[celli];
216

217
        forAll(curCellFaces, cfI)
218
219
220
221
222
223
224
225
226
        {
            newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
            masterFaceLookup[curCellFaces[cfI]] = nNewFaces;

            nNewFaces++;
        }
    }

    // Mirror internal faces
227
    for (label facei = 0; facei < nOldInternalFaces; facei++)
228
    {
229
        const face& oldFace = oldFaces[facei];
230
231
232
233
234
235
236
237
238
239
        face& nf = newFaces[nNewFaces];
        nf.setSize(oldFace.size());

        nf[0] = mirrorPointLookup[oldFace[0]];

        for (label i = 1; i < oldFace.size(); i++)
        {
            nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
        }

240
        mirrorFaceLookup[facei] = nNewFaces;
241
242
243
244
245
        nNewFaces++;
    }

    // Mirror boundary faces patch by patch

246
247

    labelList newToOldPatch(boundary().size(), -1);
248
249
250
251
    labelList newPatchSizes(boundary().size(), -1);
    labelList newPatchStarts(boundary().size(), -1);
    label nNewPatches = 0;

252
    forAll(boundaryMesh(), patchi)
253
    {
254
255
256
        const label curPatchSize = boundaryMesh()[patchi].size();
        const label curPatchStart = boundaryMesh()[patchi].start();
        const boolList& curInserted = insertedBouFace[patchi];
257
258
259
260

        newPatchStarts[nNewPatches] = nNewFaces;

        // Master side
261
        for (label facei = 0; facei < curPatchSize; facei++)
262
263
264
        {
            // Check if the face has already been added.  If not, add it and
            // insert the numbering details.
265
            if (!curInserted[facei])
266
            {
267
                newFaces[nNewFaces] = oldFaces[curPatchStart + facei];
268

269
                masterFaceLookup[curPatchStart + facei] = nNewFaces;
270
271
272
273
274
                nNewFaces++;
            }
        }

        // Mirror side
275
        for (label facei = 0; facei < curPatchSize; facei++)
276
277
278
        {
            // Check if the face has already been added.  If not, add it and
            // insert the numbering details.
279
            if (!curInserted[facei])
280
            {
281
                const face& oldFace = oldFaces[curPatchStart + facei];
282
283
284
285
286
287
288
289
290
291
                face& nf = newFaces[nNewFaces];
                nf.setSize(oldFace.size());

                nf[0] = mirrorPointLookup[oldFace[0]];

                for (label i = 1; i < oldFace.size(); i++)
                {
                    nf[i] = mirrorPointLookup[oldFace[oldFace.size() - i]];
                }

292
                mirrorFaceLookup[curPatchStart + facei] = nNewFaces;
293
294
295
296
297
                nNewFaces++;
            }
            else
            {
                // Grab the index of the master face for the mirror side
298
299
                mirrorFaceLookup[curPatchStart + facei] =
                    masterFaceLookup[curPatchStart + facei];
300
301
302
303
304
305
            }
        }

        // If patch exists, grab the name and type of the original patch
        if (nNewFaces > newPatchStarts[nNewPatches])
        {
306
            newToOldPatch[nNewPatches] = patchi;
307

308
309
310
311
312
313
314
315
316
            newPatchSizes[nNewPatches] =
                nNewFaces - newPatchStarts[nNewPatches];

            nNewPatches++;
        }
    }

    // Tidy up the lists
    newFaces.setSize(nNewFaces);
Mark Olesen's avatar
Mark Olesen committed
317
    Info<< " New faces: " << nNewFaces << endl;
318

319
    newToOldPatch.setSize(nNewPatches);
320
321
322
    newPatchSizes.setSize(nNewPatches);
    newPatchStarts.setSize(nNewPatches);

Mark Olesen's avatar
Mark Olesen committed
323
    Info<< "Mirroring patches. Old patches: " << boundary().size()
324
325
326
327
328
329
330
331
332
        << " New patches: " << nNewPatches << endl;

    Info<< "Mirroring cells.  Old cells: " << oldCells.size()
        << " New cells: " << 2*oldCells.size() << endl;

    cellList newCells(2*oldCells.size());
    label nNewCells = 0;

    // Grab the original cells.  Take care of face renumbering.
333
    forAll(oldCells, celli)
334
    {
335
        const cell& oc = oldCells[celli];
336
337
338
339

        cell& nc = newCells[nNewCells];
        nc.setSize(oc.size());

340
        forAll(oc, i)
341
342
343
344
345
346
347
348
        {
            nc[i] = masterFaceLookup[oc[i]];
        }

        nNewCells++;
    }

    // Mirror the cells
349
    forAll(oldCells, celli)
350
    {
351
        const cell& oc = oldCells[celli];
352
353
354
355

        cell& nc = newCells[nNewCells];
        nc.setSize(oc.size());

356
        forAll(oc, i)
357
358
359
360
361
362
363
364
        {
            nc[i] = mirrorFaceLookup[oc[i]];
        }

        nNewCells++;
    }

    // Mirror the cell shapes
Mark Olesen's avatar
Mark Olesen committed
365
    Info<< "Mirroring cell shapes." << endl;
366

Mark Olesen's avatar
Mark Olesen committed
367
    Info<< nl << "Creating new mesh" << endl;
368
369
370
    mirrorMeshPtr_ = new fvMesh
    (
        io,
371
372
373
        xferMove(newPoints),
        xferMove(newFaces),
        xferMove(newCells)
374
375
376
377
378
    );

    fvMesh& pMesh = *mirrorMeshPtr_;

    // Add the boundary patches
379
    List<polyPatch*> p(newPatchSizes.size());
380

381
    forAll(p, patchi)
382
    {
383
        p[patchi] = boundaryMesh()[newToOldPatch[patchi]].clone
384
        (
385
            pMesh.boundaryMesh(),
386
387
388
            patchi,
            newPatchSizes[patchi],
            newPatchStarts[patchi]
389
390
391
392
393
394
395
396
397
398
399
400
401
402
        ).ptr();
    }

    pMesh.addPatches(p);
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::mirrorFvMesh::~mirrorFvMesh()
{}


// ************************************************************************* //