[cig-commits] r13064 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Oct 15 02:08:11 PDT 2008
Author: luis
Date: 2008-10-15 02:08:11 -0700 (Wed, 15 Oct 2008)
New Revision: 13064
Added:
cs/cigma/trunk/src/fe_hex8.cpp
cs/cigma/trunk/src/fe_hex8.h
cs/cigma/trunk/src/fe_quad4.cpp
cs/cigma/trunk/src/fe_quad4.h
cs/cigma/trunk/src/fe_tet4.cpp
cs/cigma/trunk/src/fe_tet4.h
cs/cigma/trunk/src/fe_tri3.cpp
cs/cigma/trunk/src/fe_tri3.h
Modified:
cs/cigma/trunk/src/Cell.cpp
cs/cigma/trunk/src/Cell.h
Log:
Cell subclasses
Modified: cs/cigma/trunk/src/Cell.cpp
===================================================================
--- cs/cigma/trunk/src/Cell.cpp 2008-10-15 09:08:08 UTC (rev 13063)
+++ cs/cigma/trunk/src/Cell.cpp 2008-10-15 09:08:11 UTC (rev 13064)
@@ -1 +1,414 @@
+#include <cstdlib>
+#include <cassert>
#include "Cell.h"
+#include "Numeric.h"
+
+using namespace cigma;
+
+
+Cell::Cell()
+{
+ //cout << "Calling Cell()" << endl;
+ refverts = 0;
+ globverts = 0;
+}
+
+
+Cell::~Cell()
+{
+ //cout << "Calling ~Cell()" << endl;
+ if (refverts != 0)
+ {
+ delete [] refverts;
+ }
+ if (globverts != 0)
+ {
+ delete [] globverts;
+ globverts = 0;
+ }
+}
+
+
+void Cell::set_reference_vertices(double *vertices)
+{
+ int i,j;
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ const int dim = n_dim();
+
+ if (refverts != 0) delete [] refverts;
+ if (globverts != 0) delete [] globverts;
+
+ refverts = new double[nno*celldim];
+ globverts = new double[nno*dim];
+
+ // copy data from vertices
+ for (i = 0; i < nno; i++)
+ {
+ for (j = 0; j < dim; j++)
+ {
+ globverts[dim*i + j] = 0.0;
+ }
+ for (j = 0; j < celldim; j++)
+ {
+ int n = celldim*i + j;
+ refverts[n] = vertices[n];
+ globverts[dim*i+j] = vertices[n];
+ }
+ }
+
+ return;
+}
+
+
+void Cell::set_global_vertices(double *vertices)
+{
+ // Copy reference?
+ //globverts = vertices;
+
+ // Copy data instead
+ int i,j;
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ const int dim = n_dim();
+ for (i = 0; i < nno; i++)
+ {
+ for (j = 0; j < dim; j++)
+ {
+ int n = dim * i + j;
+ globverts[n] = vertices[n];
+ }
+ }
+
+ return;
+}
+
+
+double Cell::jacobian(double *point, double jac[3][3])
+{
+ const int celldim = n_celldim();
+ switch (celldim)
+ {
+ case 3: return jacobian(point[0], point[1], point[2], jac);
+ case 2: return jacobian(point[0], point[1], 0.0, jac);
+ case 1: return jacobian(point[0], 0.0, 0.0, jac);
+ }
+ assert(false);
+ return 0.0;
+}
+
+
+double Cell::jacobian(double u, double v, double w, double jac[3][3])
+{
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.0;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.0;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.0;
+
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+
+ double uvw[3] = {u,v,w};
+ double grad[nno*3];
+ double *s;
+
+ #define X(i) globverts[3*(i) + 0]
+ #define Y(i) globverts[3*(i) + 1]
+ #define Z(i) globverts[3*(i) + 2]
+
+ switch (celldim)
+ {
+ /* Cell dimension is 3 */
+ case 3:
+
+ grad_shape(1, uvw, grad);
+
+ for(int i = 0; i < n_nodes(); i++)
+ {
+ s = &grad[3*i];
+
+ jac[0][0] += X(i) * s[0];
+ jac[0][1] += Y(i) * s[0];
+ jac[0][2] += Z(i) * s[0];
+
+ jac[1][0] += X(i) * s[1];
+ jac[1][1] += Y(i) * s[1];
+ jac[1][2] += Z(i) * s[1];
+
+ jac[2][0] += X(i) * s[2];
+ jac[2][1] += Y(i) * s[2];
+ jac[2][2] += Z(i) * s[2];
+ }
+ return fabs(
+ + jac[0][0] * jac[1][1] * jac[2][2]
+ + jac[0][2] * jac[1][0] * jac[2][1]
+ + jac[0][1] * jac[1][2] * jac[2][0]
+ - jac[0][2] * jac[1][1] * jac[2][0]
+ - jac[0][0] * jac[1][2] * jac[2][1]
+ - jac[0][1] * jac[1][0] * jac[2][2]);
+
+ /* Cell dimension is 2 */
+ case 2 :
+
+ grad_shape(1, uvw, grad);
+
+ for(int i = 0; i < n_nodes(); i++)
+ {
+ s = &grad[2*i];
+
+ jac[0][0] += X(i) * s[0];
+ jac[0][1] += Y(i) * s[0];
+ jac[0][2] += Z(i) * s[0];
+
+ jac[1][0] += X(i) * s[1];
+ jac[1][1] += Y(i) * s[1];
+ jac[1][2] += Z(i) * s[1];
+ }
+ {
+ double a[3], b[3], c[3];
+
+ a[0]= X(1) - X(0);
+ a[1]= Y(1) - Y(0);
+ a[2]= Z(1) - Z(0);
+
+ b[0]= X(2) - X(0);
+ b[1]= Y(2) - Y(0);
+ b[2]= Z(2) - Z(0);
+
+ prodve(a, b, c);
+
+ jac[2][0] = c[0];
+ jac[2][1] = c[1];
+ jac[2][2] = c[2];
+ }
+ return sqrt(SQR(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
+ SQR(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
+ SQR(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));
+
+ /* Cell dimension is 1 */
+ case 1:
+
+ grad_shape(1, uvw, grad);
+
+ for(int i = 0; i < n_nodes(); i++)
+ {
+ s = &grad[i];
+
+ jac[0][0] += X(i) * s[0];
+ jac[0][1] += Y(i) * s[0];
+ jac[0][2] += Z(i) * s[0];
+ }
+ {
+ double a[3], b[3], c[3];
+ a[0]= X(1) - X(0);
+ a[1]= Y(1) - Y(0);
+ a[2]= Z(1) - Z(0);
+
+ if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
+ (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2])))
+ {
+ b[0] = a[1];
+ b[1] = -a[0];
+ b[2] = 0.0;
+ }
+ else
+ {
+ b[0] = 0.0;
+ b[1] = a[2];
+ b[2] = -a[1];
+ }
+
+ prodve(a, b, c);
+
+ jac[1][0] = b[0];
+ jac[1][1] = b[1];
+ jac[1][2] = b[2];
+
+ jac[2][0] = c[0];
+ jac[2][1] = c[1];
+ jac[2][2] = c[2];
+ }
+ return sqrt(SQR(jac[0][0])+SQR(jac[0][1])+SQR(jac[0][2]));
+ }
+ #undef X
+ #undef Y
+ #undef Z
+
+ return 0.0;
+}
+
+
+/*
+ * \tilde{\phi}(x) = \sum_{i} d_i N_i(\xi(x))
+ * @param valdim dimension of value points
+ * @param dofs dofs is a rank-2 array of dimension [ndofs x valdim]
+ * @param point point is a rank-1 array of dimension [celldim]
+ * @param value value is a rank-1 array of dimension [valdim]
+ */
+void Cell::interpolate(double *dofs, double *point, double *value, int valdim)
+{
+ const int nno = n_nodes();
+ double N[nno];
+
+ shape(1, point, N);
+
+ //
+ // value[i] = N[0]*dofs[0,i] + N[1]*dofs[1,i] + ...
+ //
+ for (int i = 0; i < valdim; i++)
+ {
+ double sum = 0.0;
+ for (int j = 0; j < nno; j++)
+ sum += dofs[i + valdim*j] * N[j];
+ value[i] = sum;
+ }
+}
+
+
+/*
+ * @param dofs is a rank-2 array of dimension [ndofs x valdim]
+ * @param point is a rank-1 array of dimension [celldim]
+ * @param value is a rank-1 array of dimension [valdim]
+ */
+void Cell::interpolate_grad(double *dofs, double *point, double *value, int stride, double invjac[3][3])
+{
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+
+ //assert(nno > 0);
+ //assert(celldim > 0);
+ //assert(ndofs > 0); <-- ndofs is (cell->nno)
+
+ int i, j, k;
+
+ double dfdu[3] = {0.0, 0.0, 0.0};
+ double grad[nno*3];
+ double *s;
+
+ grad_shape(1, point, grad);
+
+ k = 0;
+ for (i = 0; i < nno; i++)
+ {
+ s = &grad[3*i];
+ for (j = 0; j < celldim; j++)
+ {
+ dfdu[j] += dofs[k] * s[j];
+ }
+ k += stride;
+ }
+
+ if (invjac)
+ {
+ matvec(invjac, dfdu, value);
+ }
+ else
+ {
+ double jac[3][3], inv[3][3];
+ jacobian(point, jac);
+ inv3x3(jac, inv);
+ matvec(inv, dfdu, value);
+ }
+}
+
+
+/* Inverse reference map */
+void cigma::Cell::xyz2uvw(double xyz[3], double uvw[3])
+{
+ #define X(i) globverts[3*(i) + 0]
+ #define Y(i) globverts[3*(i) + 1]
+ #define Z(i) globverts[3*(i) + 2]
+
+ // general newton routine for the nonlinear case
+ uvw[0] = uvw[1] = uvw[2] = 0.0;
+
+ const int nno = n_nodes();
+ int iter = 1, maxiter = 20;
+ double error = 1.0, tol = 1.0e-6;
+
+ while ((error > tol) && (iter < maxiter))
+ {
+ double jac[3][3];
+ if (!jacobian(uvw[0], uvw[1], uvw[2], jac))
+ break;
+
+ double s[nno];
+ shape(1, uvw, s);
+
+ double xn = 0.0, yn = 0.0, zn = 0.0;
+ for (int i = 0; i < nno; i++)
+ {
+ xn += X(i) * s[i];
+ yn += Y(i) * s[i];
+ zn += Z(i) * s[i];
+ }
+
+ double inv[3][3];
+ inv3x3(jac, inv);
+
+ double un = uvw[0] + inv[0][0] * (xyz[0] - xn)
+ + inv[1][0] * (xyz[1] - yn)
+ + inv[2][0] * (xyz[2] - zn);
+
+ double vn = uvw[1] + inv[0][1] * (xyz[0] - xn)
+ + inv[1][1] * (xyz[1] - yn)
+ + inv[2][1] * (xyz[2] - zn);
+
+ double wn = uvw[2] + inv[0][2] * (xyz[0] - xn)
+ + inv[1][2] * (xyz[1] - yn)
+ + inv[2][2] * (xyz[2] - zn);
+
+ error = sqrt(SQR(un - uvw[0]) +
+ SQR(vn - uvw[1]) +
+ SQR(wn - uvw[2]));
+
+ uvw[0] = un;
+ uvw[1] = vn;
+ uvw[2] = wn;
+
+ iter++;
+ }
+
+ /*
+ if (error > tol)
+ {
+ std::cerr << "Warning: Newton did not converge in xyz2uvw" << std::endl;
+ } // */
+
+ #undef X
+ #undef Y
+ #undef Z
+}
+
+
+/* isoparametric reference map */
+void cigma::Cell::uvw2xyz(double uvw[3], double xyz[3])
+{
+ assert(globverts != 0);
+ interpolate(globverts, uvw, xyz, 3);
+}
+
+
+void cigma::Cell::bbox(double *min, double *max)
+{
+ const int nno = n_nodes();
+ const int nsd = n_dim();
+ minmax(globverts, nno, nsd, min, max);
+}
+
+
+void cigma::Cell::centroid(double c[3])
+{
+ const int nno = n_nodes();
+ const int nsd = n_dim();
+ cigma::centroid(globverts, nno, nsd, c);
+}
+
+
+bool cigma::Cell::global_interior(double xyz[3])
+{
+ double uvw[3];
+ xyz2uvw(xyz,uvw);
+ return interior(uvw[0], uvw[1], uvw[2]);
+}
+
+
Modified: cs/cigma/trunk/src/Cell.h
===================================================================
--- cs/cigma/trunk/src/Cell.h 2008-10-15 09:08:08 UTC (rev 13063)
+++ cs/cigma/trunk/src/Cell.h 2008-10-15 09:08:11 UTC (rev 13064)
@@ -1,6 +1,9 @@
#ifndef __CIGMA_CELL_H__
#define __CIGMA_CELL_H__
+#include "Quadrature.h"
+#include <boost/shared_ptr.hpp>
+
namespace cigma
{
class Cell;
@@ -10,7 +13,50 @@
{
public:
Cell();
- ~Cell();
+ virtual ~Cell();
+
+ virtual int n_nodes() const = 0;
+ virtual int n_celldim() const = 0;
+ virtual int n_dim() const = 0;
+
+ virtual void shape(int num, double *points, double *values) = 0;
+ virtual void grad_shape(int num, double *points, double *values) = 0;
+
+ double jacobian(double *point, double jac[3][3]);
+ double jacobian(double u, double v, double w, double jac[3][3]);
+
+ void interpolate(double *dofs, double *point, double *value, int valdim);
+ void interpolate_grad(double *dofs, double *point, double *value, int stride=1, double invjac[3][3]=0);
+
+ virtual void xyz2uvw(double xyz[3], double uvw[3]);
+ void uvw2xyz(double uvw[3], double xyz[3]);
+
+public:
+ typedef enum {
+ POINT = 0,
+ TRIANGLE,
+ QUADRANGLE,
+ TETRAHEDRON,
+ HEXAHEDRON
+ } Geometry;
+
+ virtual Geometry geometry() = 0;
+ virtual boost::shared_ptr<Quadrature> default_quadrature() = 0;
+
+ void set_reference_vertices(double *vertices);
+ void set_global_vertices(double *vertices);
+
+ virtual double volume() = 0;
+ virtual bool interior(double u, double v, double w) = 0;
+ virtual bool global_interior(double xyz[3]);
+
+ void bbox(double *min, double *max);
+ void centroid(double c[3]);
+
+public:
+ double *refverts; // [nno x celldim]
+ double *globverts; // [nno x nsd]
+
};
#endif
Added: cs/cigma/trunk/src/fe_hex8.cpp
===================================================================
--- cs/cigma/trunk/src/fe_hex8.cpp (rev 0)
+++ cs/cigma/trunk/src/fe_hex8.cpp 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,373 @@
+//#include <iostream>
+#include <cmath>
+#include "fe_hex8.h"
+
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+static void hex_shape(double u, double v, double w, double s[8])
+{
+ s[0] = (1.0 - u) * (1.0 - v) * (1.0 - w) * 0.125;
+ s[1] = (1.0 + u) * (1.0 - v) * (1.0 - w) * 0.125;
+ s[2] = (1.0 + u) * (1.0 + v) * (1.0 - w) * 0.125;
+ s[3] = (1.0 - u) * (1.0 + v) * (1.0 - w) * 0.125;
+ s[4] = (1.0 - u) * (1.0 - v) * (1.0 + w) * 0.125;
+ s[5] = (1.0 + u) * (1.0 - v) * (1.0 + w) * 0.125;
+ s[6] = (1.0 + u) * (1.0 + v) * (1.0 + w) * 0.125;
+ s[7] = (1.0 - u) * (1.0 + v) * (1.0 + w) * 0.125;
+}
+
+static void hex_grad_shape(double u, double v, double w, double s[8*3])
+{
+ #define S(i,j) s[3*(i) + (j)]
+
+ S(0,0) = -0.125 * (1.0 - v) * (1.0 - w);
+ S(0,1) = -0.125 * (1.0 - u) * (1.0 - w);
+ S(0,2) = -0.125 * (1.0 - u) * (1.0 - v);
+
+ S(1,0) = 0.125 * (1.0 - v) * (1.0 - w);
+ S(1,1) = -0.125 * (1.0 + u) * (1.0 - w);
+ S(1,2) = -0.125 * (1.0 + u) * (1.0 - v);
+
+ S(2,0) = 0.125 * (1.0 + v) * (1.0 - w);
+ S(2,1) = 0.125 * (1.0 + u) * (1.0 - w);
+ S(2,2) = -0.125 * (1.0 + u) * (1.0 + v);
+
+ S(3,0) = -0.125 * (1.0 + v) * (1.0 - w);
+ S(3,1) = 0.125 * (1.0 - u) * (1.0 - w);
+ S(3,2) = -0.125 * (1.0 - u) * (1.0 + v);
+
+ S(4,0) = -0.125 * (1.0 - v) * (1.0 + w);
+ S(4,1) = -0.125 * (1.0 - u) * (1.0 + w);
+ S(4,2) = 0.125 * (1.0 - u) * (1.0 - v);
+
+ S(5,0) = 0.125 * (1.0 - v) * (1.0 + w);
+ S(5,1) = -0.125 * (1.0 + u) * (1.0 + w);
+ S(5,2) = 0.125 * (1.0 + u) * (1.0 - v);
+
+ S(6,0) = 0.125 * (1.0 + v) * (1.0 + w);
+ S(6,1) = 0.125 * (1.0 + u) * (1.0 + w);
+ S(6,2) = 0.125 * (1.0 + u) * (1.0 + v);
+
+ S(7,0) = -0.125 * (1.0 + v) * (1.0 + w);
+ S(7,1) = 0.125 * (1.0 - u) * (1.0 + w);
+ S(7,2) = 0.125 * (1.0 - u) * (1.0 + v);
+
+ #undef S
+}
+
+
+// ---------------------------------------------------------------------------
+
+hex8::hex8()
+{
+ const int nno = 8;
+ const int celldim = 3;
+ double verts[nno*celldim] = {
+ -1.0, -1.0, -1.0,
+ +1.0, -1.0, -1.0,
+ +1.0, +1.0, -1.0,
+ -1.0, +1.0, -1.0,
+ -1.0, -1.0, +1.0,
+ +1.0, -1.0, +1.0,
+ +1.0, +1.0, +1.0,
+ -1.0, +1.0, +1.0
+ };
+ set_reference_vertices(verts);
+}
+
+hex8::~hex8()
+{
+}
+
+
+boost::shared_ptr<Quadrature> hex8::default_quadrature()
+{
+ // hex8_qr(3)
+ const int nno = 8;
+ const int celldim = 3;
+ double qpts[nno * celldim] = {
+ -0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, 0.57735027, 0.57735027,
+ -0.57735027, 0.57735027, 0.57735027
+ };
+ double qwts[nno] = { 1., 1., 1., 1., 1., 1., 1., 1. };
+
+ boost::shared_ptr<Quadrature> Q(new Quadrature());
+ Q->setData(nno, celldim, qpts, qwts);
+ return Q;
+}
+
+/*
+void hex8::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // hex_qr(3)
+ const int hex_nno = 8;
+ const int hex_celldim = 3;
+ double hex_qpts[hex_nno * hex_celldim] = {
+ -0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, 0.57735027, 0.57735027,
+ -0.57735027, 0.57735027, 0.57735027
+ };
+ double hex_qwts[hex_nno] = { 1., 1., 1., 1., 1., 1., 1., 1. };
+
+ int i,j,n;
+ *npts = hex_nno;
+ *ndim = hex_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < hex_nno; i++)
+ {
+ (*wts)[i] = hex_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = hex_qpts[n];
+ }
+ }
+}
+*/
+
+
+// ---------------------------------------------------------------------------
+
+/*
+ * @param points points is an [num x celldim] array (in)
+ * @param values values is an [num x ndofs] array (out)
+ *
+ */
+void hex8::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];
+ hex_shape(u, v, w, &values[nno*i]);
+ }
+}
+
+/*
+ *
+ * @param points points is an [num x celldim] array (in)
+ * @param values values is an [num x ndofs x celldim] array (out)
+ */
+void hex8::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];
+ hex_grad_shape(u, v, w, &values[stride*i]);
+ }
+}
+
+/*
+ *
+ */
+#define ZERO (-1.0e-6)
+#define ONE (1-ZERO)
+
+bool hex8::interior(double u, double v, double w)
+{
+
+ if ((u < -ONE) || (v < -ONE) || (w < -ONE) ||
+ (u > +ONE) || (v > +ONE) || (w > +ONE))
+ {
+ return false;
+ }
+ return true;
+
+}
+
+#undef ONE
+#undef ZERO
+
+
+double hex8::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
+}
+
Added: cs/cigma/trunk/src/fe_hex8.h
===================================================================
--- cs/cigma/trunk/src/fe_hex8.h (rev 0)
+++ cs/cigma/trunk/src/fe_hex8.h 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,33 @@
+#ifndef __FE_HEX8_H__
+#define __FE_HEX8_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+ class hex8;
+}
+
+class cigma::hex8 : public Cell
+{
+public:
+ hex8();
+ ~hex8();
+
+public:
+ int n_nodes() const { return 8; }
+ int n_celldim() const { return 3; }
+ int n_dim() const { return 3; }
+ Geometry geometry() { return HEXAHEDRON; }
+ boost::shared_ptr<Quadrature> default_quadrature();
+
+public:
+ 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
Added: cs/cigma/trunk/src/fe_quad4.cpp
===================================================================
--- cs/cigma/trunk/src/fe_quad4.cpp (rev 0)
+++ cs/cigma/trunk/src/fe_quad4.cpp 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,202 @@
+#include "fe_quad4.h"
+
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+static void quad_shape(double u, double v, double s[4])
+{
+ s[0] = 0.25 * (1.0 - u) * (1.0 - v);
+ s[1] = 0.25 * (1.0 + u) * (1.0 - v);
+ s[2] = 0.25 * (1.0 + u) * (1.0 + v);
+ s[3] = 0.25 * (1.0 - u) * (1.0 + v);
+}
+
+static void quad_grad_shape(double u, double v, double s[4*2])
+{
+ #define S(i,j) s[2*(i) + (j)]
+
+ S(0,0) = -0.25 * (1.0 - v);
+ S(0,1) = -0.25 * (1.0 - u);
+
+ S(1,0) = +0.25 * (1.0 - v);
+ S(1,1) = -0.25 * (1.0 + u);
+
+ S(2,0) = +0.25 * (1.0 + v);
+ S(2,1) = +0.25 * (1.0 + u);
+
+ S(3,0) = -0.25 * (1.0 + v);
+ S(3,1) = +0.25 * (1.0 - u);
+
+ #undef S
+}
+
+// ---------------------------------------------------------------------------
+
+quad4::quad4()
+{
+ const int nno = 4;
+ const int celldim = 2;
+ double verts[nno*celldim] = {
+ -1.0, -1.0,
+ +1.0, -1.0,
+ +1.0, +1.0,
+ -1.0, +1.0
+ };
+ set_reference_vertices(verts);
+}
+
+quad4::~quad4()
+{
+}
+
+boost::shared_ptr<Quadrature> quad4::default_quadrature()
+{
+ // quad4_qr(7)
+ const int nno = 16;
+ const int celldim = 2;
+ double qpts[nno*celldim] = {
+ -0.86113631, -0.86113631,
+ -0.33998104, -0.86113631,
+ 0.33998104, -0.86113631,
+ 0.86113631, -0.86113631,
+ 0.86113631, -0.33998104,
+ 0.86113631, 0.33998104,
+ 0.86113631, 0.86113631,
+ 0.33998104, 0.86113631,
+ -0.33998104, 0.86113631,
+ -0.86113631, 0.86113631,
+ -0.86113631, 0.33998104,
+ -0.86113631, -0.33998104,
+ -0.33998104, -0.33998104,
+ 0.33998104, -0.33998104,
+ -0.33998104, 0.33998104,
+ 0.33998104, 0.33998104
+ };
+ double qwts[nno] = {
+ 0.12100299, 0.22685185, 0.22685185, 0.12100299, 0.22685185,
+ 0.22685185, 0.12100299, 0.22685185, 0.22685185, 0.12100299,
+ 0.22685185, 0.22685185, 0.4252933 , 0.4252933 , 0.4252933 ,
+ 0.4252933
+ };
+
+ boost::shared_ptr<Quadrature> Q(new Quadrature());
+ Q->setData(nno, celldim, qpts, qwts);
+ return Q;
+}
+
+/*
+void quad4::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // quad_qr(7)
+ const int quad_nno = 16;
+ const int quad_celldim = 2;
+ double quad_qpts[quad_nno * quad_celldim] = {
+ -0.86113631, -0.86113631,
+ -0.33998104, -0.86113631,
+ 0.33998104, -0.86113631,
+ 0.86113631, -0.86113631,
+ 0.86113631, -0.33998104,
+ 0.86113631, 0.33998104,
+ 0.86113631, 0.86113631,
+ 0.33998104, 0.86113631,
+ -0.33998104, 0.86113631,
+ -0.86113631, 0.86113631,
+ -0.86113631, 0.33998104,
+ -0.86113631, -0.33998104,
+ -0.33998104, -0.33998104,
+ 0.33998104, -0.33998104,
+ -0.33998104, 0.33998104,
+ 0.33998104, 0.33998104
+ };
+ double quad_qwts[quad_nno] = {
+ 0.12100299, 0.22685185, 0.22685185, 0.12100299, 0.22685185,
+ 0.22685185, 0.12100299, 0.22685185, 0.22685185, 0.12100299,
+ 0.22685185, 0.22685185, 0.4252933 , 0.4252933 , 0.4252933 ,
+ 0.4252933
+ };
+
+ int i,j,n;
+ *npts = quad_nno;
+ *ndim = quad_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < quad_nno; i++)
+ {
+ (*wts)[i] = quad_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = quad_qpts[n];
+ }
+ }
+}
+*/
+
+// ---------------------------------------------------------------------------
+
+/*
+ *
+ * @param points points is an [num x celldim] array (in)
+ * @param values values is an [num x ndofs] array (out)
+ */
+void quad4::shape(int num, double *points, double *values)
+{
+ const int nno = n_nodes();
+ for (int i = 0; i < num; i++)
+ {
+ double u = points[2*i + 0];
+ double v = points[2*i + 1];
+ quad_shape(u, v, &values[nno*i]);
+ }
+}
+
+/*
+ *
+ * @param points points is an [num x celldim] array (in)
+ * @param values values is an [num x ndofs x celldim] array (out)
+ */
+void quad4::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[2*i + 0];
+ double v = points[2*i + 1];
+ quad_grad_shape(u, v, &values[stride*i]);
+ }
+}
+
+
+/*
+ *
+ *
+ */
+#define ZERO (-1.0e-6)
+#define ONE (1-ZERO)
+
+bool quad4::interior(double u, double v, double w)
+{
+ if ((u < -ONE) || (v < -ONE) ||
+ (u > +ONE) || (v > +ONE))
+ {
+ return false;
+ }
+ return true;
+}
+
+#undef ONE
+#undef ZERO
+
+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));
+ #undef X
+ #undef Y
+}
+
+// ---------------------------------------------------------------------------
Added: cs/cigma/trunk/src/fe_quad4.h
===================================================================
--- cs/cigma/trunk/src/fe_quad4.h (rev 0)
+++ cs/cigma/trunk/src/fe_quad4.h 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,31 @@
+#ifndef __FE_QUAD4_H__
+#define __FE_QUAD4_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+ class quad4;
+}
+
+class cigma::quad4 : public Cell
+{
+public:
+ quad4();
+ ~quad4();
+
+public:
+ int n_nodes() const { return 4; }
+ int n_celldim() const { return 2; }
+ int n_dim() const { return 3; }
+ Geometry geometry() { return QUADRANGLE; }
+ boost::shared_ptr<Quadrature> default_quadrature();
+
+public:
+ 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
Added: cs/cigma/trunk/src/fe_tet4.cpp
===================================================================
--- cs/cigma/trunk/src/fe_tet4.cpp (rev 0)
+++ cs/cigma/trunk/src/fe_tet4.cpp 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,271 @@
+#include <cmath>
+#include "fe_tet4.h"
+#include "Numeric.h"
+
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+static void tet_shape(double u, double v, double w, double s[4])
+{
+ s[0] = 1.0 - u - v - w;
+ s[1] = u ;
+ s[2] = v ;
+ s[3] = w;
+}
+
+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) = -1.0;
+ S(0,1) = -1.0;
+ S(0,2) = -1.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;
+
+ S(3,0) = 0.0;
+ S(3,1) = 0.0;
+ S(3,2) = +1.0;
+
+ #undef S
+}
+
+
+// ---------------------------------------------------------------------------
+
+tet4::tet4()
+{
+ const int nno = 4;
+ const int celldim = 3;
+ double verts[nno*celldim] = {
+ 0.0, 0.0, 0.0,
+ 1.0, 0.0, 0.0,
+ 0.0, 1.0, 0.0,
+ 0.0, 0.0, 1.0
+ };
+ set_reference_vertices(verts);
+}
+
+tet4::~tet4()
+{
+}
+
+boost::shared_ptr<Quadrature> tet4::default_quadrature()
+{
+ // tet4_qr(3)
+ const int nno = 8;
+ const int celldim = 3;
+ double qpts[nno * celldim] = {
+ -0.68663473, -0.72789005, -0.75497035,
+ -0.83720867, -0.85864055, 0.08830369,
+ -0.86832263, 0.13186633, -0.75497035,
+ -0.93159441, -0.4120024 , 0.08830369,
+ 0.16949513, -0.72789005, -0.75497035,
+ -0.39245447, -0.85864055, 0.08830369,
+ -0.50857335, 0.13186633, -0.75497035,
+ -0.74470688, -0.4120024 , 0.08830369
+ };
+ double qwts[nno] = {
+ 0.29583885, 0.12821632, 0.16925605, 0.07335544, 0.29583885,
+ 0.12821632, 0.16925605, 0.07335544
+ };
+
+ for (int i = 0; i < nno; i++)
+ {
+ for (int j = 0; j < celldim; j++)
+ {
+ int n = celldim * i + j;
+ qpts[n] += 1;
+ qpts[n] *= 0.5;
+ }
+ // don't forget to adjust the integration weights
+ // to account for the reference domain transformation!
+ qwts[i] *= 0.125;
+ }
+
+ boost::shared_ptr<Quadrature> Q(new Quadrature());
+ Q->setData(nno, celldim, qpts, qwts);
+ return Q;
+}
+
+/*
+void tet4::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // tet_qr(3)
+ const int tet_nno = 8;
+ const int tet_celldim = 3;
+ double tet_qpts[tet_nno * tet_celldim] = {
+ -0.68663473, -0.72789005, -0.75497035,
+ -0.83720867, -0.85864055, 0.08830369,
+ -0.86832263, 0.13186633, -0.75497035,
+ -0.93159441, -0.4120024 , 0.08830369,
+ 0.16949513, -0.72789005, -0.75497035,
+ -0.39245447, -0.85864055, 0.08830369,
+ -0.50857335, 0.13186633, -0.75497035,
+ -0.74470688, -0.4120024 , 0.08830369
+ };
+ double tet_qwts[tet_nno] = {
+ 0.29583885, 0.12821632, 0.16925605, 0.07335544, 0.29583885,
+ 0.12821632, 0.16925605, 0.07335544
+ };
+
+ int i,j,n;
+
+ for (i = 0; i < tet_nno; i++)
+ {
+ for (j = 0; j < tet_celldim; j++)
+ {
+ n = tet_celldim * i + j;
+ tet_qpts[n] += 1;
+ tet_qpts[n] *= 0.5;
+ }
+ // don't forget to adjust the integration weights
+ // to account for the reference domain transformation!
+ tet_qwts[i] *= 0.125;
+ }
+
+ *npts = tet_nno;
+ *ndim = tet_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < tet_nno; i++)
+ {
+ (*wts)[i] = tet_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = tet_qpts[n];
+ }
+ }
+}
+*/
+
+void tet4::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];
+ tet_shape(u, v, w, &values[nno*i]);
+ }
+}
+
+void tet4::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];
+ tet_grad_shape(u, v, w, &values[stride*i]);
+ }
+}
+
+void tet4::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 tet4::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 tet4::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_tet4.h
===================================================================
--- cs/cigma/trunk/src/fe_tet4.h (rev 0)
+++ cs/cigma/trunk/src/fe_tet4.h 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,34 @@
+#ifndef __FE_TET4_H__
+#define __FE_TET4_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+ class tet4;
+}
+
+class cigma::tet4 : public Cell
+{
+public:
+ tet4();
+ ~tet4();
+ void qr_default(double **wts, double **pts, int *npts, int *ndim);
+
+public:
+ int n_nodes() const { return 4; }
+ int n_celldim() const { return 3; }
+ int n_dim() const { return 3; }
+ Geometry geometry() { return TETRAHEDRON; }
+ boost::shared_ptr<Quadrature> default_quadrature();
+
+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);
+ double volume();
+};
+
+
+#endif
Added: cs/cigma/trunk/src/fe_tri3.cpp
===================================================================
--- cs/cigma/trunk/src/fe_tri3.cpp (rev 0)
+++ cs/cigma/trunk/src/fe_tri3.cpp 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,182 @@
+#include "fe_tri3.h"
+
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+static void tri_shape(double u, double v, 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 s[3*2])
+{
+ #define S(i,j) s[2*(i) + (j)]
+
+ S(0,0) = -1.0;
+ S(0,1) = -1.0;
+
+ S(1,0) = +1.0;
+ S(1,1) = 0.0;
+
+ S(2,0) = 0.0;
+ S(2,1) = +1.0;
+
+ #undef S
+}
+
+
+// ---------------------------------------------------------------------------
+
+tri3::tri3()
+{
+ const int nno = 3;
+ const int celldim = 2;
+ double verts[nno*celldim] = {
+ 0.0, 0.0,
+ 1.0, 0.0,
+ 0.0, 1.0
+ };
+ set_reference_vertices(verts);
+}
+
+tri3::~tri3()
+{
+}
+
+boost::shared_ptr<Quadrature> tri3::default_quadrature()
+{
+ // tri3_qr(5)
+ const int nno = 9;
+ const int celldim = 2;
+ double qpts[nno*celldim] = {
+ -0.79456469, -0.82282408,
+ -0.86689186, -0.18106627,
+ -0.95213774, 0.57531892,
+ -0.08858796, -0.82282408,
+ -0.40946686, -0.18106627,
+ -0.78765946, 0.57531892,
+ 0.61738877, -0.82282408,
+ 0.04795814, -0.18106627,
+ -0.62318119, 0.57531892
+ };
+ double qwts[nno] = {
+ 0.22325768, 0.25471234, 0.07758553, 0.35721229, 0.40753974,
+ 0.12413685, 0.22325768, 0.25471234, 0.07758553
+ };
+
+ boost::shared_ptr<Quadrature> Q(new Quadrature());
+ Q->setData(nno, celldim, qpts, qwts);
+ return Q;
+}
+
+/*
+void tri3::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // tri_qr(5)
+ const int tri_nno = 9;
+ const int tri_celldim = 2;
+ double tri_qpts[tri_nno * tri_celldim] = {
+ -0.79456469, -0.82282408,
+ -0.86689186, -0.18106627,
+ -0.95213774, 0.57531892,
+ -0.08858796, -0.82282408,
+ -0.40946686, -0.18106627,
+ -0.78765946, 0.57531892,
+ 0.61738877, -0.82282408,
+ 0.04795814, -0.18106627,
+ -0.62318119, 0.57531892
+ };
+ double tri_qwts[tri_nno] = {
+ 0.22325768, 0.25471234, 0.07758553, 0.35721229, 0.40753974,
+ 0.12413685, 0.22325768, 0.25471234, 0.07758553
+ };
+
+ int i,j,n;
+
+ for (i = 0; i < tri_nno; i++)
+ {
+ for (j = 0; j < tri_celldim; j++)
+ {
+ n = tri_celldim * i + j;
+ tri_qpts[n] += 1;
+ tri_qpts[n] *= 0.5;
+ }
+ }
+
+ *npts = tri_nno;
+ *ndim = tri_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < tri_nno; i++)
+ {
+ (*wts)[i] = tri_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = tri_qpts[n];
+ }
+ }
+}
+*/
+
+
+// ---------------------------------------------------------------------------
+
+void tri3::shape(int num, double *points, double *values)
+{
+ const int nno = n_nodes();
+ for (int i = 0; i < num; i++)
+ {
+ double u = points[2*i + 0];
+ double v = points[2*i + 1];
+ tri_shape(u, v, &values[nno*i]);
+ }
+}
+
+void tri3::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[2*i + 0];
+ double v = points[2*i + 1];
+ tri_grad_shape(u, v, &values[stride*i]);
+ }
+}
+
+
+//void tri3::xyz2uvw(double xyz[3], double uvw[3]) {}
+
+
+#define ZERO (-1.0e-6)
+#define ONE (1 - ZERO)
+
+bool tri3::interior(double u, double v, double w)
+{
+
+ if ((u < ZERO) || (v < ZERO) || (u > (ONE - v)))
+ {
+ return false;
+ }
+ return true;
+}
+
+#undef ZERO
+#undef ONE
+
+
+double tri3::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
+}
+
+// ---------------------------------------------------------------------------
Added: cs/cigma/trunk/src/fe_tri3.h
===================================================================
--- cs/cigma/trunk/src/fe_tri3.h (rev 0)
+++ cs/cigma/trunk/src/fe_tri3.h 2008-10-15 09:08:11 UTC (rev 13064)
@@ -0,0 +1,32 @@
+#ifndef __FE_TRI3_H__
+#define __FE_TRI3_H__
+
+#include "Cell.h"
+
+namespace cigma
+{
+ class tri3;
+}
+
+class cigma::tri3 : public Cell
+{
+public:
+ tri3();
+ ~tri3();
+
+public:
+ int n_nodes() const { return 3; }
+ int n_celldim() const { return 2; }
+ int n_dim() const { return 3; }
+ Geometry geometry() { return TRIANGLE; }
+ boost::shared_ptr<Quadrature> default_quadrature();
+
+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);
+ double volume();
+};
+
+#endif
More information about the cig-commits
mailing list