[cig-commits] r13992 - in cs/cigma/trunk: . src
luis at geodynamics.org
luis at geodynamics.org
Wed Jan 28 09:05:26 PST 2009
Author: luis
Date: 2009-01-28 09:05:26 -0800 (Wed, 28 Jan 2009)
New Revision: 13992
Added:
cs/cigma/trunk/src/fe_tet10.cpp
cs/cigma/trunk/src/fe_tet10.h
Modified:
cs/cigma/trunk/Makefile.am
cs/cigma/trunk/src/Cell.cpp
cs/cigma/trunk/src/Cell.h
cs/cigma/trunk/src/cli_info_cmd.cpp
Log:
Add quadratic tet10 element. This should work fine (modulo node ordering) on the appropriate mesh files.
Modified: cs/cigma/trunk/Makefile.am
===================================================================
--- cs/cigma/trunk/Makefile.am 2009-01-28 17:05:24 UTC (rev 13991)
+++ cs/cigma/trunk/Makefile.am 2009-01-28 17:05:26 UTC (rev 13992)
@@ -160,6 +160,8 @@
src/Cell.cpp \
src/fe_tet4.h \
src/fe_tet4.cpp \
+ src/fe_tet10.h \
+ src/fe_tet10.cpp \
src/fe_hex8.h \
src/fe_hex8.cpp \
src/fe_quad4.h \
Modified: cs/cigma/trunk/src/Cell.cpp
===================================================================
--- cs/cigma/trunk/src/Cell.cpp 2009-01-28 17:05:24 UTC (rev 13991)
+++ cs/cigma/trunk/src/Cell.cpp 2009-01-28 17:05:26 UTC (rev 13992)
@@ -7,6 +7,7 @@
#include "Cell.h"
#include "fe_hex8.h"
#include "fe_tet4.h"
+#include "fe_tet10.h"
#include "fe_quad4.h"
#include "fe_tri3.h"
#include "fe_tri6.h"
@@ -22,6 +23,7 @@
CellTypeMapEntry cellTypeMapEntries[] = {
CellTypeMapEntry("hex8", Cell::HEX8),
CellTypeMapEntry("tet4", Cell::TET4),
+ CellTypeMapEntry("tet10", Cell::TET10),
CellTypeMapEntry("quad4", Cell::QUAD4),
CellTypeMapEntry("tri3", Cell::TRI3),
CellTypeMapEntry("tri6", Cell::TRI6)
@@ -53,6 +55,7 @@
{
case HEX8: return "hex8";
case TET4: return "tet4";
+ case TET10: return "tet10";
case QUAD4: return "quad4";
case TRI3: return "tri3";
case TRI6: return "tri6";
@@ -67,6 +70,7 @@
{
case HEX8: tmp.reset(new hex8); break;
case TET4: tmp.reset(new tet4); break;
+ case TET10: tmp.reset(new tet10); break;
case QUAD4: tmp.reset(new quad4); break;
case TRI3: tmp.reset(new tri3); break;
case TRI6: tmp.reset(new tri6); break;
Modified: cs/cigma/trunk/src/Cell.h
===================================================================
--- cs/cigma/trunk/src/Cell.h 2009-01-28 17:05:24 UTC (rev 13991)
+++ cs/cigma/trunk/src/Cell.h 2009-01-28 17:05:26 UTC (rev 13992)
@@ -41,6 +41,7 @@
TRI6,
QUAD4,
TET4,
+ TET10,
HEX8
} type;
Modified: cs/cigma/trunk/src/cli_info_cmd.cpp
===================================================================
--- cs/cigma/trunk/src/cli_info_cmd.cpp 2009-01-28 17:05:24 UTC (rev 13991)
+++ cs/cigma/trunk/src/cli_info_cmd.cpp 2009-01-28 17:05:26 UTC (rev 13992)
@@ -22,8 +22,10 @@
//
const char *indent = " ";
cout << "List of elements available for use in a Field:" << endl;
- cout << indent << "tet4, hex8" << endl;
- cout << indent << "tri3, tri6, quad4" << endl;
+ cout << indent << "tet4, tet10" << endl;
+ cout << indent << "hex8" << endl;
+ cout << indent << "tri3, tri6" << endl;
+ cout << indent << "quad4" << endl;
// XXX: restore these one by one.
//cout << indent << "hex8, hex20, hex27" << endl;
Added: cs/cigma/trunk/src/fe_tet10.cpp
===================================================================
--- cs/cigma/trunk/src/fe_tet10.cpp (rev 0)
+++ cs/cigma/trunk/src/fe_tet10.cpp 2009-01-28 17:05:26 UTC (rev 13992)
@@ -0,0 +1,217 @@
+#include <cmath>
+#include "fe_tet10.h"
+#include "Numeric.h"
+
+using namespace cigma;
+
+static void tet_shape(double u, double v, double w, double s[4])
+{
+ s[0] = 1 - 3*u - 3*v + 4*u*v + 4*u*w + 4*v*w + 2*u*u + 2*v*v + 2*w*w;
+ s[1] = 2*u*u - u;
+ s[2] = 2*v*v - v;
+ s[3] = 2*w*w - w;
+ s[4] = 4*u*v;
+ s[5] = 4*v*w;
+ s[6] = 4*u*w;
+ s[7] = 4*(w - u*w - v*w - w*w);
+ s[8] = 4*(v - u*v - v*w - v*v);
+ s[9] = 4*(u - u*v - u*w - u*u);
+}
+
+static void tet_grad_shape(double u, double v, double w, double s[4*3])
+{
+ #define S(i,j) s[3*(i) + (j)]
+
+ S(0,0) = -3 + 4*v + 4*w + 4*u;
+ S(0,1) = -3 + 4*v + 4*w + 4*u;
+ S(0,2) = -3 + 4*v + 4*w + 4*u;
+
+ S(1,0) = 4*u - 1;
+ S(1,1) = 0;
+ S(1,2) = 0;
+
+ S(2,0) = 0;
+ S(2,1) = 4*v - 1;
+ S(2,2) = 0;
+
+ S(3,0) = 0;
+ S(3,1) = 0;
+ S(3,2) = 4*w - 1;
+
+ S(4,0) = 4*v;
+ S(4,1) = 4*u;
+ S(4,2) = 0;
+
+ S(5,0) = 0;
+ S(5,1) = 4*w;
+ S(5,2) = 4*v;
+
+ S(6,0) = 4*w;
+ S(6,1) = 0;
+ S(6,2) = 4*u;
+
+ S(7,0) = -4*w;
+ S(7,1) = -4*w;
+ S(7,2) = 4 - 4*u - 4*v - 8*w;
+
+ S(8,0) = -4*v;
+ S(8,1) = 4 - 4*u - 4*w - 8*v;
+ S(8,2) = -4*v;
+
+ S(9,0) = 4 - 4*v - 4*w - 8*u;
+ S(9,1) = -4*u;
+ S(9,2) = -4*u;
+
+ #undef S
+}
+
+double tet10::refverts[10*3] = {
+ 0.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 0.0, 1.0, 0.0,
+ 0.0, 0.0, 1.0,
+ 0.5, 0.5, 0.0,
+ 0.0, 0.5, 0.5,
+ 0.5, 0.0, 0.5,
+ 0.0, 0.0, 0.5,
+ 0.0, 0.5, 0.0,
+ 0.5, 0.0, 0.0
+};
+
+tet10::tet10()
+{
+ Cell::init();
+}
+
+tet10::tet10(const tet10& t)
+{
+ Cell::reinit(t);
+}
+
+tet10::~tet10()
+{
+}
+
+void tet10::shape(int num, double *points, double *values)
+{
+ const int nno = n_nodes();
+ for (int i = 0; i < num; i++)
+ {
+ const double u = points[3*i + 0];
+ const double v = points[3*i + 1];
+ const double w = points[3*i + 2];
+ tet_shape(u, v, w, &values[nno*i]);
+ }
+}
+
+void tet10::grad_shape(int num, double *points, double *values)
+{
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ const int stride = nno * celldim;
+ for (int i = 0; i < num; i++)
+ {
+ const double u = points[3*i + 0];
+ const double v = points[3*i + 1];
+ const double w = points[3*i + 2];
+ tet_grad_shape(u, v, w, &values[stride*i]);
+ }
+}
+
+void tet10::xyz2uvw(double xyz[3], double uvw[3])
+{
+ double det;
+ double mat[3][3], b[3];
+
+ /*
+ * vertices = {(x_0,y_0,z_0),
+ * (x_1,y_1,z_1),
+ * (x_2,y_2,z_2),
+ * (x_3,y_3,z_3)}
+ *
+ * [ (x[1]-x[0]) (x[2]-x[0]) (x[3]-x[0]) ]
+ * mat = [ (y[1]-y[0]) (y[2]-y[0]) (y[3]-y[0]) ]
+ * [ (z[1]-z[0]) (z[2]-z[0]) (z[3]-z[0]) ]
+ *
+ * x[i] = globverts[3*i+0]
+ * y[i] = globverts[3*i+1]
+ * z[i] = globverts[3*i+2]
+ *
+ */
+
+ #define X(i) globverts[3*(i)+0]
+ #define Y(i) globverts[3*(i)+1]
+ #define Z(i) globverts[3*(i)+2]
+
+ mat[0][0] = X(1) - X(0);
+ mat[0][1] = X(2) - X(0);
+ mat[0][2] = X(3) - X(0);
+
+ mat[1][0] = Y(1) - Y(0);
+ mat[1][1] = Y(2) - Y(0);
+ mat[1][2] = Y(3) - Y(0);
+
+ mat[2][0] = Z(1) - Z(0);
+ mat[2][1] = Z(2) - Z(0);
+ mat[2][2] = Z(3) - Z(0);
+
+ b[0] = xyz[0] - X(0);
+ b[1] = xyz[1] - Y(0);
+ b[2] = xyz[2] - Z(0);
+
+ #undef X
+ #undef Y
+ #undef Z
+
+ sys3x3(mat, b, uvw, &det);
+
+}
+
+bool tet10::interior(double u, double v, double w)
+{
+ #define ZERO (-1.0e-6)
+ #define ONE (1 - ZERO)
+
+ if ((u < ZERO) || (v < ZERO) || (w < ZERO) || (u > (ONE - v - w)))
+ return false;
+ return true;
+
+ #undef ONE
+ #undef ZERO
+}
+
+double tet10::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;
+}
+
Added: cs/cigma/trunk/src/fe_tet10.h
===================================================================
--- cs/cigma/trunk/src/fe_tet10.h (rev 0)
+++ cs/cigma/trunk/src/fe_tet10.h 2009-01-28 17:05:26 UTC (rev 13992)
@@ -0,0 +1,38 @@
+#ifndef __FE_TET10_H__
+#define __FE_TET10_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+ class tet10;
+}
+
+/**
+ * Subparametric quadratic tetrahedron.
+ */
+class cigma::tet10 : public Cell
+{
+public:
+ tet10();
+ tet10(const tet10& t);
+ ~tet10();
+
+ int n_nodes() const { return 10; }
+ int n_celldim() const { return 3; }
+
+ Cell::type cell_type() const { return TET10; }
+
+ void shape(int num, double *points, double *values);
+ 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();
+
+public:
+ static double refverts[10*3];
+
+};
+
+
+#endif
More information about the CIG-COMMITS
mailing list