[cig-commits] r6892 - in cs/cigma/trunk: include/cigma src
luis at geodynamics.org
luis at geodynamics.org
Wed May 16 10:12:39 PDT 2007
Author: luis
Date: 2007-05-16 10:12:39 -0700 (Wed, 16 May 2007)
New Revision: 6892
Added:
cs/cigma/trunk/include/cigma/cube.h
cs/cigma/trunk/src/cube.c
Log:
Moved cube subroutines into the cigma library
Added: cs/cigma/trunk/include/cigma/cube.h
===================================================================
--- cs/cigma/trunk/include/cigma/cube.h 2007-05-16 17:11:39 UTC (rev 6891)
+++ cs/cigma/trunk/include/cigma/cube.h 2007-05-16 17:12:39 UTC (rev 6892)
@@ -0,0 +1,28 @@
+#ifndef __CIGMA_CUBE_H__
+#define __CIGMA_CUBE_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+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);
+
+int cube_partition(cube_t *cube, int L, int M, int N);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif
Added: cs/cigma/trunk/src/cube.c
===================================================================
--- cs/cigma/trunk/src/cube.c 2007-05-16 17:11:39 UTC (rev 6891)
+++ cs/cigma/trunk/src/cube.c 2007-05-16 17:12:39 UTC (rev 6892)
@@ -0,0 +1,292 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <cigma/cube.h>
+
+
+static int cube_hex_connectivity(cube_t *cube);
+static int 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);
+}
+
+
+
+// --------------------------------------------------------------------------
+int 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;
+
+ int ierr;
+
+ 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");
+ return -1;
+ }
+
+ dx = 1.0 / L;
+ dy = 1.0 / M;
+ dz = 1.0 / N;
+
+ 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;
+ }
+ }
+ }
+
+ ierr = cube_hex_connectivity(cube);
+ if (ierr < 0)
+ {
+ return ierr;
+ }
+
+ ierr = cube_tet_connectivity(cube);
+ if (ierr < 0)
+ {
+ return ierr;
+ }
+
+ return 0;
+}
+
+static int 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");
+ return -1;
+ }
+
+ // 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 0;
+}
+
+static int 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");
+ return -1;
+ }
+
+ // 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 0;
+}
More information about the cig-commits
mailing list