mirrorFvMesh.C 11.7 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
     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
7
8
9
10
-------------------------------------------------------------------------------
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
            IOobject::NO_WRITE
        )
    ),
46
    mirrorMeshPtr_(nullptr)
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
{
    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
        nNewFaces++;
    }

    // Mirror boundary faces patch by patch
245
    labelList newPatchSizes(boundary().size(), 0);
246
247
    labelList newPatchStarts(boundary().size(), -1);

248
    forAll(boundaryMesh(), patchi)
249
    {
250
251
252
        const label curPatchSize = boundaryMesh()[patchi].size();
        const label curPatchStart = boundaryMesh()[patchi].start();
        const boolList& curInserted = insertedBouFace[patchi];
253

254
        newPatchStarts[patchi] = nNewFaces;
255
256

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

265
                masterFaceLookup[curPatchStart + facei] = nNewFaces;
266
267
268
269
270
                nNewFaces++;
            }
        }

        // Mirror side
271
        for (label facei = 0; facei < curPatchSize; facei++)
272
273
274
        {
            // Check if the face has already been added.  If not, add it and
            // insert the numbering details.
275
            if (!curInserted[facei])
276
            {
277
                const face& oldFace = oldFaces[curPatchStart + facei];
278
279
280
281
282
283
284
285
286
287
                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]];
                }

288
                mirrorFaceLookup[curPatchStart + facei] = nNewFaces;
289
290
291
292
293
                nNewFaces++;
            }
            else
            {
                // Grab the index of the master face for the mirror side
294
295
                mirrorFaceLookup[curPatchStart + facei] =
                    masterFaceLookup[curPatchStart + facei];
296
297
298
299
            }
        }

        // If patch exists, grab the name and type of the original patch
300
        if (nNewFaces > newPatchStarts[patchi])
301
        {
302
303
            newPatchSizes[patchi] =
                nNewFaces - newPatchStarts[patchi];
304
305
306
307
308
        }
    }

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

Mark Olesen's avatar
Mark Olesen committed
311
    Info<< "Mirroring patches. Old patches: " << boundary().size()
312
        << " New patches: " << boundary().size() << endl;
313
314
315
316
317
318
319
320

    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.
321
    forAll(oldCells, celli)
322
    {
323
        const cell& oc = oldCells[celli];
324
325
326
327

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

328
        forAll(oc, i)
329
330
331
332
333
334
335
336
        {
            nc[i] = masterFaceLookup[oc[i]];
        }

        nNewCells++;
    }

    // Mirror the cells
337
    forAll(oldCells, celli)
338
    {
339
        const cell& oc = oldCells[celli];
340
341
342
343

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

344
        forAll(oc, i)
345
346
347
348
349
350
351
352
        {
            nc[i] = mirrorFaceLookup[oc[i]];
        }

        nNewCells++;
    }

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

Mark Olesen's avatar
Mark Olesen committed
355
    Info<< nl << "Creating new mesh" << endl;
356
357
358
    mirrorMeshPtr_ = new fvMesh
    (
        io,
359
360
361
        xferMove(newPoints),
        xferMove(newFaces),
        xferMove(newCells)
362
363
364
365
366
    );

    fvMesh& pMesh = *mirrorMeshPtr_;

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

369
    forAll(p, patchi)
370
    {
371
        p[patchi] = boundaryMesh()[patchi].clone
372
        (
373
            pMesh.boundaryMesh(),
374
375
376
            patchi,
            newPatchSizes[patchi],
            newPatchStarts[patchi]
377
378
379
380
381
382
383
384
385
386
387
388
389
390
        ).ptr();
    }

    pMesh.addPatches(p);
}


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

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


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