[cig-commits] r6671 - cs/cigma/trunk/sandbox/c
luis at geodynamics.org
luis at geodynamics.org
Tue Apr 24 15:08:22 PDT 2007
Author: luis
Date: 2007-04-24 15:08:22 -0700 (Tue, 24 Apr 2007)
New Revision: 6671
Added:
cs/cigma/trunk/sandbox/c/tet4.c
cs/cigma/trunk/sandbox/c/tet4.h
Log:
Tabulation of shape functions for tet4 element
Added: cs/cigma/trunk/sandbox/c/tet4.c
===================================================================
--- cs/cigma/trunk/sandbox/c/tet4.c 2007-04-24 22:05:03 UTC (rev 6670)
+++ cs/cigma/trunk/sandbox/c/tet4.c 2007-04-24 22:08:22 UTC (rev 6671)
@@ -0,0 +1,290 @@
+#include <stdlib.h>
+#include <assert.h>
+#include "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] = poly(-0.5,-0.5,-0.5,-0.5)
+// TN[1] = poly(+0.5, 0.5, 0.0, 0.0)
+// TN[2] = poly(+0.5, 0.0, 0.5, 0.0)
+// TN[3] = poly(+0.5, 0.0, 0.0, 0.5)
+//
+typedef double (*tet4_basis_t)(double xi[3]);
+
+static tet4_basis_t TN[4] = {TN0,TN1,TN2,TN3};
+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
+// polynomials
+//
+// dTN[i,j] = \frac{\partial{TN_i}}{\partial{\xi_j}}
+//
+// TN[0]_xi = poly(-0.5, 0, 0, 0)
+// TN[0]_eta = poly(-0.5, 0, 0, 0)
+// TN[0]_zeta = poly(-0.5, 0, 0, 0)
+//
+// TN[1]_xi = poly(+0.5, 0, 0, 0)
+// TN[1]_eta = poly( 0.0, 0, 0, 0)
+// TN[1]_zeta = poly( 0.0, 0, 0, 0)
+//
+// TN[2]_xi = poly( 0.0, 0, 0, 0)
+// TN[2]_eta = poly(+0.5, 0, 0, 0)
+// TN[2]_zeta = poly( 0.0, 0, 0, 0)
+//
+// TN[3]_xi = poly( 0.0, 0, 0, 0)
+// TN[3]_eta = poly( 0.0, 0, 0, 0)
+// TN[3]_zeta = poly(+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_jacobian(tet4_t *t, double xi[3], double *J)
+{
+ //
+ // 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;
+ double *x = t->x;
+ 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_2(tet4_t *t, double *J)
+{
+ //
+ // [ 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]
+ //
+ double *x = t->x;
+
+ 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(tet4_t *t, double *J)
+{
+ /* XXX: from notes.. */
+ J[3*0+0] = 0;
+ J[3*0+1] = 0;
+ J[3*0+2] = 0;
+
+ J[3*1+0] = 0;
+ J[3*1+1] = 0;
+ J[3*1+2] = 0;
+
+ J[3*2+0] = 0;
+ J[3*2+1] = 0;
+ J[3*2+2] = 0;
+}
+
+
+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;
+ 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], double xi[3], double *val)
+{
+ *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], double xi[3], double *val)
+{
+ 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], 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];
+ }
+}
+
Added: cs/cigma/trunk/sandbox/c/tet4.h
===================================================================
--- cs/cigma/trunk/sandbox/c/tet4.h 2007-04-24 22:05:03 UTC (rev 6670)
+++ cs/cigma/trunk/sandbox/c/tet4.h 2007-04-24 22:08:22 UTC (rev 6671)
@@ -0,0 +1,31 @@
+#ifndef __CIGMA_TET4_H__
+#define __CIGMA_TET4_H__
+
+extern const int tet4_nsd; // 3
+extern const int tet4_ndof; // 4
+
+typedef struct {
+ 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_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_eval(double *dof, int dims, double *xi, double *val);
+void tet4_eval1(double dof[4], double xi[3], double *val);
+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_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_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_eval3(double dof[4*3], double *tab, int npts, double *vals);
+
+#endif
More information about the cig-commits
mailing list