[cig-commits] r6771 - cs/cigma/trunk/sandbox/c
luis at geodynamics.org
luis at geodynamics.org
Thu May 3 16:15:14 PDT 2007
Author: luis
Date: 2007-05-03 16:15:14 -0700 (Thu, 03 May 2007)
New Revision: 6771
Added:
cs/cigma/trunk/sandbox/c/main.c
Log:
Example main loop
Added: cs/cigma/trunk/sandbox/c/main.c
===================================================================
--- cs/cigma/trunk/sandbox/c/main.c 2007-05-03 23:14:54 UTC (rev 6770)
+++ cs/cigma/trunk/sandbox/c/main.c 2007-05-03 23:15:14 UTC (rev 6771)
@@ -0,0 +1,370 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+
+#include "hdf5.h"
+#include "tet4.h"
+#include "common.h"
+#include "mesh.h"
+
+#define NQ 14
+double cools_weights[NQ] = {
+ 0.025396825396825397,
+ 0.025396825396825397,
+ 0.025396825396825397,
+ 0.025396825396825397,
+ 0.025396825396825397,
+ 0.025396825396825397,
+ 0.11811976632397428,
+ 0.11811976632397428,
+ 0.11811976632397428,
+ 0.11811976632397428,
+ 0.17711832891412096,
+ 0.17711832891412096,
+ 0.17711832891412096,
+ 0.17711832891412096,
+};
+
+double cools_points[NQ*3] = {
+ -1.0, 0.0, 0.0,
+ 0.0, 0.0, -1.0,
+ -1.0, 0.0, -1.0,
+ -1.0, -1.0, 0.0,
+ 0.0, -1.0, -1.0,
+ 0.0, -1.0, 0.0,
+ -0.79894646954959103, -0.79894646954959103, -0.79894646954959103,
+ 0.39683940864877321, -0.79894646954959103, -0.79894646954959103,
+ -0.79894646954959103, 0.39683940864877321, -0.79894646954959103,
+ -0.79894646954959103, -0.79894646954959103, 0.39683940864877321,
+ -0.37125425301361559, -0.37125425301361559, -0.88623724095915324,
+ -0.37125425301361559, -0.37125425301361559, -0.37125425301361559,
+ -0.88623724095915324, -0.37125425301361559, -0.37125425301361559,
+ -0.37125425301361559, -0.88623724095915324, -0.37125425301361559
+};
+
+typedef struct
+{
+ int nno;
+ int nsd;
+ int nel;
+ int ndof;
+ int ndim;
+ mesh_t *mesh;
+ tet4_t t;
+ double *dofs; // [nno x ndim] matrix
+} field_t;
+
+void field_init(field_t *field, char *filename, char *path)
+{
+ hid_t file_id;
+ hid_t dataset;
+ hid_t dataspace;
+
+ int rank;
+ hsize_t dims[2];
+ herr_t status;
+
+ /* open the file in read-only mode */
+ file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+
+ /* read dataset */
+ dataset = H5Dopen(file_id, path);
+
+ dataspace = H5Dget_space(dataset);
+ rank = H5Sget_simple_extent_ndims(dataspace); assert(rank == 2);
+ status = H5Sget_simple_extent_dims(dataspace, dims, NULL);
+ status = H5Sclose(dataspace);
+
+ assert(field->nno == (int)dims[0]);
+ assert(field->ndim == (int)dims[1]);
+
+ field->dofs = (double *)malloc((field->nno * field->ndim) * sizeof(double));
+
+ status = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, field->dofs);
+
+ status = H5Dclose(dataset);
+ status = H5Fclose(file_id);
+}
+
+void field_fini(field_t *field)
+{
+ free(field->dofs);
+}
+
+void field_dofs(field_t *field, int e, double *dof)
+{
+ int i;
+ int n;
+ int *elt = &(field->mesh->connect->elements[4*e]);
+ for (i = 0; i < 4; i++)
+ {
+ n = elt[i];
+ dof[3*i + 0] = field->dofs[3*n + 0];
+ dof[3*i + 1] = field->dofs[3*n + 1];
+ dof[3*i + 2] = field->dofs[3*n + 2];
+ }
+}
+
+void field_eval(field_t *B, double *x, int nq, double *g, int dims)
+{
+ int q;
+ int e;
+ int nel = B->mesh->connect->nel;
+ int ierr;
+
+ double vols[4];
+ double verts[4*3];
+ double dof[4*3];
+
+ for (q = 0; q < nq; q++)
+ {
+ for (e = 0; e < nel; e++)
+ {
+ mesh_coords(B->mesh, e, verts);
+ ierr = tet_volumes(verts, &x[3*q], vols);
+ if (ierr == 0)
+ {
+ // found the element...evaluate and break element loop
+ field_dofs(B, e, dof);
+ tet4_eval3(dof, &x[3*q], &g[dims * q]);
+ break;
+ }
+ }
+ }
+}
+
+double main_loop(field_t *A, field_t *B, double *errs)
+{
+ int e; // loop counter for elements
+ int q; // loop counter for quadrature points
+ int i,n;
+
+ int nno;
+ int nel;
+ int nq;
+
+ int dims;
+ int ndof;
+
+ //fe_space_t *fe;
+
+ mesh_t *mesh;
+ double *coords; // coordinates of dofs
+
+ double dof[4*3];
+
+ double *x;
+ double *xi;
+ double *w;
+
+ double *tab;
+
+ double J[3*3];
+ double *detJ;
+
+ double *f; // from primary mesh
+ double *g;
+
+ double global_acc;
+ double local_acc;
+ double delta;
+
+ //
+ // Set up environment
+ //
+
+ //fe = A->fe;
+ mesh = A->mesh;
+
+ //dims = fe->ndim;
+ //ndof = fe->ndof;
+ dims = A->ndim;
+ ndof = A->ndof;
+
+ //xi = qrule->points;
+ //w = qrule->weights;
+ //nq = qrule->nq;
+ xi = cools_points;
+ w = cools_weights;
+ nq = NQ;
+
+ nno = mesh->coords->nno;
+ nel = mesh->connect->nel;
+
+ coords = (double *)malloc((ndof * 3) * sizeof(double));
+
+ x = (double *)malloc((nq * 3) * sizeof(double));
+
+ tab = (double *)malloc((nq * ndof) * sizeof(double));
+
+ detJ = (double *)malloc(nq * sizeof(double *));
+
+ f = (double *)malloc((nq * dims) * sizeof(double));
+ g = (double *)malloc((nq * dims) * sizeof(double));
+
+
+ //
+ // Main loop
+ //
+ global_acc = 0.0;
+ for (e = 0; e < nel; e++)
+ {
+ mesh_coords(mesh, e, coords);
+ tet4_set1(&(A->t), coords);
+
+ //fe_tabulate(fe, xi, nq, tab);
+ //fe_batch_eval(fe, coords, tab, nq, x);
+ //fe_batch_eval(fe, field->dofs[e], tab, nq, f);
+
+ tet4_tabulate(xi, nq, tab);
+
+ tet4_batch_eval3(coords, tab, nq, f);
+
+ field_dofs(A, e, dof);
+ tet4_batch_eval3(dof, tab, nq, f);
+
+ field_eval(B, x, nq, g, dims);
+
+ for (q = 0; q < nq; q++)
+ {
+ //fe_jacobian(fe, xi, J);
+ tet4_jacobian(&(A->t), &xi[3*q], J);
+ detJ[q] = det3x3(J);
+ }
+
+ // apply integration rule
+ local_acc = 0.0;
+ for (q = 0; q < nq; q++)
+ {
+ for (n = 0; n < dims; n++)
+ {
+ i = dims*q + n;
+ delta = f[i] - g[i];
+ local_acc += w[q] * detJ[q] * delta * delta;
+ }
+ }
+
+ // accumulate local contribution
+ global_acc += local_acc;
+
+ // report!
+ if (e % 100 == 0)
+ printf("e=%d\t(%g, %g)\n", e, local_acc, global_acc);
+ }
+ global_acc = sqrt(global_acc);
+
+ //
+ // Clean up
+ //
+ free(f);
+ free(g);
+ free(x);
+ free(tab);
+ free(coords);
+ free(detJ);
+
+ return global_acc;
+}
+
+
+void save_vtk(double *coor, int nno,
+ int *conn, int nel,
+ double *err)
+{
+ int n;
+ int e;
+ FILE *fp = fopen("local.vtk", "w");
+
+ fprintf(fp, "# vtk DataFile Version 3.0\n");
+ fprintf(fp, "This line is a comment\n");
+ fprintf(fp, "ASCII\n");
+ fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
+
+ fprintf(fp, "POINTS %d float\n", nno);
+ for (n = 0; n < nno; n++)
+ fprintf(fp, "%g %g %g\n", coor[3*n], coor[3*n+1], coor[3*n+2]);
+
+ fprintf(fp, "CELLS %d %d\n", nel, 4*nel);
+ for (e = 0; e < nel; e++)
+ fprintf(fp, "4 %d %d %d %d\n", conn[4*e], conn[4*e+1], conn[4*e+2], conn[4*e+3]);
+
+ fprintf(fp, "CELL_TYPES %d\n", nel);
+ for (e = 0; e < nel; e++)
+ fprintf(fp, "10\n");
+
+ fprintf(fp, "CELL_DATA %d\n", nel);
+ fprintf(fp, "SCALARS local_error float 1\n");
+ fprintf(fp, "LOOKUP_TABLE default\n");
+ for (e = 0; e < nel; e++)
+ fprintf(fp, "%g\n", err[e]);
+
+ fclose(fp);
+}
+
+
+int main(int argc, char *argv[])
+{
+ char fileA[] = "/home/luis/benchmark/bmssnog-250m-lintet.h5";
+ char fileB[] = "/home/luis/benchmark/bmssnog-1000m-lintet.h5";
+ mesh_t meshA;
+ mesh_t meshB;
+
+ field_t fieldA;
+ field_t fieldB;
+
+ double *local_err;
+
+ mesh_init(&meshA, fileA, "/model/geometry/coordinates", "/model/topology/simplices");
+ printf("Reading meshA\n");
+
+ mesh_init(&meshB, fileB, "/model/geometry/coordinates", "/model/topology/simplices");
+ printf("Reading meshB\n");
+
+ fieldA.nsd = 3;
+ fieldB.nsd = 3;
+
+ fieldA.nno = meshA.coords->nno;
+ fieldB.nno = meshB.coords->nno;
+
+ fieldA.nel = meshA.connect->nel;
+ fieldB.nel = meshA.connect->nel;
+
+ fieldA.ndof = meshA.connect->n;
+ fieldB.ndof = meshB.connect->n;
+
+ fieldA.ndim = 3;
+ fieldB.ndim = 3;
+
+ fieldA.mesh = &meshA;
+ fieldB.mesh = &meshB;
+
+ field_init(&fieldA, fileA, "/model/solution/snapshot0/displacement");
+ field_init(&fieldB, fileB, "/model/solution/snapshot0/displacement");
+
+ local_err = (double *)malloc(fieldA.nel * sizeof(double));
+
+
+
+ printf("Entering main_loop\n");
+ main_loop(&fieldA, &fieldB, local_err);
+
+
+ printf("Saving local errors to file\n");
+ save_vtk(fieldA.mesh->coords->coords, fieldA.nno,
+ fieldA.mesh->connect->elements, fieldA.nel,
+ local_err);
+
+
+
+ free(local_err);
+
+ field_fini(&fieldA);
+ field_fini(&fieldB);
+
+ mesh_fini(&meshA);
+ mesh_fini(&meshB);
+
+ return EXIT_SUCCESS;
+}
More information about the cig-commits
mailing list