[cig-commits] r7504 - in cs/cigma/branches/cigma-0.9: bin build src
luis at geodynamics.org
luis at geodynamics.org
Tue Jun 26 09:51:43 PDT 2007
Author: luis
Date: 2007-06-26 09:51:42 -0700 (Tue, 26 Jun 2007)
New Revision: 7504
Added:
cs/cigma/branches/cigma-0.9/src/elements.h
cs/cigma/branches/cigma-0.9/src/hex8.c
cs/cigma/branches/cigma-0.9/src/hex8.h
cs/cigma/branches/cigma-0.9/src/tet4.c
cs/cigma/branches/cigma-0.9/src/tet4.h
Modified:
cs/cigma/branches/cigma-0.9/bin/cigma.c
cs/cigma/branches/cigma-0.9/build/Makefile
Log:
darcs patches:
* Changed include in tet4.c
* Changed include in tet4.c
* Fixed typo (argv) in cigma.c
* Fixed makefile targets for tet4.o and hex8.o
* Added disloc3d solution to sandbox
Modified: cs/cigma/branches/cigma-0.9/bin/cigma.c
===================================================================
--- cs/cigma/branches/cigma-0.9/bin/cigma.c 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/bin/cigma.c 2007-06-26 16:51:42 UTC (rev 7504)
@@ -1,5 +1,5 @@
-int main(int argc, char *argc)
+int main(int argc, char *argv[])
{
return 0;
}
Modified: cs/cigma/branches/cigma-0.9/build/Makefile
===================================================================
--- cs/cigma/branches/cigma-0.9/build/Makefile 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/build/Makefile 2007-06-26 16:51:42 UTC (rev 7504)
@@ -48,8 +48,8 @@
$(CIGMA)/src/split.o
ELTOBJS = \
- $(CIGMA)/tet4.o \
- $(CIGMA)/hex8.o
+ $(CIGMA)/src/tet4.o \
+ $(CIGMA)/src/hex8.o
###############################################################################
Added: cs/cigma/branches/cigma-0.9/src/elements.h
===================================================================
--- cs/cigma/branches/cigma-0.9/src/elements.h 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/src/elements.h 2007-06-26 16:51:42 UTC (rev 7504)
@@ -0,0 +1,7 @@
+#ifndef __CIGMA_ELEMENTS_H__
+#define __CIGMA_ELEMENTS_H__
+
+#include "tet4.h"
+#include "hex8.h"
+
+#endif
Added: cs/cigma/branches/cigma-0.9/src/hex8.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/hex8.c 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/src/hex8.c 2007-06-26 16:51:42 UTC (rev 7504)
@@ -0,0 +1,138 @@
+#include <stdlib.h>
+#include <assert.h>
+#include "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]
+ }
+ }
+}
Added: cs/cigma/branches/cigma-0.9/src/hex8.h
===================================================================
--- cs/cigma/branches/cigma-0.9/src/hex8.h 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/src/hex8.h 2007-06-26 16:51:42 UTC (rev 7504)
@@ -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/branches/cigma-0.9/src/tet4.c
===================================================================
--- cs/cigma/branches/cigma-0.9/src/tet4.c 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/src/tet4.c 2007-06-26 16:51:42 UTC (rev 7504)
@@ -0,0 +1,411 @@
+#include <stdlib.h>
+#include <assert.h>
+#include "det.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] = (-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/branches/cigma-0.9/src/tet4.h
===================================================================
--- cs/cigma/branches/cigma-0.9/src/tet4.h 2007-06-26 16:48:09 UTC (rev 7503)
+++ cs/cigma/branches/cigma-0.9/src/tet4.h 2007-06-26 16:51:42 UTC (rev 7504)
@@ -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
More information about the cig-commits
mailing list