[cig-commits] r11699 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Tue Apr 1 01:40:23 PDT 2008
Author: luis
Date: 2008-04-01 01:40:22 -0700 (Tue, 01 Apr 2008)
New Revision: 11699
Modified:
cs/benchmark/cigma/trunk/src/Hex.cpp
cs/benchmark/cigma/trunk/src/Hex.h
cs/benchmark/cigma/trunk/src/Quad.cpp
cs/benchmark/cigma/trunk/src/Quad.h
cs/benchmark/cigma/trunk/src/Tet.cpp
cs/benchmark/cigma/trunk/src/Tet.h
cs/benchmark/cigma/trunk/src/Tri.cpp
cs/benchmark/cigma/trunk/src/Tri.h
Log:
Implemented volume method on each cell type
Modified: cs/benchmark/cigma/trunk/src/Hex.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.cpp 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Hex.cpp 2008-04-01 08:40:22 UTC (rev 11699)
@@ -1,3 +1,5 @@
+//#include <iostream>
+#include <cmath>
#include "Hex.h"
using namespace cigma;
@@ -137,3 +139,174 @@
#undef ONE
#undef ZERO
+
+double Hex::volume()
+{
+ //using namespace std;
+
+
+ /*
+ //
+ // 6V = [(x7-x0),(x1-x0),(x3-x5)]
+ // + [(x7-x0),(x4-x0),(x5-x6)]
+ // + [(x7-x0),(x2-x0),(x6-x3)]
+ //
+ // where
+ //
+ // xn = (globverts[3*n+0],globverts[3*n+1],globverts[3*n+2])
+ //
+ // | A_x B_x C_x |
+ // [A,B,C] = | A_y B_y C_y |
+ // | A_z B_z C_z |
+ //
+ double x70[3],x10[3],x35[3],x40[3],x56[3],x20[3],x63[3];
+ #define X(j) globverts[3*(j)+i]
+ for (int i = 0; i < 3; i++)
+ {
+ x70[i] = X(7) - X(0);
+ x10[i] = X(1) - X(0);
+ x35[i] = X(3) - X(5);
+ x40[i] = X(4) - X(0);
+ x56[i] = X(5) - X(6);
+ x20[i] = X(2) - X(0);
+ x63[i] = X(6) - X(3);
+ }
+ #undef X
+ return (
+ x70[0] * x10[1] * x35[2]
+ - x70[0] * x35[1] * x10[2]
+ - x70[1] * x10[0] * x35[2]
+ + x70[1] * x35[0] * x10[2]
+ + x70[2] * x10[0] * x35[1]
+ - x70[2] * x35[0] * x10[1]
+ + x70[0] * x40[1] * x56[2]
+ - x70[0] * x56[1] * x40[2]
+ - x70[1] * x40[0] * x56[2]
+ + x70[1] * x56[0] * x40[2]
+ + x70[2] * x40[0] * x56[1]
+ - x70[2] * x56[0] * x40[1]
+ + x70[0] * x20[1] * x63[2]
+ - x70[0] * x63[1] * x20[2]
+ - x70[1] * x20[0] * x63[2]
+ + x70[1] * x63[0] * x20[2]
+ + x70[2] * x20[0] * x63[1]
+ - x70[2] * x63[0] * x20[1]
+ ) / 6.0;*/
+
+
+ /*
+ #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];
+
+ vol = 0.0;
+
+ for (i = 0; i < 3; i++)
+ {
+ A[i] = (X(7) - X(0));
+ B[i] = (X(1) + X(6)) - (X(4) + X(5));
+ C[i] = -(X(1) - X(6)) + (X(4) - X(5));
+ }
+
+ vol += TRIPLE;
+
+ for (i = 0; i < 3; i++)
+ {
+ B[i] = (X(1) + X(6)) - (X(3) + X(2));
+ C[i] = (X(1) - X(6)) + (X(3) - X(2));
+ }
+
+ vol += TRIPLE;
+
+ 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));
+ B[i] = (X(7) - X(2));
+ C[i] = (X(3) - X(0));
+ }
+ vol += (TRIPLE);
+
+ for (i = 0; i < 3; i++)
+ {
+ A[i] = (X(6) - X(0));
+ B[i] = (X(7) - X(2)) + (X(5) - X(0));
+ C[i] = (X(7) - X(4));
+ }
+ vol += (TRIPLE);
+
+ for (i = 0; i < 3; i++)
+ {
+ A[i] = (X(7) - X(1));
+ B[i] = (X(5) - X(0));
+ C[i] = (X(7) - X(4)) + (X(3) - X(0));
+ }
+ vol += (TRIPLE);
+
+ return vol / 12.0;
+
+ #undef TRIPLE
+ #undef X
+ */
+
+
+ /*
+ * 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));
+ B[i] = (X(6) - X(3));
+ C[i] = (X(2) - X(0));
+ }
+ vol += (TRIPLE);
+
+ for (i = 0; i < 3; i++)
+ {
+ A[i] = (X(7) - X(0));
+ B[i] = (X(6) - X(3)) + (X(5) - X(0));
+ C[i] = (X(6) - X(4));
+ }
+ vol += (TRIPLE);
+
+ for (i = 0; i < 3; i++)
+ {
+ A[i] = (X(6) - X(1));
+ B[i] = (X(5) - X(0));
+ C[i] = (X(6) - X(4)) + (X(2) - X(0));
+ }
+ vol += (TRIPLE);
+
+ return vol / 12.0;
+
+ #undef TRIPLE
+ #undef X
+}
Modified: cs/benchmark/cigma/trunk/src/Hex.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.h 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Hex.h 2008-04-01 08:40:22 UTC (rev 11699)
@@ -24,6 +24,7 @@
void shape(int num, double *points, double *values);
void grad_shape(int num, double *points, double *values);
bool interior(double u, double v, double w);
+ double volume();
};
Modified: cs/benchmark/cigma/trunk/src/Quad.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Quad.cpp 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Quad.cpp 2008-04-01 08:40:22 UTC (rev 11699)
@@ -108,4 +108,13 @@
#undef ONE
#undef ZERO
+double Quad::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));
+ #undef X
+ #undef Y
+}
+
// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/Quad.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Quad.h 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Quad.h 2008-04-01 08:40:22 UTC (rev 11699)
@@ -24,6 +24,7 @@
void shape(int num, double *points, double *values);
void grad_shape(int num, double *points, double *values);
bool interior(double u, double v, double w);
+ double volume();
};
#endif
Modified: cs/benchmark/cigma/trunk/src/Tet.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.cpp 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Tet.cpp 2008-04-01 08:40:22 UTC (rev 11699)
@@ -1,3 +1,5 @@
+#include <cmath>
+
#include "Tet.h"
#include "Numeric.h"
@@ -144,3 +146,37 @@
#undef ZERO
}
+double Tet::volume()
+{
+ //
+ // | x0 y0 z0 1 |
+ // 6V = | x1 y1 z1 1 |
+ // | x2 y2 z2 1 |
+ // | x3 y3 z3 1 |
+ //
+ double *a = globverts;
+ return fabs(a[0] * a[4] * a[8]
+ - a[0] * a[4] * a[11]
+ + a[0] * a[7] * a[11]
+ - a[0] * a[7] * a[5]
+ + a[0] * a[10] * a[5]
+ - a[0] * a[10] * a[8]
+ - a[3] * a[1] * a[8]
+ + a[3] * a[1] * a[11]
+ - a[3] * a[7] * a[11]
+ + a[3] * a[7] * a[2]
+ - a[3] * a[10] * a[2]
+ + a[3] * a[10] * a[8]
+ + a[6] * a[1] * a[5]
+ - a[6] * a[1] * a[11]
+ + a[6] * a[4] * a[11]
+ - a[6] * a[4] * a[2]
+ + a[6] * a[10] * a[2]
+ - a[6] * a[10] * a[5]
+ - a[9] * a[1] * a[5]
+ + a[9] * a[1] * a[8]
+ - a[9] * a[4] * a[8]
+ + a[9] * a[4] * a[2]
+ - a[9] * a[7] * a[2]
+ + a[9] * a[7] * a[5]) / 6.0;
+}
Modified: cs/benchmark/cigma/trunk/src/Tet.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.h 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Tet.h 2008-04-01 08:40:22 UTC (rev 11699)
@@ -25,6 +25,7 @@
void grad_shape(int num, double *points, double *values);
void xyz2uvw(double xyz[3], double uvw[3]);
bool interior(double u, double v, double w);
+ double volume();
};
Modified: cs/benchmark/cigma/trunk/src/Tri.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.cpp 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Tri.cpp 2008-04-01 08:40:22 UTC (rev 11699)
@@ -93,4 +93,13 @@
#undef ONE
+double Tri::volume()
+{
+ #define X(n) globverts[3*(n)+0]
+ #define Y(n) globverts[3*(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
+}
+
// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/Tri.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.h 2008-04-01 08:40:21 UTC (rev 11698)
+++ cs/benchmark/cigma/trunk/src/Tri.h 2008-04-01 08:40:22 UTC (rev 11699)
@@ -25,6 +25,7 @@
void grad_shape(int num, double *points, double *values);
//void xyz2uvw(double xyz[3], double uvw[3]);
bool interior(double u, double v, double w);
+ double volume();
};
#endif
More information about the cig-commits
mailing list