[cig-commits] r6819 - in cs/cigma/trunk: include/cigma src tests
luis at geodynamics.org
luis at geodynamics.org
Wed May 9 17:02:46 PDT 2007
Author: luis
Date: 2007-05-09 17:02:46 -0700 (Wed, 09 May 2007)
New Revision: 6819
Added:
cs/cigma/trunk/include/cigma/qr.h
cs/cigma/trunk/src/qr.c
cs/cigma/trunk/tests/qr.h5
cs/cigma/trunk/tests/test_qr.c
Log:
Added functions related to quadrature rules. Also, included
test data in qr.h5
Added: cs/cigma/trunk/include/cigma/qr.h
===================================================================
--- cs/cigma/trunk/include/cigma/qr.h 2007-05-10 00:00:09 UTC (rev 6818)
+++ cs/cigma/trunk/include/cigma/qr.h 2007-05-10 00:02:46 UTC (rev 6819)
@@ -0,0 +1,32 @@
+#ifndef __CIGMA_QR_H__
+#define __CIGMA_QR_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include <hdf5.h>
+
+typedef struct
+{
+ int nq;
+ int nsd;
+ double *weights;
+ double *points;
+} rule_t;
+
+typedef void (*generator_t)(int, int, double *, double *);
+typedef double (*functor_t)(double *);
+
+int qr_init(rule_t *qr, char *filename, char *path);
+int qr_init_h5(rule_t *qr, hid_t loc_id, char *path_weights, char *path_points);
+int qr_init_txt(rule_t *qr, char *filename);
+int qr_init_gen(rule_t *qr, generator_t gen);
+int qr_free(rule_t *qr);
+
+void qr_apply(rule_t *qr, functor_t f, double *sum);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif
Added: cs/cigma/trunk/src/qr.c
===================================================================
--- cs/cigma/trunk/src/qr.c 2007-05-10 00:00:09 UTC (rev 6818)
+++ cs/cigma/trunk/src/qr.c 2007-05-10 00:02:46 UTC (rev 6819)
@@ -0,0 +1,139 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <hdf5.h>
+#include <cigma/dataset.h>
+#include <cigma/qr.h>
+
+int qr_init(rule_t *qr, char *filename, 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-only mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+ 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 = qr_init_h5(qr, group_id, "weights", "points");
+ 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 qr_init_h5(rule_t *qr, hid_t loc_id, char *path_weights, char *path_points)
+{
+ int nq;
+ int ierr;
+
+ ierr = dataset_read1(loc_id, path_weights, H5T_NATIVE_DOUBLE,
+ (void **)&(qr->weights), &nq);
+ if (ierr < 0)
+ {
+ return -1;
+ }
+
+ ierr = dataset_read2(loc_id, path_points, H5T_NATIVE_DOUBLE,
+ (void **)&(qr->points), &(qr->nq), &(qr->nsd));
+ if (ierr < 0)
+ {
+ return -2;
+ }
+ if (nq != qr->nq)
+ {
+ return -3;
+ }
+
+ return 0;
+}
+
+int qr_init_txt(rule_t *qr, char *filename)
+{
+ int q, nq;
+ int i, nsd;
+
+ FILE *fp;
+ int ret;
+
+ fp = fopen(filename, "r");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ ret = fscanf(fp, "%d", &nq);
+ ret = fscanf(fp, "%d", &nsd);
+
+ qr->nq = nq;
+ qr->nsd = nsd;
+ qr->weights = (double *)malloc(nq * sizeof(double));
+ qr->points = (double *)malloc(nq * nsd * sizeof(double));
+
+ 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 qr_init_gen(rule_t *qr, generator_t gen)
+{
+ gen(qr->nq, qr->nsd, qr->weights, qr->points);
+ return 0;
+}
+
+int qr_free(rule_t *qr)
+{
+ free(qr->weights);
+ free(qr->points);
+ return 0;
+}
+
+
+void qr_apply(rule_t *qr, double (*f)(double *), double *sum)
+{
+ int q; // for looping over quadrature points
+ int nq = qr->nq; // number of quadrature points
+ int nsd = qr->nsd; // dimension of points
+
+ double *x = qr->points;
+ double *w = qr->weights;
+ double acc = 0.0;
+
+ for (q = 0; q < nq; q++)
+ acc += w[q] * f(&x[nsd*q]);
+
+ *sum = acc;
+}
+
Added: cs/cigma/trunk/tests/qr.h5
===================================================================
(Binary files differ)
Property changes on: cs/cigma/trunk/tests/qr.h5
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: cs/cigma/trunk/tests/test_qr.c
===================================================================
--- cs/cigma/trunk/tests/test_qr.c 2007-05-10 00:00:09 UTC (rev 6818)
+++ cs/cigma/trunk/tests/test_qr.c 2007-05-10 00:02:46 UTC (rev 6819)
@@ -0,0 +1,43 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <cigma/qr.h>
+
+int main(int argc, char *argv[])
+{
+ int n, q;
+ int ierr;
+ rule_t qr;
+
+ char *rules[] = {
+ "/FIAT/cools/tet_1",
+ "/FIAT/cools/tet_2",
+ "/FIAT/cools/tet_3",
+ "/FIAT/cools/tet_4",
+ };
+
+
+ for (n = 0; n < 4; n++)
+ {
+ ierr = qr_init(&qr, "qr.h5", rules[n]);
+
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Could not find qr.h5");
+ return EXIT_FAILURE;
+ }
+
+ printf("Reading rule %s with %d points on %d dimensions\n",
+ rules[n], qr.nq, qr.nsd);
+
+ for (q = 0; q < qr.nq; q++)
+ {
+ printf("w[%2d] = %lf; x[%2d] = (%+lf, %+lf, %+lf)\n",
+ q, qr.weights[q],
+ q, qr.points[3*q+0], qr.points[3*q+1], qr.points[3*q+2]);
+ }
+ printf("\n");
+
+ qr_free(&qr);
+ }
+ return EXIT_SUCCESS;
+}
More information about the cig-commits
mailing list