[cig-commits] r6828 - in cs/cigma/trunk: include/cigma src tests
luis at geodynamics.org
luis at geodynamics.org
Wed May 9 17:26:22 PDT 2007
Author: luis
Date: 2007-05-09 17:26:22 -0700 (Wed, 09 May 2007)
New Revision: 6828
Added:
cs/cigma/trunk/include/cigma/tet4.h
cs/cigma/trunk/src/tet4.c
cs/cigma/trunk/tests/test_tet4.c
Log:
Default finite element space for tetrahedral elements
Added: cs/cigma/trunk/include/cigma/tet4.h
===================================================================
--- cs/cigma/trunk/include/cigma/tet4.h 2007-05-10 00:24:12 UTC (rev 6827)
+++ cs/cigma/trunk/include/cigma/tet4.h 2007-05-10 00:26:22 UTC (rev 6828)
@@ -0,0 +1,40 @@
+#ifndef __CIGMA_TET4_H__
+#define __CIGMA_TET4_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+extern const int tet4_nsd; // 3
+extern const int tet4_ndof; // 4
+
+void tet4_set(double x[4*3], double a[3], double b[3], double c[3], double d[3]);
+void tet4_get(double x[4*3], double a[3], double b[3], double c[3], double d[3]);
+
+void tet4_jacobian(double x[4*3], double xi[3], double J[3*3]);
+void tet4_jacobian_1(double x[4*3], double xi[3], double J[3*3]);
+void tet4_jacobian_2(double x[4*3], double J[3*3]);
+void tet4_jacobian_3(double x[4*3], double J[3*3]);
+
+void tet4_eval(double *dof, int dims, double *xi, 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*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 *tab);
+void tet4_batch_eval(double *dof, int dims, 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);
+
+double tet4_volume(double a[3], double b[3], double c[3], double d[3]);
+int tet4_test(double x[4*3], double pt[3]);
+int tet4_test2(double x[4*3], double pt[3], double vols[3]);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif
Added: cs/cigma/trunk/src/tet4.c
===================================================================
--- cs/cigma/trunk/src/tet4.c 2007-05-10 00:24:12 UTC (rev 6827)
+++ cs/cigma/trunk/src/tet4.c 2007-05-10 00:26:22 UTC (rev 6828)
@@ -0,0 +1,411 @@
+#include <stdlib.h>
+#include <assert.h>
+#include <cigma/det.h>
+#include <cigma/tet4.h>
+
+// Constants for this finite element space
+const int tet4_nsd = 3;
+const int tet4_ndof = 4;
+
+static double TN0(double xi[3]) {return 0.5*(-1.0 - xi[0] - xi[1] - xi[2]);}
+static double TN1(double xi[3]) {return 0.5*( 1.0 + xi[0] );}
+static double TN2(double xi[3]) {return 0.5*( 1.0 + xi[1] );}
+static double TN3(double xi[3]) {return 0.5*( 1.0 + xi[2] );}
+
+static double TN00(double xi[3]) {return 0.5*(-1.0);}
+static double TN10(double xi[3]) {return 0.5*(+1.0);}
+static double TN20(double xi[3]) {return 0.5*( 0.0);}
+static double TN30(double xi[3]) {return 0.5*( 0.0);}
+
+static double TN01(double xi[3]) {return 0.5*(-1.0);}
+static double TN11(double xi[3]) {return 0.5*( 0.0);}
+static double TN21(double xi[3]) {return 0.5*(+1.0);}
+static double TN31(double xi[3]) {return 0.5*( 0.0);}
+
+static double TN02(double xi[3]) {return 0.5*(-1.0);}
+static double TN12(double xi[3]) {return 0.5*( 0.0);}
+static double TN22(double xi[3]) {return 0.5*( 0.0);}
+static double TN32(double xi[3]) {return 0.5*(+1.0);}
+
+//
+// note that under the basis (1,xi,eta,zeta)
+// these functions are represented by:
+//
+// TN[0] = (-0.5,-0.5,-0.5,-0.5)
+// TN[1] = (+0.5, 0.5, 0.0, 0.0)
+// TN[2] = (+0.5, 0.0, 0.5, 0.0)
+// TN[3] = (+0.5, 0.0, 0.0, 0.5)
+//
+typedef double (*tet4_shape_fn)(double xi[3]);
+
+static tet4_shape_fn TN[4] = {TN0,TN1,TN2,TN3};
+
+static tet4_shape_fn dTN[4][3] = {
+ {TN00,TN01,TN02},
+ {TN10,TN11,TN12},
+ {TN20,TN21,TN22},
+ {TN30,TN31,TN32}
+};
+
+
+//
+// Represent derivatives as an array of constants
+// althought technically, it should be an array of
+// polynomials
+//
+// dTN[i,j] = \frac{\partial{TN_i}}{\partial{\xi_j}}
+//
+// TN[0]_xi = (-0.5, 0, 0, 0)
+// TN[0]_eta = (-0.5, 0, 0, 0)
+// TN[0]_zeta = (-0.5, 0, 0, 0)
+//
+// TN[1]_xi = (+0.5, 0, 0, 0)
+// TN[1]_eta = ( 0.0, 0, 0, 0)
+// TN[1]_zeta = ( 0.0, 0, 0, 0)
+//
+// TN[2]_xi = ( 0.0, 0, 0, 0)
+// TN[2]_eta = (+0.5, 0, 0, 0)
+// TN[2]_zeta = ( 0.0, 0, 0, 0)
+//
+// TN[3]_xi = ( 0.0, 0, 0, 0)
+// TN[3]_eta = ( 0.0, 0, 0, 0)
+// TN[3]_zeta = (+0.5, 0, 0, 0)
+//
+static double DTN[4*3] = {
+ -0.5, -0.5, -0.5,
+ 0.5, 0.0, 0.0,
+ 0.0, 0.5, 0.0,
+ 0.0, 0.0, 0.5
+};
+
+
+// ---------------------------------------------------------------------------
+
+
+void tet4_set(double x[4*3], double a[3], double b[3], double c[3], double d[3])
+{
+ int j;
+ for (j = 0; j < 3; j++)
+ {
+ x[3*0 + j] = a[j];
+ x[3*1 + j] = b[j];
+ x[3*2 + j] = c[j];
+ x[3*3 + j] = d[j];
+ }
+}
+
+void tet4_get(double x[4*3], double a[3], double b[3], double c[3], double d[3])
+{
+ int j;
+ for (j = 0; j < 3; j++)
+ {
+ a[j] = x[3*0 + j];
+ b[j] = x[3*1 + j];
+ c[j] = x[3*2 + j];
+ d[j] = x[3*3 + j];
+ }
+}
+
+// ---------------------------------------------------------------------------
+
+void tet4_jacobian(double x[4*3], double xi[3], double J[3*3])
+{
+ //
+ // J[i,j] = \frac{\partial{x_i}}{\partial{\xi_j}}
+ //
+ // where x_i = \sum_k (TN_k(\xi) x_{i,k})
+ //
+ // J[i,j] = \sum_k \frac{\partial{TN_k}}{\partial{\xi_j}} x_{i,k}
+ //
+ //
+ int i,j,k;
+ for (i = 0; i < 3; i++)
+ {
+ for (j = 0; j < 3; j++)
+ {
+ J[3*i + j] = 0.0;
+ for (k = 0; k < 4; k++)
+ J[3*i + j] += DTN[3*k + j] * x[3*i + k];
+ }
+ }
+}
+
+void tet4_jacobian_1(double x[4*3], double xi[3], double J[3*3])
+{
+ int i,j,k;
+ for (i = 0; i < 3; i++)
+ {
+ for (j = 0; j < 3; j++)
+ {
+ J[3*i + j] = 0.0;
+ for (k = 0; k < 4; k++)
+ J[3*i + j] += dTN[k][j](xi) * x[3*i + k];
+ }
+ }
+}
+
+void tet4_jacobian_2(double x[4*3], 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] ]
+ // J = [ dN[0]_eta dN[1]_eta dN[2]_eta dN[3]_eta ] [ x[1][0] x[1][1] x[1][2] ]
+ // [ dN[0]_zeta dN[1]_zeta dN[2]_zeta dN[3]_zeta ] [ x[2][0] x[2][1] x[2][2] ]
+ // [ x[3][0] x[3][1] x[3][2] ]
+ //
+ // Component-wise: J[i,j] = \sum_k dTN[k,i] x[k,j]
+ //
+ J[3*0+0] = DTN[3*0+0]*x[3*0+0] + DTN[3*1+0]*x[3*1+0] + DTN[3*2+0]*x[3*2+0] + DTN[3*3+0]*x[3*3+0];
+ J[3*0+1] = DTN[3*0+0]*x[3*0+1] + DTN[3*1+0]*x[3*1+1] + DTN[3*2+0]*x[3*2+1] + DTN[3*3+0]*x[3*3+1];
+ J[3*0+2] = DTN[3*0+0]*x[3*0+2] + DTN[3*1+0]*x[3*1+2] + DTN[3*2+0]*x[3*2+2] + DTN[3*3+0]*x[3*3+2];
+
+ J[3*1+0] = DTN[3*0+1]*x[3*0+0] + DTN[3*1+1]*x[3*1+0] + DTN[3*2+1]*x[3*2+0] + DTN[3*3+1]*x[3*3+0];
+ J[3*1+1] = DTN[3*0+1]*x[3*0+1] + DTN[3*1+1]*x[3*1+1] + DTN[3*2+1]*x[3*2+1] + DTN[3*3+1]*x[3*3+1];
+ J[3*1+2] = DTN[3*0+1]*x[3*0+2] + DTN[3*1+1]*x[3*1+2] + DTN[3*2+1]*x[3*2+2] + DTN[3*3+1]*x[3*3+2];
+
+ J[3*2+0] = DTN[3*0+2]*x[3*0+0] + DTN[3*1+2]*x[3*1+0] + DTN[3*2+2]*x[3*2+0] + DTN[3*3+2]*x[3*3+0];
+ J[3*2+1] = DTN[3*0+2]*x[3*0+1] + DTN[3*1+2]*x[3*1+1] + DTN[3*2+2]*x[3*2+1] + DTN[3*3+2]*x[3*3+1];
+ 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(double x[4*3], double J[3*3])
+{
+ //
+ // x = [[x0,y0,z0], //0,1,2
+ // [x1,y1,z1], //3,4,5
+ // [x2,y2,z2], //6,7,8
+ // [x3,y3,z3]] //9,10,11
+ //
+ // J = [[(x1-x0)/2, (x2-x0)/2, (x3-x0)/2], //0,1,2
+ // [(y1-x0)/2, (y2-x0)/2, (y3-x0)/2], //3,4,5
+ // [(z1-x0)/2, (z2-x0)/2, (z3-x0)/2]] //6,7,8
+ //
+
+ // first row
+ J[0] = 0.5 * (x[3] - x[0]);
+ J[1] = 0.5 * (x[6] - x[0]);
+ J[2] = 0.5 * (x[9] - x[0]);
+
+ // second row
+ J[3] = 0.5 * (x[4] - x[1]);
+ J[4] = 0.5 * (x[7] - x[1]);
+ J[5] = 0.5 * (x[10] - x[1]);
+
+ // third row
+ J[6] = 0.5 * (x[5] - x[2]);
+ J[7] = 0.5 * (x[8] - x[2]);
+ J[8] = 0.5 * (x[11] - x[2]);
+}
+
+// ---------------------------------------------------------------------------
+
+void tet4_eval(double *dof, int dims, double *xi, double *val)
+{
+ int i, j;
+ for (i = 0; i < dims; i++)
+ {
+ //
+ // val_i(xi) = \sum_j (dof[i,j] * TN_j(xi)
+ //
+ val[i] = 0.0;
+ for (j = 0; j < 4; j++)
+ {
+ val[i] += dof[4*i + j] * TN[j](xi);
+ }
+ }
+}
+
+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);
+}
+
+void tet4_eval3(double dof[4*3], double xi[3], double val[3])
+{
+ val[0] = dof[4*0+0]*TN0(xi) + dof[4*0+1]*TN1(xi) + dof[4*0+2]*TN2(xi) + dof[4*0+3]*TN3(xi);
+ val[1] = dof[4*1+0]*TN0(xi) + dof[4*1+1]*TN1(xi) + dof[4*1+2]*TN2(xi) + dof[4*1+3]*TN3(xi);
+ val[2] = dof[4*2+0]*TN0(xi) + dof[4*2+1]*TN1(xi) + dof[4*2+2]*TN2(xi) + dof[4*2+3]*TN3(xi);
+}
+
+void tet4_eval6(double dof[4*6], double xi[3], double val[6])
+{
+ int i;
+ for (i = 0; i < 6; i++)
+ {
+ val[i] = 0.0;
+ val[i] += dof[4*i+0] * TN0(xi);
+ val[i] += dof[4*i+1] * TN1(xi);
+ val[i] += dof[4*i+2] * TN2(xi);
+ val[i] += dof[4*i+3] * TN3(xi);
+ }
+}
+
+void tet4_eval1_2(double dof[4*1], double xi[3], double val[1])
+{
+ int j;
+ *val = 0.0;
+ for (j = 0; j < 4; j++)
+ *val += dof[j] * TN[j](xi);
+}
+
+void tet4_eval3_2(double dof[4*3], double xi[3], double val[3])
+{
+ int i,j;
+ for (i = 0; i < 3; i++)
+ {
+ val[i] = 0.0;
+ for (j = 0; j < 4; j++)
+ val[i] += dof[4*i + j] * TN[j](xi);
+ }
+}
+
+void tet4_eval6_2(double dof[4*6], double xi[3], double val[6])
+{
+ int i,j;
+ for (i = 0; i < 6; i++)
+ {
+ val[i] = 0.0;
+ for (j = 0; j < 4; j++)
+ val[i] += dof[4*i + j] * TN[j](xi);
+ }
+}
+
+
+/*
+ * Input : xs -- an [npts x 3] matrix
+ * Output: tab -- an [npts x 4] matrix
+ */
+void tet4_tabulate(double *xs, int npts, double *tab)
+{
+ int i;
+ double *x;
+ for (i = 0; i < npts; i++)
+ {
+ x = &(xs[3*i]);
+ tab[4*i + 0] = TN0(x);
+ tab[4*i + 1] = TN1(x);
+ tab[4*i + 2] = TN2(x);
+ tab[4*i + 3] = TN3(x);
+ }
+}
+
+
+/*
+ * Inputs:
+ * tab - an [npts x 4] matrix
+ * dof - an [4 x dims] matrix
+ *
+ * Output:
+ * vals - an [npts x dims] matrix
+ */
+void tet4_batch_eval(double *dof, int dims, double *tab, int npts, double *vals)
+{
+ int i,j,k;
+ for (i = 0; i < npts; i++)
+ {
+ for (j = 0; j < dims; j++)
+ {
+ //
+ // vals[i,j] = \sum_k (tab[i,k] * dof[k,j])
+ //
+ vals[dims*i+j] = 0.0;
+ for (k = 0; k < 4; k++)
+ vals[dims*i+j] += tab[4*i + k] * dof[dims*k + j]; // tab[i,k] * dof[k,j]
+ }
+ }
+}
+
+void tet4_batch_eval1(double dof[4*1], double *tab, int npts, double *vals)
+{
+ int i;
+ for (i = 0; i < npts; i++)
+ {
+ vals[i] = tab[4*i+0]*dof[0] + tab[4*i+1]*dof[1] + tab[4*i+2]*dof[2] + tab[4*i+3]*dof[3];
+ }
+}
+
+void tet4_batch_eval3(double dof[4*3], double *tab, int npts, double *vals)
+{
+ int i;
+ for (i = 0; i < npts; i++)
+ {
+ vals[3*i+0] = tab[4*i+0] * dof[3*0+0]
+ + tab[4*i+1] * dof[3*1+0]
+ + tab[4*i+2] * dof[3*2+0]
+ + tab[4*i+3] * dof[3*3+0];
+
+ vals[3*i+1] = tab[4*i+0] * dof[3*0+1]
+ + tab[4*i+1] * dof[3*1+1]
+ + tab[4*i+2] * dof[3*2+1]
+ + tab[4*i+3] * dof[3*3+1];
+
+ vals[3*i+2] = tab[4*i+0] * dof[3*0+2]
+ + tab[4*i+1] * dof[3*1+2]
+ + tab[4*i+2] * dof[3*2+2]
+ + tab[4*i+3] * dof[3*3+2];
+ }
+}
+
+// ---------------------------------------------------------------------------
+
+double tet4_volume(double a[3], double b[3], double c[3], double d[3])
+{
+ // cf. http://mathworld.wolfram.com/Tetrahedron.html (32)
+ double V[4*4];
+ int j;
+
+ 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 tet4_test(double x[4*3], double pt[3])
+{
+ double vols[4];
+ double *a = &x[3*0];
+ double *b = &x[3*1];
+ double *c = &x[3*2];
+ double *d = &x[3*3];
+
+ vols[0] = tet4_volume(a,b,c,pt);
+ vols[1] = tet4_volume(b,d,c,pt);
+ vols[2] = tet4_volume(a,c,d,pt);
+ vols[3] = tet4_volume(a,d,b,pt);
+
+ if ((vols[0] < 0) || (vols[1] < 0) || (vols[2] < 0) || (vols[3] < 0))
+ {
+ return -1;
+ }
+
+ return 0;
+}
+
+int tet4_test2(double x[4*3], double pt[3], double vols[4])
+{
+ double *a = &x[3*0];
+ double *b = &x[3*1];
+ double *c = &x[3*2];
+ double *d = &x[3*3];
+
+ vols[0] = tet4_volume(a,b,c,pt);
+ vols[1] = tet4_volume(b,d,c,pt);
+ vols[2] = tet4_volume(a,c,d,pt);
+ vols[3] = tet4_volume(a,d,b,pt);
+
+ if ((vols[0] < 0) || (vols[1] < 0) || (vols[2] < 0) || (vols[3] < 0))
+ {
+ return -1;
+ }
+
+ return 0;
+}
Added: cs/cigma/trunk/tests/test_tet4.c
===================================================================
--- cs/cigma/trunk/tests/test_tet4.c 2007-05-10 00:24:12 UTC (rev 6827)
+++ cs/cigma/trunk/tests/test_tet4.c 2007-05-10 00:26:22 UTC (rev 6828)
@@ -0,0 +1,49 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <cigma/tet4.h>
+#include "test_data.h"
+
+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]);
+ }
+}
+
+int main(int argc, char *argv[])
+{
+ int n;
+ double T[4*3];
+ double J[9];
+ double tab[NQ*4];
+ 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];
+
+ // calculate jacobian matrix
+ tet4_jacobian_2(T, J);
+
+ // 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");
+ }
+
+ return 0;
+}
+
+
More information about the cig-commits
mailing list