Implement ExodusII in GridIn.

Implement ExodusII in GridIn.

The tests were all generated by hand. For future reference I included the code here.

Here is the basic 2D test (five quadrilaterals):

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Create a disk mesh with five quadrilaterals.
//
//  ____________________________
// |\                          /|
// | \                        / |
// |  \           3          /  |
// |   \                    /   |
// |    \                  /    |
// |     \                /     |
// |      \______________/      |
// |      |              |      |
// |      |              |      |
// |      |              |      |
// |   4  |      5       |   2  |
// |      |              |      |
// |      |______________|      |
// |      /              \      |
// |     /                \     |
// |    /                  \    |
// |   /                    \   |
// |  /          1           \  |
// | /                        \ |
// |/                          \|
// -----------------------------
//
// element 1 has external side sets 0, 1
// element 2 has external side sets 0, 1, 10
// element 3 has no external side sets
// element 4 has no external side set 0
//
// The inner element is in block 10 - the rest are in block 1.
int main()
{
    const int n_elems = 5;
    const int n_nodes = 8;
    int elements_1[n_elems - 1][4] = {
    {1, 2, 3, 4},
    {3, 2, 8, 7},
    {5, 7, 8, 6},
    {1, 4, 5, 6}};
    int elements_10[1][4] = {{4, 3, 7, 5}};

    double xs[n_nodes] = {0, 7, 5, 2, 2, 0, 5, 7};
    double ys[n_nodes] = {0, 0, 2, 2, 5, 7, 5, 7};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("five.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 2;
    const int n_blocks = 2;
    const int n_node_sets = 0;
    const int n_side_sets = 3;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, nullptr);
    assert(ierr == 0);

    // element blocks:
    ierr = ex_put_elem_block(exid, 1, "quad", 4, 4, 0);
    assert(ierr == 0);
    ierr = ex_put_elem_block(exid, 10, "quad", 1, 4, 0);
    assert(ierr == 0);

    ierr = ex_put_elem_conn(exid, 10, elements_10);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 0, 3, 0);
    assert(ierr == 0);
    int sideset_0_elems[] = {1, 2, 4};
    int sideset_0_sides[] = {1, 2, 4};
    ierr = ex_put_side_set(exid, 0, sideset_0_elems, sideset_0_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 1, 2, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 2};
    int sideset_1_sides[] = {1, 2};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 10, 1, 0);
    assert(ierr == 0);
    int sideset_10_elems[] = {2};
    int sideset_10_sides[] = {2};
    ierr = ex_put_side_set(exid, 10, sideset_10_elems, sideset_10_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

A 3D mixed mesh:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

int main()
{
    const int n_elems = 4;
    const int n_nodes = 13;
    int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][5] = {{5, 6, 7, 8, 9}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 0;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex", 1, 8, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid", 1, 5, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Check that we can truncate higher-order elements correctly (like PYRAMID14 ->
our five-point pyramid)

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

int main()
{
    const int n_elems = 4;
    const int n_nodes = 29;
    int elements_1[1][20] = {{1, 2, 3, 4, 5, 6, 7, 8, 17, 14, 15, 16, 21, 18, 19, 20, 25, 22, 23, 24}};
    // int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][13] = {{5, 6, 7, 8, 9, 25, 22, 23, 24, 26, 27, 28, 29}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8, 3, 1.5, 0, 1.5, 3, 3, 0, 0, 3, 1.5, 0, 1.5, 0.5, 2, 2, 0.5};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1, 1.5, 3, 1.5, 0, 0, 3, 3, 0, 1.5, 3, 1.5, 0, 0.5, 0.5, 2, 2};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5, 0, 0, 0, 0, 1.5, 1.5, 1.5, 1.5, 3, 3, 3, 3, 4, 4, 4, 4};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four-p.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 0;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex20", 1, 20, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid13", 1, 13, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Get things working with a nonzero codimension:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Create a disk mesh with five quadrilaterals.
//
//  ____________________________
// |\                          /|
// | \                        / |
// |  \           3          /  |
// |   \                    /   |
// |    \                  /    |
// |     \                /     |
// |      \______________/      |
// |      |              |      |
// |      |              |      |
// |      |              |      |
// |   4  |      5       |   2  |
// |      |              |      |
// |      |______________|      |
// |      /              \      |
// |     /                \     |
// |    /                  \    |
// |   /                    \   |
// |  /          1           \  |
// | /                        \ |
// |/                          \|
// -----------------------------
//
// element 1 has external side sets 0, 1
// element 2 has external side sets 0, 1, 10
// element 3 has no external side sets
// element 4 has no external side set 0
//
// The inner element is in block 10 - the rest are in block 1.
int main()
{
    const int n_elems = 5;
    const int n_nodes = 8;
    int elements_1[n_elems - 1][4] = {
    {1, 2, 3, 4},
    {1, 4, 5, 6},
    {2, 8, 7, 3},
    {6, 5, 7, 8}};
    int elements_10[1][4] = {{4, 3, 7, 5}};

    double xs[n_nodes] = {0, 7, 5, 2, 2, 0, 5, 7};
    double ys[n_nodes] = {0, 0, 2, 2, 5, 7, 5, 7};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("five-3d.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = 2;
    const int n_node_sets = 0;
    const int n_side_sets = 3;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, xs);
    assert(ierr == 0);

    // element blocks:
    ierr = ex_put_elem_block(exid, 1, "quad", 4, 4, 0);
    assert(ierr == 0);
    ierr = ex_put_elem_block(exid, 10, "quad", 1, 4, 0);
    assert(ierr == 0);

    ierr = ex_put_elem_conn(exid, 10, elements_10);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 0, 3, 0);
    assert(ierr == 0);
    int sideset_0_elems[] = {1, 2, 4};
    int sideset_0_sides[] = {3, 4, 6};
    ierr = ex_put_side_set(exid, 0, sideset_0_elems, sideset_0_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 1, 2, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 2};
    int sideset_1_sides[] = {3, 4};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 10, 1, 0);
    assert(ierr == 0);
    int sideset_10_elems[] = {2};
    int sideset_10_sides[] = {4};
    ierr = ex_put_side_set(exid, 10, sideset_10_elems, sideset_10_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}

Check sidesets in 3D:

#include <exodusII.h>

#include <algorithm>
#include <cassert>
#include <cstring>

// Like writer-2 but with side sets.

int main()
{
    const int n_elems = 4;
    const int n_nodes = 13;
    int elements_1[1][8] = {{1, 2, 3, 4, 5, 6, 7, 8}};
    int elements_2[1][5] = {{5, 6, 7, 8, 9}};
    int elements_3[1][4] = {{6, 10, 7, 9}};
    int elements_4[1][6] = {{10, 6, 7, 13, 11, 12}};

    double xs[n_nodes] = {0, 3, 3, 0, 0, 3, 3, 0, 1, 3, 8, 8, 8};
    double ys[n_nodes] = {0, 0, 3, 3, 0, 0, 3, 3, 1, 1, 0, 3, 1};
    double zs[n_nodes] = {0, 0, 0, 0, 3, 3, 3, 3, 5, 5, 3, 3, 5};

    int sizeof_data = sizeof(double);
    int sizeof_geom = sizeof(double);
    const int exid = ex_create("four-sidesets.e", EX_CLOBBER, &sizeof_data, &sizeof_geom);

    const int dim = 3;
    const int n_blocks = n_elems;
    const int n_node_sets = 0;
    const int n_side_sets = 5;
    int ierr = ex_put_init(exid, "mytitle", dim, n_nodes, n_elems, n_blocks,
                           n_node_sets, n_side_sets);
    assert(ierr == 0);

    // nodes:
    ierr = ex_put_coord(exid, xs, ys, zs);
    assert(ierr == 0);

    // block 0 (hex):
    ierr = ex_put_elem_block(exid, 1, "hex", 1, 8, 0);
    assert(ierr == 0);
    // block 1 (pyramid):
    ierr = ex_put_elem_block(exid, 10, "pyramid", 1, 5, 0);
    assert(ierr == 0);
    // block 2 (tet):
    ierr = ex_put_elem_block(exid, 100, "tet", 1, 4, 0);
    assert(ierr == 0);
    // block 3 (wedge):
    ierr = ex_put_elem_block(exid, 1000, "wedge", 1, 6, 0);
    assert(ierr == 0);

    // connectivity:
    ierr = ex_put_elem_conn(exid, 1, elements_1);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 10, elements_2);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 100, elements_3);
    assert(ierr == 0);
    ierr = ex_put_elem_conn(exid, 1000, elements_4);
    assert(ierr == 0);

    // side sets:
    ierr = ex_put_side_set_param(exid, 1, 5, 0);
    assert(ierr == 0);
    int sideset_1_elems[] = {1, 1, 2, 3, 4};
    int sideset_1_sides[] = {1, 5, 1, 1, 1};
    ierr = ex_put_side_set(exid, 1, sideset_1_elems, sideset_1_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 2, 3, 0);
    assert(ierr == 0);
    int sideset_2_elems[] = {1, 3, 4};
    int sideset_2_sides[] = {2, 2, 2};
    ierr = ex_put_side_set(exid, 2, sideset_2_elems, sideset_2_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 3, 3, 0);
    assert(ierr == 0);
    int sideset_3_elems[] = {1, 2, 4};
    int sideset_3_sides[] = {3, 3, 3};
    ierr = ex_put_side_set(exid, 3, sideset_3_elems, sideset_3_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 4, 2, 0);
    assert(ierr == 0);
    int sideset_4_elems[] = {1, 2};
    int sideset_4_sides[] = {4, 4};
    ierr = ex_put_side_set(exid, 4, sideset_4_elems, sideset_4_sides);
    assert(ierr == 0);

    ierr = ex_put_side_set_param(exid, 5, 4, 0);
    assert(ierr == 0);
    int sideset_5_elems[] = {1, 1, 2, 4};
    int sideset_5_sides[] = {4, 5, 4, 5};
    ierr = ex_put_side_set(exid, 5, sideset_5_elems, sideset_5_sides);
    assert(ierr == 0);

    ierr = ex_close(exid);
    assert(ierr == 0);
}