[cig-commits] r6897 - cs/cigma/trunk/tools
luis at geodynamics.org
luis at geodynamics.org
Wed May 16 10:19:21 PDT 2007
Author: luis
Date: 2007-05-16 10:19:21 -0700 (Wed, 16 May 2007)
New Revision: 6897
Added:
cs/cigma/trunk/tools/extract-cmd.c
Log:
Added extract command
Added: cs/cigma/trunk/tools/extract-cmd.c
===================================================================
--- cs/cigma/trunk/tools/extract-cmd.c 2007-05-16 17:19:07 UTC (rev 6896)
+++ cs/cigma/trunk/tools/extract-cmd.c 2007-05-16 17:19:21 UTC (rev 6897)
@@ -0,0 +1,264 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <argtable2.h>
+
+#include <cigma/qr.h>
+#include <cigma/mesh.h>
+#include <cigma/field.h>
+
+//#include <cigma/points.h>
+#include <cigma/array.h>
+#include <cigma/dataset.h>
+
+#include <cigma/det.h>
+#include <cigma/tet4.h>
+
+
+
+int extract_command(const char *rulefile, const char *rulepath,
+ const char *modelfile, const char *coorpath, const char *connpath, const char *fieldpath,
+ const char *outfile, const char *outpath)
+{
+ rule_t *qr;
+ mesh_t *mesh;
+ field_t *field;
+
+ //points_t *points;
+ array_t *points;
+ hid_t points_file;
+ hid_t points_dset;
+ herr_t status;
+
+ int e,q,d;
+ int nel,nq;
+ int nsd, ndim;
+ int estride, qstride;
+
+ double J[3*3];
+ double detJ;
+
+ double x[4*3];
+ double dofs[4*3];
+ double value[3];
+
+ double *qx;
+ double *qw;
+ double *data;
+
+ int ierr;
+ int exitcode = EXIT_SUCCESS;
+
+
+ //
+ // TODO: distinguish between possible errors in code below
+ //
+
+
+ /* Read quadrature info */
+ qr = (rule_t *)malloc(sizeof(rule_t));
+ ierr = qr_init(qr, rulefile, rulepath);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Error: could not access quadrature rules from '%s'\n", rulefile);
+ exitcode = EXIT_FAILURE;
+ goto exit_qr;
+ }
+
+
+ /* Read mesh info */
+ mesh = (mesh_t *)malloc(sizeof(mesh_t));
+ ierr = mesh_init(mesh, modelfile, coorpath, connpath);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Error: could not access mesh from '%s'\n", modelfile);
+ exitcode = EXIT_FAILURE;
+ goto exit_mesh;
+ }
+ assert(mesh->nsd == qr->nsd);
+
+
+ /* Read field info */
+ field = (field_t *)malloc(sizeof(field_t));
+ ierr = field_init(field, modelfile, fieldpath);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Error: could not access field %s from '%s'\n", fieldpath, modelfile);
+ exitcode = EXIT_FAILURE;
+ goto exit_field;
+ }
+ field->mesh = mesh; //TODO: fix explicit assignment
+
+
+ /* Create array for quadrature points */
+ //points = (points_t *)malloc(sizeof(points_t));
+ //points_init(points);
+ points = (array_t *)malloc(sizeof(array_t));
+ array_init(points, H5T_NATIVE_DOUBLE, 3, NULL, NULL);
+ points->shape[0] = mesh->nel;
+ points->shape[1] = qr->nq;
+ points->shape[2] = qr->nsd;
+ points->data = malloc(array_npoints(points) * sizeof(double));
+
+
+ /*
+ * Map quadrature points
+ */
+
+ nel = mesh->nel;
+ nq = qr->nq;
+ nsd = qr->nsd;
+ ndim = 3;
+
+ qx = qr->points;
+ qw = qr->weights;
+ data = (double *) points->data;
+
+ estride = ndim * nq;
+ qstride = ndim;
+
+ #ifdef DEBUG
+ printf("nno = %d\n", mesh->nno);
+ printf("nel = %d\n", nel);
+ printf("nq = %d\n", nq);
+ printf("nsd = %d\n", nsd);
+ printf("ndim = %d\n", ndim);
+ printf("npoints = %d\n", array_npoints(points));
+ #endif
+
+
+ for (e = 0; e < nel; e++)
+ {
+ mesh_coords(mesh, e, x);
+ field_dofs(field, e, dofs);
+
+ tet4_jacobian_3(x, J);
+ detJ = det3x3(J);
+
+ for (q = 0; q < nq; q++)
+ {
+ tet4_eval(dofs, ndim, &qx[3*q], value);
+ for (d = 0; d < ndim; d++)
+ {
+ data[e*estride + q*qstride + d] = value[d];
+ }
+ }
+ }
+
+
+ /*
+ * Write point data
+ */
+
+ points_file = H5Fopen(outfile, H5F_ACC_RDWR, H5P_DEFAULT);
+ points_dset = array_create(points, points_file, outpath);
+ if (points_dset < 0)
+ {
+ fprintf(stderr, "Error: could not create dataset %s on '%s'\n", outpath, outfile);
+ exitcode = EXIT_FAILURE;
+ goto exit_points;
+ }
+ ierr = array_slice_write(points, points_dset, 0, nel);
+ if (ierr < 0)
+ {
+ fprintf(stderr, "Error: could not write to '%s'\n", outfile);
+ exitcode = EXIT_FAILURE;
+ }
+ status = H5Fclose(points_file);
+
+
+ /*
+ * Cleanup in reverse so that labels can be effective.
+ */
+
+exit_points:
+ array_free(points);
+ free(points);
+
+exit_field:
+ field_free(field);
+ free(field);
+
+exit_mesh:
+ mesh_free(mesh);
+ free(mesh);
+
+exit_qr:
+ qr_free(qr);
+ free(qr);
+
+ return exitcode;
+}
+
+
+
+int extract_main(int argc, char *argv[])
+{
+ const char *progname = "extract";
+
+ struct arg_lit *help = arg_lit0 ("h", "help", "display this help and exit");
+ struct arg_file *rulefile = arg_file1("r", "rules", "<rules.h5>", "quadrature rules file");
+ struct arg_str *rulepath = arg_str1 (NULL, "rule", NULL, "path to quadrature rule");
+ struct arg_file *modelfile = arg_file1("m", "model", "<model.h5>", "mesh file");
+ struct arg_str *coorpath = arg_str1 (NULL, "coord", NULL, "path to coords dataset");
+ struct arg_str *connpath = arg_str1 (NULL, "conn", NULL, "path to connectivity dataset");
+ struct arg_str *fieldpath = arg_str1 (NULL, "field", NULL, "path to field dataset");
+ struct arg_file *outfile = arg_file1("o", NULL, "<output.h5>", "output file");
+ struct arg_str *outpath = arg_str1 (NULL, "path", NULL, "output path for points dataset");
+ struct arg_end *end = arg_end(20);
+
+ void *argtable[] = {
+ rulefile, rulepath,
+ modelfile, coorpath, connpath, fieldpath,
+ outfile, outpath,
+ help,
+ end
+ };
+
+ int nerrors;
+ int exitcode = 0;
+
+ /* verify the argtable entries were allocated successfully */
+ if (arg_nullcheck(argtable) != 0)
+ {
+ fprintf(stderr, "%s: failed to allocate argtable\n", progname);
+ exitcode = 1;
+ goto exit;
+ }
+
+ nerrors = arg_parse(argc, argv, argtable);
+
+ /* special case: '--help' takes precedence over error reporting */
+ if (help->count > 0)
+ {
+ printf("Usage: %s", progname);
+ arg_print_syntax(stdout, argtable, "\n\n");
+ arg_print_glossary(stdout, argtable, " %-25s %s\n");
+ exitcode = 0;
+ goto exit;
+ }
+
+ /* If the parser returned any errors then display them and exit */
+ if (nerrors > 0)
+ {
+ /* Display the error details contained in the arg_end struct. */
+ arg_print_errors(stdout, end, progname);
+ exitcode = 1;
+ goto exit;
+ }
+
+ /* normal case: take the command line options */
+ exitcode = extract_command(rulefile->filename[0],
+ rulepath->sval[0],
+ modelfile->filename[0],
+ coorpath->sval[0],
+ connpath->sval[0],
+ fieldpath->sval[0],
+ outfile->filename[0],
+ outpath->sval[0]);
+
+exit:
+ arg_freetable(argtable, sizeof(argtable)/sizeof(argtable[0]));
+
+ return exitcode;
+}
More information about the cig-commits
mailing list