[cig-commits] r9071 - in cs/benchmark/cigma/trunk/src: . tests

[luis at geodynamics.org]

Author: luis
Date: 2008-01-16 11:32:22 -0800 (Wed, 16 Jan 2008)
New Revision: 9071

Added:
   cs/benchmark/cigma/trunk/src/Tri.cpp
   cs/benchmark/cigma/trunk/src/Tri.h
   cs/benchmark/cigma/trunk/src/tests/TestTri.cpp
Log:
Add 2D triangle element

Added: cs/benchmark/cigma/trunk/src/Tri.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.cpp	                        (rev 0)
+++ cs/benchmark/cigma/trunk/src/Tri.cpp	2008-01-16 19:32:22 UTC (rev 9071)
@@ -0,0 +1,98 @@
+#include "Tri.h"
+#include <cassert>
+
+// ---------------------------------------------------------------------------
+
+static void tri_shape(double u, double v, double w, double s[3])
+{
+    s[0] = 1.0 - u - v ;
+    s[1] =       u     ;
+    s[2] =           v ;
+}
+
+static void tri_grad_shape(double u, double v, double w, double s[3*3])
+{
+    #define S(i,j) s[3*(i) + (j)]
+
+    S(0,0) = -1.0;
+    S(0,1) = -1.0;
+    S(0,2) =  0.0;
+
+    S(1,0) = +1.0;
+    S(1,1) =  0.0;
+    S(1,2) =  0.0;
+
+    S(2,0) =  0.0;
+    S(2,1) = +1.0;
+    S(2,2) =  0.0;
+
+    #undef S
+}
+
+
+// ---------------------------------------------------------------------------
+
+cigma::Tri::Tri()
+{
+    const int tri_nno = 3;
+    const int tri_celldim = 3; // XXX
+    double verts[tri_nno * tri_celldim] = {
+        0.0, 0.0, 0.0,
+        1.0, 0.0, 0.0,
+        0.0, 1.0, 0.0,
+    };
+    set_reference_vertices(verts, tri_nno);
+}
+
+cigma::Tri::~Tri()
+{
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::Tri::shape(int num, double *points, double *values)
+{
+    const int nno = n_nodes();
+    for (int i = 0; i < num; i++)
+    {
+        double u = points[3*i + 0];
+        double v = points[3*i + 1];
+        double w = points[3*i + 2];
+        tri_shape(u, v, w, &values[nno*i]);
+    }
+}
+
+void cigma::Tri::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++)
+    {
+        double u = points[3*i + 0];
+        double v = points[3*i + 1];
+        double w = points[3*i + 2];
+        tri_grad_shape(u, v, w, &values[stride*i]);
+    }
+}
+
+
+//void cigma::Tri::xyz2uvw(double xyz[3], double uvw[3]) {}
+
+
+bool cigma::Tri::interior(double u, double v, double w)
+{
+    #define ZERO    (-1.0e-6)
+    #define ONE     (1 - ZERO)
+
+    if ((u < ZERO) || (v < ZERO) || (u > (ONE - v)))
+    {
+        return false;
+    }
+    return true;
+
+    #undef ZERO
+    #undef ONE
+}
+
+// ---------------------------------------------------------------------------

Added: cs/benchmark/cigma/trunk/src/Tri.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.h	                        (rev 0)
+++ cs/benchmark/cigma/trunk/src/Tri.h	2008-01-16 19:32:22 UTC (rev 9071)
@@ -0,0 +1,29 @@
+#ifndef __TRI_H__
+#define __TRI_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+    class Tri;
+}
+
+class cigma::Tri : public Cell
+{
+public:
+    Tri();
+    ~Tri();
+
+public:
+    int n_nodes() { return 3; }
+    int n_celldim() { return 3; } // XXX: change to 2
+    int n_dim() { return 3; }
+
+public:
+    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);
+};
+
+#endif

Added: cs/benchmark/cigma/trunk/src/tests/TestTri.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/tests/TestTri.cpp	                        (rev 0)
+++ cs/benchmark/cigma/trunk/src/tests/TestTri.cpp	2008-01-16 19:32:22 UTC (rev 9071)
@@ -0,0 +1,81 @@
+#include <iostream>
+#include <cstdlib>
+#include <cassert>
+
+//#include "../TextWriter.h"
+//#include "../VtkUgReader.h"
+#include "../Tri.h"
+
+using namespace cigma;
+using namespace std;
+
+int main(void)
+{
+    //TextWriter writer;
+    //writer.fp = stdout;
+
+    Tri tricell;
+
+    int ndofs = tricell.n_nodes();
+    const int npts = 4;
+    const int nsd = 3;
+
+    double refpts[npts*nsd] = {
+        0.0, 0.0, 0.0,
+        1.0, 0.0, 0.0,
+        0.0, 1.0, 0.0,
+        0.5, 0.5, 0.0,
+    };
+
+    double tab[npts*ndofs];
+    double tab_jet[npts*ndofs*nsd];
+    tricell.shape(npts, refpts, tab);
+    tricell.grad_shape(npts, refpts, tab_jet);
+
+    for (int i = 0; i < npts; i++)
+    {
+        double *pt = &refpts[nsd*i];
+        double u = pt[0];
+        double v = pt[1];
+        double w = pt[2];
+
+        cout << i << " : ";
+
+        cout << "x = ("
+             << u << " "
+             << v << " "
+             << w << ")  ";
+
+        bool inside = tricell.interior(u,v,w);
+        cout << "in = " << inside << "  ";
+
+        double *phi = &tab[ndofs*i];
+        cout << "phi = ("
+             << phi[0] << " "
+             << phi[1] << " "
+             << phi[2] << ")  ";
+
+        double *grad_phi = &tab_jet[i*ndofs*nsd];
+        cout << "grad_phi = ("
+             << "(" << grad_phi[3*0 + 0] << " "
+                    << grad_phi[3*0 + 1] << " "
+                    << grad_phi[3*0 + 2] << "),"
+
+             << "(" << grad_phi[3*1 + 0] << " "
+                    << grad_phi[3*1 + 1] << " "
+                    << grad_phi[3*1 + 2] << "),"
+
+             << "(" << grad_phi[3*2 + 0] << " "
+                    << grad_phi[3*2 + 1] << " "
+                    << grad_phi[3*2 + 2] << ")"
+             << ")  ";
+
+        double J[3][3];
+        double detJ = tricell.jacobian(u,v,w,J);
+        cout << "detJ = " << detJ;
+
+        cout << endl;
+    }
+
+    return 0;
+}