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

[luis at geodynamics.org]

Author: luis
Date: 2007-05-09 17:38:10 -0700 (Wed, 09 May 2007)
New Revision: 6830

Added:
   cs/cigma/trunk/include/cigma/fe.h
   cs/cigma/trunk/src/fe.c
Log:
Quick generalization of fe_space


Added: cs/cigma/trunk/include/cigma/fe.h
===================================================================
--- cs/cigma/trunk/include/cigma/fe.h	2007-05-10 00:27:52 UTC (rev 6829)
+++ cs/cigma/trunk/include/cigma/fe.h	2007-05-10 00:38:10 UTC (rev 6830)
@@ -0,0 +1,27 @@
+#ifndef __CIGMA_FE_H__
+#define __CIGMA_FE_H__
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+typedef double (*shape_fn)(double xi[3]);
+
+typedef struct {
+    int nsd;    /* number of spatial dimensions (1d, 2d, 3d) */
+    int ndim;   /* number of components (scalar, vector, tensor) */
+    int ndof;   /* number of degrees of freedom */
+    shape_fn *N;    /* array of basis functions */
+    shape_fn *dN;   /* [nsd x ndof] matrix of derivatives */
+} fe_space_t;
+
+void fe_jacobian(fe_space_t *fe, double *x, double xi[3], double J[3*3]);
+void fe_tabulate(fe_space_t *fe, double *xs, int npts, double *tab);
+void fe_eval(fe_space_t *fe, double *dof, double *xi, double *val);
+void fe_batch_eval(fe_space_t *fe, double *dof, double *tab, int npts, double *vals);
+
+
+#ifdef __cplusplus
+}
+#endif
+#endif

Added: cs/cigma/trunk/src/fe.c
===================================================================
--- cs/cigma/trunk/src/fe.c	2007-05-10 00:27:52 UTC (rev 6829)
+++ cs/cigma/trunk/src/fe.c	2007-05-10 00:38:10 UTC (rev 6830)
@@ -0,0 +1,86 @@
+#include <cigma/fe.h>
+
+void fe_jacobian(fe_space_t *fe, double *x, double xi[3], double J[3*3])
+{
+    int i,j,k;
+    int nsd = fe->nsd;
+    //int dims = fe->ndim;
+    int ndof = fe->ndof;
+    shape_fn *dN = fe->dN;
+    for (i = 0; i < nsd; i++)
+    {
+        for (j = 0; j < nsd; j++)
+        {
+            //
+            // J[i,j] = dx[i]/dx_j
+            //        = \sum_k  x[i,k] (dN[j]/dx_i)
+            //
+            J[nsd*i + j] = 0.0;
+            for (k = 0; k < ndof; k++)
+                J[nsd*i + j] += dN[nsd*k + j](xi) * x[nsd*i + k];
+        }
+    }
+}
+
+/*
+ * Input  : xs  - an [npts x nsd] matrix
+ * Output : tab - an [npts x ndof] matrix
+ */
+void fe_tabulate(fe_space_t *fe, double *xs, int npts, double *tab)
+{
+    int i,j;
+    double *x;
+    int nsd  = fe->nsd;
+    int ndof = fe->ndof;
+    shape_fn *N = fe->N;
+    for (i = 0; i < npts; i++)
+    {
+        x = &(xs[nsd*i]);
+        for (j = 0; j < ndof; j++)
+            tab[ndof*i + j] = N[j](x);
+    }
+}
+
+void fe_eval(fe_space_t *fe, double *dof, double *xi, double *val)
+{
+    int i,j;
+    int dims = fe->ndim;
+    int ndof = fe->ndof;
+    shape_fn *N = fe->N;
+    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 < ndof; j++)
+            val[i] += dof[ndof*i + j] * N[j](xi);
+    }
+}
+
+/*
+ * Inputs:
+ *  tab - an [npts x ndof] matrix
+ *  dof - an [ndof x dims] matrix
+ *
+ * Output:
+ *  vals - an [npts x dims] matrix
+ */
+void fe_batch_eval(fe_space_t *fe, double *dof, double *tab, int npts, double *vals)
+{
+    int i,j,k;
+    int dims = fe->ndim;
+    int ndof = fe->ndof;
+    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 < ndof; k++)
+                vals[dims*i + j] += tab[ndof*i + k] * dof[dims*k + j]; // tab[i,k] * dof[k,j]
+        }
+    }
+}