[cig-commits] r12856 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2008-09-10 12:51:12 -0700 (Wed, 10 Sep 2008)
New Revision: 12856

Modified:
   cs/cigma/trunk/src/Cell.cpp
Log:
Generic cell methods, using shape() and grad_shape() from children classes




Modified: cs/cigma/trunk/src/Cell.cpp
===================================================================
--- cs/cigma/trunk/src/Cell.cpp	2008-09-10 19:51:04 UTC (rev 12855)
+++ cs/cigma/trunk/src/Cell.cpp	2008-09-10 19:51:12 UTC (rev 12856)
@@ -128,17 +128,130 @@
     jac[1][0] = jac[1][1] = jac[1][2] = 0.0;
     jac[2][0] = jac[2][1] = jac[2][2] = 0.0;
 
-    /*
-    double s[3];
-    switch (celldim)
+    #define X(i)  globverts[3*(i) + 0]
+    #define Y(i)  globverts[3*(i) + 1]
+    #define Z(i)  globverts[3*(i) + 2]
+
+    double uvw[3] = {u,v,w};
+    double *grad = new double[ndofs*3];
+    double *s;
+
+    switch (nsd)
     {
-    case 3:
-    case 2:
-    case 1:
-    default:
-        assert(false);
-    };*/
+        case 3 :
 
+            grad_shape(1, uvw, grad);
+
+            for(int i = 0; i < n_nodes(); i++)
+            {
+                s = &grad[3*i];
+
+                jac[0][0] += X(i) * s[0];
+                jac[0][1] += Y(i) * s[0];
+                jac[0][2] += Z(i) * s[0];
+
+                jac[1][0] += X(i) * s[1];
+                jac[1][1] += Y(i) * s[1];
+                jac[1][2] += Z(i) * s[1];
+
+                jac[2][0] += X(i) * s[2];
+                jac[2][1] += Y(i) * s[2];
+                jac[2][2] += Z(i) * s[2];
+            }
+            return fabs(
+                    + jac[0][0] * jac[1][1] * jac[2][2]
+                    + jac[0][2] * jac[1][0] * jac[2][1]
+                    + jac[0][1] * jac[1][2] * jac[2][0]
+                    - jac[0][2] * jac[1][1] * jac[2][0]
+                    - jac[0][0] * jac[1][2] * jac[2][1]
+                    - jac[0][1] * jac[1][0] * jac[2][2]);
+
+        case 2 :
+
+            grad_shape(1, uvw, grad);
+
+            for(int i = 0; i < n_nodes(); i++)
+            {
+                s = &grad[2*i];
+
+                jac[0][0] += X(i) * s[0];
+                jac[0][1] += Y(i) * s[0];
+                jac[0][2] += Z(i) * s[0];
+
+                jac[1][0] += X(i) * s[1];
+                jac[1][1] += Y(i) * s[1];
+                jac[1][2] += Z(i) * s[1];
+            }
+            {
+                double a[3], b[3], c[3];
+
+                a[0]= X(1) - X(0);
+                a[1]= Y(1) - Y(0);
+                a[2]= Z(1) - Z(0);
+
+                b[0]= X(2) - X(0);
+                b[1]= Y(2) - Y(0);
+                b[2]= Z(2) - Z(0);
+
+                prodve(a, b, c);
+
+                jac[2][0] = c[0];
+                jac[2][1] = c[1];
+                jac[2][2] = c[2];
+          }
+          return sqrt(SQR(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+                      SQR(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
+                      SQR(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
+
+        case 1:
+
+            grad_shape(1, uvw, grad);
+
+            for(int i = 0; i < n_nodes(); i++)
+            {
+                s = &grad[i];
+
+                jac[0][0] += X(i) * s[0];
+                jac[0][1] += Y(i) * s[0];
+                jac[0][2] += Z(i) * s[0];
+            }
+            {
+                double a[3], b[3], c[3];
+                a[0]= X[1] - X(0);
+                a[1]= Y[1] - Y(0);
+                a[2]= Z[1] - Z(0);
+
+                if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
+                   (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2])))
+                {
+                    b[0] =  a[1];
+                    b[1] = -a[0];
+                    b[2] =  0.0;
+                }
+                else
+                {
+                    b[0] =  0.0;
+                    b[1] =  a[2];
+                    b[2] = -a[1];
+                }
+
+                prodve(a, b, c);
+
+                jac[1][0] = b[0];
+                jac[1][1] = b[1];
+                jac[1][2] = b[2];
+
+                jac[2][0] = c[0];
+                jac[2][1] = c[1];
+                jac[2][2] = c[2];
+          }
+          return sqrt(SQR(jac[0][0])+SQR(jac[0][1])+SQR(jac[0][2]));
+    }
+
+    #undef X
+    #undef Y
+    #undef Z
+    
     return 0.0;
 }
 
@@ -156,7 +269,9 @@
 
     shape(1, point, N);
 
+    // 
     // value[i] = N[0]*dofs[0,i] + N[1]*dofs[1,i] + ... 
+    //
     for (int i = 0; i < valdim; i++)
     {
         double sum = 0.0;
@@ -167,6 +282,51 @@
 }
 
 
+/*
+ * @param dofs is a rank-2 array of dimension [ndofs x valdim]
+ * @param point is a rank-1 array of dimension [celldim]
+ * @param value is a rank-1 array of dimension [valdim]
+ */
+void cigma::Cell::interpolate_grad(double *dofs, double *point, double *value, int stride=1, double invjac[3][3]=NULL)
+{
+    assert(celldim > 0);
+    assert(nno > 0);
+    assert(ndofs > 0);
+    
+    int i, j, k;
+
+    double dfdu[3] = {0.0, 0.0, 0.0};
+    double grad[ndofs*3];
+    double *s;
+
+    grad_shape(1, point, grad);
+
+    k = 0;
+    for (int i = 0; i < ndofs; i++)
+    {
+        s = &grad[3*i];
+        for (j = 0; j < celldim; j++)
+        {
+            dfdu[j] += dofs[k] * s[j]
+        }
+        k += stride;
+    }
+
+    if (invjac)
+    {
+        matvec(invjac, dNdu, value);
+    }
+    else
+    {
+        double jac[3][3], inv[3][3];
+        jacobian(point, jac);
+        inv3x3(jac, inv);
+        matvec(inv, dfdu, value);
+    }
+}
+
+
+
 /* isoparametric reference map */
 void cigma::Cell::uvw2xyz(double uvw[3], double xyz[3])
 {