[cig-commits] r8914 - cs/benchmark/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2007-12-19 12:02:06 -0800 (Wed, 19 Dec 2007)
New Revision: 8914

Modified:
   cs/benchmark/cigma/trunk/src/Cell.cpp
   cs/benchmark/cigma/trunk/src/Cell.h
Log:
Changes to Cell.cpp
 * Method Cell::xyz2uvw() shouldn't be a pure virtual method.
 * Also, Cell::set_reference_vertices() should reset globverts.

Modified: cs/benchmark/cigma/trunk/src/Cell.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Cell.cpp	2007-12-19 20:01:58 UTC (rev 8913)
+++ cs/benchmark/cigma/trunk/src/Cell.cpp	2007-12-19 20:02:06 UTC (rev 8914)
@@ -63,6 +63,7 @@
             refverts[n] = vertices[n];
         }
     }
+    globverts = refverts;
 }
 
 
@@ -324,7 +325,73 @@
     }
 }
 
+
+/*
+ * Inverse reference map
+ */
+void cigma::Cell::xyz2uvw(double xyz[3], double uvw[3])
+{
+    // general newton routine for the nonlinear case
+    uvw[0] = uvw[1] = uvw[2] = 0.0;
 
+    int iter = 1, maxiter = 20;
+    double error = 1.0, tol = 1.0e-6;
+
+    while ((error > tol) && (iter < maxiter))
+    {
+        double jac[3][3];
+        if (!jacobian(uvw[0], uvw[1], uvw[2], jac))
+            break;
+
+        double xn = 0.0, yn = 0.0, zn = 0.0;
+
+
+        double s[ndofs];
+
+        shape(1, uvw, s);
+
+        #define X(i)  globverts[3*(i) + 0]
+        #define Y(i)  globverts[3*(i) + 1]
+        #define Z(i)  globverts[3*(i) + 2]
+
+        for (int i = 0; i < nno; i++)
+        {
+            xn += X(i) * s[i];
+            yn += Y(i) * s[i];
+            zn += Z(i) * s[i];
+        }
+
+        #undef X
+        #undef Y
+        #undef Z
+
+        double inv[3][3];
+        inv3x3(jac, inv);
+
+        double un = uvw[0] + inv[0][0] * (xyz[0] - xn)
+                           + inv[1][0] * (xyz[1] - yn)
+                           + inv[2][0] * (xyz[2] - zn);
+
+        double vn = uvw[1] + inv[0][1] * (xyz[0] - xn)
+                           + inv[1][1] * (xyz[1] - yn)
+                           + inv[2][1] * (xyz[2] - zn);
+
+        double wn = uvw[2] + inv[0][2] * (xyz[0] - xn)
+                           + inv[1][2] * (xyz[1] - yn)
+                           + inv[2][2] * (xyz[2] - zn);
+
+        error = sqrt(SQR(un - uvw[0]) +
+                     SQR(vn - uvw[1]) +
+                     SQR(wn - uvw[2]));
+
+        uvw[0] = un;
+        uvw[1] = vn;
+        uvw[2] = wn;
+
+        iter++;
+    }
+}
+
 
 /* isoparametric reference map */
 void cigma::Cell::uvw2xyz(double uvw[3], double xyz[3])

Modified: cs/benchmark/cigma/trunk/src/Cell.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Cell.h	2007-12-19 20:01:58 UTC (rev 8913)
+++ cs/benchmark/cigma/trunk/src/Cell.h	2007-12-19 20:02:06 UTC (rev 8914)
@@ -55,7 +55,7 @@
     void interpolate(double *dofs, double *point, double *value, int valdim);
     void interpolate_grad(double *dofs, double *point, double *value, int stride=1, double invjac[3][3]=0);
 
-    virtual void xyz2uvw(double xyz[3], double uvw[3]) = 0;
+    virtual void xyz2uvw(double xyz[3], double uvw[3]);
     void uvw2xyz(double uvw[3], double xyz[3]);