This commit was manufactured by cvs2svn to create tag 'gmsh_2_2_1'.

Common/ShapeFunctions.h

+#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

DataStr/List.cpp

+// $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;
+}

DataStr/List.h

+#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
+

benchmarks/2d/Square-03.geo

+/******************************       
+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};

examples/cube.geo

-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

variables.msvc

+# 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