[cig-commits] r6824 - in cs/cigma/trunk: include/cigma src tests
luis at geodynamics.org
luis at geodynamics.org
Wed May 9 17:17:10 PDT 2007
Author: luis
Date: 2007-05-09 17:17:10 -0700 (Wed, 09 May 2007)
New Revision: 6824
Added:
cs/cigma/trunk/include/cigma/field.h
cs/cigma/trunk/src/field.c
cs/cigma/trunk/tests/test_field.c
Log:
Functions for reading and interpreting a solution field
Added: cs/cigma/trunk/include/cigma/field.h
===================================================================
--- cs/cigma/trunk/include/cigma/field.h 2007-05-10 00:14:32 UTC (rev 6823)
+++ cs/cigma/trunk/include/cigma/field.h 2007-05-10 00:17:10 UTC (rev 6824)
@@ -0,0 +1,36 @@
+#ifndef __CIGMA_FIELD_H__
+#define __CIGMA_FIELD_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+#include <hdf5.h>
+#include <cigma/fe.h>
+#include <cigma/mesh.h>
+
+
+typedef int (*extern_field)(int npts, double *x, int nsd, double *y, int dims);
+
+
+typedef struct {
+ int nno; // number of nodes
+ int ndim; // number of components in field: scalar=1 vector=2,3 tensor=6
+ mesh_t *mesh; // pointer to parent mesh object
+ fe_space_t *fe; // finite element space -- metadata for interpreting dof[]
+ extern_field eval; // function pointer for externally defined fields
+ double *dofs; // [nno x ndim] matrix
+} field_t;
+
+
+int field_init(field_t *field, char *filename, char *path);
+int field_init_h5(field_t *field, hid_t loc_id, const char *name);
+int field_free(field_t *field);
+
+void field_dofs(field_t *field, int e, double *dof);
+void field_eval(field_t *B, double *x, int nq, double *g, int dims);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif
Added: cs/cigma/trunk/src/field.c
===================================================================
--- cs/cigma/trunk/src/field.c 2007-05-10 00:14:32 UTC (rev 6823)
+++ cs/cigma/trunk/src/field.c 2007-05-10 00:17:10 UTC (rev 6824)
@@ -0,0 +1,107 @@
+#include <stdlib.h>
+#include <hdf5.h>
+#include <cigma/array.h>
+#include <cigma/field.h>
+#include <cigma/tet4.h>
+
+
+int field_init(field_t *field, char *filename, char *path)
+{
+ hid_t file_id;
+ herr_t status;
+ int ierr;
+
+ /* open the file in read-only mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+ if (file_id < 0)
+ {
+ return -1;
+ }
+
+ ierr = field_init_h5(field, file_id, path);
+ if (ierr < 0)
+ {
+ H5Fclose(file_id);
+ }
+
+ status = H5Fclose(file_id);
+ return 0;
+}
+
+int field_init_h5(field_t *field, hid_t loc_id, const char *name)
+{
+ array_t a;
+ int ierr;
+
+ array_init(&a, H5T_NATIVE_DOUBLE, 2, NULL, NULL);
+
+ ierr = array_read(&a, loc_id, name);
+ if (ierr < 0)
+ {
+ array_free_dims(&a);
+ return -1;
+ }
+
+ field->dofs = (double *) a.data;
+ array_dims2(&a, &(field->nno), &(field->ndim));
+ array_free_dims(&a);
+
+ return 0;
+}
+
+int field_free(field_t *field)
+{
+ if (field->dofs != NULL)
+ free(field->dofs);
+
+ field->nno = 0;
+ field->ndim = 0;
+ field->mesh = NULL;
+
+ return 0;
+}
+
+
+
+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;
+ }
+ }
+ }
+}
Added: cs/cigma/trunk/tests/test_field.c
===================================================================
--- cs/cigma/trunk/tests/test_field.c 2007-05-10 00:14:32 UTC (rev 6823)
+++ cs/cigma/trunk/tests/test_field.c 2007-05-10 00:17:10 UTC (rev 6824)
@@ -0,0 +1,99 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <cigma/mesh.h>
+#include <cigma/field.h>
+
+int main(int argc, char *argv[])
+{
+ mesh_t mesh;
+ field_t displacement;
+ field_t velocity;
+
+ int e,n;
+ int i,j;
+ double *coords;
+
+ int ierr;
+
+ if (argc < 2)
+ {
+ printf("Usage: %s model.h5\n", argv[0]);
+ return EXIT_FAILURE;
+ }
+
+ ierr = mesh_init(&mesh, argv[1], "/model/geometry/coordinates", "/model/topology/simplices");
+ if (ierr < 0)
+ {
+ printf("Could not read mesh from '%s'!\n", argv[1]);
+ return EXIT_FAILURE;
+ }
+
+ ierr = field_init(&displacement, argv[1], "/model/solution/snapshot0/displacement");
+ if (ierr < 0)
+ {
+ printf("Could not read displacement field from '%s'!\n", argv[1]);
+ return EXIT_FAILURE;
+ }
+
+ ierr = field_init(&velocity, argv[1], "/model/solution/snapshot0/velocity");
+ if (ierr < 0)
+ {
+ printf("Could not read velocity field from '%s'!\n", argv[1]);
+ return EXIT_FAILURE;
+ }
+
+ printf("Reading file %s\n", argv[1]);
+ printf("Found %d nodes and %d elements\n", mesh.nno, mesh.nel);
+
+ coords = (double *)malloc((mesh.ndof * 3) * sizeof(double));
+ for (e = 100; e < 105; e++)
+ {
+ mesh_coords(&mesh, e, coords);
+
+ printf("coords(e=%d) = {", e);
+ for (i = 0; i < mesh.ndof; i++)
+ {
+ printf("{");
+ for (j = 0; j < 3; j++)
+ {
+ printf("%g", coords[3*i + j]);
+ if (j != 2)
+ printf(",");
+ }
+ printf("}");
+ if (i != (mesh.ndof - 1))
+ printf(",");
+ }
+ printf("};\n");
+ }
+ free(coords);
+
+ for (n = 500; n < 515; n++)
+ {
+ printf("displacement(n=%d) = {", n);
+ for (i = 0; i < displacement.ndim; i++)
+ {
+ printf("%g", displacement.dofs[displacement.ndim * n + i]);
+ if (i != 2)
+ printf(",");
+ }
+ printf("};\n");
+ }
+
+ for (n = 500; n < 515; n++)
+ {
+ printf("velocity(n=%d) = {", n);
+ for (i = 0; i < velocity.ndim; i++)
+ {
+ printf("%g", velocity.dofs[velocity.ndim * n + i]);
+ if (i != 2)
+ printf(",");
+ }
+ printf("};\n");
+ }
+
+ mesh_free(&mesh);
+ field_free(&displacement);
+ field_free(&velocity);
+ return EXIT_SUCCESS;
+}
More information about the cig-commits
mailing list