[cig-commits] r6770 - in cs/cigma/trunk/sandbox/c: . tests
luis at geodynamics.org
luis at geodynamics.org
Thu May 3 16:14:54 PDT 2007
Author: luis
Date: 2007-05-03 16:14:54 -0700 (Thu, 03 May 2007)
New Revision: 6770
Modified:
cs/cigma/trunk/sandbox/c/Makefile
cs/cigma/trunk/sandbox/c/common.c
cs/cigma/trunk/sandbox/c/common.h
cs/cigma/trunk/sandbox/c/mesh.c
cs/cigma/trunk/sandbox/c/mesh.h
cs/cigma/trunk/sandbox/c/tests/Makefile
cs/cigma/trunk/sandbox/c/tests/test_mesh.c
cs/cigma/trunk/sandbox/c/tests/test_tet4.c
cs/cigma/trunk/sandbox/c/tet4.c
cs/cigma/trunk/sandbox/c/tet4.h
Log:
Updated sandbox/c directory
Modified: cs/cigma/trunk/sandbox/c/Makefile
===================================================================
--- cs/cigma/trunk/sandbox/c/Makefile 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/Makefile 2007-05-03 23:14:54 UTC (rev 6770)
@@ -14,7 +14,7 @@
qr.o \
common.o
-TARGETS = libcigma.a libcigma.so
+TARGETS = libcigma.a
all: $(TARGETS)
@@ -24,6 +24,9 @@
libcigma.so: $(OBJFILES)
$(LD) -shared $(LIBS) $^ -o $@
+main: main.c
+ gcc -Wall -g -DNDEBUG main.c -L. -lcigma -lhdf5 -o main
+
.c.o:
$(CC) $(CFLAGS) -c $<
Modified: cs/cigma/trunk/sandbox/c/common.c
===================================================================
--- cs/cigma/trunk/sandbox/c/common.c 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/common.c 2007-05-03 23:14:54 UTC (rev 6770)
@@ -104,4 +104,39 @@
return m[0]*c[0] - m[1]*c[1] + m[2]*c[2] - m[3]*c[3];
}
+double tet_volume(double a[3], double b[3], double c[3], double d[3])
+{
+ int j;
+ double V[4*4];
+ V[4*0 + 0] = 1.0;
+ V[4*1 + 0] = 1.0;
+ V[4*2 + 0] = 1.0;
+ V[4*3 + 0] = 1.0;
+
+ for (j = 1; j < 4; j++)
+ {
+ V[4*0 + j] = a[j-1];
+ V[4*1 + j] = b[j-1];
+ V[4*2 + j] = c[j-1];
+ V[4*3 + j] = d[j-1];
+ }
+ return det4x4(V)/6.0;
+}
+
+int tet_volumes(double vertices[4*3], double x[3], double volumes[4])
+{
+ double *a = &vertices[3*0];
+ double *b = &vertices[3*1];
+ double *c = &vertices[3*2];
+ double *d = &vertices[3*3];
+ volumes[0] = tet_volume(a,b,c,x);
+ volumes[1] = tet_volume(b,d,c,x);
+ volumes[2] = tet_volume(a,c,d,x);
+ volumes[3] = tet_volume(a,d,b,x);
+ if (volumes[0] < 0) return -1;
+ if (volumes[1] < 0) return -1;
+ if (volumes[2] < 0) return -1;
+ if (volumes[3] < 0) return -1;
+ return 0;
+}
Modified: cs/cigma/trunk/sandbox/c/common.h
===================================================================
--- cs/cigma/trunk/sandbox/c/common.h 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/common.h 2007-05-03 23:14:54 UTC (rev 6770)
@@ -4,4 +4,7 @@
double det3x3(double *m);
double det4x4(double *m);
+double tet_volume(double a[3], double b[3], double c[3], double d[3]);
+int tet_volumes(double vertices[4*3], double x[3], double volumes[4]);
+
#endif
Modified: cs/cigma/trunk/sandbox/c/mesh.c
===================================================================
--- cs/cigma/trunk/sandbox/c/mesh.c 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/mesh.c 2007-05-03 23:14:54 UTC (rev 6770)
@@ -9,7 +9,31 @@
static int mesh_coords_fini(coordinates_t *nodes);
static int mesh_connect_fini(connectivity_t *conn);
+/*
+ * Accessor functions
+ */
+//
+// Input : mesh object, element number e
+// Output: dof coordinates -- [ndof x nsd] matrix
+//
+void mesh_coords(mesh_t *m, int e, double *coords)
+{
+
+ double *coor = m->coords->coords;
+ int *conn = m->connect->elements;
+ int ndof = m->connect->n;
+
+ int i,j;
+ for (i = 0; i < ndof; i++)
+ {
+ j = conn[ndof*e + i];
+ coords[3*i + 0] = coor[3*j + 0];
+ coords[3*i + 1] = coor[3*j + 1];
+ coords[3*i + 2] = coor[3*j + 2];
+ }
+}
+
/*
* Initialize mesh object
Modified: cs/cigma/trunk/sandbox/c/mesh.h
===================================================================
--- cs/cigma/trunk/sandbox/c/mesh.h 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/mesh.h 2007-05-03 23:14:54 UTC (rev 6770)
@@ -1,8 +1,7 @@
#ifndef __CIGMA_MESH_H__
#define __CIGMA_MESH_H__
-typedef struct
-{
+typedef struct {
int n; // number of components per node (always 3)
int nno; // number of nodes
double *coords;
@@ -14,8 +13,7 @@
int *elements;
} connectivity_t;
-typedef struct
-{
+typedef struct {
coordinates_t *coords;
connectivity_t *connect;
} mesh_t;
@@ -28,4 +26,6 @@
int mesh_init(mesh_t *m, char *filename, char *path_coords, char *path_connect);
int mesh_fini(mesh_t *m);
+void mesh_coords(mesh_t *m, int e, double *coords);
+
#endif
Modified: cs/cigma/trunk/sandbox/c/tests/Makefile
===================================================================
--- cs/cigma/trunk/sandbox/c/tests/Makefile 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/tests/Makefile 2007-05-03 23:14:54 UTC (rev 6770)
@@ -10,7 +10,8 @@
TARGETS = test_mesh \
test_qr \
test_tet4 \
- test_common
+ test_common \
+ test_volume
all: $(TARGETS)
@@ -26,6 +27,9 @@
test_common: test_common.c
$(CC) $(CFLAGS) $< -o $@ $(LIBS)
+test_volume: test_volume.c
+ $(CC) $(CFLAGS) $< -o $@ $(LIBS)
+
clean:
rm -f $(TARGETS)
Modified: cs/cigma/trunk/sandbox/c/tests/test_mesh.c
===================================================================
--- cs/cigma/trunk/sandbox/c/tests/test_mesh.c 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/tests/test_mesh.c 2007-05-03 23:14:54 UTC (rev 6770)
@@ -5,6 +5,12 @@
int main(int argc, char *argv[])
{
mesh_t mesh;
+ int ndof;
+ int i,j;
+
+ int e;
+ double *coords;
+
int ierr;
if (argc < 2)
@@ -23,6 +29,29 @@
printf("Reading %s\n", argv[1]);
printf("Found %d nodes and %d elements\n", mesh.coords->nno, mesh.connect->nel);
+ ndof = mesh.connect->n;
+ coords = (double *)malloc((ndof * 3) * sizeof(double));
+
+ e = 100;
+ mesh_coords(&mesh, e, coords);
+ printf("coords(%d) = {\n", e);
+ for (i = 0; i < ndof; i++)
+ {
+ printf("\t{");
+ for (j = 0; j < 3; j++)
+ {
+ printf("%g", coords[3*i + j]);
+ if (j != 2)
+ printf(",");
+ }
+ printf("}");
+ if (i != ndof-1)
+ printf(",");
+ printf("\n");
+ }
+ printf("};\n");
+ free(coords);
+
mesh_fini(&mesh);
return EXIT_SUCCESS;
}
Modified: cs/cigma/trunk/sandbox/c/tests/test_tet4.c
===================================================================
--- cs/cigma/trunk/sandbox/c/tests/test_tet4.c 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/tests/test_tet4.c 2007-05-03 23:14:54 UTC (rev 6770)
@@ -1,68 +1,18 @@
#include <stdio.h>
#include <stdlib.h>
#include "tet4.h"
+#include "test_data.h"
-#define NTET 21
-double tets[NTET][12] = {
- {-1, -1, -1, +1, -1, -1, -1, +1, -1, -1, -1, +1},
- {11750,250,-22250,11750,0,-22500,11500,0,-22500,11500,250,-22500},
- {1750,16000,-7250,1500,15750,-7500,1500,16000,-7250,1500,16000,-7500},
- {1000,1250,-24000,750,1250,-23750,750,1000,-23750,1000,1000,-23750},
- {6250,250,-23500,6250,0,-23250,6000,250,-23250,6000,0,-23250},
- {5750,500,-22750,5750,250,-22500,5500,500,-22500,5500,250,-22500},
- {4500,500,-24000,4500,250,-23750,4750,250,-24000,4750,500,-24000},
- {2000,8250,-8000,2000,8250,-8250,2000,8000,-8000,1750,8250,-8000},
- {10750,14500,-13000,11000,14750,-13250,11000,14500,-13000,11000,14500,-13250},
- {22250,16500,-13750,22250,16500,-13500,22000,16500,-13500,22250,16250,-13500},
- {7500,10500,-14750,7750,10500,-14750,7750,10250,-14750,7750,10500,-15000},
- {7500,3750,-15750,7500,4000,-15750,7500,4000,-15500,7750,4000,-15500},
- {500,750,-23000,500,500,-22750,250,750,-22750,250,500,-22750},
- {500,1750,-23750,750,1500,-23750,750,1750,-23500,500,1500,-23750},
- {16000,750,-24000,16250,750,-23750,16000,500,-24000,16250,500,-24000},
- {1500,3000,-24000,1500,3250,-24000,1750,3250,-23750,1750,3000,-24000},
- {14750,500,-23250,15000,250,-23250,14750,250,-23250,15000,500,-23500},
- {11000,18000,-4500,10750,17750,-4750,10750,18000,-4500,10750,18000,-4750},
- {250,0,-24000,250,250,-24000,250,0,-23750,500,0,-24000},
- {20000,24000,0,20250,24000,0,20250,23750,0,20250,24000,-250},
- {20250,24000,0,20500,24000,0,20250,23750,0,20250,24000,-250}
-};
+void display_points(double *qpts, double nq)
+{
+ int i;
+ for (i = 0; i < nq; i++)
+ {
+ printf("\tx[%2d] = (%lf, %lf, %lf)\n",
+ i, qpts[3*i], qpts[3*i+1], qpts[3*i+2]);
+ }
+}
-#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
-};
-
-void display_points(double *qpts, double nq);
-
int main(int argc, char *argv[])
{
int n;
@@ -72,6 +22,9 @@
double qpts[NQ*3];
double *verts;
+ double refzero[3] = {0.0, 0.0, 0.0};
+ double zero[3];
+
for (n = 0; n < NTET; n++)
{
verts = (double *) &tets[n];
@@ -83,8 +36,10 @@
// extract quadrature points
tet4_tabulate(cools_points, NQ, tab);
tet4_batch_eval3(verts, tab, NQ, qpts);
+ tet4_eval3(verts, refzero, zero);
printf("tet[%d]\n", n);
+ printf("\t zero = (%lf, %lf, %lf)\n", zero[0], zero[1], zero[2]);
display_points(qpts, NQ);
printf("\n");
}
@@ -93,12 +48,3 @@
}
-void display_points(double *qpts, double nq)
-{
- int i;
- for (i = 0; i < nq; i++)
- {
- printf("\tx[%2d] = (%lf, %lf, %lf)\n",
- i, qpts[3*i], qpts[3*i+1], qpts[3*i+2]);
- }
-}
Modified: cs/cigma/trunk/sandbox/c/tet4.c
===================================================================
--- cs/cigma/trunk/sandbox/c/tet4.c 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/tet4.c 2007-05-03 23:14:54 UTC (rev 6770)
@@ -41,6 +41,7 @@
static tet4_basis_t dN[3][4] = {{TN00,TN10,TN20,TN30},
{TN01,TN11,TN21,TN31},
{TN02,TN12,TN22,TN32}};
+
//
// Represent derivatives as an array of constants
// althought technically, it should be an array of
@@ -71,8 +72,34 @@
0.0, 0.0, 0.5
};
-void tet4_jacobian(tet4_t *t, double xi[3], double *J)
+
+// ---------------------------------------------------------------------------
+
+
+void tet4_set1(tet4_t *t, double x[4*3])
{
+ int i;
+ for (i = 0; i < 4*3; i++)
+ {
+ t->x[i] = x[i];
+ }
+}
+
+void tet4_set4(tet4_t *t, double a[3], double b[3], double c[3], double d[3])
+{
+ int i;
+ for (i = 0; i < 3; i++)
+ {
+ t->x[3*0 + i] = a[i];
+ t->x[3*1 + i] = b[i];
+ t->x[3*2 + i] = c[i];
+ t->x[3*3 + i] = d[i];
+ }
+}
+
+
+void tet4_jacobian(tet4_t *t, double xi[3], double J[3*3])
+{
//
// J[i,j] = \frac{\partial{x_i}}{\partial{\xi_j}}
//
@@ -93,7 +120,7 @@
}
}
}
-void tet4_jacobian_2(tet4_t *t, double *J)
+void tet4_jacobian_2(tet4_t *t, double J[3*3])
{
//
// [ dN[0]_xi dN[1]_xi dN[2]_xi dN[3]_xi ] [ x[0][0] x[0][1] x[0][2] ]
@@ -118,7 +145,7 @@
J[3*2+2] = dTN[3*0+2]*x[3*0+2] + dTN[3*1+2]*x[3*1+2] + dTN[3*2+2]*x[3*2+2] + dTN[3*3+2]*x[3*3+2];
}
-void tet4_jacobian_3(tet4_t *t, double *J)
+void tet4_jacobian_3(tet4_t *t, double J[3*3])
{
/* XXX: from notes.. */
J[3*0+0] = 0;
@@ -135,27 +162,6 @@
}
-void tet4_set1(tet4_t *t, double *x)
-{
- int i;
- for (i = 0; i < 4*3; i++)
- {
- t->x[i] = x[i];
- }
-}
-
-void tet4_set4(tet4_t *t, double *a, double *b, double *c, double *d)
-{
- int i;
- for (i = 0; i < 3; i++)
- {
- t->x[3*0 + i] = a[i];
- t->x[3*1 + i] = b[i];
- t->x[3*2 + i] = c[i];
- t->x[3*3 + i] = d[i];
- }
-}
-
void tet4_eval(double *dof, int dims, double *xi, double *val)
{
int i, j;
@@ -172,7 +178,7 @@
}
}
-void tet4_eval1(double dof[4], double xi[3], double *val)
+void tet4_eval1(double dof[4*1], double xi[3], double val[1])
{
*val = dof[0]*TN0(xi) + dof[1]*TN1(xi) + dof[2]*TN2(xi) + dof[3]*TN3(xi);
}
@@ -197,7 +203,7 @@
}
}
-void tet4_eval1_2(double dof[4], double xi[3], double *val)
+void tet4_eval1_2(double dof[4*1], double xi[3], double val[1])
{
int j;
*val = 0.0;
@@ -269,7 +275,7 @@
}
}
}
-void tet4_batch_eval1(double dof[4], double *tab, int npts, double *vals)
+void tet4_batch_eval1(double dof[4*1], double *tab, int npts, double *vals)
{
int i;
for (i = 0; i < npts; i++)
Modified: cs/cigma/trunk/sandbox/c/tet4.h
===================================================================
--- cs/cigma/trunk/sandbox/c/tet4.h 2007-05-03 23:13:10 UTC (rev 6769)
+++ cs/cigma/trunk/sandbox/c/tet4.h 2007-05-03 23:14:54 UTC (rev 6770)
@@ -8,24 +8,24 @@
double x[4*3]; // [4 x 3] matrix - coordinates of degrees of freedom
} tet4_t;
-void tet4_set1(tet4_t *t, double *x);
-void tet4_set4(tet4_t *t, double *a, double *b, double *c, double *d);
+void tet4_set1(tet4_t *t, double x[4*3]);
+void tet4_set4(tet4_t *t, double a[3], double b[3], double c[3], double d[3]);
-void tet4_jacobian(tet4_t *t, double xi[3], double *J);
-void tet4_jacobian_2(tet4_t *t, double *J);
-void tet4_jacobian_3(tet4_t *t, double *J);
+void tet4_jacobian(tet4_t *t, double xi[3], double J[3*3]);
+void tet4_jacobian_2(tet4_t *t, double J[3*3]);
+void tet4_jacobian_3(tet4_t *t, double J[3*3]);
void tet4_eval(double *dof, int dims, double *xi, double *val);
-void tet4_eval1(double dof[4], double xi[3], double *val);
+void tet4_eval1(double dof[4*1], double xi[3], double val[1]);
void tet4_eval3(double dof[4*3], double xi[3], double val[3]);
void tet4_eval6(double dof[4*6], double xi[3], double val[6]);
-void tet4_eval1_2(double dof[4], double xi[3], double *val);
+void tet4_eval1_2(double dof[4*1], double xi[3], double val[1]);
void tet4_eval3_2(double dof[4*3], double xi[3], double val[3]);
void tet4_eval6_2(double dof[4*6], double xi[3], double val[6]);
-void tet4_tabulate(double *xs, int npts, double *ys);
+void tet4_tabulate(double *xs, int npts, double *tab);
void tet4_batch_eval(double *dof, int dims, double *tab, int npts, double *vals);
-void tet4_batch_eval1(double dof[4], double *tab, int npts, double *vals);
+void tet4_batch_eval1(double dof[4*1], double *tab, int npts, double *vals);
void tet4_batch_eval3(double dof[4*3], double *tab, int npts, double *vals);
#endif
More information about the cig-commits
mailing list