This commit was generated by cvs2svn to compensate for changes in r2,

which included commits to RCS files with non-trunk default branches.

Adapt/Adapt.h

+#ifndef _ADAPT_H_
+#define _ADAPT_H_
+
+#define ADAPT_P1 1
+#define ADAPT_P2 2
+#define ADAPT_H1 3
+#define ADAPT_H2 4
+
+double optimesh (
+		 int N,        /* Nombre d'elements a traiter        */
+		 int method,   /* H1 , P1 , H2 ou P2                 */
+		 int dim,      /* 2 pour 2D et 3 pour 3D             */
+		 double *err,  /* erreurs elementaires               */
+		 double *h,    /* tailles de mailles elementaires    */
+		 double *p,    /* exposante elementaires             */
+		 double e0,    /* erreur prescrite par l'utilisateur */
+		 double N0     /* nbre d'elements ds le maillage opt */
+		 );
+
+#endif

Adapt/Makefile

+#
+# Makefile for "libAdapt.a"
+#
+
+.IGNORE:
+
+CC        = c++
+C_FLAGS   = -g
+
+OS_FLAGS  = -D_UNIX
+
+RM        = rm
+RMFLAGS   = -f
+RANLIB    = ranlib
+
+LIB       = ../lib/libAdapt.a
+INCLUDE   = -I../Common -I../DataStr
+
+CFLAGS    = $(C_FLAGS) $(OS_FLAGS) $(INCLUDE) 
+
+SRC = brent.cpp \
+      mnbrak.cpp \
+      nrutil.cpp \
+      optimesh.cpp \
+      dsvdcmp.cpp \
+      newt.cpp \
+      fmin.cpp \
+      fdjac.cpp \
+      lnsrch.cpp \
+      lubksb.cpp \
+      ludcmp.cpp
+
+OBJ = $(SRC:.cpp=.o)
+
+.SUFFIXES: .o .cpp
+
+$(LIB): $(OBJ)
+	ar ruvs $(LIB) $(OBJ)
+	$(RANLIB) $(LIB)
+
+.cpp.o:
+	$(CC) $(CFLAGS) -c $<
+
+clean:
+	$(RM) $(RMFLAGS) *.o 
+
+lint:
+	$(LINT) $(CFLAGS) $(SRC)
+
+depend:
+	(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+	$(CC) -MM $(CFLAGS) ${SRC} \
+	) >Makefile.new
+	cp Makefile Makefile.bak
+	cp Makefile.new Makefile
+	$(RM) $(RMFLAGS) Makefile.new
+
+# DO NOT DELETE THIS LINE
+brent.o: brent.cpp nrutil.h ../Common/Const.h
+mnbrak.o: mnbrak.cpp nrutil.h ../Common/Const.h
+nrutil.o: nrutil.cpp
+optimesh.o: optimesh.cpp ../Common/Gmsh.h ../Common/Message.h \
+  ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
+  ../DataStr/avl.h ../DataStr/Tools.h Adapt.h nrutil.h ../Common/Const.h
+dsvdcmp.o: dsvdcmp.cpp nrutil.h ../Common/Const.h
+newt.o: newt.cpp nrutil.h ../Common/Const.h
+fmin.o: fmin.cpp nrutil.h ../Common/Const.h
+fdjac.o: fdjac.cpp nrutil.h ../Common/Const.h
+lnsrch.o: lnsrch.cpp nrutil.h ../Common/Const.h
+lubksb.o: lubksb.cpp
+ludcmp.o: ludcmp.cpp nrutil.h ../Common/Const.h

Adapt/brent.cpp

+
+
+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define ITMAX 100
+#define CGOLD 0.3819660
+#define ZEPS 1.0e-10
+#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
+
+double 
+brent (double ax, double bx, double cx, double (*f) (double), double tol,
+       double *xmin)
+{
+  int iter;
+  double a, b, d=0.0, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w,
+    x, xm;
+  double e = 0.0;
+
+  a = (ax < cx ? ax : cx);
+  b = (ax > cx ? ax : cx);
+  x = w = v = bx;
+  fw = fv = fx = (*f) (x);
+  for (iter = 1; iter <= ITMAX; iter++)
+    {
+      xm = 0.5 * (a + b);
+      tol2 = 2.0 * (tol1 = tol * fabs (x) + ZEPS);
+      if (fabs (x - xm) <= (tol2 - 0.5 * (b - a)))
+	{
+	  *xmin = x;
+	  return fx;
+	}
+      if (fabs (e) > tol1)
+	{
+	  r = (x - w) * (fx - fv);
+	  q = (x - v) * (fx - fw);
+	  p = (x - v) * q - (x - w) * r;
+	  q = 2.0 * (q - r);
+	  if (q > 0.0)
+	    p = -p;
+	  q = fabs (q);
+	  etemp = e;
+	  e = d;
+	  if (fabs (p) >= fabs (0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x))
+	    d = CGOLD * (e = (x >= xm ? a - x : b - x));
+	  else
+	    {
+	      d = p / q;
+	      u = x + d;
+	      if (u - a < tol2 || b - u < tol2)
+		d = SIGN (tol1, xm - x);
+	    }
+	}
+      else
+	{
+	  d = CGOLD * (e = (x >= xm ? a - x : b - x));
+	}
+      u = (fabs (d) >= tol1 ? x + d : x + SIGN (tol1, d));
+      fu = (*f) (u);
+      if (fu <= fx)
+	{
+	  if (u >= x)
+	    a = x;
+	  else
+	    b = x;
+	  SHFT (v, w, x, u)
+	    SHFT (fv, fw, fx, fu)
+	}
+      else
+	{
+	  if (u < x)
+	    a = u;
+	  else
+	    b = u;
+	  if (fu <= fw || w == x)
+	    {
+	      v = w;
+	      w = u;
+	      fv = fw;
+	      fw = fu;
+	    }
+	  else if (fu <= fv || v == x || v == w)
+	    {
+	      v = u;
+	      fv = fu;
+	    }
+	}
+    }
+  nrerror ("Too many iterations in brent");
+  *xmin = x;
+  return fx;
+}
+#undef ITMAX
+#undef CGOLD
+#undef ZEPS
+#undef SHFT
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/dsvdcmp.cpp

+#include <math.h>
+
+#include "nrutil.h"
+
+double dpythag(double a, double b)
+{
+  double absa,absb;
+  absa=fabs(a);
+  absb=fabs(b);
+  if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
+  else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
+}
+
+void dsvdcmp(double **a, int m, int n, double w[], double **v)
+{
+  double dpythag(double a, double b);
+  int flag,i,its,j,jj,k,l,nm;
+  double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
+
+  rv1=dvector(1,n);
+  g=scale=anorm=0.0;
+  for (i=1;i<=n;i++) {
+    l=i+1;
+    rv1[i]=scale*g;
+    g=s=scale=0.0;
+    if (i <= m) {
+      for (k=i;k<=m;k++) scale += fabs(a[k][i]);
+      if (scale) {
+	for (k=i;k<=m;k++) {
+	  a[k][i] /= scale;
+	  s += a[k][i]*a[k][i];
+	}
+	f=a[i][i];
+	g = -SIGN(sqrt(s),f);
+	h=f*g-s;
+	a[i][i]=f-g;
+	for (j=l;j<=n;j++) {
+	  for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
+	  f=s/h;
+	  for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
+	}
+	for (k=i;k<=m;k++) a[k][i] *= scale;
+      }
+    }
+    w[i]=scale *g;
+    g=s=scale=0.0;
+    if (i <= m && i != n) {
+      for (k=l;k<=n;k++) scale += fabs(a[i][k]);
+      if (scale) {
+	for (k=l;k<=n;k++) {
+	  a[i][k] /= scale;
+	  s += a[i][k]*a[i][k];
+	}
+	f=a[i][l];
+	g = -SIGN(sqrt(s),f);
+	h=f*g-s;
+	a[i][l]=f-g;
+	for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
+	for (j=l;j<=m;j++) {
+	  for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
+	  for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
+	}
+	for (k=l;k<=n;k++) a[i][k] *= scale;
+      }
+    }
+    anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
+  }
+  for (i=n;i>=1;i--) {
+    if (i < n) {
+      if (g) {
+	for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
+	for (j=l;j<=n;j++) {
+	  for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
+	  for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
+	}
+      }
+      for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
+    }
+    v[i][i]=1.0;
+    g=rv1[i];
+    l=i;
+  }
+  for (i=IMIN(m,n);i>=1;i--) {
+    l=i+1;
+    g=w[i];
+    for (j=l;j<=n;j++) a[i][j]=0.0;
+    if (g) {
+      g=1.0/g;
+      for (j=l;j<=n;j++) {
+	for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
+	f=(s/a[i][i])*g;
+	for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
+      }
+      for (j=i;j<=m;j++) a[j][i] *= g;
+    } else for (j=i;j<=m;j++) a[j][i]=0.0;
+    ++a[i][i];
+  }
+  for (k=n;k>=1;k--) {
+    for (its=1;its<=30;its++) {
+      flag=1;
+      for (l=k;l>=1;l--) {
+	nm=l-1;
+	if ((double)(fabs(rv1[l])+anorm) == anorm) {
+	  flag=0;
+	  break;
+	}
+	if ((double)(fabs(w[nm])+anorm) == anorm) break;
+      }
+      if (flag) {
+	c=0.0;
+	s=1.0;
+	for (i=l;i<=k;i++) {
+	  f=s*rv1[i];
+	  rv1[i]=c*rv1[i];
+	  if ((double)(fabs(f)+anorm) == anorm) break;
+	  g=w[i];
+	  h=dpythag(f,g);
+	  w[i]=h;
+	  h=1.0/h;
+	  c=g*h;
+	  s = -f*h;
+	  for (j=1;j<=m;j++) {
+	    y=a[j][nm];
+	    z=a[j][i];
+	    a[j][nm]=y*c+z*s;
+	    a[j][i]=z*c-y*s;
+	  }
+	}
+      }
+      z=w[k];
+      if (l == k) {
+	if (z < 0.0) {
+	  w[k] = -z;
+	  for (j=1;j<=n;j++) v[j][k] = -v[j][k];
+	}
+	break;
+      }
+      if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations");
+      x=w[l];
+      nm=k-1;
+      y=w[nm];
+      g=rv1[nm];
+      h=rv1[k];
+      f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
+      g=dpythag(f,1.0);
+      f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
+      c=s=1.0;
+      for (j=l;j<=nm;j++) {
+	i=j+1;
+	g=rv1[i];
+	y=w[i];
+	h=s*g;
+	g=c*g;
+	z=dpythag(f,h);
+	rv1[j]=z;
+	c=f/z;
+	s=h/z;
+	f=x*c+g*s;
+	g = g*c-x*s;
+	h=y*s;
+	y *= c;
+	for (jj=1;jj<=n;jj++) {
+	  x=v[jj][j];
+	  z=v[jj][i];
+	  v[jj][j]=x*c+z*s;
+	  v[jj][i]=z*c-x*s;
+	}
+	z=dpythag(f,h);
+	w[j]=z;
+	if (z) {
+	  z=1.0/z;
+	  c=f*z;
+	  s=h*z;
+	}
+	f=c*g+s*y;
+	x=c*y-s*g;
+	for (jj=1;jj<=m;jj++) {
+	  y=a[jj][j];
+	  z=a[jj][i];
+	  a[jj][j]=y*c+z*s;
+	  a[jj][i]=z*c-y*s;
+	}
+      }
+      rv1[l]=0.0;
+      rv1[k]=f;
+      w[k]=x;
+    }
+  }
+  free_dvector(rv1,1,n);
+}
+
+
+/* cf. Numerical Recipes in C, p. 62 */
+
+#define PREC   1.e-16
+
+void invert_singular_matrix(double **M, int n, double **I){
+  double  **V, **T, *W;
+  int     i, j, k; 
+
+  V = dmatrix(1,n,1,n);
+  T = dmatrix(1,n,1,n);
+  W = dvector(1,n);
+
+  dsvdcmp(M, n, n, W, V);
+
+  for(i=1 ; i<=n ; i++){
+    for(j=1 ; j<=n ; j++){
+      I[i][j] = 0.0 ;
+      T[i][j] = 0.0 ;
+    }
+  }
+
+  for(i=1 ; i<=n ; i++){
+    for(j=1 ; j<=n ; j++){
+      if(fabs(W[i]) > PREC){
+	T[i][j] += M[j][i] / W[i] ;
+      }
+      /*
+      else{
+	T[i][j] += 0.0 ;
+      }
+      */
+    }
+  }
+  for(i=1 ; i<=n ; i++){
+    for(j=1 ; j<=n ; j++){
+      for(k=1 ; k<=n ; k++){
+	I[i][j] += V[i][k] * T[k][j] ;
+      }
+    }
+  }
+
+  free_dmatrix(V,1,n,1,n);
+  free_dmatrix(T,1,n,1,n);
+  free_dvector(W,1,n);
+}
+
+#undef PREC 

Adapt/fdjac.cpp

+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define EPS 1.0e-4
+
+void 
+fdjac (int n, float x[], float fvec[], float **df,
+       void (*vecfunc) (int, float[], float[]))
+{
+  int i, j;
+  float h, temp, *f;
+
+  f = vector (1, n);
+  for (j = 1; j <= n; j++)
+    {
+      temp = x[j];
+      h = EPS * fabs (temp);
+      if (h == 0.0)
+	h = EPS;
+      x[j] = temp + h;
+      h = x[j] - temp;
+      (*vecfunc) (n, x, f);
+      x[j] = temp;
+      for (i = 1; i <= n; i++)
+	df[i][j] = (f[i] - fvec[i]) / h;
+    }
+  free_vector (f, 1, n);
+}
+#undef EPS
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/fmin.cpp

+#define NRANSI
+#include "nrutil.h"
+
+extern int nn;
+extern float *fvec;
+extern void (*nrfuncv) (int n, float v[], float f[]);
+
+float 
+fmin (float x[])
+{
+  int i;
+  float sum;
+
+  (*nrfuncv) (nn, x, fvec);
+  for (sum = 0.0, i = 1; i <= nn; i++)
+    sum += SQR (fvec[i]);
+  return 0.5 * sum;
+}
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/lnsrch.cpp

+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define ALF 1.0e-4
+#define TOLX 1.0e-7
+
+void 
+lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
+	float *f, float stpmax, int *check, float (*func) (float[]))
+{
+  int i;
+  float a, alam, alam2, alamin, b, disc, f2, fold2, rhs1, rhs2, slope,
+    sum, temp, test, tmplam;
+
+  *check = 0;
+  for (sum = 0.0, i = 1; i <= n; i++)
+    sum += p[i] * p[i];
+  sum = sqrt (sum);
+  if (sum > stpmax)
+    for (i = 1; i <= n; i++)
+      p[i] *= stpmax / sum;
+  for (slope = 0.0, i = 1; i <= n; i++)
+    slope += g[i] * p[i];
+  test = 0.0;
+  for (i = 1; i <= n; i++)
+    {
+      temp = fabs (p[i]) / FMAX (fabs (xold[i]), 1.0);
+      if (temp > test)
+	test = temp;
+    }
+  alamin = TOLX / test;
+  alam = 1.0;
+  for (;;)
+    {
+      for (i = 1; i <= n; i++)
+	x[i] = xold[i] + alam * p[i];
+      *f = (*func) (x);
+      if (alam < alamin)
+	{
+	  for (i = 1; i <= n; i++)
+	    x[i] = xold[i];
+	  *check = 1;
+	  return;
+	}
+      else if (*f <= fold + ALF * alam * slope)
+	return;
+      else
+	{
+	  if (alam == 1.0)
+	    tmplam = -slope / (2.0 * (*f - fold - slope));
+	  else
+	    {
+	      rhs1 = *f - fold - alam * slope;
+	      rhs2 = f2 - fold2 - alam2 * slope;
+	      a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
+	      b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
+	      if (a == 0.0)
+		tmplam = -slope / (2.0 * b);
+	      else
+		{
+		  disc = b * b - 3.0 * a * slope;
+		  if (disc < 0.0)
+		    nrerror ("Roundoff problem in lnsrch.");
+		  else
+		    tmplam = (-b + sqrt (disc)) / (3.0 * a);
+		}
+	      if (tmplam > 0.5 * alam)
+		tmplam = 0.5 * alam;
+	    }
+	}
+      alam2 = alam;
+      f2 = *f;
+      fold2 = fold;
+      alam = FMAX (tmplam, 0.1 * alam);
+    }
+}
+#undef ALF
+#undef TOLX
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/lubksb.cpp

+void 
+lubksb (float **a, int n, int *indx, float b[])
+{
+  int i, ii = 0, ip, j;
+  float sum;
+
+  for (i = 1; i <= n; i++)
+    {
+      ip = indx[i];
+      sum = b[ip];
+      b[ip] = b[i];
+      if (ii)
+	for (j = ii; j <= i - 1; j++)
+	  sum -= a[i][j] * b[j];
+      else if (sum)
+	ii = i;
+      b[i] = sum;
+    }
+  for (i = n; i >= 1; i--)
+    {
+      sum = b[i];
+      for (j = i + 1; j <= n; j++)
+	sum -= a[i][j] * b[j];
+      b[i] = sum / a[i][i];
+    }
+}
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/ludcmp.cpp

+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define TINY 1.0e-20;
+
+void 
+ludcmp (float **a, int n, int *indx, float *d)
+{
+  int i, imax, j, k;
+  float big, dum, sum, temp;
+  float *vv;
+
+  vv = vector (1, n);
+  *d = 1.0;
+  for (i = 1; i <= n; i++)
+    {
+      big = 0.0;
+      for (j = 1; j <= n; j++)
+	if ((temp = fabs (a[i][j])) > big)
+	  big = temp;
+      if (big == 0.0)
+	nrerror ("Singular matrix in routine ludcmp");
+      vv[i] = 1.0 / big;
+    }
+  for (j = 1; j <= n; j++)
+    {
+      for (i = 1; i < j; i++)
+	{
+	  sum = a[i][j];
+	  for (k = 1; k < i; k++)
+	    sum -= a[i][k] * a[k][j];
+	  a[i][j] = sum;
+	}
+      big = 0.0;
+      for (i = j; i <= n; i++)
+	{
+	  sum = a[i][j];
+	  for (k = 1; k < j; k++)
+	    sum -= a[i][k] * a[k][j];
+	  a[i][j] = sum;
+	  if ((dum = vv[i] * fabs (sum)) >= big)
+	    {
+	      big = dum;
+	      imax = i;
+	    }
+	}
+      if (j != imax)
+	{
+	  for (k = 1; k <= n; k++)
+	    {
+	      dum = a[imax][k];
+	      a[imax][k] = a[j][k];
+	      a[j][k] = dum;
+	    }
+	  *d = -(*d);
+	  vv[imax] = vv[j];
+	}
+      indx[j] = imax;
+      if (a[j][j] == 0.0)
+	a[j][j] = TINY;
+      if (j != n)
+	{
+	  dum = 1.0 / (a[j][j]);
+	  for (i = j + 1; i <= n; i++)
+	    a[i][j] *= dum;
+	}
+    }
+  free_vector (vv, 1, n);
+}
+#undef TINY
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/mnbrak.cpp

+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define GOLD 1.618034
+#define GLIMIT 100.0
+#define TINY 1.0e-20
+#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
+
+void 
+mnbrak (double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
+	double (*func) (double))
+{
+  double ulim, u, r, q, fu, dum;
+
+  *fa = (*func) (*ax);
+  *fb = (*func) (*bx);
+  if (*fb > *fa)
+    {
+      SHFT (dum, *ax, *bx, dum)
+	SHFT (dum, *fb, *fa, dum)
+    }
+  *cx = (*bx) + GOLD * (*bx - *ax);
+  *fc = (*func) (*cx);
+  while (*fb > *fc)
+    {
+      r = (*bx - *ax) * (*fb - *fc);
+      q = (*bx - *cx) * (*fb - *fa);
+      u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
+	(2.0 * SIGN (FMAX (fabs (q - r), TINY), q - r));
+      ulim = (*bx) + GLIMIT * (*cx - *bx);
+      if ((*bx - u) * (u - *cx) > 0.0)
+	{
+	  fu = (*func) (u);
+	  if (fu < *fc)
+	    {
+	      *ax = (*bx);
+	      *bx = u;
+	      *fa = (*fb);
+	      *fb = fu;
+	      return;
+	    }
+	  else if (fu > *fb)
+	    {
+	      *cx = u;
+	      *fc = fu;
+	      return;
+	    }
+	  u = (*cx) + GOLD * (*cx - *bx);
+	  fu = (*func) (u);
+	}
+      else if ((*cx - u) * (u - ulim) > 0.0)
+	{
+	  fu = (*func) (u);
+	  if (fu < *fc)
+	    {
+	      SHFT (*bx, *cx, u, *cx + GOLD * (*cx - *bx))
+		SHFT (*fb, *fc, fu, (*func) (u))
+	    }
+	}
+      else if ((u - ulim) * (ulim - *cx) >= 0.0)
+	{
+	  u = ulim;
+	  fu = (*func) (u);
+	}
+      else
+	{
+	  u = (*cx) + GOLD * (*cx - *bx);
+	  fu = (*func) (u);
+	}
+      SHFT (*ax, *bx, *cx, u)
+	SHFT (*fa, *fb, *fc, fu)
+    }
+}
+#undef GOLD
+#undef GLIMIT
+#undef TINY
+#undef SHFT
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/newt.cpp

+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define MAXITS 200
+#define TOLF 1.0e-4
+#define TOLMIN 1.0e-6
+#define TOLX 1.0e-7
+#define STPMX 100.0
+
+int nn;
+float *fvec;
+void (*nrfuncv) (int n, float v[], float f[]);
+#define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\
+	free_vector(p,1,n);free_vector(g,1,n);free_matrix(fjac,1,n,1,n);\
+	free_ivector(indx,1,n);return;}
+
+void 
+newt (float x[], int n, int *check,
+      void (*vecfunc) (int, float[], float[]))
+{
+  void fdjac (int n, float x[], float fvec[], float **df,
+	      void (*vecfunc) (int, float[], float[]));
+  float fmin (float x[]);
+  void lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
+	       float *f, float stpmax, int *check, float (*func) (float[]));
+  void lubksb (float **a, int n, int *indx, float b[]);
+  void ludcmp (float **a, int n, int *indx, float *d);
+  int i, its, j, *indx;
+  float d, den, f, fold, stpmax, sum, temp, test, **fjac, *g, *p, *xold;
+
+  indx = ivector (1, n);
+  fjac = matrix (1, n, 1, n);
+  g = vector (1, n);
+  p = vector (1, n);
+  xold = vector (1, n);
+  fvec = vector (1, n);
+  nn = n;
+  nrfuncv = vecfunc;
+  f = fmin (x);
+  test = 0.0;
+  for (i = 1; i <= n; i++)
+    if (fabs (fvec[i]) > test)
+      test = fabs (fvec[i]);
+  if (test < 0.01 * TOLF)
+    FREERETURN
+      for (sum = 0.0, i = 1; i <= n; i++)
+      sum += SQR (x[i]);
+  stpmax = STPMX * FMAX (sqrt (sum), (float) n);
+  for (its = 1; its <= MAXITS; its++)
+    {
+      fdjac (n, x, fvec, fjac, vecfunc);
+      for (i = 1; i <= n; i++)
+	{
+	  for (sum = 0.0, j = 1; j <= n; j++)
+	    sum += fjac[j][i] * fvec[j];
+	  g[i] = sum;
+	}
+      for (i = 1; i <= n; i++)
+	xold[i] = x[i];
+      fold = f;
+      for (i = 1; i <= n; i++)
+	p[i] = -fvec[i];
+      ludcmp (fjac, n, indx, &d);
+      lubksb (fjac, n, indx, p);
+      lnsrch (n, xold, fold, g, p, x, &f, stpmax, check, fmin);
+      test = 0.0;
+      for (i = 1; i <= n; i++)
+	if (fabs (fvec[i]) > test)
+	  test = fabs (fvec[i]);
+      if (test < TOLF)
+	{
+	  *check = 0;
+	  FREERETURN
+	}
+      if (*check)
+	{
+	  test = 0.0;
+	  den = FMAX (f, 0.5 * n);
+	  for (i = 1; i <= n; i++)
+	    {
+	      temp = fabs (g[i]) * FMAX (fabs (x[i]), 1.0) / den;
+	      if (temp > test)
+		test = temp;
+	    }
+	  *check = (test < TOLMIN ? 1 : 0);
+	  FREERETURN
+	}
+      test = 0.0;
+      for (i = 1; i <= n; i++)
+	{
+	  temp = (fabs (x[i] - xold[i])) / FMAX (fabs (x[i]), 1.0);
+	  if (temp > test)
+	    test = temp;
+	}
+      if (test < TOLX)
+	FREERETURN
+	}
+	nrerror ("MAXITS exceeded in newt");
+    }
+#undef MAXITS
+#undef TOLF
+#undef TOLMIN
+#undef TOLX
+#undef STPMX
+#undef FREERETURN
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */

Adapt/nrutil.cpp

+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#include <malloc.h>
+
+
+#define NR_END 1
+#define FREE_ARG char*
+
+void nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+	fprintf(stderr,"Numerical Recipes run-time error...\n");
+	fprintf(stderr,"%s\n",error_text);
+	fprintf(stderr,"...now exiting to system...\n");
+	exit(1);
+}
+
+float *vector(long nl, long nh)
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+	float *v;
+
+	v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
+	if (!v) nrerror("allocation failure in vector()");
+	return v-nl+NR_END;
+}
+
+int *ivector(long nl, long nh)
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+	int *v;
+
+	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
+	if (!v) nrerror("allocation failure in ivector()");
+	return v-nl+NR_END;
+}
+
+unsigned char *cvector(long nl, long nh)
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+	unsigned char *v;
+
+	v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+	if (!v) nrerror("allocation failure in cvector()");
+	return v-nl+NR_END;
+}
+
+unsigned long *lvector(long nl, long nh)
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+	unsigned long *v;
+
+	v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
+	if (!v) nrerror("allocation failure in lvector()");
+	return v-nl+NR_END;
+}
+
+double *dvector(long nl, long nh)
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+	double *v;
+
+	v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
+	if (!v) nrerror("allocation failure in dvector()");
+	return v-nl+NR_END;
+}
+
+float **matrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+double **dmatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	double **m;
+
+	/* allocate pointers to rows */
+	m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+int **imatrix(long nrl, long nrh, long ncl, long nch)
+/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
+{
+	long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	int **m;
+
+	/* allocate pointers to rows */
+	m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
+	if (!m) nrerror("allocation failure 1 in matrix()");
+	m += NR_END;
+	m -= nrl;
+
+
+	/* allocate rows and set pointers to them */
+	m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
+	if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
+	m[nrl] += NR_END;
+	m[nrl] -= ncl;
+
+	for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
+	long newrl, long newcl)
+/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
+{
+	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
+	float **m;
+
+	/* allocate array of pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in submatrix()");
+	m += NR_END;
+	m -= newrl;
+
+	/* set pointers to rows */
+	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
+
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
+/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
+declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
+and ncol=nch-ncl+1. The routine should be called with the address
+&a[0][0] as the first argument. */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
+	float **m;
+
+	/* allocate pointers to rows */
+	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
+	if (!m) nrerror("allocation failure in convert_matrix()");
+	m += NR_END;
+	m -= nrl;
+
+	/* set pointers to rows */
+	m[nrl]=a-ncl;
+	for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
+	/* return pointer to array of pointers to rows */
+	return m;
+}
+
+float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
+/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
+{
+	long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
+	float ***t;
+
+	/* allocate pointers to pointers to rows */
+	t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
+	if (!t) nrerror("allocation failure 1 in f3tensor()");
+	t += NR_END;
+	t -= nrl;
+
+	/* allocate pointers to rows and set pointers to them */
+	t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
+	if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
+	t[nrl] += NR_END;
+	t[nrl] -= ncl;
+
+	/* allocate rows and set pointers to them */
+	t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
+	if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
+	t[nrl][ncl] += NR_END;
+	t[nrl][ncl] -= ndl;
+
+	for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
+	for(i=nrl+1;i<=nrh;i++) {
+		t[i]=t[i-1]+ncol;
+		t[i][ncl]=t[i-1][ncl]+ncol*ndep;
+		for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
+	}
+
+	/* return pointer to array of pointers to rows */
+	return t;
+}
+
+void free_vector(float *v, long nl, long nh)
+/* free a float vector allocated with vector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(int *v, long nl, long nh)
+/* free an int vector allocated with ivector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_cvector(unsigned char *v, long nl, long nh)
+/* free an unsigned char vector allocated with cvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_lvector(unsigned long *v, long nl, long nh)
+/* free an unsigned long vector allocated with lvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_dvector(double *v, long nl, long nh)
+/* free a double vector allocated with dvector() */
+{
+	free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
+/* free a float matrix allocated by matrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
+/* free a double matrix allocated by dmatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
+/* free an int matrix allocated by imatrix() */
+{
+	free((FREE_ARG) (m[nrl]+ncl-NR_END));
+	free((FREE_ARG) (m+nrl-NR_END));
+}
+
+void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a submatrix allocated by submatrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
+/* free a matrix allocated by convert_matrix() */
+{
+	free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
+	long ndl, long ndh)
+/* free a float f3tensor allocated by f3tensor() */
+{
+	free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
+	free((FREE_ARG) (t[nrl]+ncl-NR_END));
+	free((FREE_ARG) (t+nrl-NR_END));
+}
+

Adapt/nrutil.h

+#ifndef _NR_UTILS_H_
+#define _NR_UTILS_H_
+
+#include "Const.h"
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+
+void nrerror (char error_text[]);
+float *vector (long nl, long nh);
+int *ivector (long nl, long nh);
+unsigned char *cvector (long nl, long nh);
+unsigned long *lvector (long nl, long nh);
+double *dvector (long nl, long nh);
+float **matrix (long nrl, long nrh, long ncl, long nch);
+double **dmatrix (long nrl, long nrh, long ncl, long nch);
+int **imatrix (long nrl, long nrh, long ncl, long nch);
+float **submatrix (float **a, long oldrl, long oldrh, long oldcl, long oldch,
+		   long newrl, long newcl);
+float **convert_matrix (float *a, long nrl, long nrh, long ncl, long nch);
+float ***f3tensor (long nrl, long nrh, long ncl, long nch, long ndl, long ndh);
+void free_vector (float *v, long nl, long nh);
+void free_ivector (int *v, long nl, long nh);
+void free_cvector (unsigned char *v, long nl, long nh);
+void free_lvector (unsigned long *v, long nl, long nh);
+void free_dvector (double *v, long nl, long nh);
+void free_matrix (float **m, long nrl, long nrh, long ncl, long nch);
+void free_dmatrix (double **m, long nrl, long nrh, long ncl, long nch);
+void free_imatrix (int **m, long nrl, long nrh, long ncl, long nch);
+void free_submatrix (float **b, long nrl, long nrh, long ncl, long nch);
+void free_convert_matrix (float **b, long nrl, long nrh, long ncl, long nch);
+void free_f3tensor (float ***t, long nrl, long nrh, long ncl, long nch,
+		    long ndl, long ndh);
+
+
+#endif /* _NR_UTILS_H_ */

Adapt/optimesh.cpp

+
+#include "Gmsh.h"
+#include "Adapt.h"
+#include "nrutil.h"
+
+#define TOL 1.e-08
+#define MAXDEG 999
+
+void frprmn(double p[], int n, double ftol, int *iter, double *fret,
+	double (*func)(double []), void (*dfunc)(double [], double []));
+
+extern void SetError(char *, char *);
+
+static int NN,METHOD;
+static double MINH , *ERR , *HH , *PP , E0, DIM ;
+
+/* METODE H VERSION 1 : MINIMISER LE NOMBRE D'ELEMENTS
+   TOUT EN GARDANT UNE ERREUR GLOBALE DONNEE.  ON MODIFIE
+   ICI LE VECTEUR TAILLE DE MAILLE
+*/
+
+double fH1 ( double l ){
+
+  int i;
+  double val1,val2;
+
+  val1 = 0.0;
+
+  for(i=1;i<=NN;i++){
+    val1 += pow(2.*l*DSQR(ERR[i])*PP[i]/DIM,(DIM/(2.*PP[i]+DIM)));
+  }
+  val2 = 0.0;
+  for(i=1;i<=NN;i++){
+    val2 += DSQR(ERR[i])*pow(2.*l*DSQR(ERR[i])*PP[i]/DIM,(-2.*PP[i]/(2.*PP[i]+DIM)));
+  }
+  return -(val1 + l * ( val2 - DSQR(E0)));
+}
+
+/* METODE H VERSION 2 : MINIMISER L'ERREUR
+   TOUT EN GARDANT UN NOMBRE D'ELEMENTS DONNE.  ON MODIFIE
+   ICI LE VECTEUR TAILLE DE MAILLE
+*/
+
+double fH2 ( double l ){
+
+  int i;
+  double val1,val2,qi;
+
+  val1 = val2 = 0.0;
+  for(i=1;i<=NN;i++){
+    qi = pow((DIM*l)/(2.*PP[i] * DSQR(ERR[i])) 
+	     ,-DIM/(DIM+2.*PP[i]));
+    val1 += DSQR(ERR[i]) * pow(qi,-2.*PP[i]/DIM);
+    val2 += qi;
+  }
+/*
+  printf( "%12.5e %12.5e\n",l,val1 + l * ( val2 - E0));
+*/
+  return -(val1 + l * ( val2 - E0));
+}
+
+
+/* METODE P VERSION 1 : MINIMISER LE NOMBRE D'ELEMENTS
+   TOUT EN GARDANT UNE ERREUR GLOBALE DONNEE.  ON MODIFIE
+   ICI LE VECTEUR DEGRE D'INTERPOLATION
+*/
+
+double fP1 ( double l ){
+
+  int i;
+  double val1,val2,qi,e;
+
+  val1 = val2 = 0.0;
+  for(i=1;i<=NN;i++){
+    e = ERR[i];
+    if(e==0.0)e=1.e-12;
+    qi = - log (2. * l * log (HH[i]/MINH) * DSQR(e)) / log (HH[i]/MINH);
+    val1 -= .5 * qi;
+    val2 += pow(HH[i]/MINH,qi) * DSQR(e);
+  }
+/*
+  printf( "%12.5e %12.5e\n",l,val1 + l * ( val2 - DSQR(E0)));
+*/
+  return -(val1 + l * ( val2 - DSQR(E0)));
+}
+
+
+double min1d ( double (*funct)(double), double *xmin){
+  double brent(double ax, double bx, double cx,
+	      double (*f)(double), double tol, double *xmin);
+  void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
+		double *fc, double (*func)(double));
+  double xx,fx,fb,fa,bx,ax;
+
+  switch(METHOD){
+  case ADAPT_H1:
+  case ADAPT_P1:
+    ax=1.e-12;
+    xx=1.e2;
+    break;
+  default:
+    ax=1.e-15;
+    xx=1.e-12;
+    break;
+  }    
+  mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,funct);
+  return brent(ax,xx,bx,funct,TOL,xmin);
+}
+  
+
+/* optimesh renvoie la contrainte (N0 ou e0) pour le probleme optimise */
+
+double optimesh (int N,        /* Nombre d'elements a traiter        */
+		 int method,   /* ADAPT_H1 , ADAPT_P1 , ADAPT_H2 ou ADAPT_P2 */
+		 int dim,      /* 2 pour 2D et 3 pour 3D             */
+		 double *err,  /* erreurs elementaires               */
+		 double *h,    /* tailles de mailles elementaires    */
+		 double *p,    /* exposante elementaires             */
+		 double e0,    /* erreur prescrite par l'utilisateur */
+		 double N0     /* nbre d'elements ds le maillage opt */
+		 )
+{
+
+  int i;
+  double contr,pivrai,lambda,minf,qi,ri,pi,obj,obj2,minri,maxri;
+  double errmin,errmax;
+
+
+  Msg(INFOS, "N=%d Meth=%d dim=%d err[1]=%g err[2]=%g p[1]=%g p[2]=%g prescr=%g",
+      N,method,dim,err[1],err[2],p[1],p[2],e0);
+
+  METHOD = method;
+
+  h[N+1] = 1.0;
+  p[N+1] = 1.0;
+
+  NN  = N;
+  ERR = err;
+  HH  = h;
+  PP = p;
+  NN = N;
+  E0 = e0;
+  DIM = (double)dim;
+
+  
+  for(i=1;i<=N;i++){
+    if(i==1){
+      errmin=errmax=err[i];
+    }
+    else{
+      errmin = DMIN(errmin,err[i]);
+      errmax = DMAX(errmax,err[i]);
+    }
+  }
+
+  switch (method) {
+  case ADAPT_H1 :
+
+    minf = min1d (fH1,&lambda);
+    obj = 0.0;
+    for(i=1;i<=N;i++){
+      qi = pow(2.*lambda*DSQR(err[i])*p[i]/DIM,DIM/(2.*p[i]+DIM));
+      ri = pow(qi,1./DIM);
+      if(i==1){
+		minri=maxri=ri;
+      }
+      if(err[i]==0.0)ri = .5;
+
+      minri = DMIN(minri,ri);
+      maxri = DMAX(maxri,ri);
+      obj += DSQR(err[i]) * pow(ri,-2.*p[i]) ; 
+      h[i-1] = sqrt(2.) * h[i]/ri;
+      p[i-1] = ri;
+    }
+
+    contr = fabs(minf);
+
+    Msg(INFOS, "Constraint  : asked %g <==> obtained %g",e0,sqrt(obj));
+    Msg(INFOS, "Objective function (Nb. of elements) : %g",-minf);
+    Msg(INFOS, "Minimum reduction factor : %g maximum : %g",minri,maxri);
+    break;
+
+  case ADAPT_P1 :
+    MINH=h[1];
+    for(i=1;i<=N;i++){
+      MINH =DMIN(h[i],MINH);
+    }
+    MINH/=2.;
+
+    minf = min1d (fP1,&lambda);
+    obj = obj2 = 0.0;
+    for(i=1;i<=N;i++){
+      qi = -log(2.*lambda*DSQR(err[i])*log(h[i]/MINH))/log(h[i]/MINH);
+      pi = p[i] - .5 * qi;
+      pivrai = DMIN(DMAX(1.,(double)(int)(pi+.99)),MAXDEG);
+      obj2 += pow(h[i]/MINH,2.*(p[i]-pivrai))*DSQR(err[i]);
+      obj += DSQR(err[i]) * pow(h[i]/MINH,qi) ; 
+      h[i-1] = h[i];
+      p[i-1] = pi;
+    }
+    Msg(INFOS, "Constraint : %g = %g ==> %g",e0,sqrt(obj),sqrt(obj2));
+    Msg(INFOS, "Objective function : %g",minf);
+    contr = fabs(minf);
+    break;
+
+  case ADAPT_H2 :
+    minf = min1d (fH2,&lambda);
+    obj = 0.0;
+    for(i=1;i<=N;i++){
+      qi = pow((DIM*lambda)/(2.*DSQR(err[i])*p[i])
+	       ,-DIM/(DIM+2.*p[i]));
+      ri = pow(qi,1./DIM);
+      if(i==1){
+	minri=maxri=ri;
+      }
+      minri = DMIN(minri,ri);
+      maxri = DMAX(maxri,ri);
+      obj += pow(ri,DIM) ; 
+      h[i-1] = h[i]/ri;
+      p[i-1] = p[i];
+    }
+
+    Msg(INFOS, "Constraint : %g = %g",e0,obj);
+    Msg(INFOS, "Objective function (Error in %%) : %g", 100. * sqrt(fabs(minf)));
+    Msg(INFOS, "Minri : %g maximum %g",minri,maxri);
+    contr = sqrt(fabs(minf));
+    break;
+
+  case ADAPT_P2 :
+    minf = min1d (fH1,&lambda);
+    break;
+
+  default :
+    Msg(WARNING, "Unknown mesh optimisation method");
+  }
+
+  return contr;
+}

Box/Box.cpp

+
+#include <signal.h>
+
+#include "Gmsh.h"
+#include "Const.h"
+#include "Geo.h"
+#include "Mesh.h"
+#include "Views.h"
+#include "Parser.h"
+#include "Context.h"
+#include "Main.h"
+#include "MinMax.h"
+
+#include "Static.h"
+
+int VERBOSE = 0 ;
+
+char progname[]  = "This is Gmsh (non-interactive)" ;
+char copyright[] = "Copyright (C) 1997-2000 C. Geuzaine, J.-F. Remacle" ;
+char clargs[]    = 
+  "Usage: %s [options] [files]\n"
+  "Mesh options:\n"
+  "  -0                    parse input and exit\n"
+  "  -1, -2, -3            batch 1-, 2- or 3-dimensional mesh\n"
+  "  -smooth int           mesh smoothing (default: 3)\n"
+  "  -degree int           mesh degree (default: 1)\n"
+  "  -format msh|unv|gref  mesh format (default: msh)\n"
+  "  -algo iso|aniso       mesh algorithm (default: iso)\n"
+  "  -scale float          scaling factor (default: 1.0)\n"
+  "  -recombine            recombine extruded meshes\n"
+  "  -bgm file             load backround mesh from file\n"
+  "Other options:\n"	  
+  "  -v                    print debug information\n"
+  "  -path string          set path for included files\n"
+  "  -version              show version number\n"
+  "  -info                 show detailed version information\n"
+  "  -help                 show this message\n"
+  ;
+
+
+/* dummy defs for link purposes */
+
+void color_table_init_param (int number, ColorTable * ct, int rgb_flag, int alpha_flag){;}
+void color_table_recompute (ColorTable * ct, int rgb_flag, int alpha_flag){;}
+void ZeroHighlight(Mesh *){;}
+void AddView(int, char *, int){;}
+void draw_polygon_2d (double, double, double, int, double *, double *, double *){;}
+
+/* ------------------------------------------------------------------------ */
+/*  p a r s e                                                               */
+/* ------------------------------------------------------------------------ */
+
+void ParseFile(char *f){
+  strncpy(yyname,f,NAME_STR_L);
+  yyerrorstate=0;
+  yylineno=1;
+  if(!(yyin = fopen(yyname,"r"))){
+    Msg(INFO, "File '%s' dos not exist", f);
+    return;
+  }
+  while(!feof(yyin)) yyparse();
+  fclose(yyin);
+}
+
+void MergeProblem(char *name){
+  Msg(INFOS, "Merging %s",name); 
+
+  ParseFile(name);  
+  if (yyerrorstate) return;
+}
+
+void OpenProblem(char *name){
+  char ext[6];
+  
+  InitSymbols();
+  Init_Mesh(&M, 1);
+  BD_EXISTS = 1;
+
+  strncpy(TheFileName,name,NAME_STR_L);
+  strncpy(TheBaseFileName,name,NAME_STR_L);
+
+  strcpy(ext,name+(strlen(name)-4));
+  if(!strcmp(ext,".GEO") || 
+     !strcmp(ext,".geo") || 
+     !strcmp(ext,".msh") || 
+     !strcmp(ext,".pos")){
+    TheBaseFileName[strlen(name)-4] = '\0';
+  }
+  else{
+    strcat(TheFileName,".geo");
+  }
+
+  strncpy(THEM->name, TheBaseFileName,NAME_STR_L);
+
+  Msg(INFOS, "Opening %s", TheFileName); 
+
+  ParseFile(TheFileName);  
+
+  mai3d(THEM,0);  
+  
+  Maillage_Dimension_0(&M);
+  ZeroHighlight(&M); 
+  CalculateMinMax(THEM->Points);  
+}
+
+/* ------------------------------------------------------------------------ */
+/*  G e t _ O p t i o n s                                                   */
+/* ------------------------------------------------------------------------ */
+
+void Get_Options (int argc, char *argv[]) {
+  int i=1;
+
+  if(argc < 2) Info(0,argv[0]);
+
+  strncpy(TheFileNameTab[0], "unnamed.geo",NAME_STR_L);
+  
+  while (i < argc) {
+    
+    if (argv[i][0] == '-') {
+      
+      if(!strcmp(argv[i]+1, "0")){ 
+	CTX.interactive = -1; i++;
+      }
+      else if(!strcmp(argv[i]+1, "1")){ 
+	CTX.interactive = 1; i++;
+      }
+      else if(!strcmp(argv[i]+1, "2")){ 
+	CTX.interactive = 2; i++;
+      }
+      else if(!strcmp(argv[i]+1, "3")){ 
+	CTX.interactive = 3; i++;
+      }
+      else if(!strcmp(argv[i]+1, "v")){ 
+	VERBOSE = 1; i++;
+      }
+      else if(!strcmp(argv[i]+1, "path")){ 
+	i++;
+	if(argv[i] != NULL){
+	  strncpy(ThePathForIncludes,argv[i++],NAME_STR_L);
+	}
+      }
+      else if(!strcmp(argv[i]+1, "bgm")){ 
+	i++;
+	if(argv[i] != NULL){
+	  strncpy(TheBgmFileName,argv[i++],NAME_STR_L);
+	  INITIALBGMESH = ONFILE;
+	}
+      }
+      else if(!strcmp(argv[i]+1, "smooth")){ 
+	i++;
+	LISSAGE = atoi(argv[i]); i++;
+      }
+      else if(!strcmp(argv[i]+1, "scale")){
+	i++;
+	GLOBALSCALINGFACTOR = atof(argv[i]); i++;
+      }
+      else if(!strcmp(argv[i]+1, "degree")){  
+	i++;
+	if(argv[i]!=NULL){
+	  CTX.mesh.degree = atoi(argv[i]); i++;
+	  if(CTX.mesh.degree != 1 || CTX.mesh.degree != 2){
+	    fprintf(stderr, "Error: Wrong degree\n");
+	    exit(1);
+	  }
+	}
+	else {	  
+	  fprintf(stderr, "Error: Missing Number\n");
+	  exit(1);
+	}
+      }
+      else if(!strcmp(argv[i]+1, "format")){  
+	i++;
+	if(argv[i]!=NULL){
+	  if(!strcmp(argv[i],"msh"))
+	    CTX.mesh.format = FORMAT_MSH ;
+	  else if(!strcmp(argv[i],"unv"))
+	    CTX.mesh.format = FORMAT_UNV ;
+	  else if(!strcmp(argv[i],"gref"))
+	    CTX.mesh.format = FORMAT_GREF ;
+	  else{
+	    fprintf(stderr, "Error: Unknown mesh format\n");
+	    exit(1);
+	  }
+	  i++;
+	}
+	else {	  
+	  fprintf(stderr, "Error: Missing format\n");
+	  exit(1);
+	}
+      }
+      else if(!strcmp(argv[i]+1, "algo")){  
+	i++;
+	if(argv[i]!=NULL){
+	  if(!strcmp(argv[i],"iso"))
+	    CTX.mesh.algo = DELAUNAY_OLDALGO ;
+	  else if(!strcmp(argv[i],"aniso"))
+	    CTX.mesh.algo = DELAUNAY_NEWALGO ;
+	  else{
+	    fprintf(stderr, "Error: Unknown mesh algorithm\n");
+	    exit(1);
+	  }
+	  i++;
+	}
+	else {	  
+	  fprintf(stderr, "Error: Missing algorithm\n");
+	  exit(1);
+	}
+      }
+      else if(!strcmp(argv[i]+1, "info")){
+	Info(2,argv[0]); 
+      }
+      else if(!strcmp(argv[i]+1, "version")){
+	Info(1,argv[0]); 
+      }
+      else if(!strcmp(argv[i]+1, "help")){
+	Info(0,argv[0]);
+      }
+      else{
+	fprintf(stderr, "Warning: Unknown option '%s'\n", argv[i]);
+	Info(0,argv[0]);
+      }
+    }
+
+    else {
+      if(NbFileName<MAX_OPEN_FILES){
+	strncpy(TheFileNameTab[NbFileName++], argv[i++], NAME_STR_L); 
+      }
+      else{
+	fprintf(stderr, "Error: Too many input files\n");
+	exit(1);
+      }
+    }
+
+  }
+
+  strncpy(TheFileName, TheFileNameTab[0], NAME_STR_L);
+
+}
+
+/* ------------------------------------------------------------------------ */
+/*  m a i n                                                                 */
+/* ------------------------------------------------------------------------ */
+
+int main(int argc, char *argv[]){
+  int     i;
+
+  InitContext(&CTX);
+  Get_Options(argc, argv);
+
+  signal(SIGINT,  Signal); 
+  signal(SIGSEGV, Signal);
+  signal(SIGFPE,  Signal); 
+
+  OpenProblem(TheFileName);
+  if(yyerrorstate)
+    exit(1);
+  else{
+    if(NbFileName>1){
+      for(i=1;i<NbFileName;i++) MergeProblem(TheFileNameTab[i]);
+    }
+    if(INITIALBGMESH == ONFILE){
+      MergeProblem(TheBgmFileName);
+      if(List_Nbr(Post_ViewList)){
+	BGMWithView((Post_View*)List_Pointer(Post_ViewList, List_Nbr(Post_ViewList)-1));
+	TYPBGMESH = ONFILE; 
+	Create_BgMesh(TYPBGMESH,.2,THEM);
+      }
+      else{
+	fprintf(stderr, "Error: invalid BGM (no view)\n"); exit(1);
+      }
+    }
+    if(CTX.interactive > 0){
+      mai3d(THEM, CTX.interactive);
+      Print_Mesh(THEM,NULL,CTX.mesh.format);
+    }
+    exit(1);
+  }    
+
+}
+
+
+
+/* ------------------------------------------------------------------------ */
+/*  I n f o                                                                 */
+/* ------------------------------------------------------------------------ */
+
+void Info (int level, char *arg0){
+  switch(level){
+  case 0 :
+    fprintf(stderr, "%s\n", progname);
+    fprintf(stderr, "%s\n", copyright);
+    fprintf(stderr, clargs, arg0);
+    exit(1);
+  case 1:
+    fprintf(stderr, "%g\n", GMSH_VERSION);
+    exit(1) ; 
+  case 2:
+    fprintf(stderr, "Version    : %g\n", GMSH_VERSION);
+    fprintf(stderr, "OS         : %s\n", GMSH_OS);
+    fprintf(stderr, "Build Date : %s\n", GMSH_DATE);
+    fprintf(stderr, "Build Host : %s\n", GMSH_HOST);
+    fprintf(stderr, "Packager   : %s\n", GMSH_PACKAGER);
+    exit(1) ; 
+  default :
+    break;
+  }
+}
+
+
+/* ------------------------------------------------------------------------ */
+/*  S i g n a l                                                             */
+/* ------------------------------------------------------------------------ */
+
+
+void Signal (int sig_num){
+
+  switch (sig_num){
+  case SIGSEGV : Msg(ERROR, "Segmentation Violation (invalid memory reference)"); break;
+  case SIGFPE  : Msg(ERROR, "Floating point exception (division by zero?)"); break;
+  case SIGINT  : Msg(ERROR, "Interrupt (generated from terminal special char)"); break;
+  default      : Msg(ERROR, "Unknown signal"); break;
+  }
+}
+
+
+/* ------------------------------------------------------------------------ */
+/*  M s g                                                                   */
+/* ------------------------------------------------------------------------ */
+
+void Msg(int level, char *fmt, ...){
+  va_list  args;
+  int      abort=0;
+  int      nb, nbvis;
+
+  va_start (args, fmt);
+
+  switch(level){
+
+  case PARSER_ERROR :
+    fprintf(stderr, "Parse Error: "); vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
+    break ;
+
+  case PARSER_INFO :
+    if(VERBOSE){
+      fprintf(stderr, "Parse Info: "); vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
+    }
+    break ;
+
+  case ERROR :
+    fprintf(stderr, "Error: "); vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
+    abort = 1 ;
+    break ;
+
+  case WARNING :
+    fprintf(stderr, "Warning: "); vfprintf(stderr, fmt,args); fprintf(stderr, "\n");
+    break;
+
+  case INFOS :
+  case INFO :
+  case SELECT :
+  case STATUS :
+    if(VERBOSE){
+      vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
+    }
+    break;
+  }
+
+  va_end (args);
+
+  if(abort) exit(1);
+
+}
+
+/* ------------------------------------------------------------------------ */
+/*  C p u                                                                   */
+/* ------------------------------------------------------------------------ */
+
+double Cpu(void){
+  return 0.;
+}
+
+/* ------------------------------------------------------------------------ */
+/*  P r o g r e s s                                                         */
+/* ------------------------------------------------------------------------ */
+
+void Progress(int i){
+}

Box/Makefile

+#
+# Makefile for ".a"
+#
+
+.IGNORE:
+
+CC            = c++
+C_FLAGS       = -g
+
+VERSION_FLAGS = 
+OS_FLAGS      = -D_UNIX 
+
+RM       = rm
+RMFLAGS  = -f
+RANLIB   = ranlib
+
+LIB      = ../lib/libBox.a
+INCLUDE  = -I../Common -I../DataStr -I../Geo\
+           -I../Graphics -I../Mesh -I../Parser -I../Unix
+
+CFLAGS   = $(C_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE) 
+
+SRC = Box.cpp
+
+OBJ = $(SRC:.cpp=.o)
+
+.SUFFIXES: .o .cpp
+
+$(LIB): $(OBJ) 
+	ar ruvs $(LIB) $(OBJ) 
+	$(RANLIB) $(LIB)
+
+.cpp.o:
+	$(CC) $(CFLAGS) -c $<
+
+clean:
+	$(RM) $(RMFLAGS) *.o
+
+lint:
+	$(LINT) $(CFLAGS) $(SRC)
+
+depend:
+	(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+	$(CC) -MM $(CFLAGS) ${SRC} \
+	) >Makefile.new
+	cp Makefile Makefile.bak
+	cp Makefile.new Makefile
+	$(RM) $(RMFLAGS) Makefile.new
+
+# DO NOT DELETE THIS LINE
+box.o: box.cpp ../includes/msh_Include.h ../GmshDataStr/listman.h \
+ ../GmshDataStr/treeman.h ../GmshDataStr/avl.h ../GmshDataStr/ualloc.h \
+ ../GmshDataStr/outil.h ../base/base.h ../geo/geo.h \
+ ../GmshUnix/message.h ../includes/msh_Const.h ../mesh/structure.h \
+ ../includes/msh_Tree.h ../includes/msh_Proto.h \
+ ../includes/msh_Version.h ../geo/context.h ../GmshGL/Draw_Post.h \
+ ../GmshGL/colortable.h

Common/Const.h

+#ifndef _CONSTS_H_
+#define _CONSTS_H_
+
+//#define RAND_LONG    LC * ((rand()%1000)/1.E08)
+//#define RAND_LONG    LC * ((rand()%1000)/1.E08)/10.
+
+//RAND_LONG in [0, LC/1.e6]
+#define RAND_LONG    (LC/1.e6*rand()/RAND_MAX)
+
+//EPSILON_LONG in [0, LC/1.e12]
+#define EPSILON_LC   (LC/1.e12*rand()/RAND_MAX)
+
+#define RADTODEG      57.29578
+
+#define TEXT_BUFFER_SIZE       1024
+#define SELECTION_BUFFER_SIZE  1024
+#define LABEL_STR_L            16
+#define NAME_STR_L             256
+#define MAX_OPEN_FILES         20
+
+#define RacineDeDeux  1.41421356237
+#define RacineDeTrois 1.73205080757
+#define Pi            3.14159265359
+#define Deux_Pi       6.28318530718
+#define UnTiers       0.33333333333
+
+#define MIN(a,b) (((a)<(b))?(a):(b))
+#define MAX(a,b) (((a)<(b))?(b):(a))
+#define SQR(a)   (((a)==0.0)?0.0:((a)*(a)))
+
+#define IMIN MIN
+#define LMIN MIN
+#define FMIN MIN
+#define DMIN MIN
+
+#define IMAX MAX
+#define LMAX MAX
+#define FMAX MAX
+#define DMAX MAX
+
+#define DSQR SQR
+#define FSQR SQU
+
+#define THRESHOLD(a,b,c) (((a)>(c))?(c):((a)<(b)?(b):(a)))
+
+#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b)))
+
+#define sign(x)   (((x)>=0)?1:-1)
+#define Pred(x)   ((x)->prev)
+#define Succ(x)   ((x)->next)
+#define square(x) ((x)*(x))
+
+#endif

Common/Context.cpp

+
+#include "Gmsh.h"
+#include "Const.h"
+#include "Geo.h"
+#include "Mesh.h"
+#include "Draw.h"
+#include "Context.h"
+
+void InitColors(rgbacolors * col, int num){
+
+  if(col->id == num){
+    return ;
+  }
+
+  col->id = num ;
+
+  switch(num){
+
+  case 0 : /* default drawing colors: black background */
+  case 1 : /* alternative drawing colors: white background */
+    switch(num){
+    case 0 :
+      col->bg                 = PACK_COLOR(0,   0,   0,   255) ;
+      col->fg                 = PACK_COLOR(255, 255, 255, 255) ;
+      col->text               = PACK_COLOR(255, 255, 255, 255) ;
+      col->axes               = PACK_COLOR(255, 255, 0,   255) ;
+      col->little_axes        = PACK_COLOR(255, 255, 255, 255) ;
+      break;
+    case 1 :
+      col->bg                 = PACK_COLOR(255, 255, 255, 255) ;
+      col->fg                 = PACK_COLOR(0,   0,   0,   255) ;
+      col->text               = PACK_COLOR(0,   0,   0,   255) ;
+      col->axes               = PACK_COLOR(128, 128, 128, 255) ;
+      col->little_axes        = PACK_COLOR(0,   0,   0,   255) ;
+      break;
+    }
+    col->geom.point         = PACK_COLOR(178,   182, 129,   255) ;
+    col->geom.line          = PACK_COLOR(0,   0,   255, 255) ;
+    col->geom.surface       = PACK_COLOR(128, 128, 128, 255) ;
+    col->geom.volume        = PACK_COLOR(128, 128, 128, 255) ;
+    col->geom.point_sel     = PACK_COLOR(255, 0,   0,   255) ;
+    col->geom.line_sel      = PACK_COLOR(255, 0,   0,   255) ;
+    col->geom.surface_sel   = PACK_COLOR(255, 0,   0,   255) ;
+    col->geom.volume_sel    = PACK_COLOR(255, 0,   0,   255) ;
+    col->geom.point_hlt     = PACK_COLOR(0,   255, 0,   255) ;
+    col->geom.line_hlt      = PACK_COLOR(0,   0,   255, 255) ;
+    col->geom.surface_hlt   = PACK_COLOR(128, 128, 128, 255) ;
+    col->geom.volume_hlt    = PACK_COLOR(128, 128, 128, 255) ;
+    col->geom.tangents      = PACK_COLOR(255, 255, 0,   255) ;
+    col->geom.normals       = PACK_COLOR(255, 0,   0,   255) ;
+
+    col->mesh.vertex        = PACK_COLOR(0  , 123, 59 , 255) ;
+    col->mesh.vertex_supp   = PACK_COLOR(255, 0,   255, 255) ;
+    col->mesh.line          = PACK_COLOR(0,   255, 0,   255) ; 
+    col->mesh.triangle      = PACK_COLOR(0,   255, 0,   255) ;
+    col->mesh.quadrangle    = PACK_COLOR(0,   255, 0,   255) ; 
+    col->mesh.tetrahedron   = PACK_COLOR(0,   255, 0,   255) ;
+    col->mesh.hexahedron    = PACK_COLOR(128, 255, 0,   255) ;
+    col->mesh.prism         = PACK_COLOR(0,   255, 128, 255) ;
+    col->mesh.pyramid       = PACK_COLOR(128, 255, 128, 255) ;
+    col->mesh.tangents      = PACK_COLOR(128, 128, 128, 255) ;
+    col->mesh.normals       = PACK_COLOR(128, 128, 128, 255) ;
+
+    col->mesh.carousel[0]   = PACK_COLOR(0  , 82 , 138, 255) ;
+    col->mesh.carousel[1]   = PACK_COLOR(255, 0  , 0  , 255) ;
+    col->mesh.carousel[2]   = PACK_COLOR(31 , 110, 171, 255) ;
+    col->mesh.carousel[3]   = PACK_COLOR(255, 255, 0  , 255) ; 
+    col->mesh.carousel[4]   = PACK_COLOR(255, 0  , 255, 255) ; 
+    col->mesh.carousel[4]   = PACK_COLOR(0  , 255, 255, 255) ; 
+    col->mesh.carousel[5]   = PACK_COLOR(128, 128, 0  , 255) ; 
+    col->mesh.carousel[6]   = PACK_COLOR(128, 0  , 255, 255) ; 
+    col->mesh.carousel[7]   = PACK_COLOR(128, 128, 255, 255) ; 
+    col->mesh.carousel[8]   = PACK_COLOR(128, 128, 255, 255) ; 
+    col->mesh.carousel[9]   = PACK_COLOR(0  , 0  , 255, 255) ; 
+    break;
+    
+  case 2 : /* grayscale */
+    col->bg                 = PACK_COLOR(255, 255, 255, 255) ;
+    col->fg                 = PACK_COLOR(0,   0,   0,   255) ;
+    col->text               = PACK_COLOR(0,   0,   0,   255) ;
+    col->axes               = PACK_COLOR(0,   0,   0,   255) ;
+    col->little_axes        = PACK_COLOR(0,   0,   0,   255) ;
+    
+    col->geom.point         = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.line          = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.surface       = PACK_COLOR(0,   0,   0,   255) ; 
+    col->geom.volume        = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.point_sel     = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.line_sel      = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.surface_sel   = PACK_COLOR(0,   0,   0,   255) ; 
+    col->geom.volume_sel    = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.point_hlt     = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.line_hlt      = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.surface_hlt   = PACK_COLOR(0,   0,   0,   255) ; 
+    col->geom.volume_hlt    = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.tangents      = PACK_COLOR(0,   0,   0,   255) ;
+    col->geom.normals       = PACK_COLOR(0,   0,   0,   255) ;
+    
+    col->mesh.vertex        = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.vertex_supp   = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.line          = PACK_COLOR(0,   0,   0,   255) ; 
+    col->mesh.triangle      = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.quadrangle    = PACK_COLOR(0,   0,   0,   255) ; 
+    col->mesh.tetrahedron   = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.hexahedron    = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.prism         = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.pyramid       = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.tangents      = PACK_COLOR(0,   0,   0,   255) ;
+    col->mesh.normals       = PACK_COLOR(0,   0,   0,   255) ;
+
+    col->mesh.carousel[0]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[1]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[2]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[3]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[4]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[4]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[5]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[6]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[7]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[8]   = PACK_COLOR(255, 255, 255, 255) ;
+    col->mesh.carousel[9]   = PACK_COLOR(255, 255, 255, 255) ;
+
+  }
+
+}
+
+
+void InitContext(Context_T *ctx){
+
+  ctx->interactive       = 0 ;
+
+  ctx->r[0]              = 0.0 ;
+  ctx->r[1]              = 0.0 ;
+  ctx->r[2]              = 0.0 ;
+  ctx->t[0]              = 0.0 ;
+  ctx->t[1]              = 0.0 ;
+  ctx->t[2]              = 0.0 ;
+  ctx->s[0]              = 1.0 ;
+  ctx->s[1]              = 1.0 ;
+  ctx->s[2]              = 1.0 ;
+
+  ctx->min[0]          = 0.0 ;
+  ctx->min[1]          = 0.0 ;
+  ctx->min[2]          = 0.0 ;
+  ctx->max[0]          = 0.0 ;
+  ctx->max[1]          = 0.0 ;
+  ctx->max[2]          = 0.0 ;
+  ctx->range[0]          = 0.0 ;
+  ctx->range[1]          = 0.0 ;
+  ctx->range[2]          = 0.0 ;
+
+  ctx->viewport[0]       = 0 ;
+  ctx->viewport[1]       = 0 ;
+  ctx->viewport[2]       = 1 ;
+  ctx->viewport[3]       = 1 ;
+  ctx->render_mode       = GMSH_RENDER ;
+  ctx->pixel_equiv_x     = 0. ;
+  ctx->pixel_equiv_y     = 0. ; 
+
+#ifdef _UNIX
+  ctx->font_string          = "-*-helvetica-medium-r-*-*-*-*-*-*-*-*-*-*";
+  ctx->colorbar_font_string = "fixed";
+#else
+  ctx->font_string          = "dummy";
+  ctx->colorbar_font_string = "dummy";
+#endif
+
+  ctx->light0[0]         = 0.5 ;
+  ctx->light0[1]         = 0.3 ;
+  ctx->light0[2]         = 1.0 ;
+  ctx->light0[3]         = 0.0 ;
+  ctx->shine             = 0.4 ;
+  ctx->alpha             = 0 ; /* disable alpha blending by default */
+  ctx->flash             = 0 ;
+  ctx->same_visual       = 0 ;
+
+  ctx->db                = 1 ;
+  ctx->overlay           = 1 ;
+  ctx->stream            = TO_SCREEN ;
+  ctx->axes              = 1 ;
+  ctx->little_axes       = 1 ;
+  ctx->ortho             = 1 ;
+  ctx->fast              = 1 ;
+  ctx->display_lists     = 0 ; 
+  ctx->command_win       = 0 ;
+  ctx->threads           = 0 ; // bugge avec Open3D sur DEC ...
+  ctx->threads_lock      = 0 ;
+
+  ctx->geom.vis_type     = 0 ;
+  ctx->geom.points       = 1 ;
+  ctx->geom.lines        = 1 ;
+  ctx->geom.surfaces     = 0 ;
+  ctx->geom.volumes      = 0 ;
+  ctx->geom.points_num   = 0 ;
+  ctx->geom.lines_num    = 0 ;
+  ctx->geom.surfaces_num = 0 ;
+  ctx->geom.volumes_num  = 0 ;
+  ctx->geom.level        = ELEMENTARY ;
+  ctx->geom.normals      = 0.0 ;
+  ctx->geom.tangents     = 0.0 ;
+  ctx->geom.highlight    = 1 ;
+  ctx->geom.hidden       = 0 ;
+  ctx->geom.shade        = 0 ;
+
+  ctx->mesh.vis_type     = 0 ;
+  ctx->mesh.draw         = 1 ;  
+  ctx->mesh.points       = 1 ;
+  ctx->mesh.lines        = 1 ;
+  ctx->mesh.surfaces     = 1 ;
+  ctx->mesh.volumes      = 1 ;
+  ctx->mesh.points_num   = 0 ;
+  ctx->mesh.lines_num    = 0 ;
+  ctx->mesh.surfaces_num = 0 ;
+  ctx->mesh.volumes_num  = 0 ;
+  ctx->mesh.normals      = 0.0 ;
+  ctx->mesh.tangents     = 0.0 ;  
+  ctx->mesh.explode      = 1.0 ;
+  ctx->mesh.hidden       = 0 ;  
+  ctx->mesh.shade        = 0 ;  
+  ctx->mesh.format       = FORMAT_MSH ;
+  ctx->mesh.nb_smoothing = 0 ;
+  ctx->mesh.algo         = DELAUNAY_OLDALGO ;
+  ctx->mesh.degree       = 1 ;
+  ctx->mesh.is_limit_gamma = 0 ;
+  ctx->mesh.limit_gamma  = 0.1 ;
+  ctx->mesh.dual         = 0 ;
+  ctx->mesh.reco_extrude = 0 ;
+
+  ctx->post.draw         = 1 ;
+  ctx->post.scales       = 1 ;
+  ctx->post.link         = 0 ;
+  ctx->post.font         = "Courier" ;
+  ctx->post.fontsize     = 12 ;
+  ctx->post.initial_visibility = 1 ;
+  ctx->post.initial_intervals = DRAW_POST_ISO ;
+  ctx->post.initial_nbiso = 15 ;
+  ctx->post.anim_delay    = 0 ;
+
+#ifdef _UNIX
+  ctx->print.type        = GLPRPAINTER ;
+  ctx->print.format      = FORMAT_EPS ;
+#else
+  ctx->print.type        = -1 ;
+  ctx->print.format      = -1 ;
+#endif
+
+  ctx->color.id = -1;
+  InitColors(&ctx->color, 0) ;
+
+}
+

Common/Context.h

+#ifndef _CONTEXT_H_
+#define _CONTEXT_H_
+
+/* 
+   Interface-independant context
+ */
+
+/* How RGBA values are packed and unpacked into/from a 4-byte integer */
+
+#  ifdef _LITTLE
+#    define PACK_COLOR(R,G,B,A)   ( (A)<<24 | (B)<<16 | (G)<<8 | (R) )
+#    define UNPACK_RED(X)         ( (X) & 0xff )
+#    define UNPACK_GREEN(X)       ( ( (X) >> 8 ) & 0xff )
+#    define UNPACK_BLUE(X)        ( ( (X) >> 16 ) & 0xff )
+#    define UNPACK_ALPHA(X)       ( ( (X) >> 24 ) & 0xff )
+#  else
+#    define PACK_COLOR(R,G,B,A)   ( (R)<<24 | (G)<<16 | (B)<<8 | (A) )
+#    define UNPACK_RED(X)         ( ( (X) >> 24 ) & 0xff )
+#    define UNPACK_GREEN(X)       ( ( (X) >> 16 ) & 0xff )
+#    define UNPACK_BLUE(X)        ( ( (X) >> 8 ) & 0xff )
+#    define UNPACK_ALPHA(X)       ( (X) & 0xff )
+#  endif
+
+typedef struct{
+  int id;			/* the current rgbacolors id */
+  
+  /* general colors */
+  unsigned int bg, fg, text, axes, little_axes;
+  
+  /* geometry colors */
+  struct{
+    unsigned int point, line, surface, volume;
+    unsigned int point_sel, line_sel, surface_sel, volume_sel;
+    unsigned int point_hlt, line_hlt, surface_hlt, volume_hlt;
+    unsigned int tangents, normals;
+  }
+  geom;
+  
+  /* mesh colors */
+  struct{
+    unsigned int vertex, vertex_supp, line, triangle, quadrangle;
+    unsigned int tetrahedron, hexahedron, prism, pyramid;
+    unsigned int carousel[10];
+    unsigned int tangents, normals;
+  }
+  mesh;
+  
+}rgbacolors;
+
+typedef struct {
+
+  int interactive;		/* 0=full gfx; -1=just parse; 1,2,3=batch mesh */
+
+  double r[3], t[3], s[3];	/* current rotation, translation and scale */
+  int rlock[3], tlock[3], slock[3];	/* locks for r, t and s */
+  
+  double min[3];		/* x, y and z min for the current geometry */
+  double max[3];		/* x, y and z max for the current geometry */
+  double range[3];		/* maximum range in the three directions */
+  
+  int db;			/* double buffer? */
+  int overlay;	   	        /* overlay graphic window? */
+  int stream;			/* output stream: TO_SCREEN or TO_FILE */
+  int ortho;			/* orthogonal projection? */
+  int fast;			/* inhibit mesh and postpro drawing when changing r,s,t */
+  int command_win;		/* command window? */
+  int display_lists;		/* use display lists? */
+  int font_base;		/* display list indice for the font */
+  int axes, little_axes;	/* draw axes? */
+  int threads, threads_lock;	/* threads?, lock (should be a mutex...) */
+  int alpha;                    /* enable alpha blending */
+  int flash;                    /* authorize colormap flashing (beek) */
+  int same_visual;              /* force same visual for GUI and Graphics */
+  
+  char *font_string;          /* main font */
+  char *colorbar_font_string; /* font for colorbar */
+
+  /* OpenGL stuff */
+  int viewport[4];
+  float light0[4];            /* light source position */
+  float shine;                /* specular value */
+  int render_mode;		/* RENDER, SELECT, FEEDBACK */
+
+  double pixel_equiv_x, pixel_equiv_y ; /* approximative equivalent model lenght of a pixel */
+  
+  /* all colors except postpro colormaps */
+  rgbacolors color;
+  
+  /* geometry options */
+  struct{
+    int vis_type;
+    int points, lines, surfaces, volumes;
+    int points_num, lines_num, surfaces_num, volumes_num;
+    int hidden, shade;
+    int highlight;
+    int level;
+    double normals, tangents;
+  } geom;
+
+  /* mesh options */
+  struct {
+    int vis_type;
+    int draw;
+    int points, lines, surfaces, volumes;
+    int points_num, lines_num, surfaces_num, volumes_num;
+    int is_limit_gamma, dual;
+    double limit_gamma;
+    int hidden, shade;
+    int format, nb_smoothing, algo, degree;
+    int reco_extrude;
+    double normals, tangents, explode;
+  } mesh;
+
+  /* post processing options */
+  struct{
+    int draw, scales, link ;
+    char *font;
+    int  fontsize;
+    int  initial_visibility, initial_nbiso, initial_intervals ;
+    long anim_delay ;
+  }post;
+
+  /* print options */
+  struct{
+    int format, type;
+  } print;
+
+} Context_T;
+
+typedef struct {
+  char *string ; 
+  int int1, int2, int3, int4 ;
+} StringX4Int;
+
+typedef struct {
+  char *string ; 
+  void *Pointer ;
+} StringXPointer ;
+
+void InitContext (Context_T * ctx);
+void InitColors (rgbacolors * col, int num);
+
+#endif