[cig-commits] r6834 - cs/cigma/trunk/tools
luis at geodynamics.org
luis at geodynamics.org
Wed May 9 17:42:52 PDT 2007
Author: luis
Date: 2007-05-09 17:42:52 -0700 (Wed, 09 May 2007)
New Revision: 6834
Added:
cs/cigma/trunk/tools/cube.c
Log:
Command for generating a discretization (tet or hex) of the unit cube.
Added: cs/cigma/trunk/tools/cube.c
===================================================================
--- cs/cigma/trunk/tools/cube.c 2007-05-10 00:41:28 UTC (rev 6833)
+++ cs/cigma/trunk/tools/cube.c 2007-05-10 00:42:52 UTC (rev 6834)
@@ -0,0 +1,433 @@
+/* cube.c - Discretization of unit cube
+ *
+ */
+#include <hdf5.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <limits.h>
+#include <cigma/dataset.h>
+
+
+// --------------------------------------------------------------------------
+
+typedef struct {
+ int L, M, N;
+ int nno;
+ int hexnel;
+ int tetnel;
+ double *coords;
+ int *hexconn;
+ int *tetconn;
+} cube_t;
+
+
+void cube_init(cube_t *cube);
+void cube_free(cube_t *cube);
+void cube_partition(cube_t *cube, int L, int M, int N);
+void cube_hex_connectivity(cube_t *cube);
+void cube_tet_connectivity(cube_t *cube);
+
+
+// --------------------------------------------------------------------------
+
+void cube_init(cube_t *cube)
+{
+ cube->L = cube->M = cube->N = 0;
+ cube->nno = 0;
+ cube->hexnel = 0;
+ cube->tetnel = 0;
+ cube->coords = NULL;
+ cube->hexconn = NULL;
+ cube->tetconn = NULL;
+}
+
+void cube_free(cube_t *cube)
+{
+ if (cube->coords != NULL) free(cube->coords);
+ if (cube->hexconn != NULL) free(cube->hexconn);
+ if (cube->tetconn != NULL) free(cube->tetconn);
+}
+
+
+void cube_partition(cube_t *cube, int L, int M, int N)
+{
+ int n;
+ int nx, ny, nz;
+ int nxstride, nystride, nzstride;
+
+ double x, y, z;
+ double dx, dy, dz;
+
+ cube->L = L;
+ cube->M = M;
+ cube->N = N;
+
+ cube->nno = (L+1)*(M+1)*(N+1);
+ cube->coords = (double *)malloc((cube->nno) * 3 * sizeof(double));
+ cube->hexconn = NULL;
+ cube->tetconn = NULL;
+
+ if (cube->coords == NULL)
+ {
+ fprintf(stderr, "Error: Could not allocate cube->coords\n");
+ }
+
+ dx = 1.0/(L+1);
+ dy = 1.0/(M+1);
+ dz = 1.0/(N+1);
+
+ nxstride = (M+1)*(N+1);
+ nystride = (N+1);
+ nzstride = 1;
+
+ for (nx = 0; nx <= L; nx++)
+ {
+ x = nx * dx;
+ for (ny = 0; ny <= M; ny++)
+ {
+ y = ny * dy;
+ for (nz = 0; nz <= N; nz++)
+ {
+ z = nz * dz;
+ n = nx * nxstride + ny * nystride + nz * nzstride;
+ cube->coords[3*n + 0] = x;
+ cube->coords[3*n + 1] = y;
+ cube->coords[3*n + 2] = z;
+ }
+ }
+ }
+
+ cube_hex_connectivity(cube);
+ cube_tet_connectivity(cube);
+
+ return;
+}
+
+void cube_hex_connectivity(cube_t *cube)
+{
+ int L,M,N;
+
+ int e;
+ int ex, ey, ez;
+ int exstride, eystride, ezstride;
+
+ int n;
+ int nx, ny, nz;
+ int nxstride, nystride, nzstride;
+
+ int hex[8];
+
+ L = cube->L;
+ M = cube->M;
+ N = cube->N;
+
+ cube->hexnel = L*M*N;
+ cube->hexconn = (int *)malloc((cube->hexnel) * 8 * sizeof(int));
+
+ if (cube->hexconn == NULL)
+ {
+ fprintf(stderr, "Error: Could not allocate cube->hexconn\n");
+ }
+
+ // strides for the connectivity array
+ exstride = M*N;
+ eystride = N;
+ ezstride = 1;
+
+ // strides for the coords array
+ nxstride = (M+1)*(N+1);
+ nystride = (N+1);
+ nzstride = 1;
+
+ for (ex = 0; ex < L; ex++)
+ {
+ for (ey = 0; ey < M; ey++)
+ {
+ for (ez = 0; ez < N; ez++)
+ {
+ //
+ // Get (nx,ny,nz) for vertex v0
+ //
+ nx = ex;
+ ny = ey;
+ nz = ez;
+
+ //
+ // Local view of hex element
+ //
+ // v7-------v6
+ // /| /|
+ // / | / |
+ // v3-------v2 |
+ // | | | |
+ // | v4----|--v5
+ // | / | /
+ // |/ | /
+ // v0-------v1
+ //
+ #define HEXGLOB(x,y,z) ((x)*nxstride + (y)*nystride + (z)*(nzstride))
+ hex[0] = HEXGLOB(nx, ny, nz);
+ hex[1] = HEXGLOB(nx+1, ny, nz);
+ hex[2] = HEXGLOB(nx+1, ny+1, nz);
+ hex[3] = HEXGLOB(nx, ny+1, nz);
+ hex[4] = HEXGLOB(nx, ny, nz+1);
+ hex[5] = HEXGLOB(nx+1, ny, nz+1);
+ hex[6] = HEXGLOB(nx+1, ny+1, nz+1);
+ hex[7] = HEXGLOB(nx, ny+1, nz+1);
+ #undef HEXGLOB
+
+ //
+ // Finally, save the connectivity for the new hex elt
+ //
+ e = ex*exstride + ey*eystride + ez*ezstride;
+ for (n = 0; n < 8; n++)
+ {
+ cube->hexconn[8*e + n] = hex[n];
+ }
+ }
+ }
+ }
+
+ return;
+}
+
+void cube_tet_connectivity(cube_t *cube)
+{
+ int L,M,N;
+
+ int e, h, n;
+ int ex, ey, ez, et;
+ int exstride, eystride, ezstride, etstride;
+
+ int hex[8];
+ int tet[5][4];
+ int localmap[5][4] = {{0,1,3,4},
+ {1,2,3,6},
+ {1,3,4,6},
+ {1,4,5,6},
+ {3,4,6,7}};
+
+ // get the cube dimensions (number of elements per side)
+ L = cube->L;
+ M = cube->M;
+ N = cube->N;
+
+ // each old hex element turns into 5 tets
+ cube->tetnel = L*M*N*5;
+ cube->tetconn = (int *)malloc((cube->tetnel) * 4 * sizeof(int));
+
+ if (cube->tetconn == NULL)
+ {
+ fprintf(stderr, "Error: Could not allocate cube->tetconn\n");
+ }
+
+ // strides for element
+ exstride = M*N*5;
+ eystride = N*5;
+ ezstride = 5;
+ etstride = 1;
+
+ for (ex = 0; ex < L; ex++)
+ {
+ for (ey = 0; ey < M; ey++)
+ {
+ for (ez = 0; ez < N; ez++)
+ {
+ //
+ // First, load <v0,...,v7> into hex array
+ //
+ h = (ex*exstride + ey*eystride + ez*ezstride)/5;
+ for (n = 0; n < 8; n++)
+ {
+ hex[n] = cube->hexconn[8*h + n];
+ }
+
+ //
+ // Tetrahedral partition of hex element
+ //
+ // v7-------v6
+ // /| /|
+ // / | / |
+ // v3-------v2 |
+ // | | | |
+ // | v4----|--v5
+ // | / | /
+ // |/ | /
+ // v0-------v1
+ //
+ // T0 = <v0, v1, v3, v4>
+ // T1 = <v1, v2, v3, v6>
+ // T2 = <v1, v3, v4, v6>
+ // T3 = <v1, v4, v5, v6>
+ // T4 = <v3, v4, v6, v7>
+ //
+ // To map from the local hex element to the local tet elements,
+ // use the following
+ //
+ // int localmap[5][4] = {{0,1,3,4},
+ // {1,2,3,6},
+ // {1,3,4,6},
+ // {1,4,5,6},
+ // {3,4,6,7}};
+ //
+ for (et = 0; et < 5; et++)
+ {
+ for (n = 0; n < 4; n++)
+ {
+ tet[et][n] = hex[localmap[et][n]];
+ }
+ }
+
+ //
+ // Finally, save the connectivity for the five new elts
+ //
+ for (et = 0; et < 5; et++)
+ {
+ e = ex*exstride + ey*eystride + ez*ezstride + et*etstride;
+ for (n = 0; n < 4; n++)
+ {
+ cube->tetconn[4*e + n] = tet[et][n];
+ }
+ }
+ }
+ }
+ }
+
+ return;
+}
+
+
+// --------------------------------------------------------------------------
+
+int main(int argc, char *argv[])
+{
+ cube_t *cube;
+ int L,M,N;
+
+ hid_t file_id;
+ herr_t status;
+ int ierr;
+
+ /*
+ * Required initializations
+ */
+ L = M = N = 9;
+ cube = (cube_t *)malloc(sizeof(cube_t));
+ cube_init(cube);
+
+
+ /*
+ * Parse the command-line arguments
+ *
+ * TODO: add ability to specify L, M, N separately, from command-line
+ *
+ */
+
+ if (argc < 2)
+ {
+ fprintf(stderr, "Usage: cigma cube output.h5 [size]\n");
+ return EXIT_FAILURE;
+ }
+
+ if (argc == 3)
+ {
+ char *endptr, *str;
+ long val;
+
+ str = argv[2];
+
+ errno = 0; // to distinguish success/failure after call
+
+ val = strtol(str, &endptr, 10);
+
+ // Check for various possible errors
+ if ((errno == ERANGE && (val == LONG_MAX || val == LONG_MIN))
+ || (errno != 0 && val == 0))
+ {
+ perror("strtol");
+ exit(EXIT_FAILURE);
+ }
+
+ if (endptr == str)
+ {
+ fprintf(stderr, "No digits were found in second argument\n");
+ exit(EXIT_FAILURE);
+ }
+
+ N = (int)val;
+ }
+
+
+ /*
+ * Open a new HDF5 file
+ */
+ file_id = H5Fcreate(argv[1], H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ fprintf(stderr, "Could not create file '%s'\n", argv[1]);
+ cube_free(cube);
+ return EXIT_FAILURE;
+ }
+
+ /*
+ * Partition the unit cube
+ */
+ L = M = N;
+ cube_partition(cube, L, M, N);
+
+ #ifndef NDEBUG
+ printf("file = \"%s\"\n", argv[1]);
+ printf("L, M, N = (%d, %d, %d)\n", L, M, N);
+ printf("cube->nno = %d\n", cube->nno);
+ printf("cube->hexnel = %d\n", cube->hexnel);
+ printf("cube->tetnel = %d\n", cube->tetnel);
+ #endif
+
+
+ /*
+ * Create the /coords dataset
+ */
+ ierr = dataset_write2(file_id, "/coords", H5T_NATIVE_DOUBLE,
+ (void *)cube->coords, cube->nno, 3);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Could not create dataset /coords in '%s'\n", argv[1]);
+ cube_free(cube);
+ return EXIT_FAILURE;
+ }
+
+ /*
+ * Save the hexahedral mesh connectivity
+ */
+ ierr = dataset_write2(file_id, "/hexconn", H5T_NATIVE_INT,
+ (void *)cube->hexconn, cube->hexnel, 8);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Could not create dataset /hexconn in '%s'\n", argv[1]);
+ cube_free(cube);
+ free(cube);
+ return EXIT_FAILURE;
+ }
+
+ /*
+ * Save the tetrahedral mesh connectivity
+ */
+ ierr = dataset_write2(file_id, "/tetconn", H5T_NATIVE_INT,
+ (void *)cube->tetconn, cube->tetnel, 4);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Could not create dataset /tetconn in '%s'\n", argv[1]);
+ cube_free(cube);
+ free(cube);
+ return EXIT_FAILURE;
+ }
+
+ /*
+ * Clean up
+ */
+ status = H5Fclose(file_id);
+ cube_free(cube);
+ free(cube);
+ return EXIT_SUCCESS;
+}
More information about the cig-commits
mailing list