InputGmsh.cpp 53.8 KB
Newer Older
1 2 3 4 5 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
////////////////////////////////////////////////////////////////////////////////
//
//  File: InputGmsh.cpp
//
//  For more information, please see: http://www.nektar.info/
//
//  The MIT License
//
//  Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
//  Department of Aeronautics, Imperial College London (UK), and Scientific
//  Computing and Imaging Institute, University of Utah (USA).
//
//  License for the specific language governing rights and limitations under
//  Permission is hereby granted, free of charge, to any person obtaining a
//  copy of this software and associated documentation files (the "Software"),
//  to deal in the Software without restriction, including without limitation
//  the rights to use, copy, modify, merge, publish, distribute, sublicense,
//  and/or sell copies of the Software, and to permit persons to whom the
//  Software is furnished to do so, subject to the following conditions:
//
//  The above copyright notice and this permission notice shall be included
//  in all copies or substantial portions of the Software.
//
//  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
//  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
//  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
//  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
//  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
//  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
//  DEALINGS IN THE SOFTWARE.
//
//  Description: GMSH converter.
//
////////////////////////////////////////////////////////////////////////////////

36 37 38
#include <string>
#include <iostream>
using namespace std;
39

40 41
#include <boost/tuple/tuple.hpp>

42 43 44 45 46 47 48 49 50
#include <NekMeshUtils/MeshElements/Element.h>
#include <NekMeshUtils/MeshElements/Point.h>
#include <NekMeshUtils/MeshElements/Line.h>
#include <NekMeshUtils/MeshElements/Triangle.h>
#include <NekMeshUtils/MeshElements/Quadrilateral.h>
#include <NekMeshUtils/MeshElements/Tetrahedron.h>
#include <NekMeshUtils/MeshElements/Pyramid.h>
#include <NekMeshUtils/MeshElements/Prism.h>
#include <NekMeshUtils/MeshElements/Hexahedron.h>
51

52
#include "InputGmsh.h"
53

54
using namespace Nektar::NekMeshUtils;
55

56 57
namespace Nektar
{
58 59 60
namespace Utilities
{

61 62
ModuleKey InputGmsh::className = GetModuleFactory().RegisterCreatorFunction(
    ModuleKey(eInputModule, "msh"), InputGmsh::create, "Reads Gmsh msh file.");
63

64
std::map<unsigned int, ElmtConfig> InputGmsh::elmMap = InputGmsh::GenElmMap();
65 66

/**
Dave Moxey's avatar
Dave Moxey committed
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
 * @brief Reorder Gmsh nodes so that they appear in a tensor product format
 * suitable for the interior of a Nektar++ quadrilateral.
 *
 * For an example, consider a second order quadrilateral. This routine will
 * produce a permutation map that applies the following permutation:
 *
 * ```
 *   3---6---2         6---7---8
 *   |       |         |       |
 *   7   8   5   --->  3   4   5
 *   |       |         |       |
 *   0---4---1         0---1---2
 *
 *     Gmsh          tensor-product
 * ```
 *
 * We assume that Gmsh uses a recursive ordering system, so that interior nodes
 * are reordered by calling this function recursively.
 *
 * @param nodes   The integer IDs of the nodes to be reordered, in Gmsh format.
 * @param n       The number of nodes in one coordinate direction. If this is
 *                zero we assume no reordering needs to be done and return the
 *                identity permutation.
 *
 * @return The nodes vector in tensor-product ordering.
92
 */
93
std::vector<int> quadTensorNodeOrdering(const std::vector<int> &nodes, int n)
94 95 96 97 98 99 100 101 102
{
    std::vector<int> nodeList;

    // Triangle
    nodeList.resize(nodes.size());

    // Vertices and edges
    nodeList[0] = nodes[0];
    if (n > 1)
103
    {
104 105 106
        nodeList[n - 1]     = nodes[1];
        nodeList[n * n - 1] = nodes[2];
        nodeList[n * (n - 1)] = nodes[3];
107
    }
108
    for (int i = 1; i < n - 1; ++i)
109
    {
110
        nodeList[i] = nodes[4 + i - 1];
111
    }
112
    for (int i = 1; i < n - 1; ++i)
113
    {
114
        nodeList[n * n - 1 - i] = nodes[4 + 2 * (n - 2) + i - 1];
115
    }
116

117 118 119 120
    // Interior (recursion)
    if (n > 2)
    {
        // Reorder interior nodes
121 122 123 124
        std::vector<int> interior((n - 2) * (n - 2));
        std::copy(
            nodes.begin() + 4 + 4 * (n - 2), nodes.end(), interior.begin());
        interior = quadTensorNodeOrdering(interior, n - 2);
125

126
        // Copy into full node list
127
        for (int j = 1; j < n - 1; ++j)
128
        {
129 130
            nodeList[j * n] = nodes[4 + 3 * (n - 2) + n - 2 - j];
            for (int i = 1; i < n - 1; ++i)
131
            {
132
                nodeList[j * n + i] = interior[(j - 1) * (n - 2) + (i - 1)];
133
            }
134
            nodeList[(j + 1) * n - 1] = nodes[4 + (n - 2) + j - 1];
135 136
        }
    }
137

138 139
    return nodeList;
}
140

Dave Moxey's avatar
Dave Moxey committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
/**
 * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
 * suitable for the interior of a Nektar++ triangle.
 *
 * For an example, consider a third-order triangle. This routine will produce a
 * permutation map that applies the following permutation:
 *
 * ```
 *   2                 9
 *   | \               | \
 *   7   6             7   7
 *   |     \           |     \
 *   8   9   5         4   5   6
 *   |         \       |         \
 *   0---3---4---1     0---1---2---3
 *
 *       Gmsh          tensor-product
 * ```
 *
 * We assume that Gmsh uses a recursive ordering system, so that interior nodes
 * are reordered by calling this function recursively.
 *
 * @param nodes   The integer IDs of the nodes to be reordered, in Gmsh format.
 * @param n       The number of nodes in one coordinate direction. If this is
 *                zero we assume no reordering needs to be done and return the
 *                identity permutation.
 *
 * @return The nodes vector in tensor-product ordering.
 */
170
std::vector<int> triTensorNodeOrdering(const std::vector<int> &nodes, int n)
171 172 173
{
    std::vector<int> nodeList;
    int cnt2;
174

175
    nodeList.resize(nodes.size());
176

177 178 179 180
    // Vertices
    nodeList[0] = nodes[0];
    if (n > 1)
    {
181 182
        nodeList[n - 1] = nodes[1];
        nodeList[n * (n + 1) / 2 - 1] = nodes[2];
183
    }
184

185 186
    // Edges
    int cnt = n;
187
    for (int i = 1; i < n - 1; ++i)
188
    {
189 190 191 192
        nodeList[i]               = nodes[3 + i - 1];
        nodeList[cnt]             = nodes[3 + 3 * (n - 2) - i];
        nodeList[cnt + n - i - 1] = nodes[3 + (n - 2) + i - 1];
        cnt += n - i;
193
    }
194

195 196 197 198
    // Interior (recursion)
    if (n > 3)
    {
        // Reorder interior nodes
199 200 201 202
        std::vector<int> interior((n - 3) * (n - 2) / 2);
        std::copy(
            nodes.begin() + 3 + 3 * (n - 2), nodes.end(), interior.begin());
        interior = triTensorNodeOrdering(interior, n - 3);
203 204

        // Copy into full node list
205
        cnt  = n;
206
        cnt2 = 0;
207
        for (int j = 1; j < n - 2; ++j)
208
        {
209
            for (int i = 0; i < n - j - 2; ++i)
210
            {
211
                nodeList[cnt + i + 1] = interior[cnt2 + i];
212
            }
213 214
            cnt += n - j;
            cnt2 += n - 2 - j;
215 216
        }
    }
217

218 219
    return nodeList;
}
220

221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
typedef boost::tuple<int, int, int> Mode;
struct cmpop
{
    bool operator()(Mode const &a, Mode const &b) const
    {
        if (a.get<0>() < b.get<0>())
        {
            return true;
        }
        if (a.get<0>() > b.get<0>())
        {
            return false;
        }
        if (a.get<1>() < b.get<1>())
        {
            return true;
        }
        if (a.get<1>() > b.get<1>())
        {
            return false;
        }
        if (a.get<2>() < b.get<2>())
        {
            return true;
        }

        return false;
    }
};

Dave Moxey's avatar
Dave Moxey committed
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
/**
 * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
 * suitable for the interior of a Nektar++ tetrahedron.
 *
 * For an example, consider a second order tetrahedron. This routine will
 * produce a permutation map that applies the following permutation:
 *
 * ```
 *               2                           5
 *             ,/|`\                       ,/|`\
 *           ,/  |  `\                   ,/  |  `\
 *         ,6    '.   `5               ,3    '.   `4
 *       ,/       8     `\           ,/       8     `\
 *     ,/         |       `\       ,/         |       `\
 *    0--------4--'.--------1     0--------1--'.--------2
 *     `\.         |      ,/       `\.         |      ,/
 *        `\.      |    ,9            `\.      |    ,7
 *           `7.   '. ,/                 `6.   '. ,/
 *              `\. |/                      `\. |/
 *                 `3                          `9
 *
 *             Gmsh                     tensor-product
 * ```
 *
 * We assume that Gmsh uses a recursive ordering system, so that interior nodes
 * are reordered by calling this function recursively.
 *
 * @param nodes   The integer IDs of the nodes to be reordered, in Gmsh format.
 * @param n       The number of nodes in one coordinate direction. If this is
 *                zero we assume no reordering needs to be done and return the
 *                identity permutation.
 *
 * @return The nodes vector in tensor-product ordering.
 */
285 286 287 288 289 290 291 292 293
std::vector<int> tetTensorNodeOrdering(const std::vector<int> &nodes, int n)
{
    std::vector<int> nodeList;
    int nTri = n*(n+1)/2;
    int nTet = n*(n+1)*(n+2)/6;

    nodeList.resize(nodes.size());
    nodeList[0] = nodes[0];

Dave Moxey's avatar
Dave Moxey committed
294
    if (n == 1)
295
    {
Dave Moxey's avatar
Dave Moxey committed
296
        return nodeList;
297 298
    }

Dave Moxey's avatar
Dave Moxey committed
299 300 301 302 303 304
    // Vertices
    nodeList[n - 1]    = nodes[1];
    nodeList[nTri - 1] = nodes[2];
    nodeList[nTet - 1] = nodes[3];

    if (n == 2)
305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
    {
        return nodeList;
    }

    // Set up a map that takes (a,b,c) -> m to help us figure out where things
    // are inside the tetrahedron.
    std::map<Mode, int, cmpop> tmp;

    for (int k = 0, cnt = 0; k < n; ++k)
    {
        for (int j = 0; j < n - k; ++j)
        {
            for (int i = 0; i < n - k - j; ++i)
            {
                tmp[Mode(i,j,k)] = cnt++;
            }
        }
    }

Dave Moxey's avatar
Dave Moxey committed
324
    // Edges
325 326 327 328 329 330 331 332 333 334 335
    for (int i = 1; i < n-1; ++i)
    {
        int eI = i-1;
        nodeList[tmp[Mode(i,0,0)]]     = nodes[4 + eI];
        nodeList[tmp[Mode(n-1-i,i,0)]] = nodes[4 + (n-2) + eI];
        nodeList[tmp[Mode(0,n-1-i,0)]] = nodes[4 + 2*(n-2) + eI];
        nodeList[tmp[Mode(0,0,n-1-i)]] = nodes[4 + 3*(n-2) + eI];
        nodeList[tmp[Mode(0,i,n-1-i)]] = nodes[4 + 4*(n-2) + eI];
        nodeList[tmp[Mode(i,0,n-1-i)]] = nodes[4 + 5*(n-2) + eI];
    }

Dave Moxey's avatar
Dave Moxey committed
336 337 338 339 340 341 342 343
    if (n == 3)
    {
        return nodeList;
    }

    // For faces, we use the triTensorNodeOrdering routine to make our lives
    // slightly easier.
    int nFacePts = (n-3)*(n-2)/2;
344

Dave Moxey's avatar
Dave Moxey committed
345 346 347 348 349 350 351 352 353 354 355 356 357 358
    // Grab face points and reorder into a tensor-product type format
    vector<vector<int> > tmpNodes(4);
    int offset = 4 + 6*(n-2);

    for (int i = 0; i < 4; ++i)
    {
        tmpNodes[i].resize(nFacePts);
        for (int j = 0; j < nFacePts; ++j)
        {
            tmpNodes[i][j] = nodes[offset++];
        }
        tmpNodes[i] = triTensorNodeOrdering(tmpNodes[i], n-3);
    }

Dave Moxey's avatar
Dave Moxey committed
359
    if (n > 4)
Dave Moxey's avatar
Dave Moxey committed
360
    {
Dave Moxey's avatar
Dave Moxey committed
361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387
        // Now align faces
        vector<int> triVertId(3), toAlign(3);
        triVertId[0] = 0;
        triVertId[1] = 1;
        triVertId[2] = 2;

        // Faces 0,2: triangle verts {0,2,1} --> {0,1,2}
        HOTriangle<int> hoTri(triVertId, tmpNodes[0]);
        toAlign[0] = 0;
        toAlign[1] = 2;
        toAlign[2] = 1;

        hoTri.Align(toAlign);
        tmpNodes[0] = hoTri.surfVerts;

        hoTri.surfVerts = tmpNodes[2];
        hoTri.Align(toAlign);
        tmpNodes[2] = hoTri.surfVerts;

        // Face 3: triangle verts {1,2,0} --> {0,1,2}
        toAlign[0] = 1;
        toAlign[1] = 2;
        toAlign[2] = 0;

        hoTri.surfVerts = tmpNodes[3];
        hoTri.Align(toAlign);
        tmpNodes[3] = hoTri.surfVerts;
Dave Moxey's avatar
Dave Moxey committed
388 389
    }

Dave Moxey's avatar
Dave Moxey committed
390 391
    // Now apply faces. Note that faces 3 and 2 are swapped between Gmsh and
    // Nektar++ order.
Dave Moxey's avatar
Dave Moxey committed
392 393 394 395
    for (int j = 1, cnt = 0; j < n-2; ++j)
    {
        for (int i = 1; i < n-j-1; ++i, ++cnt)
        {
Dave Moxey's avatar
Dave Moxey committed
396 397 398 399
            nodeList[tmp[Mode(i,j,0)]]       = tmpNodes[0][cnt];
            nodeList[tmp[Mode(i,0,j)]]       = tmpNodes[1][cnt];
            nodeList[tmp[Mode(n-1-i-j,i,j)]] = tmpNodes[3][cnt];
            nodeList[tmp[Mode(0,i,j)]]       = tmpNodes[2][cnt];
Dave Moxey's avatar
Dave Moxey committed
400 401 402
        }
    }

Dave Moxey's avatar
Dave Moxey committed
403
    if (n == 4)
Dave Moxey's avatar
Dave Moxey committed
404
    {
Dave Moxey's avatar
Dave Moxey committed
405
        return nodeList;
Dave Moxey's avatar
Dave Moxey committed
406 407
    }

Dave Moxey's avatar
Dave Moxey committed
408
    // Finally, recurse on interior volume
409
    vector<int> intNodes, tmpInt;
Dave Moxey's avatar
Dave Moxey committed
410
    for (int i = offset; i < nTet; ++i)
Dave Moxey's avatar
Dave Moxey committed
411
    {
Dave Moxey's avatar
Dave Moxey committed
412
        intNodes.push_back(nodes[i]);
Dave Moxey's avatar
Dave Moxey committed
413
    }
414
    tmpInt = tetTensorNodeOrdering(intNodes, n-4);
Dave Moxey's avatar
Dave Moxey committed
415

Dave Moxey's avatar
Dave Moxey committed
416
    for (int k = 1, cnt = 0; k < n - 2; ++k)
Dave Moxey's avatar
Dave Moxey committed
417
    {
Dave Moxey's avatar
Dave Moxey committed
418 419 420 421
        for (int j = 1; j < n - k - 1; ++j)
        {
            for (int i = 1; i < n - k - j - 1; ++i)
            {
422
                nodeList[tmp[Mode(i,j,k)]] = tmpInt[cnt++];
Dave Moxey's avatar
Dave Moxey committed
423 424
            }
        }
Dave Moxey's avatar
Dave Moxey committed
425 426
    }

427 428 429
    return nodeList;
}

Dave Moxey's avatar
Dave Moxey committed
430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463
/**
 * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
 * suitable for the interior of a Nektar++ prism. This routine is specifically
 * designed for *interior* Gmsh nodes only.
 *
 * Prisms are a bit of a special case, in that interior nodes are heterogeneous
 * in the number of points they use. As an example, for a second-order interior
 * of a fourth-order prism, this routine calculates the mapping taking us
 *
 * ```
 *               ,6                        ,8
 *            .-' | \                   .-' | \
 *         ,-3    |   \              ,-5    |   \
 *      ,-'  | \  |     \         ,-'  | \  |     \
 *     2     |   \7-------8      2     |   \6-------7
 *     | \   | ,-' \   ,-'       | \   | ,-' \   ,-'
 *     |   \,5------,0'          |   \,3------,4'
 *     |,-'  \   ,-'             |,-'  \   ,-'
 *     1-------4'                0-------1'
 *
 *             Gmsh                 tensor-product
 * ```
 *
 * @todo Make this routine work for higher orders. The Gmsh ordering seems
 *       a little different than other elements, therefore the functionality is
 *       (for now) hard-coded to orders less than 5.
 *
 * @param nodes   The integer IDs of the nodes to be reordered, in Gmsh format.
 * @param n       The number of nodes in one coordinate direction. If this is
 *                zero we assume no reordering needs to be done and return the
 *                identity permutation.
 *
 * @return The nodes vector in tensor-product ordering.
 */
464 465 466 467 468 469 470 471 472 473 474
std::vector<int> prismTensorNodeOrdering(const std::vector<int> &nodes, int n)
{
    std::vector<int> nodeList;

    if (n == 0)
    {
        return nodeList;
    }

    nodeList.resize(nodes.size());

475
    if (n == 2)
476 477 478 479 480 481
    {
        nodeList[0] = nodes[1];
        nodeList[1] = nodes[0];
        return nodeList;
    }

482 483 484 485 486 487 488
    // For some reason, this ordering is different. Whereas Gmsh usually orders
    // vertices first, followed by edges, this ordering goes VVE-VVE-VVE for the
    // three edges that contain edge-interior information, but also vertex
    // orientation is not the same as the original Gmsh prism. The below is done
    // by looking at the ordering by hand and is hence why we don't support
    // order > 4 right now.
    if (n == 3)
489
    {
490 491 492 493 494 495 496 497 498
        nodeList[0] = nodes[1];
        nodeList[1] = nodes[4];
        nodeList[2] = nodes[2];
        nodeList[3] = nodes[5];
        nodeList[4] = nodes[0];
        nodeList[5] = nodes[3];
        nodeList[6] = nodes[7];
        nodeList[7] = nodes[8];
        nodeList[8] = nodes[6];
499 500
    }

501 502 503
    ASSERTL0(n < 4, "Prism Gmsh input and output is incomplete for orders "
                    "larger than 4");

504 505 506
    return nodeList;
}

Dave Moxey's avatar
Dave Moxey committed
507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539
/**
 * @brief Reorder Gmsh nodes so that they appear in a "tensor product" format
 * suitable for the interior of a Nektar++ hexahedron.
 *
 * For an example, consider a second order hexahedron. This routine will produce
 * a permutation map that applies the following permutation:
 *
 * ```
 *    3----13----2         6----7-----8
 *    |\         |\        |\         |\
 *    |15    24  | 14      |15    16  | 18
 *    9  \ 20    11 \      3  \  4    5  \
 *    |   7----19+---6     |   24---25+---26
 *    |22 |  26  | 23|     |12 |  13  | 14|
 *    0---+-8----1   |     0---+-1----2   |
 *     \ 17    25 \  18     \ 21    22 \  23
 *     10 |  21    12|       9 |  10    11|
 *       \|         \|        \|         \|
 *        4----16----5         18---19----20
 *
 *          Gmsh            tensor-product
 * ```
 *
 * We assume that Gmsh uses a recursive ordering system, so that interior nodes
 * are reordered by calling this function recursively.
 *
 * @param nodes   The integer IDs of the nodes to be reordered, in Gmsh format.
 * @param n       The number of nodes in one coordinate direction. If this is
 *                zero we assume no reordering needs to be done and return the
 *                identity permutation.
 *
 * @return The nodes vector in tensor-product ordering.
 */
Dave Moxey's avatar
Dave Moxey committed
540 541
std::vector<int> hexTensorNodeOrdering(const std::vector<int> &nodes, int n)
{
542
    int i, j, k;
Dave Moxey's avatar
Dave Moxey committed
543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566
    std::vector<int> nodeList;

    nodeList.resize(nodes.size());
    nodeList[0] = nodes[0];

    if (n == 1)
    {
        return nodeList;
    }

    // Vertices: same order as Nektar++
    nodeList[n - 1]               = nodes[1];
    nodeList[n*n -1]              = nodes[2];
    nodeList[n*(n-1)]             = nodes[3];
    nodeList[n*n*(n-1)]           = nodes[4];
    nodeList[n - 1 + n*n*(n-1)]   = nodes[5];
    nodeList[n*n -1 + n*n*(n-1)]  = nodes[6];
    nodeList[n*(n-1) + n*n*(n-1)] = nodes[7];

    if (n == 2)
    {
        return nodeList;
    }

567
    int hexEdges[12][2] = {
568 569 570 571 572
        { 0, 1 }, { n-1, n }, { n*n-1, -1 }, { n*(n-1), -n },
        { 0, n*n }, { n-1, n*n }, { n*n - 1, n*n }, { n*(n-1), n*n },
        { n*n*(n-1), 1 }, { n*n*(n-1) + n-1, n }, { n*n*n-1, -1 },
        { n*n*(n-1) + n*(n-1), -n }
    };
573
    int hexFaces[6][3] = {
574 575 576
        { 0, 1, n }, { 0, 1, n*n }, { n-1, n, n*n },
        { n*(n-1), 1, n*n }, { 0, n, n*n }, { n*n*(n-1), 1, n }
    };
577
    int gmshToNekEdge[12] = {0, -3, 4, 1, 5, 2, 6, 7, 8, -11, 9, 10};
578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601

    // Edges
    int offset = 8;
    for (int i = 0; i < 12; ++i)
    {
        int e = abs(gmshToNekEdge[i]);

        if (gmshToNekEdge[i] >= 0)
        {
            for (int j = 1; j < n-1; ++j)
            {
                nodeList[hexEdges[e][0] + j*hexEdges[e][1]] = nodes[offset++];
            }
        }
        else
        {
            for (int j = 1; j < n-1; ++j)
            {
                nodeList[hexEdges[e][0] + (n-j-1)*hexEdges[e][1]] = nodes[offset++];
            }
        }
    }

    // Faces
602
    int gmsh2NekFace[6] = {0, 1, 4, 2, 3, 5};
603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665

    // Map which defines orientation between Gmsh and Nektar++ faces.
    StdRegions::Orientation faceOrient[6] = {
        StdRegions::eDir1FwdDir2_Dir2FwdDir1,
        StdRegions::eDir1FwdDir1_Dir2FwdDir2,
        StdRegions::eDir1FwdDir2_Dir2FwdDir1,
        StdRegions::eDir1FwdDir1_Dir2FwdDir2,
        StdRegions::eDir1BwdDir1_Dir2FwdDir2,
        StdRegions::eDir1FwdDir1_Dir2FwdDir2};

    for (i = 0; i < 6; ++i)
    {
        int n2 = (n-2)*(n-2);
        int face = gmsh2NekFace[i];
        offset   = 8 + 12 * (n-2) + i * n2;

        // Create a list of interior face nodes for this face only.
        vector<int> faceNodes(n2);
        for (j = 0; j < n2; ++j)
        {
            faceNodes[j] = nodes[offset + j];
        }

        // Now get the reordering of this face, which puts Gmsh
        // recursive ordering into Nektar++ row-by-row order.
        faceNodes = quadTensorNodeOrdering(faceNodes, n-2);
        vector<int> tmp(n2);

        // Finally reorient the face according to the geometry
        // differences.
        if (faceOrient[i] == StdRegions::eDir1FwdDir1_Dir2FwdDir2)
        {
            // Orientation is the same, just copy.
            tmp = faceNodes;
        }
        else if (faceOrient[i] == StdRegions::eDir1FwdDir2_Dir2FwdDir1)
        {
            // Tranposed faces
            for (j = 0; j < n-2; ++j)
            {
                for (k = 0; k < n-2; ++k)
                {
                    tmp[j * (n-2) + k] = faceNodes[k * (n-2) + j];
                }
            }
        }
        else if (faceOrient[i] == StdRegions::eDir1BwdDir1_Dir2FwdDir2)
        {
            for (j = 0; j < n-2; ++j)
            {
                for (k = 0; k < n-2; ++k)
                {
                    tmp[j * (n-2) + k] = faceNodes[j * (n-2) + (n - k - 3)];
                }
            }
        }

        // Now put this into the right place in the output array
        for (k = 1; k < n-1; ++k)
        {
            for (j = 1; j < n-1; ++j)
            {
                nodeList[hexFaces[face][0] + j*hexFaces[face][1] + k*hexFaces[face][2]]
666
                    = faceNodes[(k-1)*(n-2) + j-1];
667 668 669 670 671
            }
        }
    }

    // Finally, recurse on interior volume
672
    vector<int> intNodes, tmpInt;
673 674 675 676 677
    for (int i = 8 + 12 * (n-2) + 6 * (n-2) * (n-2); i < n*n*n; ++i)
    {
        intNodes.push_back(nodes[i]);
    }

678
    if (intNodes.size())
679
    {
680 681
        tmpInt = hexTensorNodeOrdering(intNodes, n-2);
        for (int k = 1, cnt = 0; k < n - 1; ++k)
682
        {
683
            for (int j = 1; j < n - 1; ++j)
684
            {
685 686 687 688
                for (int i = 1; i < n - 1; ++i)
                {
                    nodeList[i + j * n + k * n * n] = tmpInt[cnt++];
                }
689 690 691
            }
        }
    }
Dave Moxey's avatar
Dave Moxey committed
692 693 694 695 696

    return nodeList;
}


697 698 699 700 701 702 703
/**
 * @brief Set up InputGmsh object.
 *
 */
InputGmsh::InputGmsh(MeshSharedPtr m) : InputModule(m)
{
}
704

705 706 707
InputGmsh::~InputGmsh()
{
}
708

709
/**
Dave Moxey's avatar
Dave Moxey committed
710 711 712 713 714
 * Gmsh file contains a list of nodes and their coordinates, along with a list
 * of elements and those nodes which define them. We read in and store the list
 * of nodes in #m_node and store the list of elements in #m_element. Each new
 * element is supplied with a list of entries from m_node which defines the
 * element. Finally some mesh statistics are printed.
715 716 717 718 719 720 721 722
 *
 * @param   pFilename           Filename of Gmsh file to read.
 */
void InputGmsh::Process()
{
    // Open the file stream.
    OpenStream();

723
    m_mesh->m_expDim   = 0;
724 725 726 727
    m_mesh->m_spaceDim = 0;
    string line;
    int nVertices = 0;
    int nEntities = 0;
728 729 730
    int elm_type  = 0;
    int prevId    = -1;
    int maxTagId  = -1;
731 732 733 734 735 736 737 738 739 740 741 742 743

    map<unsigned int, ElmtConfig>::iterator it;

    // This map takes each element ID and maps it to a permutation map
    // that is required to take Gmsh element node orderings and map them
    // to Nektar++ orderings.
    boost::unordered_map<int, vector<int> > orderingMap;
    boost::unordered_map<int, vector<int> >::iterator oIt;

    if (m_mesh->m_verbose)
    {
        cout << "InputGmsh: Start reading file..." << endl;
    }
744

745 746 747 748 749 750
    while (!m_mshFile.eof())
    {
        getline(m_mshFile, line);
        stringstream s(line);
        string word;
        s >> word;
751

752 753
        // Process nodes.
        if (word == "$Nodes")
754
        {
755 756 757 758 759
            getline(m_mshFile, line);
            stringstream s(line);
            s >> nVertices;
            int id = 0;
            for (int i = 0; i < nVertices; ++i)
760
            {
761
                getline(m_mshFile, line);
762 763 764
                stringstream st(line);
                double x = 0, y = 0, z = 0;
                st >> id >> x >> y >> z;
765

766
                if ((x * x) > 0.000001 && m_mesh->m_spaceDim < 1)
767
                {
768
                    m_mesh->m_spaceDim = 1;
769
                }
770
                if ((y * y) > 0.000001 && m_mesh->m_spaceDim < 2)
771
                {
772 773 774 775 776
                    m_mesh->m_spaceDim = 2;
                }
                if ((z * z) > 0.000001 && m_mesh->m_spaceDim < 3)
                {
                    m_mesh->m_spaceDim = 3;
777 778
                }

779
                id -= 1; // counter starts at 0
780

781
                if (id - prevId != 1)
782 783 784 785 786
                {
                    cerr << "Gmsh vertex ids should be contiguous" << endl;
                    abort();
                }
                prevId = id;
787 788
                m_mesh->m_node.push_back(
                    boost::shared_ptr<Node>(new Node(id, x, y, z)));
789 790 791 792 793 794 795 796 797
            }
        }
        // Process elements
        else if (word == "$Elements")
        {
            getline(m_mshFile, line);
            stringstream s(line);
            s >> nEntities;
            for (int i = 0; i < nEntities; ++i)
798
            {
799 800 801
                getline(m_mshFile, line);
                stringstream st(line);
                int id = 0, num_tag = 0, num_nodes = 0;
802

803 804
                st >> id >> elm_type >> num_tag;
                id -= 1; // counter starts at 0
805

806 807
                it = elmMap.find(elm_type);
                if (it == elmMap.end())
808
                {
809 810 811
                    cerr << "Error: element type " << elm_type
                         << " not supported" << endl;
                    abort();
812 813
                }

814 815 816
                // Read element tags
                vector<int> tags;
                for (int j = 0; j < num_tag; ++j)
817
                {
818 819 820
                    int tag = 0;
                    st >> tag;
                    tags.push_back(tag);
821
                }
822 823 824 825 826 827 828 829
                tags.resize(1);

                maxTagId = max(maxTagId, tags[0]);

                // Read element node list
                vector<NodeSharedPtr> nodeList;
                num_nodes = GetNnodes(elm_type);
                for (int k = 0; k < num_nodes; ++k)
830
                {
831 832 833 834
                    int node = 0;
                    st >> node;
                    node -= 1; // counter starts at 0
                    nodeList.push_back(m_mesh->m_node[node]);
835 836
                }

837 838 839 840 841
                // Look up reordering.
                oIt = orderingMap.find(elm_type);

                // If it's not created, then create it.
                if (oIt == orderingMap.end())
842
                {
843
                    oIt = orderingMap.insert(
844
                        make_pair(elm_type, CreateReordering(elm_type))).first;
845 846
                }

847 848
                // Apply reordering map where necessary.
                if (oIt->second.size() > 0)
849
                {
850
                    vector<int> &mapping      = oIt->second;
851 852
                    vector<NodeSharedPtr> tmp = nodeList;
                    for (int i = 0; i < mapping.size(); ++i)
853
                    {
854
                        nodeList[i] = tmp[mapping[i]];
855 856 857
                    }
                }

858
                // Create element
859 860
                ElementSharedPtr E = GetElementFactory().CreateInstance(
                    it->second.m_e, it->second, nodeList, tags);
861 862

                // Determine mesh expansion dimension
863 864
                if (E->GetDim() > m_mesh->m_expDim)
                {
865 866 867 868
                    m_mesh->m_expDim = E->GetDim();
                }
                m_mesh->m_element[E->GetDim()].push_back(E);
            }
869
        }
870
    }
871
    m_mshFile.reset();
872

873 874 875 876
    // Go through element and remap tags if necessary.
    map<int, map<LibUtilities::ShapeType, int> > compMap;
    map<int, map<LibUtilities::ShapeType, int> >::iterator cIt;
    map<LibUtilities::ShapeType, int>::iterator sIt;
877

878 879
    for (int i = 0; i < m_mesh->m_element[m_mesh->m_expDim].size(); ++i)
    {
880
        ElementSharedPtr el          = m_mesh->m_element[m_mesh->m_expDim][i];
881
        LibUtilities::ShapeType type = el->GetConf().m_e;
882

883
        vector<int> tags = el->GetTagList();
884
        int tag          = tags[0];
885

886
        cIt = compMap.find(tag);
887

888 889 890 891
        if (cIt == compMap.end())
        {
            compMap[tag][type] = tag;
            continue;
892
        }
893

894 895 896 897 898 899 900 901 902 903
        // Reset tag for this element.
        sIt = cIt->second.find(type);
        if (sIt == cIt->second.end())
        {
            maxTagId++;
            cIt->second[type] = maxTagId;
            tags[0] = maxTagId;
            el->SetTagList(tags);
        }
        else if (sIt->second != tag)
904
        {
905 906 907 908
            tags[0] = sIt->second;
            el->SetTagList(tags);
        }
    }
909

910 911 912 913 914 915 916 917 918
    bool printInfo = false;
    for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
    {
        if (cIt->second.size() > 1)
        {
            printInfo = true;
            break;
        }
    }
919

920 921
    if (printInfo)
    {
922
        cout << "Multiple elements in composite detected; remapped:" << endl;
923 924 925
        for (cIt = compMap.begin(); cIt != compMap.end(); ++cIt)
        {
            if (cIt->second.size() > 1)
926
            {
927
                sIt = cIt->second.begin();
928
                cout << "- Tag " << cIt->first << " => " << sIt->second << " ("
929 930
                     << LibUtilities::ShapeTypeMap[sIt->first] << ")";
                sIt++;
931

932 933 934
                for (; sIt != cIt->second.end(); ++sIt)
                {
                    cout << ", " << sIt->second << " ("
935
                         << LibUtilities::ShapeTypeMap[sIt->first] << ")";
936
                }
937

938 939
                cout << endl;
            }
940
        }
941
    }
942

943
    // Process rest of mesh.
944 945 946 947
    ProcessVertices();
    ProcessEdges();
    ProcessFaces();
    ProcessElements();
948 949
    ProcessComposites();
}
950

951 952 953 954 955 956 957
/**
 * For a given msh ID, return the corresponding number of nodes.
 */
int InputGmsh::GetNnodes(unsigned int InputGmshEntity)
{
    int nNodes;
    map<unsigned int, ElmtConfig>::iterator it;
958

959
    it = elmMap.find(InputGmshEntity);
960

961 962 963 964 965
    if (it == elmMap.end())
    {
        cerr << "Unknown element type " << InputGmshEntity << endl;
        abort();
    }
966

967
    switch (it->second.m_e)
968 969
    {
        case LibUtilities::ePoint:
970
            nNodes = Point::GetNumNodes(it->second);
971 972
            break;
        case LibUtilities::eSegment:
973
            nNodes = Line::GetNumNodes(it->second);
974 975
            break;
        case LibUtilities::eTriangle:
976
            nNodes = Triangle::GetNumNodes(it->second);
977 978 979 980 981
            break;
        case LibUtilities::eQuadrilateral:
            nNodes = Quadrilateral::GetNumNodes(it->second);
            break;
        case LibUtilities::eTetrahedron:
982 983
            nNodes = Tetrahedron::GetNumNodes(it->second);
            it->second.m_faceCurveType = LibUtilities::eNodalTriEvenlySpaced;
984 985
            break;
        case LibUtilities::ePyramid:
986
            nNodes = Pyramid::GetNumNodes(it->second);
987 988
            break;
        case LibUtilities::ePrism:
989
            nNodes = Prism::GetNumNodes(it->second);
990 991
            break;
        case LibUtilities::eHexahedron:
992
            nNodes = Hexahedron::GetNumNodes(it->second);
993 994 995 996 997 998
            break;
        default:
            cerr << "Unknown element type!" << endl;
            abort();
            break;
    }
999

1000 1001
    return nNodes;
}
1002

1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013
/**
 * @brief Create a reordering map for a given element.
 *
 * Since Gmsh and Nektar++ have different vertex, edge and face
 * orientations, we need to reorder the nodes in a Gmsh MSH file so that
 * they work with the Nektar++ orderings, since this is what is used in
 * the elements defined in the converter.
 */
vector<int> InputGmsh::CreateReordering(unsigned int InputGmshEntity)
{
    map<unsigned int, ElmtConfig>::iterator it;
1014

1015
    it = elmMap.find(InputGmshEntity);
1016

1017 1018 1019 1020 1021
    if (it == elmMap.end())
    {
        cerr << "Unknown element type " << InputGmshEntity << endl;
        abort();
    }
1022

1023 1024
    // For specific elements, call the appropriate function to perform
    // the renumbering.
1025
    switch (it->second.m_e)
1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041
    {
        case LibUtilities::eTriangle:
            return TriReordering(it->second);
            break;
        case LibUtilities::eQuadrilateral:
            return QuadReordering(it->second);
            break;
        case LibUtilities::eTetrahedron:
            return TetReordering(it->second);
            break;
        case LibUtilities::ePrism:
            return PrismReordering(it->second);
            break;
        case LibUtilities::eHexahedron:
            return HexReordering(it->second);
            break;
Michael Turner's avatar
Michael Turner committed
1042 1043
        case LibUtilities::eSegment:
            return LineReordering(it->second);
1044 1045 1046
        default:
            break;
    }
1047

1048 1049 1050 1051
    // Default: no reordering.
    vector<int> returnVal;
    return returnVal;
}
1052

Michael Turner's avatar
Michael Turner committed
1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065
vector<int> InputGmsh::LineReordering(ElmtConfig conf)
{
    const int order = conf.m_order;

    vector<int> mapping(order+1);
    for(int i = 0; i < order+1; i++)
    {
        mapping[i] = i;
    }

    return mapping;
}

1066 1067 1068 1069 1070 1071
/**
 * @brief Create a reordering for triangles.
 */
vector<int> InputGmsh::TriReordering(ElmtConfig conf)
{
    const int order = conf.m_order;
1072
    const int n     = order - 1;
1073

1074 1075 1076 1077 1078 1079
    // Copy vertices.
    vector<int> mapping(3);
    for (int i = 0; i < 3; ++i)
    {
        mapping[i] = i;
    }
1080

1081 1082 1083 1084
    if (order == 1)
    {
        return mapping;
    }
1085

1086
    // Curvilinear edges.
1087
    mapping.resize(3 + 3 * n);
1088

1089
    for (int i = 3; i < 3 + 3 * n; ++i)
1090 1091 1092
    {
        mapping[i] = i;
    }
1093

1094 1095 1096 1097
    if (!conf.m_faceNodes)
    {
        return mapping;
    }
1098

1099
    // Interior nodes.
1100
    vector<int> interior(n * (n - 1) / 2);
1101 1102
    for (int i = 0; i < interior.size(); ++i)
    {
1103
        interior[i] = i + 3 + 3 * n;
1104
    }
1105

1106 1107
    if (interior.size() > 0)
    {
1108
        interior = triTensorNodeOrdering(interior, n - 1);
1109
    }
1110

1111 1112 1113
    mapping.insert(mapping.end(), interior.begin(), interior.end());
    return mapping;
}
1114

1115 1116 1117 1118 1119 1120
/**
 * @brief Create a reordering for quadrilaterals.
 */
vector<int> InputGmsh::QuadReordering(ElmtConfig conf)
{
    const int order = conf.m_order;
1121
    const int n     = order - 1;
1122

1123 1124 1125 1126 1127 1128
    // Copy vertices.
    vector<int> mapping(4);
    for (int i = 0; i < 4; ++i)
    {
        mapping[i] = i;
    }
1129

1130 1131 1132 1133
    if (order == 1)
    {
        return mapping;
    }
1134

1135
    // Curvilinear edges.
1136
    mapping.resize(4 + 4 * n);
1137

1138
    for (int i = 4; i < 4 + 4 * n; ++i)
1139 1140 1141
    {
        mapping[i] = i;
    }
1142

1143 1144 1145 1146
    if (!conf.m_faceNodes)
    {
        return mapping;
    }
1147

1148
    // Interior nodes.
1149
    vector<int> interior(n * n);
1150 1151
    for (int i = 0; i < interior.size(); ++i)
    {
1152
        interior[i] = i + 4 + 4 * n;
1153
    }
1154

1155 1156 1157 1158 1159 1160 1161
    if (interior.size() > 0)
    {
        interior = quadTensorNodeOrdering(interior, n);
    }
    mapping.insert(mapping.end(), interior.begin(), interior.end());
    return mapping;
}
1162

1163 1164 1165 1166 1167 1168
/**
 * @brief Create a reordering for tetrahedra.
 */
vector<int> InputGmsh::TetReordering(ElmtConfig conf)
{
    const int order = conf.m_order;
1169 1170
    const int n     = order - 1;
    const int n2    = n * (n - 1) / 2;
1171

1172 1173 1174 1175 1176 1177 1178 1179
    int i, j;
    vector<int> mapping(4);

    // Copy vertices.
    for (i = 0; i < 4; ++i)
    {
        mapping[i] = i;
    }
1180

1181 1182 1183 1184
    if (order == 1)
    {
        return mapping;
    }
1185

1186
    // Curvilinear edges.
1187
    mapping.resize(4 + 6 * n);
1188

1189
    // Curvilinear edges.
1190 1191
    static int gmshToNekEdge[6] = {0, 1, 2, 3, 5, 4};
    static int gmshToNekRev[6]  = {0, 0, 1, 1, 1, 1};
1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202

    // Reorder edges.
    int offset, cnt = 4;
    for (i = 0; i < 6; ++i)
    {
        offset = 4 + n * gmshToNekEdge[i];

        if (gmshToNekRev[i])
        {
            for (int j = 0; j < n; ++j)
            {
1203
                mapping[offset + n - j - 1] = cnt++;
1204 1205 1206 1207 1208 1209
            }
        }
        else
        {
            for (int j = 0; j < n; ++j)
            {
1210
                mapping[offset + j] = cnt++;
1211
            }
1212 1213
        }
    }
1214

1215 1216 1217 1218 1219 1220
    if (conf.m_faceNodes == false || n2 == 0)
    {
        return mapping;
    }

    // Curvilinear faces.
1221
    mapping.resize(4 + 6 * n + 4 * n2);
1222

1223
    static int gmshToNekFace[4] = {0, 1, 3, 2};
1224 1225

    vector<int> triVertId(3);
1226 1227 1228
    triVertId[0] = 0;
    triVertId[1] = 1;
    triVertId[2] = 2;
1229 1230 1231 1232 1233

    // Loop over Gmsh faces
    for (i = 0; i < 4; ++i)
    {
        int face    = gmshToNekFace[i];
1234 1235
        int offset2 = 4 + 6 * n + i * n2;
        offset      = 4 + 6 * n + face * n2;
1236 1237 1238

        // Create a list of interior face nodes for this face only.
        vector<int> faceNodes(n2);
1239
        vector<int> toAlign(3);
1240 1241 1242
        for (j = 0; j < n2; ++j)
        {
            faceNodes[j] = offset2 + j;
1243 1244
        }

1245 1246
        // Now get the reordering of this face, which puts Gmsh
        // recursive ordering into Nektar++ row-by-row order.
1247
        vector<int> tmp = triTensorNodeOrdering(faceNodes, n - 1);
1248 1249 1250 1251 1252 1253
        HOTriangle<int> hoTri(triVertId, tmp);

        // Apply reorientation
        if (i == 0 || i == 2)
        {
            // Triangle verts {0,2,1} --> {0,1,2}
1254 1255 1256
            toAlign[0] = 0;
            toAlign[1] = 2;
            toAlign[2] = 1;
1257 1258 1259
            hoTri.Align(toAlign);
        }
        else if (i == 3)
1260
        {
1261
            // Triangle verts {1,2,0} --> {0,1,2}
1262 1263 1264
            toAlign[0] = 1;
            toAlign[1] = 2;
            toAlign[2] = 0;
1265 1266
            hoTri.Align(toAlign);