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

[luis at geodynamics.org]

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

Added:
   cs/cigma/trunk/src/fe_tri6.cpp
   cs/cigma/trunk/src/fe_tri6.h
Log:
Add quadratic tri6 element. This should just work with minimal changes.

Added: cs/cigma/trunk/src/fe_tri6.cpp
===================================================================
--- cs/cigma/trunk/src/fe_tri6.cpp	                        (rev 0)
+++ cs/cigma/trunk/src/fe_tri6.cpp	2009-01-28 17:05:13 UTC (rev 13983)
@@ -0,0 +1,117 @@
+#include "fe_tri6.h"
+
+using namespace cigma;
+
+static void tri_shape(double u, double v, double s[6])
+{
+    s[0] = 1 - 3*u - 3*v + 4*u*v + 2*u*u + 2*v*v;
+    s[1] = 2*u*u - u;
+    s[2] = 2*v*v - v;
+    s[3] = 4*u*v;
+    s[4] = 4*(v - u*v - v*v);
+    s[5] = 4*(u - u*v - u*u);
+}
+
+static void tri_grad_shape(double u, double v, double s[6*2])
+{
+    #define S(i,j) s[2*(i) + (j)]
+
+    S(0,0) = -3 + 4*u + 4*v;
+    S(0,1) = -3 + 4*u + 4*v;
+
+    S(1,0) = 4*u - 1;
+    S(1,1) = 0;
+
+    S(2,0) = 0;
+    S(2,1) = 4*v - 1;
+
+    S(3,0) = 4*v;
+    S(3,1) = 4*u;
+
+    S(4,0) = -4*v;
+    S(4,1) = 4 - 4*u - 8*v;
+
+    S(5,0) = 4 - 4*v - 8*u;
+    S(5,1) = -4*u;
+
+    #undef S
+}
+
+double tri6::refverts[6*2] = {
+    0.0, 0.0,
+    1.0, 0.0,
+    0.0, 1.0,
+    0.5, 0.5,
+    0.0, 0.5,
+    0.5, 0.0
+};
+
+
+tri6::tri6()
+{
+    Cell::init();
+}
+
+tri6::tri6(const tri6& t)
+{
+    Cell::reinit(t);
+}
+
+tri6::~tri6()
+{
+}
+
+void tri6::shape(int num, double *points, double *values)
+{
+    const int nno = n_nodes();
+    for (int i = 0; i < num; i++)
+    {
+        const double u = points[2*i + 0];
+        const double v = points[2*i + 1];
+        tri_shape(u, v, &values[nno*i]);
+    }
+}
+
+void tri6::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[2*i + 0];
+        const double v = points[2*i + 1];
+        tri_grad_shape(u, v, &values[stride*i]);
+    }
+}
+
+
+//void tri6::xyz2uvw(double xyz[3], double uvw[3]) {}
+
+
+#define ZERO    (-1.0e-6)
+#define ONE     (1 - ZERO)
+
+bool tri6::interior(double u, double v, double w)
+{
+
+    if ((u < ZERO) || (v < ZERO) || (u > (ONE - v)))
+    {
+        return false;
+    }
+    return true;
+}
+
+#undef ZERO
+#undef ONE
+
+
+double tri6::volume()
+{
+    #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
+}
+

Added: cs/cigma/trunk/src/fe_tri6.h
===================================================================
--- cs/cigma/trunk/src/fe_tri6.h	                        (rev 0)
+++ cs/cigma/trunk/src/fe_tri6.h	2009-01-28 17:05:13 UTC (rev 13983)
@@ -0,0 +1,36 @@
+#ifndef __FE_TRI6_H__
+#define __FE_TRI6_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+    class tri6;
+}
+
+/**
+ * Subparametric quadratic triangle.
+ */
+class cigma::tri6 : public Cell
+{
+public:
+    tri6();
+    tri6(const tri6& t);
+    ~tri6();
+
+    int n_nodes() const { return 6; }
+    int n_celldim() const { return 2; }
+
+    Cell::type cell_type() const { return TRI6; }
+
+    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[6*2];
+};
+
+#endif