mirrorFvMesh.C 12.8 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
Mark Olesen's avatar
Mark Olesen committed
5
    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    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 2 of the License, or (at your
    option) any later version.

    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
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

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

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

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

const Foam::label Foam::mirrorFvMesh::cellRenumber[8][8] =
{
    {-1, -1, -1, -1, -1, -1, -1, -1},    // unknown
Mark Olesen's avatar
Mark Olesen committed
36
37
    {-1, -1, -1, -1, -1, -1, -1, -1},    //
    {-1, -1, -1, -1, -1, -1, -1, -1},    //
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    { 0,  3,  2,  1,  4,  7,  6,  5},    // hex
    { 2,  1,  0,  5,  4,  3,  6, -1},    // wedge
    { 0,  2,  1,  3,  5,  4, -1, -1},    // prism
    { 0,  3,  2,  1,  4, -1, -1, -1},    // pyramid
    { 2,  1,  0,  3, -1, -1, -1, -1},    // tet
};

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

Foam::mirrorFvMesh::mirrorFvMesh(const IOobject& io)
:
    fvMesh(io),
    mirrorMeshDict_
    (
        IOobject
        (
            "mirrorMeshDict",
            time().system(),
            *this,
            IOobject::MUST_READ,
            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
77
    Info<< "Mirroring points. Old points: " << oldPoints.size();
78
79
80
81
82
83
84

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

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

    // Grab the old points
85
    forAll(oldPoints, pointI)
86
87
88
89
90
    {
        newPoints[nNewPoints] = oldPoints[pointI];
        nNewPoints++;
    }

91
    forAll(oldPoints, pointI)
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
    {
        scalar alpha =
            mirrorPlane.normalIntersect
            (
                oldPoints[pointI],
                mirrorPlane.normal()
            );

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

            // remember the point correspondence
            mirrorPointLookup[pointI] = nNewPoints;
            nNewPoints++;
        }
        else
        {
            // The point is on the plane and does not get mirrored
            // Adjust plane location
            newPoints[nNewPoints] =
                oldPoints[pointI] + alpha*mirrorPlane.normal();

            mirrorPointLookup[pointI] = pointI;
        }
    }

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

Mark Olesen's avatar
Mark Olesen committed
126
    Info<< "Mirroring faces. Old faces: " << oldFaces.size();
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

    // 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());

    const unallocLabelList& oldOwnerStart = lduAddr().ownerStartAddr();

148
    forAll(newCellFaces, cellI)
149
150
151
152
153
154
155
156
    {
        labelList& curFaces = newCellFaces[cellI];

        const label s = oldOwnerStart[cellI];
        const label e = oldOwnerStart[cellI + 1];

        curFaces.setSize(e - s);

157
        forAll(curFaces, i)
158
159
160
161
162
163
164
165
166
        {
            curFaces[i] = s + i;
        }
    }

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

167
    forAll(oldPatches, patchI)
168
169
    {
        const polyPatch& curPatch = oldPatches[patchI];
170
171
172
173
174
175
176
177
178
179

        if (curPatch.coupled())
        {
            WarningIn("mirrorFvMesh::mirrorFvMesh(const IOobject&)")
                << "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;
        }

180
181
182
183
184
185
186
187
188
        boolList& curInsBouFace = insertedBouFace[patchI];

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

        // Get faceCells for face insertion
        const unallocLabelList& curFaceCells = curPatch.faceCells();
        const label curStart = curPatch.start();

189
        forAll(curPatch, faceI)
190
191
192
193
194
195
196
        {
            // 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
            const face& origFace = curPatch[faceI];

            face mirrorFace(origFace.size());
197
            forAll(mirrorFace, pointI)
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
            {
                mirrorFace[pointI] = mirrorPointLookup[origFace[pointI]];
            }

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

                newCellFaces[curFaceCells[faceI]].setSize(oldSize + 1);
                newCellFaces[curFaceCells[faceI]][oldSize] = curStart + faceI;

                curInsBouFace[faceI] = true;
            }
        }
    }

    // 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
228
    forAll(newCellFaces, cellI)
229
230
231
    {
        const labelList& curCellFaces = newCellFaces[cellI];

232
        forAll(curCellFaces, cfI)
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
        {
            newFaces[nNewFaces] = oldFaces[curCellFaces[cfI]];
            masterFaceLookup[curCellFaces[cfI]] = nNewFaces;

            nNewFaces++;
        }
    }

    // Mirror internal faces
    for (label faceI = 0; faceI < nOldInternalFaces; faceI++)
    {
        const face& oldFace = oldFaces[faceI];
        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]];
        }

        mirrorFaceLookup[faceI] = nNewFaces;
        nNewFaces++;
    }

    // Mirror boundary faces patch by patch

    wordList newPatchTypes(boundary().size());
    wordList newPatchNames(boundary().size());
    labelList newPatchSizes(boundary().size(), -1);
    labelList newPatchStarts(boundary().size(), -1);
    label nNewPatches = 0;

267
    forAll(boundaryMesh(), patchI)
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
    {
        const label curPatchSize = boundaryMesh()[patchI].size();
        const label curPatchStart = boundaryMesh()[patchI].start();
        const boolList& curInserted = insertedBouFace[patchI];

        newPatchStarts[nNewPatches] = nNewFaces;

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

                masterFaceLookup[curPatchStart + faceI] = nNewFaces;
                nNewFaces++;
            }
        }

        // Mirror side
        for (label faceI = 0; faceI < curPatchSize; faceI++)
        {
            // Check if the face has already been added.  If not, add it and
            // insert the numbering details.
            if (!curInserted[faceI])
            {
                const face& oldFace = oldFaces[curPatchStart + faceI];
                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]];
                }

                mirrorFaceLookup[curPatchStart + faceI] = nNewFaces;
                nNewFaces++;
            }
            else
            {
                // Grab the index of the master face for the mirror side
                mirrorFaceLookup[curPatchStart + faceI] =
                    masterFaceLookup[curPatchStart + faceI];
            }
        }

        // If patch exists, grab the name and type of the original patch
        if (nNewFaces > newPatchStarts[nNewPatches])
        {
            newPatchTypes[nNewPatches] = boundaryMesh()[patchI].type();
            newPatchNames[nNewPatches] = boundaryMesh()[patchI].name();
            newPatchSizes[nNewPatches] =
                nNewFaces - newPatchStarts[nNewPatches];

            nNewPatches++;
        }
    }

    // Tidy up the lists
    newFaces.setSize(nNewFaces);
Mark Olesen's avatar
Mark Olesen committed
332
    Info<< " New faces: " << nNewFaces << endl;
333
334
335
336
337
338

    newPatchTypes.setSize(nNewPatches);
    newPatchNames.setSize(nNewPatches);
    newPatchSizes.setSize(nNewPatches);
    newPatchStarts.setSize(nNewPatches);

Mark Olesen's avatar
Mark Olesen committed
339
    Info<< "Mirroring patches. Old patches: " << boundary().size()
340
341
342
343
344
345
346
347
348
        << " 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.
349
    forAll(oldCells, cellI)
350
351
352
353
354
355
    {
        const cell& oc = oldCells[cellI];

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

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

        nNewCells++;
    }

    // Mirror the cells
365
    forAll(oldCells, cellI)
366
367
368
369
370
371
    {
        const cell& oc = oldCells[cellI];

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

372
        forAll(oc, i)
373
374
375
376
377
378
379
380
        {
            nc[i] = mirrorFaceLookup[oc[i]];
        }

        nNewCells++;
    }

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

Mark Olesen's avatar
Mark Olesen committed
383
    Info<< nl << "Creating new mesh" << endl;
384
385
386
    mirrorMeshPtr_ = new fvMesh
    (
        io,
387
388
389
        xferMove(newPoints),
        xferMove(newFaces),
        xferMove(newCells)
390
391
392
393
394
395
396
    );

    fvMesh& pMesh = *mirrorMeshPtr_;

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

397
    forAll(p, patchI)
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
    {
        p[patchI] = polyPatch::New
        (
            newPatchTypes[patchI],
            newPatchNames[patchI],
            newPatchSizes[patchI],
            newPatchStarts[patchI],
            patchI,
            pMesh.boundaryMesh()
        ).ptr();
    }

    pMesh.addPatches(p);
}


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

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


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