[cig-commits] r6775 - in cs/cigma/trunk/sandbox: . qpts
luis at geodynamics.org
luis at geodynamics.org
Thu May 3 17:14:54 PDT 2007
Author: luis
Date: 2007-05-03 17:14:53 -0700 (Thu, 03 May 2007)
New Revision: 6775
Added:
cs/cigma/trunk/sandbox/qpts/
cs/cigma/trunk/sandbox/qpts/map.c
cs/cigma/trunk/sandbox/qpts/qr.h5
Log:
Utility to map quadrature points from a HDF5 mesh
Added: cs/cigma/trunk/sandbox/qpts/map.c
===================================================================
--- cs/cigma/trunk/sandbox/qpts/map.c 2007-05-03 23:33:53 UTC (rev 6774)
+++ cs/cigma/trunk/sandbox/qpts/map.c 2007-05-04 00:14:53 UTC (rev 6775)
@@ -0,0 +1,277 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include "hdf5.h"
+
+const char *filename = "bmssnog-250m-pylith-lintet.h5";
+const int nel = 5308416;
+const int nno = 912673;
+
+//
+// In this comment, the following notation applies:
+//
+// x = (a[0], b[0], c[0], d[0])
+// y = (a[1], b[1], c[1], d[1])
+// z = (a[2], b[2], c[2], d[2])
+//
+// The translation part of the reference map acts with the following
+// point as the origin:
+//
+// origin = [[(x1 + x2 + x3 - x0)/2],
+// [(y1 + y2 + y3 - y0)/2],
+// [(z1 + z2 + z3 - z0)/2]]
+//
+// The rotation & shearing part of the reference map is represented by
+// the transpose of the jacobian matrix:
+//
+// J^T = [[(x1-x0)/2, (y1-y0)/2, (x1-x0)/2],
+// [(x2-x0)/2, (y2-y0)/2, (y2-x0)/2],
+// [(x3-x0)/2, (y3-y0)/2, (z3-x0)/2]]
+//
+// The reference map takes a point (xi,eta,zeta) on a reference element
+// to the point (x,y,z) in the physical element, acting in the following
+// way:
+//
+// [x] [ xi ]
+// [y] = origin + (J^T) * [ eta ]
+// [z] [ zeta ]
+//
+
+typedef struct
+{
+ double matrix[9];
+ double origin[3];
+} refmap_t;
+
+
+void calculate_reference_map(double *a, double *b, double *c, double *d, refmap_t *X)
+{
+ X->origin[0] = 0.5*(b[0] + c[0] + d[0] - a[0]);
+ X->origin[1] = 0.5*(b[1] + c[1] + d[1] - a[1]);
+ X->origin[2] = 0.5*(b[2] + c[2] + d[2] - a[2]);
+
+ X->matrix[3*0 + 0] = 0.5*(b[0] - a[0]);
+ X->matrix[3*0 + 1] = 0.5*(b[1] - a[1]);
+ X->matrix[3*0 + 2] = 0.5*(b[2] - a[2]);
+
+ X->matrix[3*1 + 0] = 0.5*(c[0] - a[0]);
+ X->matrix[3*1 + 1] = 0.5*(c[1] - a[1]);
+ X->matrix[3*1 + 2] = 0.5*(c[2] - a[2]);
+
+ X->matrix[3*2 + 0] = 0.5*(d[0] - a[0]);
+ X->matrix[3*2 + 1] = 0.5*(d[1] - a[1]);
+ X->matrix[3*2 + 2] = 0.5*(d[2] - a[2]);
+
+ return;
+}
+
+void apply_reference_map(refmap_t *X, double *xi, double *x)
+{
+ x[0] = X->origin[0] + (X->matrix[3*0+0])*xi[0] + (X->matrix[3*0+1])*xi[1] + (X->matrix[3*0+2])*xi[2];
+ x[1] = X->origin[1] + (X->matrix[3*1+0])*xi[0] + (X->matrix[3*1+1])*xi[1] + (X->matrix[3*1+2])*xi[2];
+ x[2] = X->origin[2] + (X->matrix[3*2+0])*xi[0] + (X->matrix[3*2+1])*xi[1] + (X->matrix[3*2+2])*xi[2];
+ return;
+}
+
+
+//
+// The following functions act on a mesh structure
+//
+//
+
+typedef struct
+{
+ int nno;
+ int nel;
+ double *coords;
+ int *elements;
+} mesh_t;
+
+int mesh_initialize(mesh_t *mesh)
+{
+ hid_t file_id;
+ hid_t coords_dset_id;
+ hid_t connect_dset_id;
+ herr_t status;
+
+ mesh->nno = nno;
+ mesh->nel = nel;
+
+ /* open the file in read-only mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ fprintf(stderr, "Could not open file %s\n", filename);
+ return -1;
+ }
+
+ /* read all coordinates */
+ coords_dset_id = H5Dopen(file_id, "/model/geometry/coordinates");
+ mesh->coords = (double *)malloc((mesh->nno) * 3 * sizeof(double));
+ status = H5Dread(coords_dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, mesh->coords);
+ status = H5Dclose(coords_dset_id);
+
+ /* read all elements */
+ connect_dset_id = H5Dopen(file_id, "/model/topology/simplices");
+ mesh->elements = (int *)malloc((mesh->nel) * 4 * sizeof(int));
+ status = H5Dread(connect_dset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, mesh->elements);
+ status = H5Dclose(connect_dset_id);
+
+ /* close the file */
+ status = H5Fclose(file_id);
+
+ return 0;
+}
+
+void mesh_finalize(mesh_t *mesh)
+{
+ free(mesh->coords);
+ free(mesh->elements);
+}
+
+//
+// The following functions act on a quadrature rule
+//
+
+typedef struct
+{
+ int n;
+ double *weights;
+ double *points;
+} qrule_t;
+
+#ifndef QRULE
+#define QRULE "/FIAT/cools/tet_4"
+#endif
+int qrule_initalize(qrule_t *qrule)
+{
+ hid_t file_id;
+ hid_t qrule_group_id;
+ hid_t points_dset_id;
+ hid_t weights_dset_id;
+ hid_t weights_dspace_id;
+ herr_t status;
+
+ int nq;
+
+ /* open the file in read-only mode */
+ file_id = H5Fopen("qr.h5", H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ fprintf(stderr, "Could not open file \"qr.h5\"");
+ return -1;
+ }
+
+ /* open the desired rule */
+ qrule_group_id = H5Gopen(file_id, QRULE);
+
+ /*
+ * read quadrature weights
+ */
+
+ weights_dset_id = H5Dopen(qrule_group_id, "weights");
+
+ weights_dspace_id = H5Dget_space(weights_dset_id);
+ nq = H5Sget_simple_extent_npoints(weights_dspace_id);
+ qrule->n = nq;
+ status = H5Sclose(weights_dspace_id);
+
+ qrule->weights = (double *)malloc(nq * sizeof(double));
+ status = H5Dread(weights_dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, qrule->weights);
+ status = H5Dclose(weights_dset_id);
+
+ /*
+ * read quadrature points
+ */
+
+ points_dset_id = H5Dopen(qrule_group_id, "points");
+ qrule->points = (double *)malloc(nq * 3 * sizeof(double));
+ status = H5Dread(points_dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, qrule->points);
+ status = H5Dclose(points_dset_id);
+
+ /* finalize */
+ status = H5Gclose(qrule_group_id);
+ status = H5Fclose(file_id);
+
+ return 0;
+}
+
+void qrule_finalize(qrule_t *qrule)
+{
+ free(qrule->points);
+ free(qrule->weights);
+}
+
+//
+// Finally, just execute
+//
+
+int main(int argc, char *argv[])
+{
+ int a,b,c,d; // node ids for a given element
+ int e; // for iterating through elements
+ int q; // for iterating through quadrature points
+
+ double x[3]; // for the actual integration point
+
+ mesh_t mesh;
+ refmap_t map;
+ qrule_t qrule;
+
+ int err;
+
+ /* initialize the mesh */
+ err = mesh_initialize(&mesh);
+ if (err < 0)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* read the quadrature rule */
+ err = qrule_initalize(&qrule);
+ if (err < 0)
+ {
+ return EXIT_FAILURE;
+ }
+
+ /* iterate through all elements in the mesh */
+ for (e = 0; e < mesh.nel; e++)
+ {
+ a = mesh.elements[4*e + 0];
+ b = mesh.elements[4*e + 1];
+ c = mesh.elements[4*e + 2];
+ d = mesh.elements[4*e + 3];
+
+ calculate_reference_map(&mesh.coords[3*(a-1)],
+ &mesh.coords[3*(b-1)],
+ &mesh.coords[3*(c-1)],
+ &mesh.coords[3*(d-1)],
+ &map);
+
+ for (q = 0; q < qrule.n; q++)
+ {
+ apply_reference_map(&map, &qrule.points[3*q], x);
+
+ //*
+ printf("e%d=(%d,%d,%d,%d): ", e, a, b, c, d);
+ printf("q%d=(%lf,%lf,%lf) -> ", q, qrule.points[3*q+0], qrule.points[3*q+1], qrule.points[3*q+2]);
+ printf("(%lf,%lf,%lf)\n", x[0], x[1], x[2]);
+ // */
+
+ /*
+ printf("\t [%lf,%lf,%lf],\n", map.origin[0], map.origin[1], map.origin[2]);
+ printf("\t[[%lf,%lf,%lf],\n", map.matrix[0], map.matrix[1], map.matrix[2]);
+ printf("\t [%lf,%lf,%lf],\n", map.matrix[3], map.matrix[4], map.matrix[5]);
+ printf("\t [%lf,%lf,%lf]]\n", map.matrix[6], map.matrix[7], map.matrix[8]);
+ // */
+ }
+ }
+
+ mesh_finalize(&mesh);
+ qrule_finalize(&qrule);
+
+ return EXIT_SUCCESS;
+}
Added: cs/cigma/trunk/sandbox/qpts/qr.h5
===================================================================
(Binary files differ)
Property changes on: cs/cigma/trunk/sandbox/qpts/qr.h5
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
More information about the cig-commits
mailing list