[cig-commits] r7499 - cs/cigma/branches/cigma-0.9/src
luis at geodynamics.org
luis at geodynamics.org
Tue Jun 26 09:42:39 PDT 2007
Author: luis
Date: 2007-06-26 09:42:39 -0700 (Tue, 26 Jun 2007)
New Revision: 7499
Added:
cs/cigma/branches/cigma-0.9/src/array.c
cs/cigma/branches/cigma-0.9/src/attr.c
cs/cigma/branches/cigma-0.9/src/cube.c
cs/cigma/branches/cigma-0.9/src/dataset.c
cs/cigma/branches/cigma-0.9/src/det.c
cs/cigma/branches/cigma-0.9/src/fe.c
cs/cigma/branches/cigma-0.9/src/field.c
cs/cigma/branches/cigma-0.9/src/io.c
cs/cigma/branches/cigma-0.9/src/mesh.c
cs/cigma/branches/cigma-0.9/src/points.c
cs/cigma/branches/cigma-0.9/src/rule.c
cs/cigma/branches/cigma-0.9/src/split.c
Log:
darcs patch: Added most source files for cigma
Added: cs/cigma/branches/cigma-0.9/src/array.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/array.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/array.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,428 @@
+#include <stdlib.h>
+#include <hdf5.h>
+#include <cigma/array.h>
+#include <cigma/attr.h>
+
+
+
+int array_init(array_t *array, hid_t type_id, int rank, int *shape, void *data)
+{
+ int i;
+
+ // initialize dataspace
+ array->rank = rank;
+ array->shape = (hsize_t *)malloc(rank * sizeof(hsize_t));
+ if (shape != NULL)
+ {
+ for (i = 0; i < rank; i++)
+ array->shape[i] = shape[i];
+ }
+ else
+ {
+ for (i = 0; i < rank; i++)
+ array->shape[i] = 0;
+ }
+
+ // initialize datatype
+ array->type_id = type_id;
+
+ // initialize data
+ array->data = data;
+
+ return 0;
+}
+
+int array_free_dims(array_t *array)
+{
+ free(array->shape);
+ return 0;
+}
+
+int array_free_data(array_t *array)
+{
+ if (array->data != NULL)
+ {
+ free(array->data);
+ }
+ return 0;
+}
+
+int array_free(array_t *array)
+{
+ array_free_dims(array);
+ array_free_data(array);
+ return 0;
+}
+
+
+
+int array_npoints(array_t *array)
+{
+ int i;
+ int npoints = 1;
+ for (i = 0; i < array->rank; i++)
+ {
+ npoints *= (int) array->shape[i];
+ }
+ return npoints;
+}
+
+int array_dim(array_t *array, int i)
+{
+
+ if ((0 <= i) && (i < array->rank))
+ {
+ return (int)(array->shape[i]);
+ }
+ return 0;
+}
+
+void array_dims(array_t *array, int *rank, int **shape)
+{
+ int i;
+
+ *rank = array->rank;
+ *shape = (int *)malloc((array->rank) * sizeof(int));
+ for (i = 0; i < array->rank; i++)
+ {
+ (*shape)[i] = (int) array->shape[i];
+ }
+}
+
+void array_dims1(array_t *array, int *n)
+{
+ *n = array->shape[0];
+}
+
+void array_dims2(array_t *array, int *m, int *n)
+{
+ *m = array->shape[0];
+ *n = array->shape[1];
+}
+
+void array_dims3(array_t *array, int *l, int *m, int *n)
+{
+ *l = array->shape[0];
+ *m = array->shape[1];
+ *n = array->shape[2];
+}
+
+
+
+hid_t array_open(array_t *array, hid_t loc_id, const char *name)
+{
+ hid_t dataset_id;
+ hid_t dataspace_id;
+ herr_t status;
+
+ dataset_id = H5Dopen(loc_id, name);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ dataspace_id = H5Dget_space(dataset_id);
+ if (dataspace_id < 0)
+ {
+ H5Dclose(dataset_id);
+ return -2;
+ }
+
+ array->rank = H5Sget_simple_extent_ndims(dataspace_id);
+
+ status = H5Sget_simple_extent_dims(dataspace_id, array->shape, NULL);
+
+ status = H5Sclose(dataspace_id);
+
+ return dataset_id;
+}
+
+hid_t array_create(array_t *array, hid_t loc_id, const char *name, const char *title)
+{
+ hid_t dataspace_id;
+ hid_t dataset_id;
+ herr_t status;
+
+ dataspace_id = H5Screate_simple(array->rank, array->shape, NULL);
+ if (dataspace_id < 0)
+ {
+ return -1;
+ }
+
+ dataset_id = H5Dcreate(loc_id, name, array->type_id, dataspace_id,
+ H5P_DEFAULT);
+ if (dataset_id < 0)
+ {
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = set_str_attr(dataset_id, "TITLE", title);
+ status = set_str_attr(dataset_id, "CLASS", "ARRAY");
+ status = set_str_attr(dataset_id, "FLAVOR", "numpy");
+ status = set_str_attr(dataset_id, "VERSION", "2.3");
+
+ status = H5Sclose(dataspace_id);
+
+ return dataset_id;
+}
+
+
+
+int array_read(array_t *array, hid_t loc_id, const char *name)
+{
+ int npoints;
+ hid_t dataset_id;
+ herr_t status;
+
+ dataset_id = array_open(array, loc_id, name);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ npoints = array_npoints(array);
+
+ array->data = malloc(npoints * H5Tget_size(array->type_id));
+ if (array->data == NULL)
+ {
+ H5Dclose(dataset_id);
+ return -2;
+ }
+
+ status = H5Dread(dataset_id, array->type_id, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, array->data);
+ if (status < 0)
+ {
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ status = H5Dclose(dataset_id);
+
+ return 0;
+}
+
+int array_write(array_t *array, hid_t loc_id, const char *name)
+{
+ hid_t dataspace_id;
+ hid_t dataset_id;
+ herr_t status;
+
+ dataspace_id = H5Screate_simple(array->rank, array->shape, NULL);
+ if (dataspace_id < 0)
+ {
+ return -1;
+ }
+
+ dataset_id = H5Dcreate(loc_id, name, array->type_id, dataspace_id,
+ H5P_DEFAULT);
+ if (dataset_id < 0)
+ {
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+ status = H5Dwrite(dataset_id, array->type_id, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, array->data);
+ if (status < 0)
+ {
+ H5Sclose(dataspace_id);
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ status = H5Sclose(dataspace_id);
+ status = H5Dclose(dataset_id);
+
+ return 0;
+}
+
+
+
+static herr_t array_hyperslab_init(array_t *array,
+ hsize_t **offset, hsize_t **stride,
+ hsize_t **count, hsize_t **block,
+ int start, int end)
+{
+ int i;
+
+ *offset = NULL;
+ *stride = NULL;
+ *count = NULL;
+ *block = NULL;
+
+ if ((start < 0) || (start > end))
+ {
+ return -1;
+ }
+
+ if (start >= array->rank)
+ {
+ return -2;
+ }
+
+ if (array->rank == 0)
+ {
+ return -3;
+ }
+
+ *offset = (hsize_t *)malloc((array->rank) * sizeof(hsize_t));
+ *stride = (hsize_t *)malloc((array->rank) * sizeof(hsize_t));
+ *count = (hsize_t *)malloc((array->rank) * sizeof(hsize_t));
+ *block = (hsize_t *)malloc((array->rank) * sizeof(hsize_t));
+
+ // first select everything
+ for (i = 0; i < array->rank; i++)
+ {
+ (*offset)[i] = 0;
+ (*block )[i] = array->shape[i];
+ (*stride)[i] = 1;
+ (*count )[i] = 1;
+ }
+
+ // now, do the selection on the first dimension only
+ (*offset)[0] = start;
+ (*block )[0] = end - start;
+
+ return 0;
+}
+
+static void array_hyperslab_free(hsize_t *offset, hsize_t *stride,
+ hsize_t *count, hsize_t *block)
+{
+ if (offset != NULL) free(offset);
+ if (stride != NULL) free(stride);
+ if (count != NULL) free(count);
+ if (block != NULL) free(block);
+}
+
+int array_read_slice(array_t *array, hid_t dataset_id, int start, int end)
+{
+ hsize_t *offset, *stride, *count, *block;
+ hid_t memspace_id;
+ hid_t filespace_id;
+ herr_t status;
+
+ status = array_hyperslab_init(array,
+ &offset, &stride,
+ &count, &block,
+ start, end);
+ if (status < 0)
+ {
+ return -1;
+ }
+
+ memspace_id = H5Screate_simple(array->rank, block, NULL);
+ if (memspace_id < 0)
+ {
+ return -2;
+ }
+
+ filespace_id = H5Dget_space(dataset_id);
+ if (filespace_id < 0)
+ {
+ H5Sclose(memspace_id);
+ return -3;
+ }
+
+ if (array->data == NULL)
+ {
+ int i;
+ int npoints = 1;
+
+ for (i = 0; i < array->rank; i++)
+ {
+ npoints *= block[i];
+ }
+
+ array->data = malloc(npoints * H5Tget_size(array->type_id));
+ if (array->data == NULL)
+ {
+ H5Sclose(memspace_id);
+ H5Sclose(filespace_id);
+ return -4;
+ }
+ }
+
+ status = H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET,
+ offset, stride, count, block);
+ if (status < 0)
+ {
+ H5Sclose(memspace_id);
+ H5Sclose(filespace_id);
+ return -5;
+ }
+
+ status = H5Dread(dataset_id, array->type_id, memspace_id, filespace_id,
+ H5P_DEFAULT, array->data);
+ if (status < 0)
+ {
+ H5Sclose(memspace_id);
+ H5Sclose(filespace_id);
+ return -5;
+ }
+
+ status = H5Sclose(memspace_id);
+ status = H5Sclose(filespace_id);
+
+ array_hyperslab_free(offset, stride, count, block);
+
+ return 0;
+}
+
+int array_write_slice(array_t *array, hid_t dataset_id, int start, int end)
+{
+ hsize_t *offset, *stride, *count, *block;
+ hid_t memspace_id;
+ hid_t filespace_id;
+ herr_t status;
+
+ status = array_hyperslab_init(array,
+ &offset, &stride,
+ &count, &block,
+ start, end);
+ if (status < 0)
+ {
+ return -1;
+ }
+
+ memspace_id = H5Screate_simple(array->rank, block, NULL);
+ if (memspace_id < 0)
+ {
+ return -2;
+ }
+
+ filespace_id = H5Dget_space(dataset_id);
+ if (filespace_id < 0)
+ {
+ H5Sclose(memspace_id);
+ return -3;
+ }
+
+ status = H5Sselect_hyperslab(filespace_id, H5S_SELECT_SET,
+ offset, stride, count, block);
+ if (status < 0)
+ {
+ H5Sclose(memspace_id);
+ H5Sclose(filespace_id);
+ return -4;
+ }
+
+ status = H5Dwrite(dataset_id, array->type_id, memspace_id, filespace_id,
+ H5P_DEFAULT, array->data);
+ if (status < 0)
+ {
+ H5Sclose(memspace_id);
+ H5Sclose(filespace_id);
+ return -5;
+ }
+
+ status = H5Sclose(memspace_id);
+ status = H5Sclose(filespace_id);
+
+ array_hyperslab_free(offset, stride, count, block);
+
+ return 0;
+}
+
Added: cs/cigma/branches/cigma-0.9/src/attr.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/attr.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/attr.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,572 @@
+#include <stdlib.h>
+#include <string.h>
+#include <hdf5.h>
+#include <cigma/attr.h>
+
+/****************************************************************************
+ * Some of the following functions were based from the H5ATTR.c *
+ * source file in PyTables, which is a BSD-licensed python extension *
+ * for accessing HDF5 files. *
+ * *
+ * The copyright notice is hereby retained. *
+ * *
+ * NCSA HDF *
+ * Scientific Data Technologies *
+ * National Center for Supercomputing Applications *
+ * University of Illinois at Urbana-Champaign *
+ * 605 E. Springfield, Champaign IL 61820 *
+ * *
+ * For conditions of distribution and use, see the accompanying *
+ * hdf/COPYING file. *
+ * *
+ * Modified versions of H5LT for getting and setting attributes for open *
+ * groups and leaves. *
+ * F. Altet 2005/09/29 *
+ * *
+ ****************************************************************************/
+
+
+
+/* Function : find_attr_op
+ * Purpose : operator function used by find_attribute
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : June 21, 2001
+ */
+static herr_t find_attr_op(hid_t loc_id, const char *name, void *op_data)
+{
+ /*
+ * Define a default zero value for return. This will cause the
+ * iterator to continue if the palette attribute is not found yet.
+ */
+ int ret = 0;
+
+ char *attr_name = (char *)op_data;
+
+ /* Shut the compiler up */
+ loc_id = loc_id;
+
+ /*
+ * Define a positive value for return value if the attribute was
+ * found. This will cause the iterator to immediately return that
+ * positive value, indicating short-circuit sucesss
+ */
+
+ if (strcmp(name, attr_name) == 0)
+ {
+ ret = 1;
+ }
+
+ return ret;
+}
+
+
+
+/* Function : find_attr
+ * Purpose : Inquires if an attribute named attr_name exists attached
+ * to the object loc_id.
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : June 21, 2001
+ *
+ * Comments:
+ *
+ * The function uses H5Aiterate with the operator function find_attr
+ *
+ * Return:
+ *
+ * Success: The return value of the first operator that returns
+ * non-zero, or zero if all members were processed with no
+ * operator returning non-zero.
+ *
+ * Failure: Negative if something goes wrong within the library,
+ * or the negative value returned by one of the operators.
+ */
+herr_t find_attr(hid_t loc_id, const char *attr_name)
+{
+ unsigned int attr_num;
+ herr_t ret;
+
+ attr_num = 0;
+ ret = H5Aiterate(loc_id, &attr_num, find_attr_op, (void *)attr_name);
+
+ return ret;
+}
+
+
+
+/* Function : set_attr
+ * Purpose : Private function used by
+ * set_int_attr and set_float_attr
+ * Return : Success 0, Failure -1
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : July 25, 2001
+ */
+herr_t set_attr(hid_t obj_id, const char *attr_name, hid_t type_id, const void *data)
+{
+ hid_t space_id, attr_id;
+ herr_t status;
+
+ int has_attr;
+
+ /* Create the data space for the attribute. */
+ space_id = H5Screate(H5S_SCALAR);
+ if (space_id < 0)
+ {
+ goto out;
+ }
+
+ /* Verify if the attribute already exists. */
+ has_attr = find_attr(obj_id, attr_name);
+ if (has_attr == 1)
+ {
+ /* The attribute already exists. Delete it. */
+ status = H5Adelete(obj_id, attr_name);
+ if (status < 0)
+ {
+ goto out;
+ }
+ }
+
+ /* Create the attribute. */
+ attr_id = H5Acreate(obj_id, attr_name, type_id, space_id, H5P_DEFAULT);
+ if (attr_id < 0)
+ {
+ goto out;
+ }
+
+ /* Write the attribute data. */
+ status = H5Awrite(attr_id, type_id, data);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Close the attribute. */
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Close the data space. */
+ status = H5Sclose(space_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ return 0;
+
+out:
+ return -1;
+}
+
+
+
+/* Function : get_attr
+ * Purpose : Reads an attribute named attr_name with the memory type type_id
+ * Return : Success 0, Failure -1
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : September 19, 2002
+ */
+herr_t get_attr(hid_t obj_id, const char *attr_name, hid_t type_id, void *data)
+{
+ hid_t attr_id;
+ herr_t status;
+
+ attr_id = H5Aopen_name(obj_id, attr_name);
+ if (attr_id < 0)
+ {
+ return -1;
+ }
+
+ status = H5Aread(attr_id, type_id, data);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ return -1;
+ }
+
+ return 0;
+
+out:
+ H5Aclose(attr_id);
+ return -1;
+}
+
+
+
+/* Function : get_attr_dims
+ * Purpose : Gets the dimensionality of an attribute
+ * Return : Success 0, Failure -1
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : September 4, 2001
+ */
+herr_t get_attr_dims(hid_t obj_id,
+ const char *attr_name,
+ int *rank,
+ hsize_t *dims)
+{
+ hid_t attr_id;
+ hid_t space_id;
+ herr_t status;
+
+ /* Open the attribute */
+ attr_id = H5Aopen_name(obj_id, attr_name);
+ if (attr_id < 0)
+ {
+ return -1;
+ }
+
+ /* Get the dataspace handle */
+ space_id = H5Aget_space(attr_id);
+ if (space_id < 0)
+ {
+ goto out;
+ }
+
+ /* Get number of dimensions */
+ *rank = H5Sget_simple_extent_ndims(space_id);
+ if (*rank < 0)
+ {
+ goto out;
+ }
+
+ /* Get dimensions */
+ status = H5Sget_simple_extent_dims(space_id, dims, NULL);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Terminate access to the dataspace */
+ status = H5Sclose(space_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* End access to the attribute */
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ return 0;
+
+out:
+ H5Aclose(attr_id);
+ return -1;
+}
+
+
+
+/* Function : set_str_attr
+ *
+ * Purpose : Creates and writes a string attribute named attr_name
+ * and attaches it to the object specified by obj_id
+ *
+ * Return : Success 0, Failure -1
+ * Programmer: Pedro Vicente, pvn at ncsa.uiuc.edu
+ * Date : July 23, 2001
+ * Comments : If the attribute already exists, it is overwritten.
+ */
+herr_t set_str_attr(hid_t obj_id,
+ const char *attr_name,
+ const char *attr_data)
+{
+ hid_t attr_type;
+ hid_t attr_size;
+ hid_t attr_space_id;
+ hid_t attr_id;
+ int has_attr;
+ herr_t status;
+
+ /* Create the attribute */
+ attr_type = H5Tcopy(H5T_C_S1);
+ if (attr_type < 0)
+ {
+ goto out;
+ }
+
+ attr_size = strlen(attr_data) + 1; /* extra null term */
+
+ status = H5Tset_size(attr_type, (size_t)attr_size);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ status = H5Tset_strpad(attr_type, H5T_STR_NULLTERM);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ attr_space_id = H5Screate(H5S_SCALAR);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Verify if the attribute already exists */
+ has_attr = find_attr(obj_id, attr_name);
+
+ /* The attribute already exists, delete it */
+ if (has_attr == 1)
+ {
+ status = H5Adelete(obj_id, attr_name);
+ if (status < 0)
+ {
+ goto out;
+ }
+ }
+
+ /* Create the attribute. */
+ attr_id = H5Acreate(obj_id, attr_name, attr_type, attr_space_id,
+ H5P_DEFAULT);
+ if (attr_id < 0)
+ {
+ goto out;
+ }
+
+ /* Write the attribute data. */
+ status = H5Awrite(attr_id, attr_type, attr_data);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+ status = H5Sclose(attr_space_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ status = H5Tclose(attr_type);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ return 0;
+
+out:
+ return -1;
+}
+
+
+
+/* Function : get_str_attr
+ * Purpose : Reads a string attribute named attr_name
+ * Return : Success 0, Failure -1
+ * Programmer: Francesc Altet, faltet at carabos.com
+ * Date : February 23, 2005
+ */
+herr_t get_str_attr(hid_t obj_id, const char *attr_name, char **data)
+{
+ hid_t attr_id;
+ hid_t attr_type;
+ size_t type_size;
+ herr_t status;
+
+ *data = NULL;
+
+ attr_id = H5Aopen_name(obj_id, attr_name);
+ if (attr_id < 0)
+ {
+ return -1;
+ }
+
+ attr_type = H5Aget_type(attr_id);
+ if (attr_type < 0)
+ {
+ goto out;
+ }
+
+ /* Get the size */
+ type_size = H5Tget_size(attr_type);
+ if (type_size < 0)
+ {
+ goto out;
+ }
+
+ /* Malloc enough space for the string, plus 1 for the trailing '\0' */
+ *data = (char *)malloc((type_size + 1) * sizeof(char));
+
+ status = H5Aread(attr_id, attr_type, *data);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Set the last character to \0 in case we are dealing with
+ * space padded strings
+ */
+ (*data)[type_size] = '\0';
+
+ status = H5Tclose(attr_type);
+ if (status < 0);
+ {
+ goto out;
+ }
+
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ return -1;
+ }
+
+ return 0;
+
+out:
+ H5Tclose(attr_type);
+ H5Aclose(attr_id);
+ if (*data)
+ free(*data);
+ return -1;
+}
+
+
+
+/* Function: set_array_attr
+ * Purpose : write an array attribute
+ * Return : Success 0, Failure -1
+ * Date : July 25, 2001
+ */
+herr_t set_array_attr(hid_t obj_id,
+ const char *attr_name,
+ size_t rank,
+ hsize_t *dims,
+ hid_t type_id,
+ const void *data)
+{
+ hid_t space_id, attr_id;
+ herr_t status;
+
+ int has_attr;
+
+ /* Create the data space for the attribute. */
+ space_id = H5Screate_simple(rank, dims, NULL);
+ if (space_id < 0)
+ {
+ goto out;
+ }
+
+ /* Verify if the attribute already exists. */
+ has_attr = find_attr(obj_id, attr_name);
+ if (has_attr == 1)
+ {
+ /* The attribute already exists. Delete it. */
+ status = H5Adelete(obj_id, attr_name);
+ if (status < 0)
+ {
+ goto out;
+ }
+ }
+
+ /* Create the attribute. */
+ attr_id = H5Acreate(obj_id, attr_name, type_id, space_id, H5P_DEFAULT);
+ if (attr_id < 0)
+ {
+ goto out;
+ }
+
+ /* Write the attribute data. */
+ status = H5Awrite(attr_id, type_id, data);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Close the attribute. */
+ status = H5Aclose(attr_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ /* Close the dataspace. */
+ status = H5Sclose(space_id);
+ if (status < 0)
+ {
+ goto out;
+ }
+
+ return 0;
+
+out:
+ return -1;
+}
+
+
+
+herr_t set_float_attr(hid_t obj_id, const char *attr_name, float x)
+{
+ return set_attr(obj_id, attr_name, H5T_NATIVE_FLOAT, &x);
+}
+
+herr_t set_double_attr(hid_t obj_id, const char *attr_name, double x)
+{
+ return set_attr(obj_id, attr_name, H5T_NATIVE_DOUBLE, &x);
+}
+
+herr_t set_int_attr(hid_t obj_id, const char *attr_name, int n)
+{
+ return set_attr(obj_id, attr_name, H5T_NATIVE_INT, &n);
+}
+
+herr_t set_long_attr(hid_t obj_id, const char *attr_name, long n)
+{
+ return set_attr(obj_id, attr_name, H5T_NATIVE_LONG, &n);
+}
+
+herr_t set_llong_attr(hid_t obj_id, const char *attr_name, long long n)
+{
+ return set_attr(obj_id, attr_name, H5T_NATIVE_LLONG, &n);
+}
+
+
+herr_t set_array1_attr(hid_t obj_id,
+ const char *attr_name,
+ hsize_t dim,
+ hid_t type_id,
+ const void *data)
+{
+ return set_array_attr(obj_id, attr_name, 1, &dim, type_id, data);
+}
+
+herr_t set_int_array_attr(hid_t obj_id,
+ const char *attr_name,
+ hsize_t dim,
+ const int *data)
+{
+ return set_array_attr(obj_id, attr_name, 1, &dim, H5T_NATIVE_INT, data);
+}
+
+herr_t set_float_array_attr(hid_t obj_id,
+ const char *attr_name,
+ hsize_t dim,
+ const float *data)
+{
+ return set_array_attr(obj_id, attr_name, 1, &dim, H5T_NATIVE_FLOAT, data);
+}
+
+herr_t set_double_array_attr(hid_t obj_id,
+ const char *attr_name,
+ hsize_t dim,
+ const double *data)
+{
+ return set_array_attr(obj_id, attr_name, 1, &dim, H5T_NATIVE_DOUBLE, data);
+}
+
Added: cs/cigma/branches/cigma-0.9/src/cube.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/cube.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/cube.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,292 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <cigma/cube.h>
+
+
+static int hexahedral_connectivity(cube_t *cube);
+static int tetrahedral_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 = hexahedral_connectivity(cube);
+ if (ierr < 0)
+ {
+ return ierr;
+ }
+
+ ierr = tetrahedral_connectivity(cube);
+ if (ierr < 0)
+ {
+ return ierr;
+ }
+
+ return 0;
+}
+
+static int hexahedral_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 tetrahedral_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;
+}
Added: cs/cigma/branches/cigma-0.9/src/dataset.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/dataset.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/dataset.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,254 @@
+#include <stdlib.h>
+#include <hdf5.h>
+#include <cigma/attr.h>
+#include <cigma/dataset.h>
+
+
+
+hid_t dataset_create(hid_t loc_id, const char *name, const char *title,
+ hid_t type_id, int rank, int *shape)
+{
+ hid_t dataset_id;
+ hid_t dataspace_id;
+ herr_t status;
+
+ if (rank > 0)
+ {
+ hsize_t *dims;
+ int i;
+
+ dims = (hsize_t *)malloc(rank * sizeof(hsize_t));
+ for (i = 0; i < rank; i++)
+ {
+ dims[i] = shape[i];
+ }
+
+ dataspace_id = H5Screate_simple(rank, dims, NULL);
+ if (dataspace_id < 0)
+ {
+ free(dims);
+ return -1;
+ }
+ free(dims);
+ }
+ else
+ {
+ return -3;
+ }
+
+
+ dataset_id = H5Dcreate(loc_id, name, type_id, dataspace_id, H5P_DEFAULT);
+ if (dataset_id < 0)
+ {
+ H5Sclose(dataspace_id);
+ return -2;
+ }
+
+
+ status = set_str_attr(dataset_id, "TITLE", title);
+ status = set_str_attr(dataset_id, "CLASS", "ARRAY");
+ status = set_str_attr(dataset_id, "FLAVOR", "numpy");
+ status = set_str_attr(dataset_id, "VERSION", "2.3");
+
+ status = H5Sclose(dataspace_id);
+
+ return dataset_id;
+}
+
+
+
+hid_t dataset_open(hid_t loc_id, const char *name,
+ hid_t *type_id, int *rank, int *shape,
+ int *npoints)
+{
+ hid_t dataset_id;
+ hid_t dataspace_id;
+ herr_t status;
+
+ dataset_id = H5Dopen(loc_id, name);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ if (type_id != NULL)
+ {
+ *type_id = H5Dget_type(dataset_id);
+ }
+
+ if ((rank != NULL) || (shape != NULL) || (npoints != NULL))
+ {
+ dataspace_id = H5Dget_space(dataset_id);
+
+ if (rank != NULL)
+ {
+ *rank = H5Sget_simple_extent_ndims(dataspace_id);
+ }
+
+ if (shape != NULL)
+ {
+ int i;
+ hsize_t *dims;
+
+ dims = (hsize_t *)malloc((*rank) * sizeof(hsize_t));
+ status = H5Sget_simple_extent_dims(dataspace_id, dims, NULL);
+ for (i = 0; i < *rank; i++)
+ {
+ shape[i] = (int)dims[i];
+ }
+ free(dims);
+ }
+
+ if (npoints != NULL)
+ {
+ *npoints = H5Sget_simple_extent_npoints(dataspace_id);
+ }
+
+ status = H5Sclose(dataspace_id);
+ }
+
+ return dataset_id;
+}
+
+
+
+int dataset_read(hid_t loc_id, const char *name,
+ hid_t type_id, int rank, int *shape,
+ void **data)
+{
+ int dset_rank;
+ int dset_npts;
+ hid_t dataset_id;
+ herr_t status;
+
+ dataset_id = dataset_open(loc_id, name, NULL, &dset_rank, shape, &dset_npts);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+ if (rank != dset_rank)
+ {
+ H5Dclose(dataset_id);
+ return -2;
+ }
+
+ *data = malloc(dset_npts * H5Tget_size(type_id));
+ if (*data == NULL)
+ {
+ H5Dclose(dataset_id);
+ return -3;
+ }
+
+ status = H5Dread(dataset_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, *data);
+ if (status < 0)
+ {
+ H5Dclose(dataset_id);
+ return -4;
+ }
+
+ status = H5Dclose(dataset_id);
+
+ return 0;
+}
+
+
+
+int dataset_write(hid_t loc_id, const char *name, const char *title,
+ hid_t type_id, int rank, int *shape,
+ void *data)
+{
+ hid_t dataset_id;
+ herr_t status;
+
+ dataset_id = dataset_create(loc_id, name, title, type_id, rank, shape);
+ if (dataset_id < 0)
+ {
+ return -1;
+ }
+
+ status = H5Dwrite(dataset_id, type_id, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, data);
+ if (status < 0)
+ {
+ H5Dclose(dataset_id);
+ return -2;
+ }
+
+ H5Dclose(dataset_id);
+ return 0;
+}
+
+
+
+int dataset_read1(hid_t loc_id,
+ const char *name,
+ hid_t type_id,
+ void **data, int *n)
+{
+ return dataset_read(loc_id, name, type_id, 1, n, data);
+}
+
+int dataset_read2(hid_t loc_id,
+ const char *name,
+ hid_t type_id,
+ void **data, int *m, int *n)
+{
+ int shape[2];
+ int ret;
+
+ ret = dataset_read(loc_id, name, type_id, 2, shape, data);
+
+ *m = shape[0];
+ *n = shape[1];
+
+ return ret;
+}
+
+int dataset_read3(hid_t loc_id,
+ const char *name,
+ hid_t type_id,
+ void **data, int *l, int *m, int *n)
+{
+ int shape[3];
+ int ret;
+
+ ret = dataset_read(loc_id, name, type_id, 3, shape, data);
+
+ *l = shape[0];
+ *m = shape[1];
+ *n = shape[2];
+
+ return ret;
+}
+
+
+
+int dataset_write1(hid_t loc_id,
+ const char *name,
+ const char *title,
+ hid_t type_id,
+ void *data, int n)
+{
+ return dataset_write(loc_id, name, title, type_id, 1, &n, data);
+}
+
+int dataset_write2(hid_t loc_id,
+ const char *name,
+ const char *title,
+ hid_t type_id,
+ void *data, int m, int n)
+{
+ int dims[2] = {m,n};
+ return dataset_write(loc_id, name, title, type_id, 2, dims, data);
+}
+
+int dataset_write3(hid_t loc_id,
+ const char *name,
+ const char *title,
+ hid_t type_id,
+ void *data, int l, int m, int n)
+{
+ int dims[3] = {l,m,n};
+ return dataset_write(loc_id, name, title, type_id, 3, dims, data);
+}
+
Added: cs/cigma/branches/cigma-0.9/src/det.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/det.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/det.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,106 @@
+#include <cigma/det.h>
+
+double det3x3(double *m)
+{
+ //
+ // |m[0] m[1] m[2]|
+ // |m[3] m[4] m[5]|
+ // |m[6] m[7] m[8]|
+ //
+ // = m[0]*|m[4] m[5]|
+ // |m[7] m[8]|
+ //
+ // - m[1]*|m[3] m[5]|
+ // |m[6] m[8]|
+ //
+ // + m[2]*|m[3] m[4]|
+ // |m[6] m[7]|
+ //
+ //
+ return m[0]*(m[4]*m[8] - m[5]*m[7])
+ - m[1]*(m[3]*m[8] - m[5]*m[6])
+ + m[2]*(m[3]*m[7] - m[4]*m[6]);
+}
+
+double det4x4(double *m)
+{
+ // |m[ 0] m[ 1] m[ 2] m[ 3]|
+ // |m[ 4] m[ 5] m[ 6] m[ 7]|
+ // |m[ 8] m[ 9] m[10] m[11]|
+ // |m[12] m[13] m[14] m[15]|
+ //
+ //
+ // = m[0] * |m[ 5] m[ 6] m[ 7]|
+ // |m[ 9] m[10] m[11]|
+ // |m[13] m[14] m[15]|
+ //
+ // - m[1] * |m[ 4] m[ 6] m[ 7]|
+ // |m[ 8] m[10] m[11]|
+ // |m[12] m[14] m[15]|
+ //
+ // + m[2] * |m[ 4] m[ 5] m[ 7]|
+ // |m[ 8] m[ 9] m[11]|
+ // |m[12] m[13] m[15]|
+ //
+ // - m[3] * |m[ 4] m[ 5] m[ 6]|
+ // |m[ 8] m[ 9] m[10]|
+ // |m[12] m[13] m[14]|
+ //
+ //
+ // = m[0] * ( + m[5] * |m[10] m[11]|
+ // |m[14] m[15]|
+ //
+ // - m[6] * |m[ 9] m[11]|
+ // |m[13] m[15]|
+ //
+ // + m[7] * |m[ 9] m[10]|
+ // |m[13] m[14]| )
+ //
+ // - m[1] * ( + m[4] * |m[10] m[11]|
+ // |m[14] m[15]|
+ //
+ // - m[6] * |m[ 8] m[11]|
+ // |m[12] m[15]|
+ //
+ // + m[7] * |m[ 8] m[10]|
+ // |m[12] m[14]| )
+ //
+ // + m[2] * ( + m[4] * |m[ 9] m[11]|
+ // |m[13] m[15]|
+ //
+ // - m[5] * |m[ 8] m[11]|
+ // |m[12] m[15]|
+ //
+ // + m[7] * |m[ 8] m[ 9]|
+ // |m[12] m[13]| )
+ //
+ // - m[3] * ( + m[4] * |m[ 9] m[10]|
+ // |m[13] m[14]|
+ //
+ // - m[5] * |m[ 8] m[10]|
+ // |m[12] m[14]|
+ //
+ // + m[6] * |m[ 8] m[ 9]|
+ // |m[12] m[13]| )
+
+ double c[4];
+
+ c[0] = m[5]*(m[10]*m[15]-m[11]*m[14])
+ - m[6]*(m[ 9]*m[15]-m[11]*m[13])
+ + m[7]*(m[ 9]*m[14]-m[10]*m[13]);
+
+ c[1] = m[4]*(m[10]*m[15]-m[11]*m[14])
+ - m[6]*(m[ 8]*m[15]-m[11]*m[12])
+ + m[7]*(m[ 8]*m[14]-m[10]*m[12]);
+
+ c[2] = m[4]*(m[ 9]*m[15]-m[11]*m[13])
+ - m[5]*(m[ 8]*m[15]-m[11]*m[12])
+ + m[7]*(m[ 8]*m[13]-m[ 9]*m[12]);
+
+ c[3] = m[4]*(m[ 9]*m[14]-m[10]*m[13])
+ - m[5]*(m[ 8]*m[14]-m[10]*m[12])
+ + m[6]*(m[ 8]*m[13]-m[ 9]*m[12]);
+
+ return m[0]*c[0] - m[1]*c[1] + m[2]*c[2] - m[3]*c[3];
+}
+
Added: cs/cigma/branches/cigma-0.9/src/fe.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/fe.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/fe.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,86 @@
+#include <cigma/fe.h>
+
+void fe_jacobian(fe_t *fe, double *x, double xi[3], double J[3*3])
+{
+ int i,j,k;
+ int nsd = fe->nsd;
+ //int dims = fe->ndim;
+ int ndof = fe->ndof;
+ shape_fn *dN = fe->dN;
+ for (i = 0; i < nsd; i++)
+ {
+ for (j = 0; j < nsd; j++)
+ {
+ //
+ // J[i,j] = dx[i]/dx_j
+ // = \sum_k x[i,k] (dN[j]/dx_i)
+ //
+ J[nsd*i + j] = 0.0;
+ for (k = 0; k < ndof; k++)
+ J[nsd*i + j] += dN[nsd*k + j](xi) * x[nsd*i + k];
+ }
+ }
+}
+
+/*
+ * Input : xs - an [npts x nsd] matrix
+ * Output : tab - an [npts x ndof] matrix
+ */
+void fe_tabulate(fe_t *fe, double *xs, int npts, double *tab)
+{
+ int i,j;
+ double *x;
+ int nsd = fe->nsd;
+ int ndof = fe->ndof;
+ shape_fn *N = fe->N;
+ for (i = 0; i < npts; i++)
+ {
+ x = &(xs[nsd*i]);
+ for (j = 0; j < ndof; j++)
+ tab[ndof*i + j] = N[j](x);
+ }
+}
+
+void fe_eval(fe_t *fe, double *dof, double *xi, double *val)
+{
+ int i,j;
+ int dims = fe->ndim;
+ int ndof = fe->ndof;
+ shape_fn *N = fe->N;
+ for (i = 0; i < dims; i++)
+ {
+ //
+ // val_i(xi) = \sum_j dof[i,j] * TN_j(xi)
+ //
+ val[i] = 0.0;
+ for (j = 0; j < ndof; j++)
+ val[i] += dof[ndof*i + j] * N[j](xi);
+ }
+}
+
+/*
+ * Inputs:
+ * tab - an [npts x ndof] matrix
+ * dof - an [ndof x dims] matrix
+ *
+ * Output:
+ * vals - an [npts x dims] matrix
+ */
+void fe_batch_eval(fe_t *fe, double *dof, double *tab, int npts, double *vals)
+{
+ int i,j,k;
+ int dims = fe->ndim;
+ int ndof = fe->ndof;
+ for (i = 0; i < npts; i++)
+ {
+ for (j = 0; j < dims; j++)
+ {
+ //
+ // vals[i,j] = \sum_k tab[i,k] * dof[k,j]
+ //
+ vals[dims*i + j] = 0.0;
+ for (k = 0; k < ndof; k++)
+ vals[dims*i + j] += tab[ndof*i + k] * dof[dims*k + j]; // tab[i,k] * dof[k,j]
+ }
+ }
+}
Added: cs/cigma/branches/cigma-0.9/src/field.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/field.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/field.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,124 @@
+#include <stdlib.h>
+#include <hdf5.h>
+#include <cigma/field.h>
+#include <cigma/io.h>
+#include <cigma/fe.h>
+#include <cigma/mesh.h>
+#include <cigma/dataset.h>
+#include "elements.h"
+
+
+
+void field_dofs(field_t *field, int e, double *dof)
+{
+ int n,i,j;
+ int ndof = 4; // XXX: use fe->ndof instead
+ int ndim = field->ndim;
+ int *elt = &(field->mesh->connect[ndof*e]);
+ for (i = 0; i < ndof; i++)
+ {
+ n = elt[i];
+ for (j = 0; j < ndim; j++)
+ {
+ dof[ndim*i + j] = field->dofs[ndim*n + j];
+ }
+ }
+}
+
+void field_eval(field_t *B, double *x, int nq, double *g, int dims)
+{
+ int q,e;
+ int nel = B->mesh->nel;
+ int contained;
+
+ double verts[4*3];
+ double dof[4*3];
+
+ for (q = 0; q < nq; q++)
+ {
+ // XXX: replace this search in other routines
+ for (e = 0; e < nel; e++)
+ {
+ mesh_coords(B->mesh, e, verts);
+ contained = tet4_test(verts, &x[3*q]);
+ if (contained == 0)
+ {
+ // found the element...evaluate and break element loop
+ field_dofs(B, e, verts);
+ tet4_eval3(dof, &x[3*q], &g[dims*q]);
+ break;
+ }
+ }
+ }
+}
+
+
+
+int field_init(field_t *field, mesh_t *mesh)
+{
+ field->mesh = mesh;
+ field->nno = mesh->nno;
+ field->ndim = 0;
+ return 0;
+}
+
+int field_init2(field_t *field, mesh_t *mesh, int nno, int ndim, double *dofs)
+{
+ field->mesh = mesh;
+ field->nno = nno;
+ field->ndim = ndim;
+ field->dofs = dofs;
+ return 0;
+}
+
+int field_free(field_t *field)
+{
+ if (field->dofs != NULL)
+ free(field->dofs);
+
+ field->nno = 0;
+ field->ndim = 0;
+ field->mesh = NULL;
+ field->dofs = NULL;
+
+ return 0;
+}
+
+
+
+int field_open(field_t *field, const char *filename, const char *path)
+{
+ hid_t file_id;
+ herr_t status;
+ int ierr;
+
+ /* open the file in read-write mode */
+ file_id = file_open(filename, "rw");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ ierr = field_read(field, file_id, path);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ status = H5Fclose(file_id);
+ return 0;
+}
+
+int field_read(field_t *field, hid_t loc_id, const char *name)
+{
+ return dataset_read2(loc_id, name, H5T_NATIVE_DOUBLE,
+ (void **)&(field->dofs), &(field->nno), &(field->ndim));
+}
+
+int field_write(field_t *field, hid_t loc_id, const char *name, const char *title)
+{
+ return dataset_write2(loc_id, name, title, H5T_NATIVE_DOUBLE,
+ field->dofs, field->nno, field->ndim);
+}
+
Added: cs/cigma/branches/cigma-0.9/src/io.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/io.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/io.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,275 @@
+#include <string.h>
+#include <hdf5.h>
+#include <cigma/io.h>
+#include <cigma/attr.h>
+#include <cigma/split.h>
+
+
+
+hid_t file_create(const char *filename, const char *mode)
+{
+ hid_t file_id;
+ hid_t root;
+ herr_t status;
+
+ if (strcmp(mode, "w") == 0)
+ {
+ /* Create file by truncating (i.e. overwriting previous file) */
+ file_id = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ }
+ else if (strcmp(mode, "x") == 0)
+ {
+ /* Create file exclusively (i.e. fail if it already exists) */
+ file_id = H5Fcreate(filename, H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT);
+ }
+ else
+ {
+ return -1;
+ }
+
+ if (file_id < 0)
+ {
+ return -2;
+ }
+
+ root = H5Gopen(file_id, "/");
+ status = set_str_attr(root, "TITLE", "CIGMA File");
+ status = set_str_attr(root, "CLASS", "GROUP");
+ status = set_str_attr(root, "VERSION", "1.0");
+ status = set_str_attr(root, "PYTABLES_FORMAT_VERSION", "1.5");
+ status = H5Gclose(root);
+
+ return file_id;
+}
+
+
+
+hid_t file_open(const char *filename, const char *mode)
+{
+ /*
+ * Open file for reading. Fail if file doesn't exist.
+ */
+ if (strcmp(mode, "r") == 0)
+ {
+ hid_t file_id;
+
+ /* Open file in read-only mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+
+ /* Check for failure */
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ return file_id;
+ }
+
+ /*
+ * Open file for writing. If file exists, it is truncated.
+ */
+ if (strcmp(mode, "w") == 0)
+ {
+ hid_t file_id;
+
+ file_id = file_create(filename, "w");
+
+ if (file_id < 0)
+ {
+ return -2;
+ }
+
+ return file_id;
+ }
+
+ /*
+ * Open file for reading and writing. Fail if file doesn't exist.
+ */
+ if (strcmp(mode, "rw") == 0)
+ {
+ hid_t file_id;
+
+ /* Open file in read-write mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
+
+ /* Check for failure */
+ if (file_id < 0)
+ {
+ return -3;
+ }
+
+ return file_id;
+ }
+
+ /*
+ * Open file for reading and writing. Create the file if necessary.
+ */
+ if (strcmp(mode, "rw+") == 0)
+ {
+ hid_t file_id;
+
+ /* See http://hdf.ncsa.uiuc.edu/HDF5/doc/Errors.html */
+ herr_t (*old_func)(void *);
+ void *old_client_data;
+
+ /* Save old error handler */
+ H5Eget_auto(&old_func, &old_client_data);
+
+ /* Turn off error handling */
+ H5Eset_auto(NULL, NULL);
+
+ /* Open file in read-write mode -- errors suppressed */
+ file_id = H5Fopen(filename, H5F_ACC_RDWR, H5P_DEFAULT);
+
+ /* Restore error handler */
+ H5Eset_auto(old_func, old_client_data);
+
+ /* If opening the file failed, try to create it */
+ if (file_id < 0)
+ {
+ file_id = file_create(filename, "w");
+
+ if (file_id < 0)
+ {
+ return -4;
+ }
+ }
+
+ return file_id;
+ }
+
+ /*
+ * Exclusively open file for writing. Fail if file already exists.
+ */
+ if (strcmp(mode, "x") == 0)
+ {
+ hid_t file_id;
+
+ file_id = file_create(filename, "x");
+
+ if (file_id < 0)
+ {
+ return -5;
+ }
+ }
+
+
+ /*
+ * Invalid mode
+ */
+ return -6;
+}
+
+
+
+hid_t group_create(hid_t loc_id, const char *name)
+{
+ hid_t group_id;
+ herr_t status;
+
+ group_id = H5Gcreate(loc_id, name, 0);
+ if (group_id < 0)
+ {
+ return -1;
+ }
+
+ status = set_str_attr(group_id, "TITLE", "");
+ status = set_str_attr(group_id, "CLASS", "GROUP");
+ status = set_str_attr(group_id, "VERSION", "1.0");
+ status = set_str_attr(group_id, "PYTABLES_FORMAT_VERSION", "1.5");
+
+ return group_id;
+}
+
+
+
+hid_t group_open(hid_t loc_id, const char *name)
+{
+ return H5Gopen(loc_id, name);
+}
+
+
+
+hid_t touch_group(hid_t loc_id, const char *name)
+{
+ hid_t group_id;
+
+ group_id = H5Gopen(loc_id, name);
+
+ if (group_id < 0)
+ {
+ group_id = group_create(loc_id, name);
+ if (group_id < 0)
+ {
+ return -1;
+ }
+ }
+
+ return group_id;
+}
+
+
+
+hid_t touch_path(hid_t loc_id, const char *path)
+{
+ hid_t group_id;
+
+ /* Temporarily disabling error messages:
+ * http://hdf.ncsa.uiuc.edu/HDF5/doc/Errors.html
+ */
+ herr_t (*old_func)(void *);
+ void *old_client_data;
+
+ /* Save old error handler */
+ H5Eget_auto(&old_func, &old_client_data);
+
+ /* Turn off error handling */
+ H5Eset_auto(NULL, NULL);
+
+ /* Attempt to open the group */
+ group_id = H5Gopen(loc_id, path);
+
+ if (group_id < 0)
+ {
+ herr_t status;
+ hid_t parent_id;
+ hid_t child_id;
+ char **names;
+ int i,n;
+
+ split(path, strlen(path), &names, &n, '/');
+
+ // first parent
+ parent_id = touch_group(loc_id, names[0]);
+ if (parent_id < 0)
+ {
+ H5Eset_auto(old_func, old_client_data);
+ split_free(names, n);
+ return -1;
+ }
+
+ for (i = 1; i < n; i++)
+ {
+ // get child id
+ child_id = touch_group(parent_id, names[i]);
+ if (child_id < 0)
+ {
+ H5Eset_auto(old_func, old_client_data);
+ split_free(names, n);
+ return -2;
+ }
+ // move to next parent
+ status = H5Gclose(parent_id);
+ parent_id = child_id;
+ }
+
+ // return last group
+ group_id = parent_id;
+ }
+
+ /* Restore previous error handler */
+ H5Eset_auto(old_func, old_client_data);
+
+ return group_id;
+}
+
Added: cs/cigma/branches/cigma-0.9/src/mesh.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/mesh.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/mesh.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,162 @@
+#include <stdlib.h>
+#include <assert.h>
+#include <hdf5.h>
+#include <cigma/io.h>
+#include <cigma/mesh.h>
+#include <cigma/dataset.h>
+
+
+
+/*
+ * Input : mesh object, element number e
+ * Output: dof coordinates -- [ndof x nsd] matrix
+ */
+void mesh_coords(mesh_t *mesh, int e, double *dof_coords)
+{
+
+ double *coord = mesh->coords;
+ int *conn = mesh->connect;
+ int ndof = mesh->ndof;
+
+ int i,j;
+ for (i = 0; i < ndof; i++)
+ {
+ j = conn[ndof*e + i];
+ dof_coords[3*i + 0] = coord[3*j + 0];
+ dof_coords[3*i + 1] = coord[3*j + 1];
+ dof_coords[3*i + 2] = coord[3*j + 2];
+ }
+}
+
+
+
+int mesh_init_coords(mesh_t *mesh, int nno, int nsd)
+{
+ mesh->nno = nno;
+ mesh->nsd = nsd;
+ mesh->coords = (double *)malloc(nno * nsd * sizeof(double));
+ return 0;
+}
+
+int mesh_init_connect(mesh_t *mesh, int nel, int ndof)
+{
+ mesh->nel = nel;
+ mesh->ndof = ndof;
+ mesh->connect = (int *)malloc(nel * ndof * sizeof(int));
+ return 0;
+}
+
+int mesh_free(mesh_t *mesh)
+{
+ if (mesh->connect != NULL)
+ free(mesh->connect);
+
+ if (mesh->coords != NULL)
+ free(mesh->coords);
+
+ return 0;
+}
+
+
+
+int mesh_open(mesh_t *mesh,
+ const char *filename,
+ const char *path)
+{
+ hid_t file_id;
+ hid_t model_id;
+ int ierr;
+
+ /* open the file in read-write mode */
+ file_id = file_open(filename, "rw");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ model_id = H5Gopen(file_id, path);
+ if (model_id < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ //
+ // TODO: read attributes in model_id to determine the following paths
+ //
+
+ /* read all coords */
+ ierr = mesh_read_coords(mesh, model_id, "geometry/coords");
+ if (ierr < 0)
+ {
+ H5Gclose(model_id);
+ H5Fclose(file_id);
+ return -3;
+ }
+ assert(mesh->nsd == 3);
+
+ /* read all elements */
+ ierr = mesh_read_connect(mesh, model_id, "topology/simplices");
+ if (ierr < 0)
+ {
+ H5Gclose(model_id);
+ H5Fclose(file_id);
+ return -4;
+ }
+
+ H5Gclose(model_id);
+ H5Fclose(file_id);
+ return 0;
+}
+
+int mesh_open2(mesh_t *mesh,
+ const char *filename,
+ const char *coords_path,
+ const char *connect_path)
+{
+ hid_t file_id;
+ int ierr;
+
+ /* open the file in read-write mode */
+ file_id = file_open(filename, "rw");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ /* read all coords */
+ ierr = mesh_read_coords(mesh, file_id, coords_path);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+ assert(mesh->nsd == 3);
+
+ /* read all elements */
+ ierr = mesh_read_connect(mesh, file_id, connect_path);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ return -3;
+ }
+
+ H5Fclose(file_id);
+ return 0;
+}
+
+
+
+int mesh_read_coords(mesh_t *mesh, hid_t loc_id, const char *coords_path)
+{
+ return dataset_read2(loc_id, coords_path, H5T_NATIVE_DOUBLE,
+ (void **)&(mesh->coords), &(mesh->nno), &(mesh->nsd));
+}
+
+int mesh_read_connect(mesh_t *mesh, hid_t loc_id, const char *connect_path)
+{
+
+ return dataset_read2(loc_id, connect_path, H5T_NATIVE_INT,
+ (void **)&(mesh->connect), &(mesh->nel), &(mesh->ndof));
+}
+
Added: cs/cigma/branches/cigma-0.9/src/points.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/points.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/points.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,225 @@
+#include <stdlib.h>
+#include <hdf5.h>
+#include <cigma/points.h>
+#include <cigma/array.h>
+#include <cigma/io.h>
+
+
+
+int points_init(points_t *points, int nel, int npts, int ndim, double *x)
+{
+ points->nel = nel;
+ points->npts = npts;
+ points->ndim = ndim;
+
+ points->start = 0;
+ points->end = 0;
+ points->dset = -1;
+
+ points->x = x;
+
+ if (x == NULL)
+ {
+ points->x = (double *)malloc(nel * npts * ndim * sizeof(double));
+ }
+
+ return 0;
+}
+
+int points_free_dset(points_t *points)
+{
+ if (points->dset >= 0)
+ {
+ H5Dclose(points->dset);
+ points->dset = -1;
+ }
+ return 0;
+}
+
+int points_free_data(points_t *points)
+{
+ if (points->x != NULL)
+ {
+ free(points->x);
+ points->x = NULL;
+ }
+ return 0;
+}
+
+int points_free(points_t *points)
+{
+ points_free_dset(points);
+ points_free_data(points);
+
+ points->nel = 0;
+ points->npts = 0;
+ points->ndim = 0;
+
+ points->start = 0;
+ points->end = 0;
+
+ return 0;
+}
+
+
+
+int points_open(points_t *points, const char *filename, const char *path)
+{
+ hid_t file_id;
+ int ierr;
+
+ file_id = file_open(filename, "rw");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ ierr = points_open_dset(points, file_id, path);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ ierr = points_read(points);
+ if (ierr < 0)
+ {
+ points_free_dset(points);
+ H5Fclose(file_id);
+ return -3;
+ }
+
+ H5Fclose(file_id);
+ return 0;
+}
+
+int points_create(points_t *points, const char *filename, const char *path, const char *title)
+{
+ hid_t file_id;
+ int ierr;
+
+ file_id = file_create(filename, "x");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ ierr = points_create_dset(points, file_id, path, title);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ ierr = points_write(points);
+ if (ierr < 0)
+ {
+ points_free_dset(points);
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ H5Fclose(file_id);
+ return 0;
+}
+
+
+
+int points_open_dset(points_t *points, hid_t loc_id, const char *name)
+{
+ array_t a;
+ array_init(&a, H5T_NATIVE_DOUBLE, 3, NULL, NULL);
+ points->dset = array_open(&a, loc_id, name);
+ if (points->dset < 0)
+ {
+ array_free_dims(&a);
+ return -1;
+ }
+ array_dims3(&a, &(points->nel), &(points->npts), &(points->ndim));
+ array_free_dims(&a);
+ return 0;
+}
+
+int points_create_dset(points_t *points, hid_t loc_id, const char *name, const char *title)
+{
+ array_t a;
+ int shape[3] = {points->nel, points->npts, points->ndim};
+ array_init(&a, H5T_NATIVE_DOUBLE, 3, shape, points->x);
+ points->dset = array_create(&a, loc_id, name, title);
+ if (points->dset < 0)
+ {
+ array_free_dims(&a);
+ return -1;
+ }
+ array_free_dims(&a);
+ return 0;
+}
+
+
+
+int points_read(points_t *points)
+{
+ return points_read_slice(points, 0, points->nel);
+}
+
+int points_write(points_t *points)
+{
+ return points_write_slice(points, 0, points->nel);
+}
+
+
+
+int points_read_slice(points_t *points, int start, int end)
+{
+ array_t b;
+ int shape[3] = {points->nel, points->npts, points->ndim};
+ int ierr;
+
+ ierr = array_init(&b, H5T_NATIVE_DOUBLE, 3, shape, NULL);
+ ierr = array_read_slice(&b, points->dset, start, end);
+ if (ierr < 0)
+ {
+ array_free_dims(&b);
+ return -1;
+ }
+
+ if (b.data == NULL)
+ {
+ array_free_dims(&b);
+ return -2;
+ }
+
+ points->x = (double *) b.data;
+ points->start = start;
+ points->end = end;
+
+ array_free_dims(&b);
+ return 0;
+}
+
+int points_write_slice(points_t *points, int start, int end)
+{
+ array_t b;
+ int shape[3] = {points->nel, points->npts, points->ndim};
+ int ierr;
+
+ if (points->x == NULL)
+ {
+ return -1;
+ }
+
+ array_init(&b, H5T_NATIVE_DOUBLE, 3, shape, (void *)points->x);
+ ierr = array_write_slice(&b, points->dset, start, end);
+ if (ierr < 0)
+ {
+ array_free_dims(&b);
+ return -2;
+ }
+
+ points->start = start;
+ points->end = end;
+
+ array_free_dims(&b);
+ return 0;
+}
+
Added: cs/cigma/branches/cigma-0.9/src/rule.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/rule.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/rule.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,225 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <hdf5.h>
+#include <cigma/io.h>
+#include <cigma/rule.h>
+#include <cigma/dataset.h>
+
+
+
+int rule_init(rule_t *qr, int nq, int nsd)
+{
+ qr->nq = nq;
+ qr->nsd = nsd;
+ qr->weights = (double *)malloc(nq * sizeof(double));
+ qr->points = (double *)malloc(nq * nsd * sizeof(double));
+ return 0;
+}
+
+int rule_free(rule_t *qr)
+{
+ free(qr->weights);
+ free(qr->points);
+ return 0;
+}
+
+
+
+int rule_open(rule_t *qr, const char *filename, const char *path)
+{
+ hid_t file_id;
+ hid_t group_id;
+ herr_t status;
+ int ierr;
+
+ /* for now, zero out nq & nsd */
+ qr->nq = 0;
+ qr->nsd = 0;
+
+ /* open the file in read-write mode */
+ file_id = file_open(filename, "rw");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ /* open group at specified path */
+ group_id = H5Gopen(file_id, path);
+ if (group_id < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ /* read quadrature weights & points */
+ ierr = rule_read(qr, group_id);
+ if (ierr < 0)
+ {
+ status = H5Gclose(group_id);
+ status = H5Fclose(file_id);
+ return -3;
+ }
+
+ /* finalize */
+ status = H5Gclose(group_id);
+ status = H5Fclose(file_id);
+
+ return 0;
+}
+
+int rule_open_txt(rule_t *qr, const char *filename)
+{
+ FILE *fp;
+ int ret;
+
+ int q, nq;
+ int i, nsd;
+
+ fp = fopen(filename, "r");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ ret = fscanf(fp, "%d", &nq);
+ ret = fscanf(fp, "%d", &nsd);
+
+ rule_init(qr, nq, nsd);
+
+ for (q = 0; q < nq; q++)
+ {
+ ret = fscanf(fp, "%lf", &(qr->weights[q]));
+ for (i = 0; i < nsd; i++)
+ {
+ ret = fscanf(fp, "%lf", &(qr->points[nsd*q + i]));
+ }
+ }
+
+ fclose(fp);
+ return 0;
+}
+
+
+
+int rule_create(rule_t *qr, const char *filename, const char *path)
+{
+ hid_t file_id;
+ hid_t group_id;
+ int ierr;
+
+ file_id = file_create(filename, "x");
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ group_id = touch_path(file_id, path);
+ if (group_id < 0)
+ {
+ H5Fclose(file_id);
+ return -2;
+ }
+
+ ierr = rule_write(qr, group_id);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ H5Gclose(group_id);
+ return -3;
+ }
+
+ H5Fclose(file_id);
+ H5Gclose(group_id);
+ return 0;
+}
+
+int rule_create_txt(rule_t *qr, const char *filename)
+{
+ FILE *fp;
+ int q, nq;
+ int i, nsd;
+
+ fp = fopen(filename, "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ nsd = qr->nsd;
+ nq = qr->nq;
+
+ fprintf(fp, "%d %d\n", nq, nsd);
+
+ for (q = 0; q < nq; q++)
+ {
+ fprintf(fp, "%lf", qr->weights[q]);
+ for (i = 0; i < nsd; i++)
+ {
+ fprintf(fp, " %lf", qr->points[nsd*q + i]);
+ }
+ }
+
+ fclose(fp);
+ return 0;
+}
+
+
+
+int rule_read(rule_t *qr, hid_t loc_id)
+{
+ return rule_read2(qr, loc_id, "weights", "points");
+}
+
+int rule_write(rule_t *qr, hid_t loc_id)
+{
+ return rule_write2(qr, loc_id, "weights", "points");
+}
+
+
+int rule_read2(rule_t *qr, hid_t loc_id, const char *weights_path, const char *points_path)
+{
+ int nq;
+ int ierr;
+
+ ierr = dataset_read1(loc_id, weights_path, H5T_NATIVE_DOUBLE,
+ (void **)&(qr->weights), &nq);
+ if (ierr < 0)
+ {
+ return -1;
+ }
+
+ ierr = dataset_read2(loc_id, points_path, H5T_NATIVE_DOUBLE,
+ (void **)&(qr->points), &(qr->nq), &(qr->nsd));
+ if (ierr < 0)
+ {
+ return -2;
+ }
+ if (nq != qr->nq)
+ {
+ return -3;
+ }
+
+ return 0;
+}
+
+int rule_write2(rule_t *qr, hid_t loc_id, const char *weights_path, const char *points_path)
+{
+ int ierr;
+
+ ierr = dataset_write1(loc_id, weights_path, "quadrature weights",
+ H5T_NATIVE_DOUBLE, qr->weights, qr->nq);
+ if (ierr < 0)
+ {
+ return -1;
+ }
+
+ ierr = dataset_write2(loc_id, points_path, "quadrature points",
+ H5T_NATIVE_DOUBLE, qr->points, qr->nq, qr->nsd);
+ if (ierr < 0)
+ {
+ return -2;
+ }
+
+ return 0;
+}
+
Added: cs/cigma/branches/cigma-0.9/src/split.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/split.c 2007-06-26 16:41:33 UTC (rev 7498)
+++ cs/cigma/branches/cigma-0.9/src/split.c 2007-06-26 16:42:39 UTC (rev 7499)
@@ -0,0 +1,184 @@
+#include <stdlib.h>
+#include <string.h>
+#include <cigma/split.h>
+
+
+/*
+ * Data structure for our basic string tokenizer. Basically, the strategy
+ * is to copy the desired string into our buffer, where we can zero out
+ * the separator characters in-place. Using a struct allows us to avoid
+ * defining global variables like in the version of split.c from the project
+ * http://www.nongnu.org/uri/ on which this code is inspired.
+ *
+ */
+
+typedef struct
+{
+ char **tokens; // tokens array
+ int token_count; // real size of tokens array
+ int max_token_count; // size of tokens array
+
+ char *buffer; // buffer string
+ int buffer_length; // length of buffer string
+
+} strtok_t;
+
+
+
+/*
+ * In this section, we define the methods for our string tokenizer
+ * object. Namely, a constructor, a destructor, and then the actual
+ * routine to split the string.
+ *
+ */
+
+static void strtok_init(strtok_t *st)
+{
+ /* starting up with 0 tokens */
+ st->token_count = 0;
+
+ /* setup the initial array sizes */
+ st->max_token_count = 16; // 16 tokens in array
+ st->buffer_length = 512; // 512 chars in buffer
+
+ /* allocate enough tokens and initialize buffer */
+ st->buffer = (char *)malloc(st->buffer_length * sizeof(char));
+ st->tokens = (char **)malloc(st->max_token_count * sizeof(char *));
+}
+
+static void strtok_free(strtok_t *st)
+{
+ if (st != NULL)
+ {
+ free(st->buffer);
+ free(st->tokens);
+ }
+}
+
+static void strtok_split(strtok_t *st, const char *s, int len, char separator)
+{
+ /*
+ * First, check whether our buffer is large enough to manipulate
+ * the string s, and if not, reallocate enough memory.
+ */
+ if (st->buffer_length < len)
+ {
+ st->buffer_length = (len < 512) ? 512 : len+1;
+ st->buffer = (char *)realloc(st->buffer, st->buffer_length * sizeof(char));
+ }
+
+ /*
+ * Next, copy the string s into our buffer and tokenize it in-place.
+ * Essentially, zero out the locations where we find the separator
+ * character, while remembering the beginning of each string.
+ */
+ memcpy(st->buffer, s, len);
+ st->buffer[len] = '\0';
+ {
+
+ char *first, *p;
+ int index, last;
+
+ first = st->buffer;
+ last = st->buffer_length - 1;
+
+ /* remove trailing separators */
+ while (last >= 0 && st->buffer[last] == separator)
+ {
+ st->buffer[last] = '\0';
+ last--;
+ }
+
+ /* remove leading separators */
+ while (*first == separator)
+ {
+ first++;
+ }
+
+ /* store first token */
+ index = 0;
+ st->tokens[index++] = first;
+
+ /* keep tokenizing the buffer */
+ for (p = strchr(first, separator);
+ p != NULL;
+ p = strchr(p, separator))
+ {
+ /* separator found -- zero it out */
+ *p = '\0';
+
+ /* make p point to next char */
+ p++;
+
+ /* store the next token */
+ if ((*p != separator) && (*p != '\0'))
+ {
+ st->tokens[index++] = p;
+
+ /* check whether we need to expand our tokens array,
+ * to make room for the next batch of tokens
+ */
+ if (index >= (st->max_token_count))
+ {
+ st->max_token_count += 16;
+ st->tokens = (char **)realloc(st->tokens,
+ st->max_token_count
+ * sizeof(char *));
+ }
+ }
+ }
+
+ /* store the final count */
+ st->token_count = index;
+ }
+
+ return;
+}
+
+
+
+/*
+ * Finally, we provide a procedural interface to our string tokenizer.
+ * The caller subsumes the responsibility of freeing the newly allocated
+ * list, as well as each individual string in that list.
+ */
+void split(const char *str, int len,
+ char ***split_list, int *split_count,
+ char sep)
+{
+ int i;
+ strtok_t tok;
+
+ if (sep == '\0')
+ {
+ *split_list = NULL;
+ *split_count = 0;
+ return;
+ }
+
+ strtok_init(&tok);
+ strtok_split(&tok, str, len, sep);
+
+ *split_list = (char **)malloc(tok.token_count * sizeof(char *));
+ *split_count = tok.token_count;
+
+ for (i = 0; i < tok.token_count; i++)
+ {
+ (*split_list)[i] = strdup(tok.tokens[i]);
+ }
+
+ strtok_free(&tok);
+}
+
+void split_free(char **split_list, int split_count)
+{
+ int i;
+ if (split_list != NULL)
+ {
+ for (i = 0; i < split_count; i++)
+ {
+ free(split_list[i]);
+ }
+ free(split_list);
+ }
+}
More information about the cig-commits
mailing list