#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#if defined(_WIN32) && !defined(__NUTC__)
# include <io.h>
#else
# include <unistd.h>
#endif
#include "cgnslib.h"

#ifndef CGNS_VERSION
# ifdef NofValidSimulationTypes
#  if NofValidModelTypes > 20
#   define CGNS_VERSION 2100
#  else
#   define CGNS_VERSION 2000
#  endif
# else
#  if NofValidElementTypes > 21
#   define CGNS_VERSION 1200 /* could be 1270 */
#  else
#   define CGNS_VERSION 1100
#  endif
# endif
#endif

#if NofValidBCTypes < 26
# define FamilySpecified BCTypeUserDefined
#endif

#if CGNS_VERSION >= 2100
# define USE_LINKS
#endif

#define TWOPI 6.2831853

/* set STRUCTURED_FACES to write the structured zone bocos for the
   mixed case as face range and face list instead of as point range/list
#define STRUCTURED_FACES
*/

/* set UNSTRUCTURED_1TO1 to write the unstructured case
   zone connectivities as 1to1 instead of abutting1to1
#define UNSTRUCTURED_1TO1
*/

/* set ABUTTING1TO1_FACES to write the unstructured case
   zone connectivites as abutting1to1 with faces(elements)
   instead of points. This also writes the mixed case zone 1
   to zone 2 connectivity with a face range
#define ABUTTING1TO1_FACES
*/

int CellDim = 3, PhyDim = 3;

int cgfile, cgbase, cgzone;
int CellDim, PhyDim, size[9];

#define NUM_SIDE 5
#define NODE_INDEX(I,J,K) ((I)+NUM_SIDE*(((J)-1)+NUM_SIDE*((K)-1)))
#define CELL_INDEX(I,J,K) ((I)+(NUM_SIDE-1)*(((J)-1)+(NUM_SIDE-1)*((K)-1)))

int num_coord;
float *xcoord, *ycoord, *zcoord;
int num_element, *elements;
int num_face, *faces, *parent;

int max_sol;
float *solution;

int npts, *pts, *d_pts;
float *interp;

char errmsg[128];

void init_data();
void write_structured(), write_unstructured();
void write_mixed(), write_mismatched();

void error_exit (char *where)
{
    fprintf (stderr, "ERROR:%s:%s\n", where, cg_get_error());
    exit (1);
}

int main ()
{
    init_data();
    unlink("test.cgns");
    if (cg_open("test.cgns", MODE_WRITE, &cgfile)) error_exit("cg_open");
    write_structured();
    write_unstructured();
    write_mixed();
    write_mismatched();
    if (cg_close(cgfile)) error_exit("cg_close");
    return 0;
}

void init_data()
{
    int n, i, j, k, nn, nf, np;

    /* compute coordinates - make it twice as big for use with cylindrical */

    num_coord = NUM_SIDE * NUM_SIDE * NUM_SIDE;
    xcoord = (float *) malloc (6 * num_coord * sizeof(float));
    if (NULL == xcoord) {
        fprintf(stderr, "malloc failed for coordinates\n");
        exit(1);
    }
    ycoord = xcoord + 2 * num_coord;
    zcoord = ycoord + 2 * num_coord;
    for (n = 0, k = 0; k < NUM_SIDE; k++) {
        for (j = 0; j < NUM_SIDE; j++) {
            for (i = 0; i < NUM_SIDE; i++, n++) {
                xcoord[n] = (float)i;
                ycoord[n] = (float)j;
                zcoord[n] = (float)k;
            }
        }
    }

    /* create solution vector large enough for grid + 1 rind plane */

    max_sol = (NUM_SIDE + 2) * (NUM_SIDE + 2) * (NUM_SIDE + 2);
    solution = (float *) malloc (max_sol * sizeof(float));
    if (NULL == solution) {
        fprintf(stderr, "malloc failed for solution\n");
        exit(1);
    }
    for (n = 0; n < max_sol; n++)
        solution[n] = (float)(n + 1);

    /* compute elements */

    num_element = (NUM_SIDE - 1) * (NUM_SIDE - 1) * (NUM_SIDE - 1);
    elements = (int *) malloc (8 * num_element * sizeof(int));
    if (NULL == elements) {
        fprintf(stderr, "malloc failed for elements");
        exit(1);
    }
    for (n = 0, k = 1; k < NUM_SIDE; k++) {
        for (j = 1; j < NUM_SIDE; j++) {
            for (i = 1; i < NUM_SIDE; i++) {
                nn = NODE_INDEX(i, j, k);
                elements[n++] = nn;
                elements[n++] = nn + 1;
                elements[n++] = nn + 1 + NUM_SIDE;
                elements[n++] = nn + NUM_SIDE;
                nn += NUM_SIDE * NUM_SIDE;
                elements[n++] = nn;
                elements[n++] = nn + 1;
                elements[n++] = nn + 1 + NUM_SIDE;
                elements[n++] = nn + NUM_SIDE;
            }
        }
    }

    /* compute outside face elements */

    num_face = 6 * (NUM_SIDE - 1) * (NUM_SIDE - 1);
    faces = (int *) malloc (4 * num_face * sizeof(int));
    parent = (int *) malloc (4 * num_face * sizeof(int));
    if (NULL == faces || NULL == parent) {
        fprintf(stderr, "malloc failed for elements");
        exit(1);
    }
    for (n = 0; n < 4*num_face; n++)
        parent[n] = 0;
    nf = np = 0;
    n = 2 * num_face;
    i = 1;
    for (k = 1; k < NUM_SIDE; k++) {
        for (j = 1; j < NUM_SIDE; j++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + NUM_SIDE * NUM_SIDE;
            faces[nf++]  = nn + NUM_SIDE * (NUM_SIDE + 1);
            faces[nf++]  = nn + NUM_SIDE;
            parent[np]   = CELL_INDEX(i, j, k);
            parent[np+n] = 5;
            np++;
        }
    }
    i = NUM_SIDE;
    for (k = 1; k < NUM_SIDE; k++) {
        for (j = 1; j < NUM_SIDE; j++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + NUM_SIDE;
            faces[nf++]  = nn + NUM_SIDE * (NUM_SIDE + 1);
            faces[nf++]  = nn + NUM_SIDE * NUM_SIDE;
            parent[np]   = CELL_INDEX(i-1, j, k);
            parent[np+n] = 3;
            np++;
        }
    }
    j = 1;
    for (k = 1; k < NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + 1;
            faces[nf++]  = nn + 1 + NUM_SIDE * NUM_SIDE;
            faces[nf++]  = nn + NUM_SIDE * NUM_SIDE;
            parent[np]   = CELL_INDEX(i, j, k);
            parent[np+n] = 2;
            np++;
        }
    }
    j = NUM_SIDE;
    for (k = 1; k < NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + NUM_SIDE * NUM_SIDE;
            faces[nf++]  = nn + 1 + NUM_SIDE * NUM_SIDE;
            faces[nf++]  = nn + 1;
            parent[np]   = CELL_INDEX(i, j-1, k);
            parent[np+n] = 4;
            np++;
        }
    }
    k = 1;
    for (j = 1; j < NUM_SIDE; j++) {
        for (i = 1; i < NUM_SIDE; i++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + NUM_SIDE;
            faces[nf++]  = nn + NUM_SIDE + 1;
            faces[nf++]  = nn + 1;
            parent[np]   = CELL_INDEX(i, j, k);
            parent[np+n] = 1;
            np++;
        }
    }
    k = NUM_SIDE;
    for (j = 1; j < NUM_SIDE; j++) {
        for (i = 1; i < NUM_SIDE; i++) {
            nn = NODE_INDEX(i, j, k);
            faces[nf++]  = nn;
            faces[nf++]  = nn + 1;
            faces[nf++]  = nn + NUM_SIDE + 1;
            faces[nf++]  = nn + NUM_SIDE;
            parent[np]   = CELL_INDEX(i, j, k-1);
            parent[np+n] = 6;
            np++;
        }
    }

    /* connectivity points - make it big enough to hold 4 surfaces */

    npts = NUM_SIDE * NUM_SIDE;
    pts = (int *) malloc (12 * npts * sizeof(int));
    if (NULL == pts) {
        fprintf(stderr, "malloc failed for connectivity points");
        exit(1);
    }
    d_pts = pts + 6 * npts;

    /* create interpolate data array */

    interp = (float *) malloc (6 * npts * sizeof(float));
    if (NULL == interp) {
        fprintf(stderr, "malloc failed for interpolate array");
        exit(1);
    }
}

void write_reference ()
{
    int n, i, ierr, dim = 1;
    float exps[5];

    static struct {
        char *name;
        int dim;
        int exps[5];
        float val;
    } state[] = {
        {"Mach", 0, {0, 0, 0, 0, 0}, 0.2},
        {"VelocitySound", 1, {0, 1, -1, 0, 0}, 330.0},
        {"VelocityMagnitude", 1, {0, 1, -1, 0, 0}, 66.0},
        {"VelocityUnitVectorX", 0, {0, 0, 0, 0, 0}, 1.0},
        {"VelocityUnitVectorY", 0, {0, 0, 0, 0, 0}, 0.0},
        {"VelocityUnitVectorZ", 0, {0, 0, 0, 0, 0}, 0.0},
        {"Reynolds", 0, {0, 0, 0, 0, 0}, 3.0e6},
        {"Temperature", 1, {0, 0, 0, 1, 0}, 300.0},
        {"Pressure", 1, {1, -1, -2, 0, 0}, 1.0e5},
        {"LengthReference", 1, {0, 1, 0, 0, 0}, 10.0}
    };

    if (cg_goto(cgfile, cgbase, "end") ||
        cg_state_write("reference state quantities") ||
        cg_goto(cgfile, cgbase, "ReferenceState_t", 1, "end") ||
        cg_dataclass_write(Dimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("reference state");

    for (n = 0; n < 10; n++) {
        if (cg_goto(cgfile, cgbase, "ReferenceState_t", 1, "end") ||
            cg_array_write(state[n].name, RealSingle, 1, &dim, &state[n].val) ||
            cg_goto(cgfile, cgbase, "ReferenceState_t", 1,
                "DataArray_t", n+1, "end")) {
            sprintf (errmsg, "reference state data %d", n+1);
            error_exit(errmsg);
        }
        if (state[n].dim) {
            for (i = 0; i < 5; i++)
                exps[i] = (float)state[n].exps[i];
            ierr = cg_exponents_write(RealSingle, exps);
        }
        else
            ierr = cg_dataclass_write(NondimensionalParameter);
        if (ierr) {
            sprintf (errmsg, "reference state data %d dataclass", n+1);
            error_exit(errmsg);
        }
    }
}

void write_equationset ()
{
    int n, diff[6], dim = 1;
    float g = 1.4;
    float R = 53.352;
    float ts = 110.6;
    float mu = 1.716e-5;
    float exp = 0.666;
    float pt = 0.9;

    for (n = 0; n < 6; n++)
        diff[n] = 0;
    diff[2] = 1;

    /* flow equation set */

    if (cg_goto(cgfile, cgbase, "end") ||
        cg_equationset_write (3) ||
        cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1, "end") ||
        cg_governing_write(NSTurbulent) ||
        cg_model_write("GasModel_t", Ideal) ||
        cg_model_write("ViscosityModel_t", SutherlandLaw) ||
        cg_model_write("ThermalConductivityModel_t", PowerLaw) ||
        cg_model_write("TurbulenceClosure_t", EddyViscosity) ||
        cg_model_write("TurbulenceModel_t", Algebraic_BaldwinLomax))
        error_exit("flow equation set");

#if CGNS_VERSION >= 2000
    /* not allowed prior to 2.0 */
    if (cg_dataclass_write(Dimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("flow equation set dataclass");
#endif

    /* diffusion model */

    if (cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "GoverningEquations_t", 1, "end") ||
        cg_diffusion_write(diff) ||
        cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "TurbulenceModel_t", 1, "end") ||
        cg_diffusion_write(diff))
        error_exit("diffusion model");

    /* gas model */

    if (cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "GasModel_t", 1, "end") ||
        cg_dataclass_write(DimensionlessConstant) ||
        cg_array_write("SpecificHeatRatio", RealSingle, 1, &dim, &g) ||
        cg_array_write("IdealGasConstant", RealSingle, 1, &dim, &R) ||
        cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "GasModel_t", 1, "DataArray_t", 2, "end") ||
        cg_dataclass_write(Dimensional) ||
        cg_units_write(Slug, Foot, Second, Rankine, Radian))
        error_exit("gas model");

    /* viscosity model */

    if (cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "ViscosityModel_t", 1, "end") ||
        cg_array_write("SutherlandLawConstant", RealSingle, 1, &dim, &ts) ||
        cg_array_write("ViscosityMolecularReference", RealSingle, 1, &dim, &mu))
        error_exit("viscosity model");

    /* thermal conductivity model */

    if (cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "ThermalConductivityModel_t", 1, "end") ||
        cg_array_write("PowerLawExponent", RealSingle, 1, &dim, &exp) ||
        cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "ThermalConductivityModel_t", 1, "DataArray_t", 1, "end") ||
        cg_dataclass_write(DimensionlessConstant))
        error_exit("thermal conductivity model");

    /* turbulence closure model */

    if (cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "TurbulenceClosure_t", 1, "end") ||
        cg_array_write("PrandtlTurbulent", RealSingle, 1, &dim, &pt) ||
        cg_goto(cgfile, cgbase, "FlowEquationSet_t", 1,
            "TurbulenceClosure_t", 1, "DataArray_t", 1, "end") ||
        cg_dataclass_write(DimensionlessConstant))
        error_exit("turbulence closure model");
}

void write_coords(int nz)
{
    int k, nn, n, nij, koff, cgcoord;

    koff = nz == 1 ? 1 - NUM_SIDE : 0;
    nij = NUM_SIDE * NUM_SIDE;
    for (n = 0, k = 0; k < NUM_SIDE; k++) {
        for (nn = 0; nn < nij; nn++)
            zcoord[n++] = (float)(k + koff);
    }

    if (cg_coord_write(cgfile, cgbase, nz, RealSingle,
                       "CoordinateX", xcoord, &cgcoord) ||
        cg_coord_write(cgfile, cgbase, nz, RealSingle,
                       "CoordinateY", ycoord, &cgcoord) ||
        cg_coord_write(cgfile, cgbase, nz, RealSingle,
                       "CoordinateZ", zcoord, &cgcoord)) {
        sprintf (errmsg, "zone %d coordinates", nz);
        error_exit(errmsg);
    }
}

void write_elements(int nz)
{
    int cgsect;

    if (cg_section_write(cgfile, cgbase, nz, "Elements", HEXA_8,
            1, num_element, 0, elements, &cgsect) ||
        cg_section_write(cgfile, cgbase, nz, "Faces", QUAD_4,
            num_element+1, num_element+num_face, 0, faces, &cgsect) ||
        cg_parent_data_write(cgfile, cgbase, nz, cgsect, parent)) {
        sprintf (errmsg, "zone %d elements", nz);
        error_exit(errmsg);
    }
}

#ifdef USE_LINKS

void write_zone_link(int nz, char *basename, char *nodename)
{
    char pathname[128];

    sprintf(pathname, "/%s/Zone%d/%s", basename, nz, nodename);
    if (cg_goto(cgfile, cgbase, "Zone_t", nz, "end") ||
        cg_link_write(nodename, "", pathname)) {
        sprintf (errmsg, "zone %d link", nz);
        error_exit(errmsg);
    }
}

#endif

/*------------------------------------------------------------------------
 * structured grid
 *------------------------------------------------------------------------*/

void write_structured()
{
    int i, j, k, n, cgconn, cgbc, cgfam, cggeo, cgpart;
    int range[6], d_range[6], transform[3];
    int cgsol, rind[6], cgfld;
    char name[33];

    printf ("writing structured base\n");
    fflush (stdout);

    if (cg_base_write(cgfile, "Structured", CellDim, PhyDim, &cgbase) ||
        cg_goto(cgfile, cgbase, "end") ||
        cg_descriptor_write("Descriptor", "Multi-block Structured Grid") ||
        cg_dataclass_write(Dimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("structured base");
#if CGNS_VERSION >= 2000
    if (cg_simulation_type_write(cgfile, cgbase, NonTimeAccurate))
        error_exit("simulation type");
#endif

    write_reference();
    write_equationset();

    /* write zones */

    for (n = 0; n < 3; n++) {
        size[n]   = NUM_SIDE;
        size[n+3] = NUM_SIDE - 1;
        size[n+6] = 0;
    }
    for (n = 1; n <= 2; n++) {
        sprintf(name, "Zone%d", n);
        if (cg_zone_write(cgfile, cgbase, name, size, Structured, &cgzone)) {
            sprintf (errmsg, "structured zone %d", n);
            error_exit(errmsg);
        }
        write_coords(n);
    }

    /* write zone 1 to zone 2 connectivity as 1to1 */

    for (n = 0; n < 3; n++) {
        range[n] = d_range[n] = 1;
        range[n+3] = d_range[n+3] = NUM_SIDE;
        transform[n] = n + 1;
    }
    range[2] = NUM_SIDE;
    d_range[5] = 1;
    if (cg_1to1_write(cgfile, cgbase, 1, "1to1 -> Zone2", "Zone2",
            range, d_range, transform, &cgconn))
        error_exit("1to1->zone2");

    /* write zone 2 to zone 1 connectivity as Abutting1to1 */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE;
    }
    range[5] = 1;
    for (n = 0, j = 1; j <= NUM_SIDE; j++) {
        for (i = 1; i <= NUM_SIDE; i++) {
            d_pts[n++] = i;
            d_pts[n++] = j;
            d_pts[n++] = 1;
        }
    }
    if (cg_conn_write(cgfile, cgbase, 2, "Abutting1to1 -> Zone1",
            Vertex, Abutting1to1, PointRange, 2, range, "Zone1",
            Structured, PointListDonor, Integer, npts, d_pts, &cgconn))
        error_exit("abutting1to1->zone1");

    /* write inlet BC (zone 1) as point range */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE;
    }
    range[5] = 1;
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            PointRange, 2, range, &cgbc))
        error_exit("inlet boco");

    /* write outlet BC (zone 2) as point list */

    for (n = 0, j = 1; j <= NUM_SIDE; j++) {
        for (i = 1; i <= NUM_SIDE; i++) {
            pts[n++] = i;
            pts[n++] = j;
            pts[n++] = NUM_SIDE;
        }
    }
    if (cg_boco_write(cgfile, cgbase, 2, "Outlet", BCOutflow,
            PointList, npts, pts, &cgbc))
        error_exit("outlet boco");

    /* write zone 1 wall BC as point ranges using a family to group them */

    if (cg_family_write(cgfile, cgbase, "WallFamily", &cgfam) ||
        cg_fambc_write(cgfile, cgbase, cgfam, "WallBC", BCWall, &cgbc))
        error_exit("wall family bc");

    /* write out some bogus geometry info for the family */

    if (cg_geo_write(cgfile, cgbase, cgfam, "Geometry",
            "geometry.file", "CADsystem", &cggeo) ||
        cg_part_write(cgfile, cgbase, cgfam, cggeo, "imin part", &cgpart) ||
        cg_part_write(cgfile, cgbase, cgfam, cggeo, "imax part", &cgpart) ||
        cg_part_write(cgfile, cgbase, cgfam, cggeo, "jmin part", &cgpart) ||
        cg_part_write(cgfile, cgbase, cgfam, cggeo, "jmax part", &cgpart))
        error_exit("wall family parts");

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE;
    }

    range[3] = 1;
    if (cg_boco_write(cgfile, cgbase, 1, "imin", FamilySpecified,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_famname_write("WallFamily"))
        error_exit("imin boco");

    range[0] = range[3] = NUM_SIDE;
    if (cg_boco_write(cgfile, cgbase, 1, "imax", FamilySpecified,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_famname_write("WallFamily"))
        error_exit("imax boco");

    range[0] = range[4] = 1;
    if (cg_boco_write(cgfile, cgbase, 1, "jmin", FamilySpecified,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_famname_write("WallFamily"))
        error_exit("jmin boco");

    range[1] = range[4] = NUM_SIDE;
    if (cg_boco_write(cgfile, cgbase, 1, "jmax", FamilySpecified,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_famname_write("WallFamily"))
        error_exit("jmax boco");

    /* write zone 2 wall BC as point list */

    for (n = 0, k = 1; k <= NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = i + 1;
            pts[n++] = 1;
            pts[n++] = k;
            pts[n++] = i;
            pts[n++] = NUM_SIDE;
            pts[n++] = k;
        }
        for (j = 1; j < NUM_SIDE; j++) {
            pts[n++] = 1;
            pts[n++] = j;
            pts[n++] = k;
            pts[n++] = NUM_SIDE;
            pts[n++] = j+1;
            pts[n++] = k;
        }
    }
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            PointList, n / 3, pts, &cgbc))
        error_exit("zone 2 walls boco");

    /* write solution for zone 1 as vertex with rind points */

    for (n = 0; n < 6; n++)
        rind[n] = 1;
    if (cg_sol_write(cgfile, cgbase, 1, "VertexSolution", Vertex, &cgsol) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "FlowSolution_t", cgsol, "end") ||
        cg_rind_write(rind) ||
        cg_field_write(cgfile, cgbase, 1, cgsol, RealSingle,
            "Density", solution, &cgfld))
        error_exit("zone 1 solution");

    /* write solution for zone 2 as cell center with rind points */

    rind[0] = rind[1] = 0;
    if (cg_sol_write(cgfile, cgbase, 2, "CellCenterSolution",
            CellCenter, &cgsol) ||
        cg_goto(cgfile, cgbase, "Zone_t", 2, "FlowSolution_t", cgsol, "end") ||
        cg_rind_write(rind) ||
        cg_field_write(cgfile, cgbase, 2, cgsol, RealSingle,
            "Density", solution, &cgfld))
        error_exit("zone 2 solution");
}

/*------------------------------------------------------------------------
 * unstructured grid
 *------------------------------------------------------------------------*/

void write_unstructured()
{
    int n, nelem, cgconn, cgbc;
    int range[2];
    int cgsol, cgfld;
    char name[33];
#ifdef UNSTRUCTURED_1TO1
    int d_range[2], transform;
#else
# ifdef ABUTTING1TO1_FACES
    GridLocation_t location;
# else
    int i, j;
# endif
#endif

    printf ("writing unstructured base\n");
    fflush (stdout);

    if (cg_base_write(cgfile, "Unstructured", CellDim, PhyDim, &cgbase) ||
        cg_goto(cgfile, cgbase, "end") ||
        cg_descriptor_write("Descriptor", "Multi-block Unstructured Grid") ||
        cg_dataclass_write(NormalizedByDimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("unstructured base");

    /* write zones */

    for (n = 0; n < 9; n++)
        size[n] = 0;
    size[0] = num_coord;
    size[1] = num_element;
    for (n = 1; n <= 2; n++) {
        sprintf(name, "Zone%d", n);
        if (cg_zone_write(cgfile, cgbase, name, size, Unstructured, &cgzone)) {
            sprintf (errmsg, "unstructured zone %d", n);
            error_exit(errmsg);
        }
        write_coords(n);
        write_elements(n);
    }
    nelem = (NUM_SIDE - 1) * (NUM_SIDE - 1);

#ifdef UNSTRUCTURED_1TO1

    /* write connectivities as 1to1 */

    range[0] = NODE_INDEX(1, 1, NUM_SIDE);
    range[1] = NODE_INDEX(NUM_SIDE, NUM_SIDE, NUM_SIDE);
    d_range[0] = NODE_INDEX(1, 1, 1);
    d_range[1] = NODE_INDEX(NUM_SIDE, NUM_SIDE, 1);
    transform = 1;
    if (cg_1to1_write(cgfile, cgbase, 1, "1to1 -> Zone2", "Zone2",
            range, d_range, &transform, &cgconn))
        error_exit("1to1->zone2");

    if (cg_1to1_write(cgfile, cgbase, 2, "1to1 -> Zone1", "Zone1",
            d_range, range, &transform, &cgconn))
        error_exit("1to1->zone1");

#else
# ifdef ABUTTING1TO1_FACES

    /* zone 1 to zone 2 connectivity as Abutting1to1 with element range */

    range[0] = num_element + num_face - nelem + 1;
    range[1] = num_element + num_face;
    for (n = 0; n < nelem; n++)
        d_pts[n] = range[0] + n - nelem;

#if CGNS_VERSION > 2100
    location = FaceCenter;

    /* this fail for version prior to 2.2 - see below */
    if (cg_conn_write(cgfile, cgbase, 1, "Abutting1to1 -> Zone2",
            location, Abutting1to1, PointRange, 2, range, "Zone2",
            Unstructured, PointListDonor, Integer, nelem, d_pts, &cgconn))
        error_exit("face center abutting1to1->zone2");
#else
    location = CellCenter;

    for (n = 0; n < nelem; n++)
        pts[n] = range[0] + n;
    if (cg_conn_write(cgfile, cgbase, 1, "Abutting1to1 -> Zone2",
            location, Abutting1to1, PointList, nelem, pts, "Zone2",
            Unstructured, PointListDonor, Integer, nelem, d_pts, &cgconn))
        error_exit("face center abutting1to1->zone2");
#endif

    /* zone 2 to zone 1 connectivity as Abutting1to1 with element list */

    for (n = 0; n < nelem; n++) {
        pts[n]   = num_element + num_face - 2 * nelem + 1 + n;
        d_pts[n] = pts[n] + nelem;
    }
    if (cg_conn_write(cgfile, cgbase, 2, "Abutting1to1 -> Zone1",
            location, Abutting1to1, PointList, nelem, pts, "Zone1",
            Unstructured, PointListDonor, Integer, nelem, d_pts, &cgconn))
        error_exit("face center abutting1to1->zone1");

# else

    /* zone 1 to zone 2 connectivity as Abutting1to1 with point range */

    range[0] = NODE_INDEX(1, 1, NUM_SIDE);
    range[1] = NODE_INDEX(NUM_SIDE, NUM_SIDE, NUM_SIDE);
    for (n = 0, j = 1; j <= NUM_SIDE; j++) {
        for (i = 1; i <= NUM_SIDE; i++)
            d_pts[n++] = NODE_INDEX(i, j, 1);
    }
    if (cg_conn_write(cgfile, cgbase, 1, "Abutting1to1 -> Zone2",
            Vertex, Abutting1to1, PointRange, 2, range, "Zone2",
            Unstructured, PointListDonor, Integer, npts, d_pts, &cgconn))
        error_exit("point range abutting1to1->zone2");

    /* zone 2 to zone 1 connectivity as Abutting1to1 with point list */

    for (n = 0, j = 1; j <= NUM_SIDE; j++) {
        for (i = 1; i <= NUM_SIDE; i++) {
            pts[n] = d_pts[n];
            d_pts[n++] = NODE_INDEX(i, j, NUM_SIDE);
        }
    }
    if (cg_conn_write(cgfile, cgbase, 2, "Abutting1to1 -> Zone1",
            Vertex, Abutting1to1, PointList, npts, pts, "Zone1",
            Unstructured, PointListDonor, Integer, npts, d_pts, &cgconn))
        error_exit("point list abutting1to1->zone1");

# endif
#endif

    /* write inlet BC (zone 1) and outlet BC (zone 2) as element lists
       and zone 1 and zone 2 walls as element range */

    range[0] = num_element + 1;
    range[1] = num_element + 4 * nelem;
    for (n = 0; n < nelem; n++) {
        pts[n]   = num_element + num_face - 2 * nelem + 1 + n;
        d_pts[n] = pts[n] + nelem;
    }

    /* in order to write an ElementList/Range for versions 1.27 -> 2.2 */
    /* you need to use PointList/Range with a goto and GridLocation */
    /* For versions prior to 1.27, GridLocation is not allowed under BC_t */

#if CGNS_VERSION >= 1270 && CGNS_VERSION <= 2200
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            PointList, nelem, pts, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter inlet boco");
    if (cg_boco_write(cgfile, cgbase, 2, "Outlet", BCOutflow,
            PointList, nelem, d_pts, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 2, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter outlet boco");

    /* a bug in versions prior to 2.1 causes this to fail since a */
    /* check is performed on the input with the vertex range and */
    /* the element range is outside of this range - this check */
    /* was commented  out in later versions of the library */
# if CGNS_VERSION >= 2100
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter zone 1 walls boco");
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 2, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter zone 2 walls boco");
# else
    /* use list instead of range */
    for (nelem = 0, n = range[0]; n <= range[1]; n++)
        pts[nelem++] = n;
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            PointList, nelem, pts, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter zone 1 walls boco");
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            PointList, nelem, pts, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 2, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit ("facecenter zone 2 walls boco");
# endif
    /* for version < 1.27 and > 2.2, just use ElementList/Range */
#else
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            ElementList, nelem, pts, &cgbc))
        error_exit ("elementlist inlet boco");
    if (cg_boco_write(cgfile, cgbase, 2, "Outlet", BCOutflow,
            ElementList, nelem, d_pts, &cgbc))
        error_exit ("elementlist outlet boco");
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            ElementRange, 2, range, &cgbc))
        error_exit ("elementrange zone 1 walls boco");
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            ElementRange, 2, range, &cgbc))
        error_exit ("elementrange zone 2 walls boco");
#endif

    /* write solution for zone 1 as vertex and zone 2 as cell center */

    if (cg_sol_write(cgfile, cgbase, 1, "VertexSolution", Vertex, &cgsol) ||
        cg_field_write(cgfile, cgbase, 1, cgsol, RealSingle,
            "Density", solution, &cgfld))
        error_exit("zone 1 solution");
    if (cg_sol_write(cgfile, cgbase, 2, "CellCenterSolution",
            CellCenter, &cgsol) ||
        cg_field_write(cgfile, cgbase, 2, cgsol, RealSingle,
            "Density", solution, &cgfld))
        error_exit("zone 2 solution");
}

/*------------------------------------------------------------------------
 * mixed grid
 *------------------------------------------------------------------------*/

void write_mixed()
{
    int i, j, k, n, nelem, cgconn, cgbc;
    int range[6];
#ifdef ABUTTING1TO1_FACES
    GridLocation_t location;
#endif

    printf ("writing mixed base\n");
    fflush (stdout);

    if (cg_base_write(cgfile, "Mixed", CellDim, PhyDim, &cgbase) ||
        cg_goto(cgfile, cgbase, "end") ||
        cg_descriptor_write("Descriptor",
            "Mixed Structured and Unstructured Grid") ||
        cg_dataclass_write(Dimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("mixed base");

    /* zone 1 is structured */

    for (n = 0; n < 3; n++) {
        size[n]   = NUM_SIDE;
        size[n+3] = NUM_SIDE - 1;
        size[n+6] = 0;
    }
    if (cg_zone_write(cgfile, cgbase, "StructuredZone", size,
        Structured, &cgzone))
        error_exit("structured zone");
#ifdef USE_LINKS
    write_zone_link(1, "Structured", "GridCoordinates");
#else
    write_coords(1);
#endif

    /* zone 2 is unstructured */

    for (n = 0; n < 9; n++)
        size[n] = 0;
    size[0] = num_coord;
    size[1] = num_element;
    if (cg_zone_write(cgfile, cgbase, "UnstructuredZone", size,
        Unstructured, &cgzone))
        error_exit("unstructured zone");
#ifdef USE_LINKS
    write_zone_link(2, "Unstructured", "GridCoordinates");
    write_zone_link(2, "Unstructured", "Elements");
    write_zone_link(2, "Unstructured", "Faces");
#else
    write_coords(2);
    write_elements(2);
#endif
    nelem = (NUM_SIDE - 1) * (NUM_SIDE - 1);

#ifdef ABUTTING1TO1_FACES

    /* zone 1 -> zone 2 connectivity as face range */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE - 1;
    }
    range[2] = range[5] = NUM_SIDE;
    for (n = 0; n < nelem; n++)
        d_pts[n] = num_element + num_face - 2 * nelem + 1 + n;

    /* this fails with 1.2->2.1, since a check is done on the */
    /* element range, and the face number falls outside this range */
    /* Use PointList instead, since no check is done */
#if CGNS_VERSION >= 1200 && CGNS_VERSION <= 2100
    for (n = 0, j = 1; j < NUM_SIDE; j++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = i;
            pts[n++] = j;
            pts[n++] = NUM_SIDE;
        }
    }
    if (cg_conn_write(cgfile, cgbase, 1, "Structured -> Unstructured",
        CellCenter, Abutting1to1, PointList, nelem, pts, "UnstructuredZone",
        Unstructured, PointListDonor, Integer, nelem, d_pts, &cgconn))
        error_exit("abutting1to1->unstructured face list");
#else
# if CGNS_VERSION < 2200
    location = CellCenter;
# elif CGNS_VERSION < 2300
    location = FaceCenter;
# else
    location = KFaceCenter;
# endif
    if (cg_conn_write(cgfile, cgbase, 1, "Structured -> Unstructured",
        location, Abutting1to1, PointRange, 2, range, "UnstructuredZone",
        Unstructured, PointListDonor, Integer, nelem, d_pts, &cgconn))
        error_exit("abutting1to1->unstructured face range");
#endif

#else

    /* zone 1 -> zone 2 connectivity as point range */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE;
    }
    range[2] = NUM_SIDE;
    for (n = 0; n < npts; n++)
        d_pts[n] = n + 1;
    if (cg_conn_write(cgfile, cgbase, 1, "Structured -> Unstructured",
            Vertex, Abutting1to1, PointRange, 2, range, "UnstructuredZone",
            Unstructured, PointListDonor, Integer, npts, d_pts, &cgconn))
        error_exit("abutting1to1->unstructured point range");

#endif

    /* zone 2 -> zone 1 connectivity as point range (k=1 surface) */

    range[0] = 1;
    range[1] = npts;
    for (n = 0, j = 1; j <= NUM_SIDE; j++) {
        for (i = 1; i <= NUM_SIDE; i++) {
            d_pts[n++] = i;
            d_pts[n++] = j;
            d_pts[n++] = NUM_SIDE;
        }
    }
    if (cg_conn_write(cgfile, cgbase, 2, "Unstructured -> Structured",
            Vertex, Abutting1to1, PointRange, 2, range, "StructuredZone",
            Structured, PointListDonor, Integer, npts, d_pts, &cgconn))
        error_exit("abutting1to1->structured point range");

#ifdef STRUCTURED_FACES

    /* write inlet BC (zone 1) as face range */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE - 1;
    }
    range[5] = 1;
#if CGNS_VERSION < 1270
    /* gridlocation not supported */
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            ElementRange, 2, range, &cgbc))
        error_exit("element range inlet boco");
#else
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            PointRange, 2, range, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(KFaceCenter))
        error_exit("kfacecenter inlet boco");
#endif

#else

    /* write inlet BC (zone 1) as point range */

    for (n = 0; n < 3; n++) {
        range[n] = 1;
        range[n+3] = NUM_SIDE;
    }
    range[5] = 1;
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            PointRange, 2, range, &cgbc))
        error_exit("point range inlet boco");

#endif

    /* write outlet BC (zone 2) as point range */

    range[0] = num_coord - npts + 1;
    range[1] = num_coord;
    if (cg_boco_write(cgfile, cgbase, 2, "Outlet", BCOutflow,
            PointRange, 2, range, &cgbc))
        error_exit("point range outlet boco");

#ifdef STRUCTURED_FACES

    /* write zone 1 walls as face list */

    for (n = 0, k = 1; k < NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = i;
            pts[n++] = 1;
            pts[n++] = k;
            pts[n++] = i;
            pts[n++] = NUM_SIDE;
            pts[n++] = k;
        }
        for (j = 1; j < NUM_SIDE; j++) {
            pts[n++] = 1;
            pts[n++] = j;
            pts[n++] = k;
            pts[n++] = NUM_SIDE;
            pts[n++] = j;
            pts[n++] = k;
        }
    }
#if CGNS_VERSION < 1270
    /* gridlocation not supported */
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            ElementList, 4 * nelem, pts, &cgbc))
        error_exit("element list zone 1 walls boco");
#else
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            PointList, 4 * nelem, pts, &cgbc) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneBC_t", 1,
            "BC_t", cgbc, "end") ||
        cg_gridlocation_write(FaceCenter))
        error_exit("face list zone 1 walls boco");
#endif

#else

    /* write zone 1 wall BC as point list */

    for (n = 0, k = 1; k <= NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = i + 1;
            pts[n++] = 1;
            pts[n++] = k;
            pts[n++] = i;
            pts[n++] = NUM_SIDE;
            pts[n++] = k;
        }
        for (j = 1; j < NUM_SIDE; j++) {
            pts[n++] = 1;
            pts[n++] = j;
            pts[n++] = k;
            pts[n++] = NUM_SIDE;
            pts[n++] = j+1;
            pts[n++] = k;
        }
    }
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            PointList, n / 3, pts, &cgbc))
        error_exit("point list zone 1 walls boco");

#endif

    /* write zone 2 walls as point list */

    for (n = 0, k = 1; k <= NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = NODE_INDEX(i + 1, 1, k);
            pts[n++] = NODE_INDEX(i, NUM_SIDE, k);
        }
        for (j = 1; j < NUM_SIDE; j++) {
            pts[n++] = NODE_INDEX(1, j, k);
            pts[n++] = NODE_INDEX(NUM_SIDE, j + 1, k);
        }
    }
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            PointList, n, pts, &cgbc))
        error_exit("point list zone 2 walls boco");
}

/*------------------------------------------------------------------------
 * mismatched zones with cylindrical coordinates
 *------------------------------------------------------------------------*/

void write_mismatched()
{
    int i, j, k, n, nj, cgcoord, cgconn, cgbc;
    int ic, jc, dims[2];
    float dx, dy, dtheta, r, t, x, y;
    int range[6], d_range[6], transform[3];
    PointSetType_t d_type;

    printf ("writing mismatched base\n");
    fflush (stdout);

    if (cg_base_write(cgfile, "Mismatched", CellDim, PhyDim, &cgbase) ||
        cg_goto(cgfile, cgbase, "end") ||
        cg_descriptor_write("Descriptor", "Mismatched Grid") ||
        cg_dataclass_write(Dimensional) ||
        cg_units_write(Kilogram, Meter, Second, Kelvin, Radian))
        error_exit("mismatched base");

    /* zone 1 is cartesian */

    for (n = 0; n < 3; n++) {
        size[n]   = NUM_SIDE;
        size[n+3] = NUM_SIDE - 1;
        size[n+6] = 0;
    }

    dx = 0.5 * (float)(NUM_SIDE - 1);
    dy = 0.5 * (float)(NUM_SIDE - 1);
    for (n = 0, k = 0; k < NUM_SIDE; k++) {
        for (j = 0; j < NUM_SIDE; j++) {
            for (i = 0; i < NUM_SIDE; i++, n++) {
                xcoord[n] -= dx;
                ycoord[n] -= dy;
            }
        }
    }
    if (cg_zone_write(cgfile, cgbase, "CartesianZone", size,
        Structured, &cgzone))
        error_exit("cartesion zone");
    write_coords(1);

    /* zone 2 is cylindrical */

    nj = 2 * NUM_SIDE;
    dtheta = (float)TWOPI / (float)(nj - 1);
    for (n = 0, k = 0; k < NUM_SIDE; k++) {
        for (j = 0; j < nj; j++) {
            for (i = 0; i < NUM_SIDE; i++, n++) {
                xcoord[n] = (float)i;
                ycoord[n] = (float)j * dtheta;
                zcoord[n] = (float)k;
            }
        }
    }
    size[1] = nj;
    size[4] = nj - 1;
    if (cg_zone_write(cgfile, cgbase, "CylindricalZone", size,
        Structured, &cgzone) ||
        cg_coord_write(cgfile, cgbase, cgzone, RealSingle,
            "CoordinateR", xcoord, &cgcoord) ||
        cg_coord_write(cgfile, cgbase, cgzone, RealSingle,
            "CoordinateTheta", ycoord, &cgcoord) ||
        cg_coord_write(cgfile, cgbase, cgzone, RealSingle,
            "CoordinateZ", zcoord, &cgcoord))
        error_exit("cylindrical zone");

    /* zone 1 -> zone 2 connectivity */

    for (n = 0, j = 0; j < NUM_SIDE; j++) {
        for (i = 0; i < NUM_SIDE; i++) {
            x = (float)i - dx;
            y = (float)j - dy;
            r = (float)sqrt (x * x + y * y);
            if (r > (float)(NUM_SIDE - 1)) continue;
            ic = (int)r;
            if (ic >= NUM_SIDE - 1) ic = NUM_SIDE - 2;
            t = (float)atan2 (y, x);
            if (t < 0.0) t += (float)TWOPI;
            jc = (int)(t / dtheta);
            if (jc >= nj - 1) jc = nj - 2;
            pts[n]      = i + 1;
            pts[n+1]    = j + 1;
            pts[n+2]    = NUM_SIDE;
            d_pts[n]    = ic + 1;
            d_pts[n+1]  = jc + 1;
            d_pts[n+2]  = 1;
            interp[n++] = r - (float)ic;
            interp[n++] = t / dtheta - (float)jc;
            interp[n++] = 0.0;
        }
    }

    /* priot to 1.27, the donor points and interpolation was */
    /* passed as a single array, which could be RealSingle, */
    /* RealDouble, or Integer. */
#if CGNS_VERSION < 1270
    for (i = 0; i < n; i++)
        interp[i] += (float)d_pts[i];
    if (cg_conn_write(cgfile, cgbase, 1, "Cartesian -> Cylindrical",
            Vertex, Abutting, PointList, n/3, pts, "CylindricalZone",
            Structured, PointListDonor, RealSingle, n/3, interp, &cgconn))
        error_exit("cartesian->cylindrical connectivity");
#else
    /* for 1.27 and later, the donor data is Integer, and the */
    /* interpolants are given as a child node data array */

    dims[0] = 3;
    dims[1] = n / 3;
    /* for 1.27, PointListDonor must be used with structured donor */
    /*       and CellListDonor with unstructured donor */
# if CGNS_VERSION < 2000
    d_type = PointListDonor;
# else
    d_type = CellListDonor;
# endif
    if (cg_conn_write(cgfile, cgbase, 1, "Cartesian -> Cylindrical",
            Vertex, Abutting, PointList, n/3, pts, "CylindricalZone",
            Structured, d_type, Integer, n/3, d_pts, &cgconn) ||
        cg_goto(cgfile, cgbase, "Zone_t", 1, "ZoneGridConnectivity_t", 1,
            "GridConnectivity_t", cgconn, "end") ||
        cg_array_write("InterpolantsDonor", RealSingle, 2, dims, interp))
        error_exit("cartesian->cylindrical connectivity");
#endif

    /* zone 2 -> zone 1 connectivity */

    for (n = 0, j = 0; j < nj; j++) {
        for (i = 0; i < NUM_SIDE; i++) {
            r = (float)i;
            t = (float)j * dtheta;
            x = r * (float)cos (t) + dx;
            y = r * (float)sin (t) + dy;
            if (x < 0.0 || x > (float)(NUM_SIDE - 1) ||
                y < 0.0 || y > (float)(NUM_SIDE - 1)) continue;
            ic = (int)x;
            if (ic >= NUM_SIDE - 1) ic = NUM_SIDE - 2;
            jc = (int)y;
            if (jc >= NUM_SIDE - 1) jc = NUM_SIDE - 2;
            pts[n]      = i + 1;
            pts[n+1]    = j + 1;
            pts[n+2]    = 1;
            d_pts[n]    = ic + 1;
            d_pts[n+1]  = jc + 1;
            d_pts[n+2]  = NUM_SIDE;
            interp[n++] = x - (float)ic;
            interp[n++] = y - (float)jc;
            interp[n++] = 0.0;
        }
    }

#if CGNS_VERSION < 1270
    for (i = 0; i < n; i++)
        interp[i] += (float)d_pts[i];
    if (cg_conn_write(cgfile, cgbase, 2, "Cylindrical -> Cartesian",
            Vertex, Abutting, PointList, n/3, pts, "CartesianZone",
            Structured, PointListDonor, RealSingle, n/3, interp, &cgconn))
        error_exit("cylindrical->cartesian connectivity");
#else
    dims[0] = 3;
    dims[1] = n / 3;
# if CGNS_VERSION < 2000
    d_type = PointListDonor;
# else
    d_type = CellListDonor;
# endif
    if (cg_conn_write(cgfile, cgbase, 2, "Cylindrical -> Cartesian",
            Vertex, Abutting, PointList, n/3, pts, "CartesianZone",
            Structured, d_type, Integer, n/3, d_pts, &cgconn) ||
        cg_goto(cgfile, cgbase, "Zone_t", 2, "ZoneGridConnectivity_t", 1,
            "GridConnectivity_t", cgconn, "end") ||
        cg_array_write("InterpolantsDonor", RealSingle, 2, dims, interp))
        error_exit("cylindrical->cartesian connectivity");
#endif

    /* periodic boundary for zone 2 */

    for (n = 0; n < 3; n++) {
        transform[n] = n + 1;
        range[n] = d_range[n] = 1;
        range[n+3] = d_range[n+3] = NUM_SIDE;
    }
    range[4] = 1;
    d_range[1] = d_range[4] = nj;
    if (cg_1to1_write(cgfile, cgbase, 2, "Periodic", "CylindricalZone",
            range, d_range, transform, &cgconn))
        error_exit("periodic 1to1");

    /* write inlet BC (zone 1) as point range */

    range[4] = NUM_SIDE;
    range[5] = 1;
    if (cg_boco_write(cgfile, cgbase, 1, "Inlet", BCInflow,
            PointRange, 2, range, &cgbc))
        error_exit("inlet boco");

    /* write outlet BC (zone 2) as point range */

    range[4] = nj;
    range[2] = range[5] = NUM_SIDE;
    if (cg_boco_write(cgfile, cgbase, 2, "Outlet", BCOutflow,
            PointRange, 2, range, &cgbc))
        error_exit("outlet boco");

    /* write zone 1 wall BC as point list */

    for (n = 0, k = 1; k <= NUM_SIDE; k++) {
        for (i = 1; i < NUM_SIDE; i++) {
            pts[n++] = i + 1;
            pts[n++] = 1;
            pts[n++] = k;
            pts[n++] = i;
            pts[n++] = NUM_SIDE;
            pts[n++] = k;
        }
        for (j = 1; j < NUM_SIDE; j++) {
            pts[n++] = 1;
            pts[n++] = j;
            pts[n++] = k;
            pts[n++] = NUM_SIDE;
            pts[n++] = j+1;
            pts[n++] = k;
        }
    }
    if (cg_boco_write(cgfile, cgbase, 1, "Walls", BCWall,
            PointList, n/3, pts, &cgbc))
        error_exit("zone 1 walls boco");

    /* write zone 2 wall BC as point list */

    for (n = 0, k = 1; k <= NUM_SIDE; k++) {
        for (j = 1; j <= nj; j++) {
            pts[n++] = NUM_SIDE;
            pts[n++] = j;
            pts[n++] = k;
        }
    }
    if (cg_boco_write(cgfile, cgbase, 2, "Walls", BCWall,
            PointList, n/3, pts, &cgbc))
        error_exit("zone 2 walls boco");
}