[cig-commits] r6829 - in cs/cigma/trunk: include/cigma src

[luis at geodynamics.org]

Author: luis
Date: 2007-05-09 17:27:52 -0700 (Wed, 09 May 2007)
New Revision: 6829

Added:
   cs/cigma/trunk/include/cigma/hex8.h
   cs/cigma/trunk/src/hex8.c
Log:
Default finite element space for hexahedral elements


Added: cs/cigma/trunk/include/cigma/hex8.h
===================================================================
--- cs/cigma/trunk/include/cigma/hex8.h	2007-05-10 00:26:22 UTC (rev 6828)
+++ cs/cigma/trunk/include/cigma/hex8.h	2007-05-10 00:27:52 UTC (rev 6829)
@@ -0,0 +1,20 @@
+#ifndef __CIGMA_HEX8_H__
+#define __CIGMA_HEX8_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+extern const int hex8_nsd;  // 3
+extern const int hex8_ndof; // 8
+
+void hex8_jacobian(double x[8*3], double xi[3], double J[3*3]);
+void hex8_eval(double *dof, int dims, double xi[3], double *val);
+void hex8_tabulate(double *xs, int npts, double *tab);
+void hex8_batch_eval(double *dof, int dims, double *tab, int npts, double *vals);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif

Added: cs/cigma/trunk/src/hex8.c
===================================================================
--- cs/cigma/trunk/src/hex8.c	2007-05-10 00:26:22 UTC (rev 6828)
+++ cs/cigma/trunk/src/hex8.c	2007-05-10 00:27:52 UTC (rev 6829)
@@ -0,0 +1,138 @@
+#include <stdlib.h>
+#include <assert.h>
+#include <cigma/hex8.h>
+
+// Constants for this finite element space
+const int hex8_nsd = 3;
+const int hex8_ndof = 8;
+
+static double HN0(double xi[3]) {return 0.125*(1-xi[0])*(1-xi[1])*(1-xi[2]);}
+static double HN1(double xi[3]) {return 0.125*(1+xi[0])*(1-xi[1])*(1-xi[2]);}
+static double HN2(double xi[3]) {return 0.125*(1-xi[0])*(1+xi[1])*(1-xi[2]);}
+static double HN3(double xi[3]) {return 0.125*(1+xi[0])*(1+xi[1])*(1-xi[2]);}
+static double HN4(double xi[3]) {return 0.125*(1-xi[0])*(1-xi[1])*(1+xi[2]);}
+static double HN5(double xi[3]) {return 0.125*(1+xi[0])*(1-xi[1])*(1+xi[2]);}
+static double HN6(double xi[3]) {return 0.125*(1-xi[0])*(1+xi[1])*(1+xi[2]);}
+static double HN7(double xi[3]) {return 0.125*(1+xi[0])*(1+xi[1])*(1+xi[2]);}
+
+static double HN00(double xi[3]) {return 0.125*(-1)*(1-xi[1])*(1-xi[2]);}
+static double HN10(double xi[3]) {return 0.125*(+1)*(1-xi[1])*(1-xi[2]);}
+static double HN20(double xi[3]) {return 0.125*(-1)*(1+xi[1])*(1-xi[2]);}
+static double HN30(double xi[3]) {return 0.125*(+1)*(1+xi[1])*(1-xi[2]);}
+static double HN40(double xi[3]) {return 0.125*(-1)*(1-xi[1])*(1+xi[2]);}
+static double HN50(double xi[3]) {return 0.125*(+1)*(1-xi[1])*(1+xi[2]);}
+static double HN60(double xi[3]) {return 0.125*(-1)*(1+xi[1])*(1+xi[2]);}
+static double HN70(double xi[3]) {return 0.125*(+1)*(1+xi[1])*(1+xi[2]);}
+
+static double HN01(double xi[3]) {return 0.125*(1-xi[0])*(-1)*(1-xi[2]);}
+static double HN11(double xi[3]) {return 0.125*(1+xi[0])*(-1)*(1-xi[2]);}
+static double HN21(double xi[3]) {return 0.125*(1-xi[0])*(+1)*(1-xi[2]);}
+static double HN31(double xi[3]) {return 0.125*(1+xi[0])*(+1)*(1-xi[2]);}
+static double HN41(double xi[3]) {return 0.125*(1-xi[0])*(-1)*(1+xi[2]);}
+static double HN51(double xi[3]) {return 0.125*(1+xi[0])*(-1)*(1+xi[2]);}
+static double HN61(double xi[3]) {return 0.125*(1-xi[0])*(+1)*(1+xi[2]);}
+static double HN71(double xi[3]) {return 0.125*(1+xi[0])*(+1)*(1+xi[2]);}
+
+static double HN02(double xi[3]) {return 0.125*(1-xi[0])*(1-xi[1])*(-1);}
+static double HN12(double xi[3]) {return 0.125*(1+xi[0])*(1-xi[1])*(-1);}
+static double HN22(double xi[3]) {return 0.125*(1-xi[0])*(1+xi[1])*(-1);}
+static double HN32(double xi[3]) {return 0.125*(1+xi[0])*(1+xi[1])*(-1);}
+static double HN42(double xi[3]) {return 0.125*(1-xi[0])*(1-xi[1])*(+1);}
+static double HN52(double xi[3]) {return 0.125*(1+xi[0])*(1-xi[1])*(+1);}
+static double HN62(double xi[3]) {return 0.125*(1-xi[0])*(1+xi[1])*(+1);}
+static double HN72(double xi[3]) {return 0.125*(1+xi[0])*(1+xi[1])*(+1);}
+
+typedef double (*hex8_shape_fn)(double xi[3]);
+
+static hex8_shape_fn HN[8] = {HN0,HN1,HN2,HN3,HN4,HN5,HN6,HN7};
+
+static hex8_shape_fn dHN[8][3] = {
+    {HN00, HN01, HN02},
+    {HN10, HN11, HN12},
+    {HN20, HN21, HN22},
+    {HN30, HN31, HN32},
+    {HN40, HN41, HN42},
+    {HN50, HN51, HN52},
+    {HN60, HN61, HN62},
+    {HN70, HN71, HN72}
+};
+
+
+void hex8_jacobian(double x[8*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 < 8; k++)
+                J[3*i + j] += dHN[k][j](xi) * x[3*i + k];
+        }
+    }
+}
+
+void hex8_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 < 8; j++)
+        {
+            val[i] += dof[8*i + j] * HN[j](xi);
+        }
+    }
+}
+
+/*
+ * Input : xs  -- an [npts x 3] matrix
+ * Output: tab -- an [npts x 8] matrix
+ *
+ */
+void hex8_tabulate(double *xs, int npts, double *tab)
+{
+    int i;
+    double *x;
+    for (i = 0; i < npts; i++)
+    {
+        x = &(xs[3*i]);
+        tab[8*i + 0] = HN0(x);
+        tab[8*i + 1] = HN1(x);
+        tab[8*i + 2] = HN2(x);
+        tab[8*i + 3] = HN3(x);
+        tab[8*i + 4] = HN4(x);
+        tab[8*i + 5] = HN5(x);
+        tab[8*i + 6] = HN6(x);
+        tab[8*i + 7] = HN7(x);
+    }
+}
+
+
+/*
+ * Inputs:
+ *  tab - an [npts x 8] matrix
+ *  dof - an [8 x dims] matrix
+ *
+ * Output:
+ *  vals - an [npts x dims] matrix
+ */
+void hex8_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 < 8; k++)
+                vals[dims*i + j] += tab[4*i + k] * dof[dims*k + j]; // tab[i,k] * dof[k,j]
+        }
+    }
+}