[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