diff --git a/Common/ShapeFunctions.h b/Common/ShapeFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..6f41b9400519cb4f0994504154c21a97c04bfccc
--- /dev/null
+++ b/Common/ShapeFunctions.h
@@ -0,0 +1,943 @@
+#ifndef _SHAPE_FUNCTIONS_H_
+#define _SHAPE_FUNCTIONS_H_
+
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include "Numeric.h"
+#include "Message.h"
+
+#define ONE (1. + 1.e-6)
+#define ZERO (-1.e-6)
+
+class element{
+protected:
+ double *_x, *_y, *_z;
+public:
+ element(double *x, double *y, double *z) : _x(x), _y(y), _z(z) {}
+ virtual ~element(){}
+ virtual int getDimension() = 0;
+ virtual int getNumNodes() = 0;
+ virtual void getNode(int num, double &u, double &v, double &w) = 0;
+ virtual int getNumEdges() = 0;
+ virtual void getEdge(int num, int &start, int &end) = 0;
+ virtual int getNumGaussPoints() = 0;
+ virtual void getGaussPoint(int num, double &u, double &v, double &w, double &weight) = 0;
+ virtual void getShapeFunction(int num, double u, double v, double w, double &s) = 0;
+ virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3]) = 0;
+ double getJacobian(double u, double v, double w, double jac[3][3])
+ {
+ jac[0][0] = jac[0][1] = jac[0][2] = 0.;
+ jac[1][0] = jac[1][1] = jac[1][2] = 0.;
+ jac[2][0] = jac[2][1] = jac[2][2] = 0.;
+ double s[3];
+ switch(getDimension()){
+ case 3 :
+ for(int i = 0; i < getNumNodes(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ 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]);
+ case 2 :
+ for(int i = 0; i < getNumNodes(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ 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]));
+ case 1:
+ for(int i = 0; i < getNumNodes(); i++) {
+ getGradShapeFunction(i, u, v, w, s);
+ 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.;
+ }
+ else {
+ b[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]));
+ default:
+ return 1.;
+ }
+ }
+ double interpolate(double val[], double u, double v, double w, int stride=1)
+ {
+ double sum = 0;
+ int j = 0;
+ for(int i = 0; i < getNumNodes(); i++){
+ double s;
+ getShapeFunction(i, u, v, w, s);
+ sum += val[j] * s;
+ j += stride;
+ }
+ return sum;
+ }
+ void interpolateGrad(double val[], double u, double v, double w, double f[3], int stride=1,
+ double invjac[3][3]=NULL)
+ {
+ double dfdu[3] = {0., 0., 0.};
+ int j = 0;
+ for(int i = 0; i < getNumNodes(); i++){
+ double s[3];
+ getGradShapeFunction(i, u, v, w, s);
+ dfdu[0] += val[j] * s[0];
+ dfdu[1] += val[j] * s[1];
+ dfdu[2] += val[j] * s[2];
+ j += stride;
+ }
+ if(invjac){
+ matvec(invjac, dfdu, f);
+ }
+ else{
+ double jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ matvec(inv, dfdu, f);
+ }
+ }
+ void interpolateCurl(double val[], double u, double v, double w, double f[3], int stride=3)
+ {
+ double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ interpolateGrad(&val[0], u, v, w, fx, stride, inv);
+ interpolateGrad(&val[1], u, v, w, fy, stride, inv);
+ interpolateGrad(&val[2], u, v, w, fz, stride, inv);
+ f[0] = fz[1] - fy[2];
+ f[1] = -(fz[0] - fx[2]);
+ f[2] = fy[0] - fx[1];
+ }
+ double interpolateDiv(double val[], double u, double v, double w, int stride=3)
+ {
+ double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
+ getJacobian(u, v, w, jac);
+ inv3x3(jac, inv);
+ interpolateGrad(&val[0], u, v, w, fx, stride, inv);
+ interpolateGrad(&val[1], u, v, w, fy, stride, inv);
+ interpolateGrad(&val[2], u, v, w, fz, stride, inv);
+ return fx[0] + fy[1] + fz[2];
+ }
+ double integrate(double val[], int stride=1)
+ {
+ double sum = 0;
+ for(int i = 0; i < getNumGaussPoints(); i++){
+ double u, v, w, weight, jac[3][3];
+ getGaussPoint(i, u, v, w, weight);
+ double det = getJacobian(u, v, w, jac);
+ double d = interpolate(val, u, v, w, stride);
+ sum += d * weight * det;
+ }
+ return sum;
+ }
+ double integrateLevelsetPositive(double val[])
+ {
+ // FIXME: explain + generalize this
+ double ones[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
+ double area = integrate(ones);
+ double sum = 0, sumabs = 0.;
+ for(int i = 0; i < getNumNodes(); i++){
+ sum += val[i];
+ sumabs += fabs(val[i]);
+ }
+ double res = 0.;
+ if(sumabs)
+ res = area * (1 + sum/sumabs) * 0.5 ;
+ return res;
+ }
+ virtual double integrateCirculation(double val[])
+ {
+ Msg::Error("integrateCirculation not available for this element");
+ return 0.;
+ }
+ virtual double integrateFlux(double val[])
+ {
+ Msg::Error("integrateFlux not available for this element");
+ return 0.;
+ }
+ virtual void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ // general newton routine for the nonlinear case (more efficient
+ // routines are implemented for simplices, where the basis
+ // functions are linear)
+ uvw[0] = uvw[1] = uvw[2] = 0.0;
+
+ int iter = 1, maxiter = 20;
+ double error = 1., tol = 1.e-6;
+
+ while (error > tol && iter < maxiter){
+ double jac[3][3];
+ if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
+
+ double xn = 0., yn = 0., zn = 0.;
+ for (int i = 0; i < getNumNodes(); i++) {
+ double s;
+ getShapeFunction(i, uvw[0], uvw[1], uvw[2], s);
+ xn += _x[i] * s;
+ yn += _y[i] * s;
+ zn += _z[i] * s;
+ }
+ 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) Msg::Warning("Newton did not converge in xyz2uvw") ;
+ }
+ virtual int isInside(double u, double v, double w) = 0;
+ double maxEdgeLength()
+ {
+ double max = 0.;
+ for(int i = 0; i < getNumEdges(); i++){
+ int n1, n2;
+ getEdge(i, n1, n2);
+ double d = sqrt(SQR(_x[n1]-_x[n2]) + SQR(_y[n1]-_y[n2]) + SQR(_z[n1]-_z[n2]));
+ if(d > max) max = d;
+ }
+ return max;
+ }
+};
+
+class point : public element{
+public:
+ point(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 0; }
+ inline int getNumNodes(){ return 1; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ u = v = w = 0.;
+ }
+ inline int getNumEdges(){ return 0; }
+ void getEdge(int num, int &start, int &end)
+ {
+ start = end = 0;
+ }
+ inline int getNumGaussPoints(){ return 1; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ u = v = w = 0.;
+ weight = 1.;
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 1.; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ s[0] = s[1] = s[2] = 0.;
+ }
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ uvw[0] = uvw[1] = uvw[2] = 0.;
+ }
+ int isInside(double u, double v, double w)
+ {
+ return 1;
+ }
+};
+
+class line : public element{
+public:
+ line(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 1; }
+ inline int getNumNodes(){ return 2; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ v = w = 0.;
+ switch(num) {
+ case 0 : u = -1.; break;
+ case 1 : u = 1.; break;
+ default: u = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 1; }
+ void getEdge(int num, int &start, int &end)
+ {
+ start = 0; end = 1;
+ }
+ inline int getNumGaussPoints(){ return 1; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ if(num < 0 || num > 0) return;
+ u = v = w = 0.;
+ weight = 2.;
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 0.5 * (1.-u); break;
+ case 1 : s = 0.5 * (1.+u); break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.5; s[1] = 0.; s[2] = 0.; break;
+ case 1 : s[0] = 0.5; s[1] = 0.; s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ double integrateCirculation(double val[])
+ {
+ double t[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
+ norme(t);
+ double v[3];
+ for(int i = 0; i < 3; i++)
+ v[i] = integrate(&val[i], 3);
+ double d;
+ prosca(t, v, &d);
+ return d;
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || u > ONE)
+ return 0;
+ return 1;
+ }
+};
+
+class triangle : public element{
+public:
+ triangle(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 2; }
+ inline int getNumNodes(){ return 3; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ w = 0.;
+ switch(num) {
+ case 0 : u = 0.; v = 0.; break;
+ case 1 : u = 1.; v = 0.; break;
+ case 2 : u = 0.; v = 1.; break;
+ default: u = 0.; v = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 3; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 1; end = 2; break;
+ case 2 : start = 2; end = 0; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 3; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u3[3] = {0.16666666666666,0.66666666666666,0.16666666666666};
+ double v3[3] = {0.16666666666666,0.16666666666666,0.66666666666666};
+ double p3[3] = {0.16666666666666,0.16666666666666,0.16666666666666};
+ if(num < 0 || num > 2) return;
+ u = u3[num];
+ v = v3[num];
+ w = 0.;
+ weight = p3[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num){
+ case 0 : s = 1.-u-v; break;
+ case 1 : s = u ; break;
+ case 2 : s = v; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -1.; s[1] = -1.; s[2] = 0.; break;
+ case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
+ case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ double integrateFlux(double val[])
+ {
+ double t1[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
+ double t2[3] = {_x[2]-_x[0], _y[2]-_y[0], _z[2]-_z[0]};
+ double n[3];
+ prodve(t1, t2, n);
+ norme(n);
+ double v[3];
+ for(int i = 0; i < 3; i++)
+ v[i] = integrate(&val[i], 3);
+ double d;
+ prosca(n, v, &d);
+ return d;
+ }
+#if 0 // faster, but only valid for triangles in the z=0 plane
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ double mat[2][2], b[2];
+ mat[0][0] = _x[1] - _x[0];
+ mat[0][1] = _x[2] - _x[0];
+ mat[1][0] = _y[1] - _y[0];
+ mat[1][1] = _y[2] - _y[0];
+ b[0] = xyz[0] - _x[0];
+ b[1] = xyz[1] - _y[0];
+ sys2x2(mat, b, uvw);
+ uvw[2] = 0.;
+ }
+#endif
+ int isInside(double u, double v, double w)
+ {
+ if(u < ZERO || v < ZERO || u > (ONE - v))
+ return 0;
+ return 1;
+ }
+};
+
+class quadrangle : public element{
+public:
+ quadrangle(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 2; }
+ inline int getNumNodes(){ return 4; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ w = 0.;
+ switch(num) {
+ case 0 : u = -1.; v = -1.; break;
+ case 1 : u = 1.; v = -1.; break;
+ case 2 : u = 1.; v = 1.; break;
+ case 3 : u = -1.; v = 1.; break;
+ default: u = 0.; v = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 4; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 1; end = 2; break;
+ case 2 : start = 2; end = 3; break;
+ case 3 : start = 3; end = 0; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 4; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u4[4] = {0.577350269189,-0.577350269189,0.577350269189,-0.577350269189};
+ double v4[4] = {0.577350269189,0.577350269189,-0.577350269189,-0.577350269189};
+ double p4[4] = {1.,1.,1.,1.};
+ if(num < 0 || num > 3) return;
+ u = u4[num];
+ v = v4[num];
+ w = 0.;
+ weight = p4[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 0.25 * (1.-u) * (1.-v); break;
+ case 1 : s = 0.25 * (1.+u) * (1.-v); break;
+ case 2 : s = 0.25 * (1.+u) * (1.+v); break;
+ case 3 : s = 0.25 * (1.-u) * (1.+v); break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.25 * (1.-v); s[1] = -0.25 * (1.-u); s[2] = 0.; break;
+ case 1 : s[0] = 0.25 * (1.-v); s[1] = -0.25 * (1.+u); s[2] = 0.; break;
+ case 2 : s[0] = 0.25 * (1.+v); s[1] = 0.25 * (1.+u); s[2] = 0.; break;
+ case 3 : s[0] = -0.25 * (1.+v); s[1] = 0.25 * (1.-u); s[2] = 0.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ double integrateFlux(double val[])
+ {
+ double t1[3] = {_x[1]-_x[0], _y[1]-_y[0], _z[1]-_z[0]};
+ double t2[3] = {_x[2]-_x[0], _y[2]-_y[0], _z[2]-_z[0]};
+ double n[3];
+ prodve(t1, t2, n);
+ norme(n);
+ double v[3];
+ for(int i = 0; i < 3; i++)
+ v[i] = integrate(&val[i], 3);
+ double d;
+ prosca(n, v, &d);
+ return d;
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || v < -ONE || u > ONE || v > ONE)
+ return 0;
+ return 1;
+ }
+};
+
+class tetrahedron : public element{
+public:
+ tetrahedron(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 3; }
+ inline int getNumNodes(){ return 4; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ switch(num) {
+ case 0 : u = 0.; v = 0.; w = 0.; break;
+ case 1 : u = 1.; v = 0.; w = 0.; break;
+ case 2 : u = 0.; v = 1.; w = 0.; break;
+ case 3 : u = 0.; v = 0.; w = 1.; break;
+ default: u = 0.; v = 0.; w = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 6; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 1; end = 2; break;
+ case 2 : start = 2; end = 0; break;
+ case 3 : start = 3; end = 0; break;
+ case 4 : start = 3; end = 2; break;
+ case 5 : start = 3; end = 1; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 4; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u4[4] = {0.138196601125,0.138196601125,0.138196601125,0.585410196625};
+ double v4[4] = {0.138196601125,0.138196601125,0.585410196625,0.138196601125};
+ double w4[4] = {0.138196601125,0.585410196625,0.138196601125,0.138196601125};
+ double p4[4] = {0.0416666666667,0.0416666666667,0.0416666666667,0.0416666666667};
+ if(num < 0 || num > 3) return;
+ u = u4[num];
+ v = v4[num];
+ w = w4[num];
+ weight = p4[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = 1.-u-v-w; break;
+ case 1 : s = u ; break;
+ case 2 : s = v ; break;
+ case 3 : s = w; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -1.; s[1] = -1.; s[2] = -1.; break;
+ case 1 : s[0] = 1.; s[1] = 0.; s[2] = 0.; break;
+ case 2 : s[0] = 0.; s[1] = 1.; s[2] = 0.; break;
+ case 3 : s[0] = 0.; s[1] = 0.; s[2] = 1.; break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ double mat[3][3], b[3];
+ 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];
+ double det;
+ sys3x3(mat, b, uvw, &det);
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < ZERO || v < ZERO || w < ZERO || u > (ONE - v - w))
+ return 0;
+ return 1;
+ }
+};
+
+class hexahedron : public element{
+public:
+ hexahedron(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 3; }
+ inline int getNumNodes(){ return 8; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ switch(num) {
+ case 0 : u = -1.; v = -1.; w = -1.; break;
+ case 1 : u = 1.; v = -1.; w = -1.; break;
+ case 2 : u = 1.; v = 1.; w = -1.; break;
+ case 3 : u = -1.; v = 1.; w = -1.; break;
+ case 4 : u = -1.; v = -1.; w = 1.; break;
+ case 5 : u = 1.; v = -1.; w = 1.; break;
+ case 6 : u = 1.; v = 1.; w = 1.; break;
+ case 7 : u = -1.; v = 1.; w = 1.; break;
+ default: u = 0.; v = 0.; w = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 12; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 0; end = 3; break;
+ case 2 : start = 0; end = 4; break;
+ case 3 : start = 1; end = 2; break;
+ case 4 : start = 1; end = 5; break;
+ case 5 : start = 2; end = 3; break;
+ case 6 : start = 2; end = 6; break;
+ case 7 : start = 3; end = 7; break;
+ case 8 : start = 4; end = 5; break;
+ case 9 : start = 4; end = 7; break;
+ case 10: start = 5; end = 6; break;
+ case 11: start = 6; end = 7; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 6; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u6[6] = { 0.40824826, 0.40824826, -0.40824826,
+ -0.40824826, -0.81649658, 0.81649658};
+ double v6[6] = { 0.70710678, -0.70710678, 0.70710678,
+ -0.70710678, 0., 0.};
+ double w6[6] = {-0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, -0.57735027, 0.57735027};
+ double p6[6] = { 1.3333333333, 1.3333333333, 1.3333333333,
+ 1.3333333333, 1.3333333333, 1.3333333333};
+ if(num < 0 || num > 5) return;
+ u = u6[num];
+ v = v6[num];
+ w = w6[num];
+ weight = p6[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = (1.-u) * (1.-v) * (1.-w) * 0.125; break;
+ case 1 : s = (1.+u) * (1.-v) * (1.-w) * 0.125; break;
+ case 2 : s = (1.+u) * (1.+v) * (1.-w) * 0.125; break;
+ case 3 : s = (1.-u) * (1.+v) * (1.-w) * 0.125; break;
+ case 4 : s = (1.-u) * (1.-v) * (1.+w) * 0.125; break;
+ case 5 : s = (1.+u) * (1.-v) * (1.+w) * 0.125; break;
+ case 6 : s = (1.+u) * (1.+v) * (1.+w) * 0.125; break;
+ case 7 : s = (1.-u) * (1.+v) * (1.+w) * 0.125; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.125 * (1.-v) * (1.-w);
+ s[1] = -0.125 * (1.-u) * (1.-w);
+ s[2] = -0.125 * (1.-u) * (1.-v); break;
+ case 1 : s[0] = 0.125 * (1.-v) * (1.-w);
+ s[1] = -0.125 * (1.+u) * (1.-w);
+ s[2] = -0.125 * (1.+u) * (1.-v); break;
+ case 2 : s[0] = 0.125 * (1.+v) * (1.-w);
+ s[1] = 0.125 * (1.+u) * (1.-w);
+ s[2] = -0.125 * (1.+u) * (1.+v); break;
+ case 3 : s[0] = -0.125 * (1.+v) * (1.-w);
+ s[1] = 0.125 * (1.-u) * (1.-w);
+ s[2] = -0.125 * (1.-u) * (1.+v); break;
+ case 4 : s[0] = -0.125 * (1.-v) * (1.+w);
+ s[1] = -0.125 * (1.-u) * (1.+w);
+ s[2] = 0.125 * (1.-u) * (1.-v); break;
+ case 5 : s[0] = 0.125 * (1.-v) * (1.+w);
+ s[1] = -0.125 * (1.+u) * (1.+w);
+ s[2] = 0.125 * (1.+u) * (1.-v); break;
+ case 6 : s[0] = 0.125 * (1.+v) * (1.+w);
+ s[1] = 0.125 * (1.+u) * (1.+w);
+ s[2] = 0.125 * (1.+u) * (1.+v); break;
+ case 7 : s[0] = -0.125 * (1.+v) * (1.+w);
+ s[1] = 0.125 * (1.-u) * (1.+w);
+ s[2] = 0.125 * (1.-u) * (1.+v); break;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < -ONE || v < -ONE || w < -ONE || u > ONE || v > ONE || w > ONE)
+ return 0;
+ return 1;
+ }
+};
+
+class prism : public element{
+public:
+ prism(double *x, double *y, double *z) : element(x, y, z) {};
+ inline int getDimension(){ return 3; }
+ inline int getNumNodes(){ return 6; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ switch(num) {
+ case 0 : u = 0.; v = 0.; w = -1.; break;
+ case 1 : u = 1.; v = 0.; w = -1.; break;
+ case 2 : u = 0.; v = 1.; w = -1.; break;
+ case 3 : u = 0.; v = 0.; w = 1.; break;
+ case 4 : u = 1.; v = 0.; w = 1.; break;
+ case 5 : u = 0.; v = 1.; w = 1.; break;
+ default: u = 0.; v = 0.; w = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 9; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 0; end = 2; break;
+ case 2 : start = 0; end = 3; break;
+ case 3 : start = 1; end = 2; break;
+ case 4 : start = 1; end = 4; break;
+ case 5 : start = 2; end = 5; break;
+ case 6 : start = 3; end = 4; break;
+ case 7 : start = 3; end = 5; break;
+ case 8 : start = 4; end = 5; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 6; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u6[6] = {0.166666666666666, 0.333333333333333, 0.166666666666666,
+ 0.166666666666666, 0.333333333333333, 0.166666666666666};
+ double v6[6] = {0.166666666666666, 0.166666666666666, 0.333333333333333,
+ 0.166666666666666, 0.166666666666666, 0.333333333333333};
+ double w6[6] = {-0.577350269189, -0.577350269189, -0.577350269189,
+ 0.577350269189, 0.577350269189, 0.577350269189};
+ double p6[6] = {0.166666666666666,0.166666666666666,0.166666666666666,
+ 0.166666666666666,0.166666666666666,0.166666666666666,};
+ if(num < 0 || num > 5) return;
+ u = u6[num];
+ v = v6[num];
+ w = w6[num];
+ weight = p6[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ switch(num) {
+ case 0 : s = (1.-u-v) * (1.-w) * 0.5; break;
+ case 1 : s = u * (1.-w) * 0.5; break;
+ case 2 : s = v * (1.-w) * 0.5; break;
+ case 3 : s = (1.-u-v) * (1.+w) * 0.5; break;
+ case 4 : s = u * (1.+w) * 0.5; break;
+ case 5 : s = v * (1.+w) * 0.5; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ switch(num) {
+ case 0 : s[0] = -0.5 * (1.-w) ;
+ s[1] = -0.5 * (1.-w) ;
+ s[2] = -0.5 * (1.-u-v) ; break ;
+ case 1 : s[0] = 0.5 * (1.-w) ;
+ s[1] = 0. ;
+ s[2] = -0.5 * u ; break ;
+ case 2 : s[0] = 0. ;
+ s[1] = 0.5 * (1.-w) ;
+ s[2] = -0.5 * v ; break ;
+ case 3 : s[0] = -0.5 * (1.+w) ;
+ s[1] = -0.5 * (1.+w) ;
+ s[2] = 0.5 * (1.-u-v) ; break ;
+ case 4 : s[0] = 0.5 * (1.+w) ;
+ s[1] = 0. ;
+ s[2] = 0.5 * u ; break ;
+ case 5 : s[0] = 0. ;
+ s[1] = 0.5 * (1.+w) ;
+ s[2] = 0.5 * v ; break ;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(w > ONE || w < -ONE || u < ZERO || v < ZERO || u > (ONE - v))
+ return 0;
+ return 1;
+ }
+};
+
+class pyramid : public element{
+public:
+ pyramid(double *x, double *y, double *z) : element(x, y, z) {}
+ inline int getDimension(){ return 3; }
+ inline int getNumNodes(){ return 5; }
+ void getNode(int num, double &u, double &v, double &w)
+ {
+ switch(num) {
+ case 0 : u = -1.; v = -1.; w = 0.; break;
+ case 1 : u = 1.; v = -1.; w = 0.; break;
+ case 2 : u = 1.; v = 1.; w = 0.; break;
+ case 3 : u = -1.; v = 1.; w = 0.; break;
+ case 4 : u = 0.; v = 0.; w = 1.; break;
+ default: u = 0.; v = 0.; w = 0.; break;
+ }
+ }
+ inline int getNumEdges(){ return 8; }
+ void getEdge(int num, int &start, int &end)
+ {
+ switch(num) {
+ case 0 : start = 0; end = 1; break;
+ case 1 : start = 0; end = 3; break;
+ case 2 : start = 0; end = 4; break;
+ case 3 : start = 1; end = 2; break;
+ case 4 : start = 1; end = 4; break;
+ case 5 : start = 2; end = 3; break;
+ case 6 : start = 2; end = 4; break;
+ case 7 : start = 3; end = 4; break;
+ default: start = end = 0; break;
+ }
+ }
+ inline int getNumGaussPoints(){ return 8; }
+ void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
+ {
+ double u8[8] = {0.3595161057791018,0.09633205020967324,
+ 0.3595161057791018,0.09633205020967324,
+ 0.6920507403468987,0.1854344369976602,
+ 0.6920507403468987,0.1854344369976602};
+ double v8[8] = {0.3595161057791018,0.3595161057791018,
+ 0.09633205020967324,0.09633205020967324,
+ 0.6920507403468987,0.6920507403468987,
+ 0.1854344369976602,0.1854344369976602};
+ double w8[8] = {0.544151844011225,0.544151844011225,
+ 0.544151844011225,0.544151844011225,
+ 0.122514822655441,0.122514822655441,
+ 0.122514822655441,0.122514822655441};
+ double p8[8] = {0.02519647051995625,0.02519647051995625,
+ 0.02519647051995625,0.02519647051995625,
+ 0.058136862813377,0.058136862813377,
+ 0.058136862813377,0.058136862813377};
+ if(num < 0 || num > 7) return;
+ u = u8[num];
+ v = v8[num];
+ w = w8[num];
+ weight = p8[num];
+ }
+ void getShapeFunction(int num, double u, double v, double w, double &s)
+ {
+ double r;
+ if(w != 1. && num != 4) r = u*v*w / (1.-w);
+ else r = 0.;
+ switch(num) {
+ case 0 : s = 0.25 * ((1.-u) * (1.-v) - w + r); break;
+ case 1 : s = 0.25 * ((1.+u) * (1.-v) - w - r); break;
+ case 2 : s = 0.25 * ((1.+u) * (1.+v) - w + r); break;
+ case 3 : s = 0.25 * ((1.-u) * (1.+v) - w - r); break;
+ case 4 : s = w; break;
+ default : s = 0.; break;
+ }
+ }
+ void getGradShapeFunction(int num, double u, double v, double w, double s[3])
+ {
+ if(w == 1. && num != 4) {
+ s[0] = 0.25;
+ s[1] = 0.25;
+ s[2] = -0.25;
+ }
+ else{
+ switch(num) {
+ case 0 : s[0] = 0.25 * ( -(1.-v) + v*w/(1.-w) ) ;
+ s[1] = 0.25 * ( -(1.-u) + u*w/(1.-w) ) ;
+ s[2] = 0.25 * ( -1. + u*v/SQR(1.-w) ) ; break ;
+ case 1 : s[0] = 0.25 * ( (1.-v) + v*w/(1.-w) ) ;
+ s[1] = 0.25 * ( -(1.+u) + u*w/(1.-w) ) ;
+ s[2] = 0.25 * ( -1. + u*v/SQR(1.-w) ) ; break ;
+ case 2 : s[0] = 0.25 * ( (1.+v) + v*w/(1.-w) ) ;
+ s[1] = 0.25 * ( (1.+u) + u*w/(1.-w) ) ;
+ s[2] = 0.25 * ( -1. + u*v/SQR(1.-w) ) ; break ;
+ case 3 : s[0] = 0.25 * ( -(1.+v) + v*w/(1.-w) ) ;
+ s[1] = 0.25 * ( (1.-u) + u*w/(1.-w) ) ;
+ s[2] = 0.25 * ( -1. + u*v/SQR(1.-w) ) ; break ;
+ case 4 : s[0] = 0. ;
+ s[1] = 0. ;
+ s[2] = 1. ; break ;
+ default : s[0] = s[1] = s[2] = 0.; break;
+ }
+ }
+ }
+ int isInside(double u, double v, double w)
+ {
+ if(u < (w - ONE) || u > (ONE - w) || v < (w - ONE) || v > (ONE - w) ||
+ w < ZERO || w > ONE)
+ return 0;
+ return 1;
+ }
+};
+
+class elementFactory{
+ public:
+ element* create(int numNodes, int dimension, double *x, double *y, double *z){
+ switch(dimension){
+ case 0: return new point(x, y, z);
+ case 1: return new line(x, y, z);
+ case 2:
+ if(numNodes == 4) return new quadrangle(x, y, z);
+ else return new triangle(x, y, z);
+ case 3:
+ if(numNodes == 8) return new hexahedron(x, y, z);
+ else if(numNodes == 6) return new prism(x, y, z);
+ else if(numNodes == 5) return new pyramid(x, y, z);
+ else return new tetrahedron(x, y, z);
+ default:
+ Msg::Error("Unknown type of element in factory");
+ return NULL;
+ }
+ }
+};
+
+#undef ONE
+#undef ZERO
+
+#endif
diff --git a/DataStr/List.cpp b/DataStr/List.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d7d4c7ca1a12bafb42369d518bd4fec96a175b3
--- /dev/null
+++ b/DataStr/List.cpp
@@ -0,0 +1,625 @@
+// $Id: List.cpp,v 1.45 2008-05-04 08:31:11 geuzaine Exp $
+//
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Marc Ume
+//
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <sys/types.h>
+
+#include "Malloc.h"
+#include "List.h"
+#include "Message.h"
+#include "SafeIO.h"
+
+static char *startptr;
+
+List_T *List_Create(int n, int incr, int size)
+{
+ List_T *liste;
+
+ if(n <= 0)
+ n = 1;
+ if(incr <= 0)
+ incr = 1;
+
+ liste = (List_T *) Malloc(sizeof(List_T));
+
+ liste->nmax = 0;
+ liste->incr = incr;
+ liste->size = size;
+ liste->n = 0;
+ liste->isorder = 0;
+ liste->array = NULL;
+
+ List_Realloc(liste, n);
+ return (liste);
+}
+
+void List_Delete(List_T * liste)
+{
+ if(!liste)
+ return;
+ Free(liste->array);
+ Free(liste);
+}
+
+void List_Realloc(List_T * liste, int n)
+{
+ if(n <= 0)
+ return;
+
+ if(liste->array == NULL) {
+ // This does not permit to allocate lists smaller that liste->incr:
+ //liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+ // So this is much better
+ liste->nmax = n;
+ liste->array = (char *)Malloc(liste->nmax * liste->size);
+ }
+ else if(n > liste->nmax) {
+ liste->nmax = ((n - 1) / liste->incr + 1) * liste->incr;
+ liste->array = (char *)Realloc(liste->array, liste->nmax * liste->size);
+ }
+}
+
+void List_Add(List_T * liste, void *data)
+{
+ liste->n++;
+
+ List_Realloc(liste, liste->n);
+ liste->isorder = 0;
+ memcpy(&liste->array[(liste->n - 1) * liste->size], data, liste->size);
+}
+
+int List_Nbr(List_T * liste)
+{
+ return (liste) ? liste->n : 0;
+}
+
+void List_Insert(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ if(List_Search(liste, data, fcmp) == 0)
+ List_Add(liste, data);
+}
+
+int List_Replace(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ void *ptr;
+
+ if(liste->isorder != 1)
+ List_Sort(liste, fcmp);
+ liste->isorder = 1;
+ ptr = (void *)bsearch(data, liste->array, liste->n, liste->size, fcmp);
+ if(ptr == NULL) {
+ List_Add(liste, data);
+ return (0);
+ }
+ else {
+ memcpy(ptr, data, liste->size);
+ return (1);
+ }
+}
+
+void List_Read(List_T * liste, int index, void *data)
+{
+ if((index < 0) || (index >= liste->n))
+ Msg::Fatal("Wrong list index (read)");
+ memcpy(data, &liste->array[index * liste->size], liste->size);
+}
+
+void List_Write(List_T * liste, int index, void *data)
+{
+ if((index < 0) || (index >= liste->n))
+ Msg::Error("Wrong list index (write)");
+ else {
+ liste->isorder = 0;
+ memcpy(&liste->array[index * liste->size], data, liste->size);
+ }
+}
+
+void List_Put(List_T * liste, int index, void *data)
+{
+ if(index < 0)
+ Msg::Error("Wrong list index (put)");
+ else {
+ if(index >= liste->n) {
+ liste->n = index + 1;
+ List_Realloc(liste, liste->n);
+ List_Write(liste, index, data);
+ }
+ else {
+ List_Write(liste, index, data);
+ }
+ }
+}
+
+void List_Pop(List_T * liste)
+{
+ if(liste->n > 0)
+ liste->n--;
+}
+
+void *List_Pointer(List_T * liste, int index)
+{
+ if((index < 0) || (index >= liste->n))
+ Msg::Fatal("Wrong list index (pointer)");
+
+ liste->isorder = 0;
+ return (&liste->array[index * liste->size]);
+}
+
+void *List_Pointer_NoChange(List_T * liste, int index)
+{
+ if((index < 0) || (index >= liste->n))
+ Msg::Fatal("Wrong list index (pointer)");
+
+ return (&liste->array[index * liste->size]);
+}
+
+void *List_Pointer_Fast(List_T * liste, int index)
+{
+ return (&liste->array[index * liste->size]);
+}
+
+void *List_Pointer_Test(List_T * liste, int index)
+{
+ if(!liste || (index < 0) || (index >= liste->n))
+ return NULL;
+
+ liste->isorder = 0;
+ return (&liste->array[index * liste->size]);
+}
+
+void List_Sort(List_T * liste, int (*fcmp) (const void *a, const void *b))
+{
+ qsort(liste->array, liste->n, liste->size, fcmp);
+}
+
+int List_Search(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ void *ptr;
+
+ if(liste->isorder != 1) {
+ List_Sort(liste, fcmp);
+ liste->isorder = 1;
+ }
+ ptr = (void *)bsearch(data, liste->array, liste->n, liste->size, fcmp);
+ if(ptr == NULL)
+ return (0);
+ return (1);
+}
+
+int List_ISearch(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ void *ptr;
+
+ if(liste->isorder != 1)
+ List_Sort(liste, fcmp);
+ liste->isorder = 1;
+ ptr = (void *)bsearch(data, liste->array, liste->n, liste->size, fcmp);
+ if(ptr == NULL)
+ return (-1);
+ return (((long)ptr - (long)liste->array) / liste->size);
+}
+
+int List_ISearchSeq(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ int i;
+
+ if(!liste)
+ return -1;
+ i = 0;
+ while((i < List_Nbr(liste)) && fcmp(data, (void *)List_Pointer(liste, i)))
+ i++;
+ if(i == List_Nbr(liste))
+ i = -1;
+ return i;
+}
+
+int List_ISearchSeqPartial(List_T * liste, void *data, int i_Start,
+ int (*fcmp) (const void *a, const void *b))
+{
+ int i;
+
+ if(!liste)
+ return -1;
+ i = i_Start;
+ while((i < List_Nbr(liste)) && fcmp(data, (void *)List_Pointer(liste, i)))
+ i++;
+ if(i == List_Nbr(liste))
+ i = -1;
+ return i;
+}
+
+int List_Query(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ void *ptr;
+
+ if(liste->isorder != 1)
+ List_Sort(liste, fcmp);
+ liste->isorder = 1;
+ ptr = (void *)bsearch(data, liste->array, liste->n, liste->size, fcmp);
+ if(ptr == NULL)
+ return (0);
+
+ memcpy(data, ptr, liste->size);
+ return (1);
+}
+
+static void *lolofind(void *data, void *array, int n, int size,
+ int (*fcmp) (const void *a, const void *b))
+{
+ char *ptr;
+ int i;
+
+ ptr = (char *)array;
+ for(i = 0; i < n; i++) {
+ if(fcmp(ptr, data) == 0)
+ break;
+ ptr += size;
+ }
+ if(i < n)
+ return (ptr);
+ return (NULL);
+}
+
+int List_LQuery(List_T *liste, void *data,
+ int (*fcmp)(const void *a, const void *b), int first)
+{
+ char *ptr;
+
+ if (first == 1) {
+ ptr = (char *) lolofind(data,liste->array,liste->n,liste->size,fcmp);
+ }
+ else {
+ if (startptr != NULL)
+ ptr = (char *) lolofind(data,startptr,
+ liste->n - (startptr-liste->array)/liste->size,
+ liste->size,fcmp);
+ else
+ return(0);
+ }
+
+ if (ptr == NULL) return(0);
+
+ startptr = ptr + liste->size;
+ if ( startptr >= ( liste->array + liste->n * liste->size))
+ startptr = NULL;
+ memcpy(data,ptr,liste->size);
+ return (1);
+}
+
+void *List_PQuery(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ void *ptr;
+
+ if(liste->isorder != 1)
+ List_Sort(liste, fcmp);
+ liste->isorder = 1;
+ ptr = (void *)bsearch(data, liste->array, liste->n, liste->size, fcmp);
+ return (ptr);
+}
+
+int List_Suppress(List_T * liste, void *data,
+ int (*fcmp) (const void *a, const void *b))
+{
+ char *ptr;
+ int len;
+
+ ptr = (char *)List_PQuery(liste, data, fcmp);
+ if(ptr == NULL)
+ return (0);
+
+ liste->n--;
+ len = liste->n - (((long)ptr - (long)liste->array) / liste->size);
+ if(len > 0)
+ memmove(ptr, ptr + liste->size, len * liste->size);
+ return (1);
+}
+
+int List_PSuppress(List_T * liste, int index)
+{
+ char *ptr;
+ int len;
+
+ ptr = (char *)List_Pointer_NoChange(liste, index);
+ if(ptr == NULL)
+ return (0);
+
+ liste->n--;
+ len = liste->n - (((long)ptr - (long)liste->array) / liste->size);
+ if(len > 0)
+ memmove(ptr, ptr + liste->size, len * liste->size);
+ return (1);
+}
+
+void List_Invert(List_T * a, List_T * b)
+{
+ int i, N;
+ N = List_Nbr(a);
+ for(i = 0; i < N; i++) {
+ List_Add(b, List_Pointer(a, N - i - 1));
+ }
+}
+
+void List_Reset(List_T * liste)
+{
+ if(!liste)
+ return;
+ liste->n = 0;
+}
+
+void List_Action(List_T * liste, void (*action) (void *data, void *dummy))
+{
+ int i, dummy;
+
+ for(i = 0; i < List_Nbr(liste); i++) {
+ (*action) (List_Pointer_NoChange(liste, i), &dummy);
+ }
+
+}
+
+void List_Action_Inverse(List_T * liste,
+ void (*action) (void *data, void *dummy))
+{
+ int i, dummy;
+
+ for(i = List_Nbr(liste); i > 0; i--) {
+ (*action) (List_Pointer_NoChange(liste, i - 1), &dummy);
+ }
+
+}
+
+void List_Copy(List_T * a, List_T * b)
+{
+ int i, N;
+ N = List_Nbr(a);
+ for(i = 0; i < N; i++) {
+ List_Add(b, List_Pointer(a, i));
+ }
+}
+
+void List_Merge(List_T * a, List_T * b)
+{
+ int i;
+
+ if(!a || !b) return;
+ for(i = 0; i < List_Nbr(a); i++) {
+ List_Add(b, List_Pointer_Fast(a, i));
+ }
+}
+
+void swap_bytes(char *array, int size, int n)
+{
+ int i, c;
+ char *x, *a;
+
+ x = (char *)Malloc(size);
+
+ for(i = 0; i < n; i++) {
+ a = &array[i * size];
+ memcpy(x, a, size);
+ for(c = 0; c < size; c++)
+ a[size - 1 - c] = x[c];
+ }
+
+ Free(x);
+}
+
+List_T *List_CreateFromFile(int n, int incr, int size, FILE * file, int format,
+ int swap)
+{
+ int i, error = 0;
+ List_T *liste;
+
+ if(!n){
+ liste = List_Create(1, incr, size);
+ return liste;
+ }
+
+ liste = List_Create(n, incr, size);
+ liste->n = n;
+ switch (format) {
+ case LIST_FORMAT_ASCII:
+ if(size == sizeof(double)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%lf", (double *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(float)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%f", (float *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(int)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%d", (int *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(char)){
+ for(i = 0; i < n; i++){
+ char c = (char)fgetc(file);
+ if(c == EOF){
+ error = 1;
+ break;
+ }
+ else{
+ liste->array[i * size] = c;
+ }
+ }
+ }
+ else{
+ Msg::Error("Bad type of data to create list from (size = %d)", size);
+ error = 1;
+ }
+ break;
+ case LIST_FORMAT_BINARY:
+ if(!fread(liste->array, size, n, file)){
+ error = 1;
+ break;
+ }
+ if(swap)
+ swap_bytes(liste->array, size, n);
+ break;
+ default:
+ Msg::Error("Unknown list format");
+ error = 1;
+ break;
+ }
+
+ if(error){
+ Msg::Error("Read error");
+ liste->n = 0;
+ }
+
+ return liste;
+}
+
+void List_WriteToFile(List_T * liste, FILE * file, int format)
+{
+ int i, n;
+
+ if(!(n = List_Nbr(liste)))
+ return;
+
+ switch (format) {
+ case LIST_FORMAT_ASCII:
+ if(liste->size == sizeof(double))
+ for(i = 0; i < n; i++)
+ fprintf(file, " %.16g", *((double *)&liste->array[i * liste->size]));
+ else if(liste->size == sizeof(float))
+ for(i = 0; i < n; i++)
+ fprintf(file, " %.16g", *((float *)&liste->array[i * liste->size]));
+ else if(liste->size == sizeof(int))
+ for(i = 0; i < n; i++)
+ fprintf(file, " %d", *((int *)&liste->array[i * liste->size]));
+ else if(liste->size == sizeof(char))
+ for(i = 0; i < n; i++)
+ fputc(*((char *)&liste->array[i * liste->size]), file);
+ else
+ Msg::Error("Bad type of data to write list to file (size = %d)",
+ liste->size);
+ break;
+ case LIST_FORMAT_BINARY:
+ safe_fwrite(liste->array, liste->size, n, file);
+ break;
+ default:
+ Msg::Error("Unknown list format");
+ break;
+ }
+}
+
+// For backward compatibility purposes:
+
+List_T *List_CreateFromFileOld(int n, int incr, int size, FILE * file, int format,
+ int swap)
+{
+ int i, error = 0;
+ List_T *liste;
+
+ if(!n){
+ liste = List_Create(1, incr, size);
+ return liste;
+ }
+
+ liste = List_Create(n, incr, size);
+ liste->n = n;
+ switch (format) {
+ case LIST_FORMAT_ASCII:
+ if(size == sizeof(double)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%lf", (double *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(float)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%f", (float *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(int)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%d", (int *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ }
+ }
+ else if(size == sizeof(char)){
+ for(i = 0; i < n; i++){
+ if(!fscanf(file, "%c", (char *)&liste->array[i * size])){
+ error = 1;
+ break;
+ }
+ if(liste->array[i * size] == '^')
+ liste->array[i * size] = '\0';
+ }
+ }
+ else {
+ Msg::Error("Bad type of data to create list from (size = %d)", size);
+ error = 1;
+ }
+ return liste;
+ case LIST_FORMAT_BINARY:
+ if(!fread(liste->array, size, n, file)){
+ error = 1;
+ break;
+ }
+ if(swap)
+ swap_bytes(liste->array, size, n);
+ return liste;
+ default:
+ Msg::Error("Unknown list format");
+ error = 1;
+ break;
+ }
+
+ if(error){
+ Msg::Error("Read error");
+ liste->n = 0;
+ }
+
+ return liste;
+}
diff --git a/DataStr/List.h b/DataStr/List.h
new file mode 100644
index 0000000000000000000000000000000000000000..a352caa8d60541fb2f53b1c3d62c5d6b5c121815
--- /dev/null
+++ b/DataStr/List.h
@@ -0,0 +1,81 @@
+#ifndef _LIST_H_
+#define _LIST_H_
+
+// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Marc Ume
+//
+
+#include <stdio.h>
+
+#define LIST_FORMAT_ASCII 0
+#define LIST_FORMAT_BINARY 1
+
+class List_T {
+public:
+ int nmax;
+ int size;
+ int incr;
+ int n;
+ int isorder;
+ char *array;
+};
+
+List_T *List_Create(int n, int incr, int size);
+void List_Delete(List_T *liste);
+void List_Realloc(List_T *liste,int n);
+void List_Add(List_T *liste, void *data);
+int List_Nbr(List_T *liste);
+void List_Insert(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_Replace(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+void List_Read(List_T *liste, int index, void *data);
+void List_Write(List_T *liste, int index, void *data);
+void List_Put(List_T *liste, int index, void *data);
+void List_Pop(List_T *liste);
+void *List_Pointer(List_T *liste, int index);
+void *List_Pointer_NoChange(List_T *liste, int index);
+void *List_Pointer_Fast(List_T *liste, int index);
+void *List_Pointer_Test(List_T *liste, int index);
+void List_Sort(List_T *liste, int (*fcmp)(const void *a, const void *b));
+int List_Search(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_ISearch(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_ISearchSeq(List_T *liste, void * data, int (*fcmp)(const void *a, const void *b));
+int List_ISearchSeqPartial(List_T *liste, void * data, int i_Start,
+ int (*fcmp)(const void *a, const void *b)) ;
+int List_Query(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_LQuery(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b), int first);
+void *List_PQuery(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_Suppress(List_T *liste, void *data, int (*fcmp)(const void *a, const void *b));
+int List_PSuppress(List_T *liste, int index);
+void List_Invert(List_T *a, List_T *b);
+void List_Reset(List_T *liste);
+void List_Action(List_T *liste, void (*action)(void *data, void *dummy));
+void List_Action_Inverse(List_T *liste, void (*action)(void *data, void *dummy));
+void List_Copy(List_T *a , List_T *b);
+void List_Merge(List_T *a , List_T *b);
+List_T *List_CreateFromFile(int n, int incr, int size, FILE *file, int format, int swap);
+void List_WriteToFile(List_T *liste, FILE *file, int format);
+
+// for backward compatibility
+List_T *List_CreateFromFileOld(int n, int incr, int size, FILE *file, int format, int swap);
+
+#endif
+
diff --git a/benchmarks/2d/Square-03.geo b/benchmarks/2d/Square-03.geo
new file mode 100644
index 0000000000000000000000000000000000000000..b303fff71b5032ba2f4dfbfbb591ec18f1e7ef78
--- /dev/null
+++ b/benchmarks/2d/Square-03.geo
@@ -0,0 +1,16 @@
+/******************************
+Square uniformly meshed
+******************************/
+lc = .49999;
+Point(1) = {0.0,0.0,0,lc};
+Point(2) = {1,0.0,0,lc};
+Point(3) = {1,1,0,lc};
+Point(4) = {0,1,0,lc};
+Point(5) = {0,2,0,lc};
+Line(1) = {3,2};
+Line(2) = {2,1};
+Line(3) = {1,4};
+Line(4) = {4,3};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+Line(7) = {4,5};
diff --git a/examples/cube.geo b/examples/cube.geo
deleted file mode 100644
index decdd4a696bb42dae6ae5bda42e91f3bb5b099e1..0000000000000000000000000000000000000000
--- a/examples/cube.geo
+++ /dev/null
@@ -1,83 +0,0 @@
-lcar1 = .0275;
-
-length = 0.25;
-height = 1.0;
-depth = 0.25;
-
-Point(newp) = {length/2,height/2,depth,lcar1}; /* Point 1 */
-Point(newp) = {length/2,height/2,0,lcar1}; /* Point 2 */
-Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point 3 */
-Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point 4 */
-Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point 5 */
-Point(newp) = {length/2,-height/2,0,lcar1}; /* Point 6 */
-Point(newp) = {-length/2,height/2,0,lcar1}; /* Point 7 */
-Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point 8 */
-Line(1) = {3,1};
-Line(2) = {3,7};
-Line(3) = {7,2};
-Line(4) = {2,1};
-Line(5) = {1,5};
-Line(6) = {5,4};
-Line(7) = {4,8};
-Line(8) = {8,6};
-Line(9) = {6,5};
-Line(10) = {6,2};
-Line(11) = {3,4};
-Line(12) = {8,7};
-Line Loop(13) = {-6,-5,-1,11};
-Plane Surface(14) = {13};
-Line Loop(15) = {4,5,-9,10};
-Plane Surface(16) = {15};
-Line Loop(17) = {-3,-12,8,10};
-Plane Surface(18) = {17};
-Line Loop(19) = {7,12,-2,11};
-Plane Surface(20) = {19};
-Line Loop(21) = {-4,-3,-2,1};
-Plane Surface(22) = {21};
-Line Loop(23) = {8,9,6,7};
-Plane Surface(24) = {23};
-
-Surface Loop(25) = {14,24,-18,22,16,-20};
-Complex Volume(26) = {25};
-
-n = 21;
-
-/*
-sizex = 4;
-sizey = 6;
-sizez = 4;
-*/
-sizex = n*length;
-sizey = n*height;
-sizez = n*depth;
-
-sizex = 9;
-sizey = 33;
-sizez = 9;
-
-
-Transfinite Line {2,4,9,7} = sizez;
-Transfinite Line {11,5,10,12} = sizey;
-Transfinite Line {3,1,6,8} = sizex;
-
-Transfinite Surface {14} = {5,4,3,1};
-Transfinite Surface {16} = {6,2,1,5};
-Transfinite Surface {18} = {6,2,7,8};
-Transfinite Surface {20} = {3,7,8,4};
-Transfinite Surface {22} = {3,7,2,1};
-Transfinite Surface {24} = {4,8,6,5};
-
-Recombine Surface {14,16,18,20,22,24};
-
-Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
-
-Physical Surface(200) = {14,16,18,20,24};
-
-Physical Volume(100) = {26};
-/*
-Mesh.use_cut_plane = 1;
-Mesh.cut_planea = 0;
-Mesh.cut_planeb = 0;
-Mesh.cut_planec = 1;
-Mesh.cut_planed = -0.125;
-*/
\ No newline at end of file
diff --git a/variables.msvc b/variables.msvc
new file mode 100644
index 0000000000000000000000000000000000000000..966725cab25d2cbbf3d005ab3a350770085f5297
--- /dev/null
+++ b/variables.msvc
@@ -0,0 +1,68 @@
+# This is a pre-filled variables file for building a blackbox version
+# of Gmsh with Microsoft Visual C++ (MSVC).
+#
+# This has been tested with MSVC 2003 and MSVC 2008. See
+# doc/README.msvc for building instructions.
+
+# OS and host
+UNAME=WIN32MSVC
+HOSTNAME=localhost
+
+# The names of the C and C++ compilers
+CC=cl
+CXX=cl /EHsc /nologo /GR /MT
+
+# Use /MLd for single-thread debug mode
+# /MTd for multi-thread debug mode
+# /MT for multi-thread release mode
+
+# append different suffix for release or debug version of library
+ifneq (,${findstring MTd,${CXX}})
+ LIBSUFFIX=_d
+else
+ LIBSUFFIX=_r
+endif
+
+# increase stack size to 16Mb to avoid stack overflows in recursive
+# tet classification for large 3D Delaunay grids
+LINKER=cl /F16777216
+
+# All compiler flags except optimization flags
+FLAGS=/DWIN32 /D_USE_MATH_DEFINES /DHAVE_NO_DLL /DHAVE_NO_VSNPRINTF /DHAVE_NO_SNPRINTF /DHAVE_NO_SOCKLEN_T /DHAVE_ANN /DHAVE_MATH_EVAL /DHAVE_TETGEN
+
+# Additional system includes ($INCLUDE is automatically defined by MSVC when
+# you launch the MSVC command prompt)
+SYSINCLUDE=/I"${INCLUDE}"
+
+# Compiler optimization flags
+OPTIM=/O2
+
+# Gmsh subdirectories
+GMSH_DIRS=Common DataStr Geo Mesh Post Numeric Parallel Parser Plugin contrib/ANN contrib/MathEval contrib/NR contrib/Tetgen
+
+# Gmsh libraries
+GMSH_LIBS=Common/Main.obj lib/*.lib
+
+# How you create a static library on this machine
+AR=LIB
+ARFLAGS=/OUT:
+RANLIB=true
+
+# The symbol used in front of compiler flags
+DASH=/
+
+# The extension to use for object files, libraries and executables
+OBJEXT=.obj
+LIBEXT=.lib
+EXEEXT=.exe
+
+# Installation directories
+prefix="S:\Lib\gmsh"
+exec_prefix=${prefix}
+bindir=${exec_prefix}/bin
+datadir=${datarootdir}
+datarootdir=${prefix}/share
+includedir=${prefix}/include
+libdir=${exec_prefix}/lib
+mandir=${datarootdir}/man
+infodir=${datarootdir}/info