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

[luis at geodynamics.org]

Author: luis
Date: 2009-01-28 09:05:02 -0800 (Wed, 28 Jan 2009)
New Revision: 13976

Modified:
   cs/cigma/trunk/src/fe_hex8.cpp
   cs/cigma/trunk/src/fe_quad4.cpp
   cs/cigma/trunk/src/fe_tri3.cpp
Log:
Fixed volume methods for hex8, quad4, tri3 elements

This was causing division by zero in some cases, since error
residuals are weighed by the corresponding element volumes

Modified: cs/cigma/trunk/src/fe_hex8.cpp
===================================================================
--- cs/cigma/trunk/src/fe_hex8.cpp	2009-01-28 17:05:01 UTC (rev 13975)
+++ cs/cigma/trunk/src/fe_hex8.cpp	2009-01-28 17:05:02 UTC (rev 13976)
@@ -1,4 +1,4 @@
-//#include <iostream>
+#include <iostream>
 #include <cmath>
 #include "fe_hex8.h"
 
@@ -141,8 +141,8 @@
 
 double hex8::volume()
 {
-    //using namespace std;
 
+    double vol = 0.0;
 
     /*
     // 
@@ -198,11 +198,8 @@
     #define TRIPLE (A[0] * B[1] * C[2] - A[0] * C[1] * B[2] - A[1] * B[0] * C[2] + A[1] * C[0] * B[2] + A[2] * B[0] * C[1] - A[2] * C[0] * B[1])
 
     int i;
-    double vol = 0;
     double A[3],B[3],C[3];
 
-    vol = 0.0;
-
     for (i = 0; i < 3; i++)
     {
         A[i] =  (X(7) - X(0));
@@ -219,24 +216,18 @@
     }
 
     vol += TRIPLE;
+    vol /= 12.0;
 
-    return vol/12.0;
+    // */
 
-    #undef X
-    #undef TRIPLE
-    */
 
-
-    /*
+    //*
     #define X(n)   (globverts[3*(n) + i])
     #define TRIPLE (A[0] * B[1] * C[2] - A[0] * C[1] * B[2] - A[1] * B[0] * C[2] + A[1] * C[0] * B[2] + A[2] * B[0] * C[1] - A[2] * C[0] * B[1])
 
     int i;
-    double vol = 0;
     double A[3], B[3], C[3];
 
-    for (int n = 0; n < 8; n++){for (i = 0; i < 3; i++){std::cerr << X(n) << " ";}std::cerr << std::endl;}
-
     for (i = 0; i < 3; i++)
     {
         A[i] = (X(7) - X(1)) + (X(6) - X(0));
@@ -260,26 +251,20 @@
         C[i] = (X(7) - X(4)) + (X(3) - X(0));
     }
     vol += (TRIPLE);
+    vol /= 12.0;
+    // */
 
-    return vol / 12.0;
 
-    #undef TRIPLE
-    #undef X
-    */
-
-
     /*
-     * Permutation of vertices: (2 3)(6 7)
-     */
+    //
+    // Permutation of vertices: (2 3)(6 7)
+    //
     #define X(n)   (globverts[3*(n) + i])
     #define TRIPLE (A[0] * B[1] * C[2] - A[0] * C[1] * B[2] - A[1] * B[0] * C[2] + A[1] * C[0] * B[2] + A[2] * B[0] * C[1] - A[2] * C[0] * B[1])
 
     int i;
-    double vol = 0;
     double A[3], B[3], C[3];
 
-    //for (int n = 0; n < 8; n++){for (i = 0; i < 3; i++){std::cerr << X(n) << " ";}std::cerr << std::endl;}
-
     for (i = 0; i < 3; i++)
     {
         A[i] = (X(6) - X(1)) + (X(7) - X(0));
@@ -303,9 +288,25 @@
         C[i] = (X(6) - X(4)) + (X(2) - X(0));
     }
     vol += (TRIPLE);
+    vol /= 12.0;
+    // */
 
-    return vol / 12.0;
 
+    /*
+    using namespace std;
+    cout << vol << "; ";
+    for (int n = 0; n < 8; n++)
+    {
+        for (i = 0; i < 3; i++)
+            cout << X(n) << " ";
+        cout << "; ";
+    }
+    cout << endl;
+    // */
+
+
+    return vol;
+
     #undef TRIPLE
     #undef X
 }

Modified: cs/cigma/trunk/src/fe_quad4.cpp
===================================================================
--- cs/cigma/trunk/src/fe_quad4.cpp	2009-01-28 17:05:01 UTC (rev 13975)
+++ cs/cigma/trunk/src/fe_quad4.cpp	2009-01-28 17:05:02 UTC (rev 13976)
@@ -1,4 +1,5 @@
 #include "fe_quad4.h"
+#include <iostream>
 
 using namespace cigma;
 
@@ -107,9 +108,25 @@
 
 double quad4::volume()
 {
-    #define X(n) globverts[3*(n)+0]
-    #define Y(n) globverts[3*(n)+1]
-    return 0.5*(X(0) * Y(1) - X(1) * Y(0) + X(1) * Y(2) - X(2) * Y(1) + X(2) * Y(3) - X(3) * Y(2) + X(3) * Y(0) - X(0) * Y(3));
+    //#define X(n) globverts[3*(n)+0]
+    //#define Y(n) globverts[3*(n)+1]
+
+    #define X(n) globverts[2*(n)+0]
+    #define Y(n) globverts[2*(n)+1]
+
+    const double vol = 0.5*(X(0) * Y(1) - X(1) * Y(0) + X(1) * Y(2) - X(2) * Y(1) + X(2) * Y(3) - X(3) * Y(2) + X(3) * Y(0) - X(0) * Y(3));
+    
+    /*
+    using namespace std;
+    cout << vol << "; ";
+    for (int n = 0; n < 4; n++)
+    {
+        cout << X(n) << " " << Y(n) << "; ";
+    }
+    cout << endl;
+    // */
+
+    return vol;
     #undef X
     #undef Y
 }

Modified: cs/cigma/trunk/src/fe_tri3.cpp
===================================================================
--- cs/cigma/trunk/src/fe_tri3.cpp	2009-01-28 17:05:01 UTC (rev 13975)
+++ cs/cigma/trunk/src/fe_tri3.cpp	2009-01-28 17:05:02 UTC (rev 13976)
@@ -93,8 +93,8 @@
 
 double tri3::volume()
 {
-    #define X(n) globverts[3*(n)+0]
-    #define Y(n) globverts[3*(n)+1]
+    #define X(n) globverts[2*(n) + 0]
+    #define Y(n) globverts[2*(n) + 1]
     return 0.5*(X(0) * Y(1) - X(0) * Y(2) - X(1) * Y(0) + X(1) * Y(2) + X(2) * Y(0) - X(2) * Y(1));
     #undef X
     #undef Y