diff --git a/Box/Main.cpp b/Box/Main.cpp
index fae497f6516904117167b35cfe9318d98e1f1d93..8477f23be580417bd3e78a715618e400cb582fec 100644
--- a/Box/Main.cpp
+++ b/Box/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.25 2003-01-24 23:13:34 geuzaine Exp $
+// $Id: Main.cpp,v 1.26 2003-02-18 05:50:04 geuzaine Exp $
//
// Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
//
@@ -113,6 +113,8 @@ int main(int argc, char *argv[]){
if(CTX.default_plugins)
GMSH_PluginManager::Instance()->RegisterDefaultPlugins();
+ check_gsl();
+
OpenProblem(CTX.filename);
if(yyerrorstate)
ParUtil::Instance()->Abort();
diff --git a/Box/Makefile b/Box/Makefile
index 0d63a0277910952795c84d4127a4b011ebf51184..1d542dac983b01f7612cab1708b5fc6bf7f2c92b 100644
--- a/Box/Makefile
+++ b/Box/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.20 2003-02-11 08:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.21 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshBox.a
INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics -I../Mesh -I../Numeric\
-I../Parser -I../Fltk -I../Plugin -I../Parallel
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Main.cpp
diff --git a/Common/Makefile b/Common/Makefile
index 76325ccd441c6a17f9eddba5a704345e281de3df..c4b94ff27a8d6cb466ace5ce61179a836802ef7f 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.37 2003-02-11 08:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.38 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshCommon.a
INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics\
-I../Mesh -I../Numeric -I../Parser -I../Fltk
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE} ${GUI_INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Context.cpp\
Views.cpp\
diff --git a/DataStr/Makefile b/DataStr/Makefile
index 9a006dcbbcde68979eaacc844a1641df7ff5abec..940aea3b489511bc9cd8764bdb64839ff525a00f 100644
--- a/DataStr/Makefile
+++ b/DataStr/Makefile
@@ -1,10 +1,10 @@
-# $Id: Makefile,v 1.18 2003-02-11 08:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.19 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshDataStr.a
INCLUDE = -I../Common
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = List.cpp \
Malloc.cpp \
diff --git a/Fltk/Main.cpp b/Fltk/Main.cpp
index cdb6cc7ceb36ca835933b5bf1bc8a2518eca1bc9..7f2740aa8043cac94c907a388ac223890a5463b3 100644
--- a/Fltk/Main.cpp
+++ b/Fltk/Main.cpp
@@ -1,4 +1,4 @@
-// $Id: Main.cpp,v 1.41 2003-01-23 20:19:19 geuzaine Exp $
+// $Id: Main.cpp,v 1.42 2003-02-18 05:50:04 geuzaine Exp $
//
// Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
//
@@ -38,6 +38,7 @@
#include "GUI.h"
#include "OpenFile.h"
#include "CommandLine.h"
+#include "Numeric.h"
char yyname[256];
int yyerrorstate;
@@ -173,6 +174,10 @@ int main(int argc, char *argv[]){
Msg(LOG_INFO, "Command line : %s", cmdline);
Msg(LOG_INFO, "-------------------------------------------------------");
+ // Check for buggy obsolete GSL versions
+
+ check_gsl();
+
// Display the GUI immediately to have a quick "a la Windows" launch time
WID->check();
diff --git a/Fltk/Makefile b/Fltk/Makefile
index d8b34f2e50fa00c4196df926816817ed53e90ae3..e02363d0621621a5993fe4d3467a2af0cd0fb9e8 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.40 2003-02-11 08:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.41 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshFltk.a
INCLUDE = -I../Common -I../DataStr -I../Graphics -I../Geo\
-I../Mesh -I../Numeric -I../Parser -I../Fltk -I../Plugin -I../utils
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE} ${GUI_INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Main.cpp \
Message.cpp \
diff --git a/Geo/CAD.cpp b/Geo/CAD.cpp
index a89e3d4e65d5bb576719947850ad27bb4073713d..cc9e8b9eac183362a1b24d3cd5d6fceab8b90e60 100644
--- a/Geo/CAD.cpp
+++ b/Geo/CAD.cpp
@@ -1,4 +1,4 @@
-// $Id: CAD.cpp,v 1.53 2003-01-23 20:19:20 geuzaine Exp $
+// $Id: CAD.cpp,v 1.54 2003-02-18 05:50:04 geuzaine Exp $
//
// Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
//
@@ -29,6 +29,13 @@
#include "CAD.h"
#include "Edge.h"
#include "Context.h"
+
+#if defined(HAVE_GSL)
+#include "gsl_newt.h"
+#include "gsl_brent.h"
+#else
+#include "NR.h"
+#endif
extern Mesh *THEM;
extern Context_T CTX;
@@ -1590,20 +1597,12 @@ void ReplaceAllDuplicates(Mesh *m){
static Curve *CURVE, *CURVE_2;
static Surface *SURFACE;
static Vertex *VERTEX;
-extern void newt(double x[], int n, int *check,
- void (*vecfunc)(int, double [], double []));
double min1d (double (*funct)(double), double *xmin){
double xx, fx, fb, fa, bx, ax;
- 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));
-
- ax=1.e-15; xx=1.e-12;
- mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,funct);
-#define TOL 1.e-08
- return( brent(ax,xx,bx,funct,TOL,xmin) );
+ ax=1.e-15; xx=1.e-12;
+ mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,funct);
+ return(brent(ax,xx,bx,funct,1.e-8,xmin));
}
static void intersectCS (int N, double x[], double res[]){
@@ -1778,7 +1777,7 @@ bool ProjectPointOnSurface (Surface *s, Vertex &p){
while(1){
newt(x,2,&check,projectPS);
vv = InterpolateSurface(s,x[1],x[2],0,0);
- if(x[1] >= UMIN && x[1] <= UMAX && x[2] >=VMIN && x[2] <= VMAX)break;
+ if(x[1] >= UMIN && x[1] <= UMAX && x[2] >=VMIN && x[2] <= VMAX) break;
x[1] = UMIN + (UMAX-UMIN)*((rand() % 10000)/10000.);
x[2] = VMIN + (VMAX-VMIN)*((rand() % 10000)/10000.);
}
@@ -1856,7 +1855,7 @@ bool ProjectPointOnSurface (Surface *s, Vertex *p,double *u, double *v){
l = sqrt(DSQR(vv.Pos.X - p->Pos.X) +
DSQR(vv.Pos.Y - p->Pos.Y) +
DSQR(vv.Pos.Z - p->Pos.Z));
- if(l < 1.e-1)break;
+ if(l < 1.e-1) break;
else {
x[1] = UMIN + (UMAX-UMIN)*((rand() % 10000)/10000.);
x[2] = VMIN + (VMAX-VMIN)*((rand() % 10000)/10000.);
@@ -1909,14 +1908,14 @@ bool IntersectCurves (Curve *c1, Curve *c2,
CURVE = c1;
x[1] = x[2] = 0.0;
newt(x,2,&check,intersectCC);
- if(check)return false;
+ if(check) return false;
v1 = InterpolateCurve(c1,x[1],0);
v2 = InterpolateCurve(c2,x[2],0);
- if(x[1] <= c1->ubeg)return false;
- if(x[1] >= c1->uend)return false;
- if(x[2] <= c2->ubeg)return false;
- if(x[2] >= c2->uend)return false;
- if(fabs(v1.Pos.Z - v2.Pos.Z) > 1.e-08 * CTX.lc)return false;
+ if(x[1] <= c1->ubeg) return false;
+ if(x[1] >= c1->uend) return false;
+ if(x[2] <= c2->ubeg) return false;
+ if(x[2] >= c2->uend) return false;
+ if(fabs(v1.Pos.Z - v2.Pos.Z) > 1.e-08 * CTX.lc) return false;
*v = Create_Vertex(NEWPOINT(), v1.Pos.X,v1.Pos.Y,v1.Pos.Z,v1.lc,x[1]);
Tree_Insert(THEM->Points,v);
DivideCurve(c1,x[1],*v,c11,c12);
diff --git a/Geo/Makefile b/Geo/Makefile
index 4068650388a2cdbbfb3938cce2db936099ff7144..a7b0ffe2b0915d86351839d685368796d50b5391 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.33 2003-02-11 08:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.34 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshGeo.a
-INCLUDE = -I../Common -I../DataStr -I../Geo -I../Mesh -I../Numeric\
+INCLUDE = -I../Common -I../DataStr -I../Geo -I../Mesh -I../Numeric -I../NR\
-I../Parser -I../Fltk
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = CAD.cpp \
DataBase.cpp \
diff --git a/Graphics/Makefile b/Graphics/Makefile
index 96c21ded3b76547aa8336b346b570fbf99e0dff8..057505ee9dc3641722adf55582f579a6cb303dd5 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.44 2003-02-12 16:12:27 geuzaine Exp $
+# $Id: Makefile,v 1.45 2003-02-18 05:50:04 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshGraphics.a
INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics\
-I../Fltk -I../Mesh -I../Numeric -I../Parser
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE} ${GUI_INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Draw.cpp \
Mesh.cpp \
diff --git a/Makefile b/Makefile
index a0818b11ef7df3dd6f77a70ac13569282d9c5baf..c249ca7eb478b1d4c130c4b66736c90ec85c52d0 100644
--- a/Makefile
+++ b/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.259 2003-02-12 22:09:50 geuzaine Exp $
+# $Id: Makefile,v 1.260 2003-02-18 05:50:04 geuzaine Exp $
include variables
@@ -47,7 +47,8 @@ source:
mkdir gmsh-${GMSH_RELEASE}
cd gmsh-${GMSH_RELEASE} && tar zxvf ../gmsh.tgz
rm -f gmsh.tgz
- cd gmsh-${GMSH_RELEASE} && rm -rf CVS */CVS */*/CVS */.globalrc ${GMSH_VERSION_FILE}
+ cd gmsh-${GMSH_RELEASE} && rm -rf CVS */CVS */*/CVS */.globalrc\
+ NR ${GMSH_VERSION_FILE}
tar zcvf gmsh-${GMSH_RELEASE}-source.tgz gmsh-${GMSH_RELEASE}
parser:
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 5d7f3cf7405818e6c926d7b80595231afcda4dbe..791b59611eb1e5b37ac04538b7d68075d459aacd 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.46 2003-02-11 08:54:56 geuzaine Exp $
+# $Id: Makefile,v 1.47 2003-02-18 05:50:05 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshMesh.a
-INCLUDE = -I../Numeric -I../Common -I../DataStr -I../Geo -I../Mesh\
+INCLUDE = -I../Numeric -I../NR -I../Common -I../DataStr -I../Geo -I../Mesh\
-I../Graphics -I../Parser -I../Fltk -I../Triangle
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = 1D_Mesh.cpp \
2D_Mesh.cpp \
@@ -247,7 +247,7 @@ Utils.o: Utils.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
../Geo/CAD.h ../Mesh/Mesh.h ../Mesh/Vertex.h ../Mesh/Simplex.h \
../Mesh/Edge.h ../Geo/ExtrudeParams.h ../Mesh/Metric.h ../Mesh/Matrix.h \
- Mesh.h Interpolation.h ../Numeric/NRUtil.h ../Common/Context.h
+ Mesh.h Interpolation.h ../Common/Context.h
Metric.o: Metric.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Numeric/Numeric.h ../Geo/Geo.h \
diff --git a/Mesh/Utils.cpp b/Mesh/Utils.cpp
index a6dffeec78d86c06a58611d311c725b417bf4852..f692103739bfc2694626839d8e038f3cd8f48ed9 100644
--- a/Mesh/Utils.cpp
+++ b/Mesh/Utils.cpp
@@ -1,4 +1,4 @@
-// $Id: Utils.cpp,v 1.14 2003-01-23 20:19:22 geuzaine Exp $
+// $Id: Utils.cpp,v 1.15 2003-02-18 05:50:05 geuzaine Exp $
//
// Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
//
@@ -25,12 +25,16 @@
#include "CAD.h"
#include "Mesh.h"
#include "Interpolation.h"
-#include "NRUtil.h"
#include "Context.h"
-extern Context_T CTX;
+#if defined(HAVE_GSL)
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_linalg.h>
+#else
+#include "NR.h"
+#endif
-void dsvdcmp(double **a, int m, int n, double w[], double **v);
+extern Context_T CTX;
void direction (Vertex * v1, Vertex * v2, double d[3]){
d[0] = v2->Pos.X - v1->Pos.X;
@@ -58,17 +62,12 @@ void Projette (Vertex * v, double mat[3][3]){
void MeanPlane(List_T *points, Surface *s){
int i, j, min, ndata, na;
- double **U, **V, *W, res[4], ex[3], t1[3], t2[3], svd[3];
+ double res[4], ex[3], t1[3], t2[3], svd[3];
Vertex *v;
double xm=0., ym=0., zm=0.;
ndata = List_Nbr(points);
na = 3;
-
- U = dmatrix(1,ndata,1,na);
- V = dmatrix(1,na,1,na);
- W = dvector(1,na);
-
for (i=0; i<ndata; i++){
List_Read(points, i, &v);
xm += v->Pos.X;
@@ -79,6 +78,36 @@ void MeanPlane(List_T *points, Surface *s){
ym/=(double)ndata;
zm/=(double)ndata;
+#if defined(HAVE_GSL)
+ gsl_matrix *U = gsl_matrix_alloc(ndata, na);
+ gsl_matrix *V = gsl_matrix_alloc(na, na);
+ gsl_vector *W = gsl_vector_alloc(na);
+ gsl_vector *TMPVEC = gsl_vector_alloc(na);
+ for (i=0; i<ndata; i++){
+ List_Read(points, i, &v);
+ gsl_matrix_set(U, i, 0, v->Pos.X-xm);
+ gsl_matrix_set(U, i, 1, v->Pos.Y-ym);
+ gsl_matrix_set(U, i, 2, v->Pos.Z-zm);
+ }
+ gsl_linalg_SV_decomp(U,V,W,TMPVEC);
+ svd[0] = gsl_vector_get(W, 0);
+ svd[1] = gsl_vector_get(W, 1);
+ svd[2] = gsl_vector_get(W, 2);
+ if(fabs(svd[0])<fabs(svd[1]) && fabs(svd[0])<fabs(svd[2])) min=0;
+ else if(fabs(svd[1])<fabs(svd[0]) && fabs(svd[1])<fabs(svd[2])) min=1;
+ else min=2;
+ res[0] = gsl_matrix_get(V, 0, min);
+ res[1] = gsl_matrix_get(V, 1, min);
+ res[2] = gsl_matrix_get(V, 2, min);
+ norme(res);
+ gsl_matrix_free(U);
+ gsl_matrix_free(V);
+ gsl_vector_free(W);
+ gsl_vector_free(TMPVEC);
+#else
+ double **U = dmatrix(1,ndata,1,na);
+ double **V = dmatrix(1,na,1,na);
+ double *W = dvector(1,na);
for (i=0; i<ndata; i++){
List_Read(points, i, &v);
U[i+1][1] = v->Pos.X-xm;
@@ -96,10 +125,10 @@ void MeanPlane(List_T *points, Surface *s){
res[1] = V[2][min];
res[2] = V[3][min];
norme(res);
-
free_dmatrix(U,1,ndata,1,na);
free_dmatrix(V,1,na,1,na);
free_dvector(W,1,na);
+#endif
// check coherence of results for non-plane surfaces
if(s->Typ != MSH_SURF_PLAN){
@@ -233,44 +262,75 @@ void find_bestuv (Surface * s, double X, double Y,
*V = minv;
}
-// Numerical Recipes in C, p. 62
-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);
+void invert_singular_matrix3x3(double _M[3][3], double _I[3][3]){
+ int i, j, k, n = 3;
+ double _T[3][3];
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
- I[i][j] = 0.0 ;
- T[i][j] = 0.0 ;
+ _I[i-1][j-1] = 0.0 ;
+ _T[i-1][j-1] = 0.0 ;
}
}
-#define PREC 1.e-16 //singular value precision
+#if defined(HAVE_GSL)
+ gsl_matrix *M = gsl_matrix_alloc(3, 3);
+ gsl_matrix *V = gsl_matrix_alloc(3, 3);
+ gsl_vector *W = gsl_vector_alloc(3);
+ gsl_vector *TMPVEC = gsl_vector_alloc(3);
+ for(i=1 ; i<=n ; i++){
+ for(j=1 ; j<=n ; j++){
+ gsl_matrix_set(M, i-1, j-1, _M[i-1][j-1]);
+ }
+ }
+ gsl_linalg_SV_decomp(M,V,W,TMPVEC);
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] ;
+ double ww = gsl_vector_get(W, i-1);
+ if(fabs(ww) > 1.e-16){ //singular value precision
+ _T[i-1][j-1] += gsl_matrix_get(M, j-1, i-1) / ww ;
}
}
}
-#undef PREC
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] ;
+ _I[i-1][j-1] += gsl_matrix_get(V, i-1, k-1) * _T[k-1][j-1] ;
}
}
}
-
+ gsl_matrix_free(M);
+ gsl_matrix_free(V);
+ gsl_vector_free(W);
+ gsl_vector_free(TMPVEC);
+#else
+ double **M = dmatrix(1,3,1,3);
+ double **V = dmatrix(1,3,1,3);
+ double *W = dvector(1,3);
+ for(i=1 ; i<=n ; i++){
+ for(j=1 ; j<=n ; j++){
+ M[i][j] = _M[i-1][j-1];
+ }
+ }
+ dsvdcmp(M, n, n, W, V);
+ for(i=1 ; i<=n ; i++){
+ for(j=1 ; j<=n ; j++){
+ if(fabs(W[i]) > 1.e-16){ //singular value precision
+ _T[i-1][j-1] += M[j][i] / W[i] ;
+ }
+ }
+ }
+ for(i=1 ; i<=n ; i++){
+ for(j=1 ; j<=n ; j++){
+ for(k=1 ; k<=n ; k++){
+ _I[i-1][j-1] += V[i][k] * _T[k-1][j-1] ;
+ }
+ }
+ }
+ free_dmatrix(M,1,n,1,n);
free_dmatrix(V,1,n,1,n);
- free_dmatrix(T,1,n,1,n);
free_dvector(W,1,n);
+#endif
}
void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
@@ -278,7 +338,7 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
double Unew,Vnew,err;
int iter;
Vertex D_u,D_v,P;
- double **mat, **jac ;
+ double mat[3][3], jac[3][3] ;
double umin, umax, vmin, vmax;
if (s->Typ == MSH_SURF_NURBS){
@@ -292,9 +352,6 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
umax = vmax = 1.0;
}
- mat = dmatrix(1,3,1,3);
- jac = dmatrix(1,3,1,3);
-
*U = *V = 0.487;
err = 1.0;
iter = 1;
@@ -304,21 +361,21 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
D_u = InterpolateSurface(s, *U, *V, 1, 1);
D_v = InterpolateSurface(s, *U, *V, 1, 2);
- mat[1][1] = D_u.Pos.X;
- mat[1][2] = D_u.Pos.Y;
- mat[1][3] = D_u.Pos.Z;
- mat[2][1] = D_v.Pos.X;
- mat[2][2] = D_v.Pos.Y;
- mat[2][3] = D_v.Pos.Z;
- mat[3][1] = 0.;
- mat[3][2] = 0.;
- mat[3][3] = 0.;
- invert_singular_matrix(mat,3,jac);
-
- Unew = *U + relax * (jac[1][1] * (X-P.Pos.X) + jac[2][1] * (Y-P.Pos.Y) +
- jac[3][1] * (Z-P.Pos.Z)) ;
- Vnew = *V + relax * (jac[1][2] * (X-P.Pos.X) + jac[2][2] * (Y-P.Pos.Y) +
- jac[3][2] * (Z-P.Pos.Z)) ;
+ mat[0][0] = D_u.Pos.X;
+ mat[0][1] = D_u.Pos.Y;
+ mat[0][2] = D_u.Pos.Z;
+ mat[1][0] = D_v.Pos.X;
+ mat[1][1] = D_v.Pos.Y;
+ mat[1][2] = D_v.Pos.Z;
+ mat[2][0] = 0.;
+ mat[2][1] = 0.;
+ mat[2][2] = 0.;
+ invert_singular_matrix3x3(mat, jac);
+
+ Unew = *U + relax *
+ (jac[0][0]*(X-P.Pos.X) + jac[1][0]*(Y-P.Pos.Y) + jac[2][0]*(Z-P.Pos.Z)) ;
+ Vnew = *V + relax *
+ (jac[0][1]*(X-P.Pos.X) + jac[1][1]*(Y-P.Pos.Y) + jac[2][1]*(Z-P.Pos.Z)) ;
err = DSQR(Unew - *U) + DSQR(Vnew - *V) ;
@@ -327,9 +384,6 @@ void XYZtoUV (Surface *s, double X, double Y, double Z, double *U, double *V,
*V = Vnew;
}
- free_dmatrix(mat,1,3,1,3);
- free_dmatrix(jac,1,3,1,3);
-
if(iter == MaxIter ||
(Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin)){
if(relax < 1.e-12)
diff --git a/NR/Makefile b/NR/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..00312ed71a94ccdb99f8006ec00b00994365b71c
--- /dev/null
+++ b/NR/Makefile
@@ -0,0 +1,58 @@
+# $Id: Makefile,v 1.1 2003-02-18 05:50:05 geuzaine Exp $
+
+include ../variables
+
+LIB = ../lib/libGmshNR.a
+INCLUDE = -I../Common -I../DataStr -I../Numeric
+# don't optimize this library: there are some problems with gcc...
+OPTIM = -O0
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
+
+SRC = brent.cpp\
+ dpythag.cpp\
+ dsvdcmp.cpp\
+ fdjac.cpp\
+ fmin.cpp\
+ lnsrch.cpp\
+ lubksb.cpp\
+ ludcmp.cpp\
+ mnbrak.cpp\
+ newt.cpp\
+ nrutil.cpp
+
+OBJ = ${SRC:.cpp=.o}
+
+.SUFFIXES: .o .cpp
+
+${LIB}: ${OBJ}
+ ${AR} ${LIB} ${OBJ}
+ ${RANLIB} ${LIB}
+
+.cpp.o:
+ ${CXX} ${CFLAGS} -c $<
+
+clean:
+ rm -f *.o
+
+depend:
+ (sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
+ ${CXX} -MM ${CFLAGS} ${SRC} \
+ ) >Makefile.new
+ cp Makefile Makefile.bak
+ cp Makefile.new Makefile
+ rm -f Makefile.new
+
+# DO NOT DELETE THIS LINE
+brent.o: brent.cpp nrutil.h ../Numeric/Numeric.h
+dpythag.o: dpythag.cpp nrutil.h ../Numeric/Numeric.h
+dsvdcmp.o: dsvdcmp.cpp nrutil.h ../Numeric/Numeric.h
+fdjac.o: fdjac.cpp nrutil.h ../Numeric/Numeric.h
+fmin.o: fmin.cpp nrutil.h ../Numeric/Numeric.h
+lnsrch.o: lnsrch.cpp nrutil.h ../Numeric/Numeric.h
+lubksb.o: lubksb.cpp
+ludcmp.o: ludcmp.cpp nrutil.h ../Numeric/Numeric.h
+mnbrak.o: mnbrak.cpp nrutil.h ../Numeric/Numeric.h
+newt.o: newt.cpp nrutil.h ../Numeric/Numeric.h
+nrutil.o: nrutil.cpp ../Common/Gmsh.h ../Common/Message.h \
+ ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
+ ../DataStr/avl.h ../DataStr/Tools.h
diff --git a/NR/NR.h b/NR/NR.h
new file mode 100644
index 0000000000000000000000000000000000000000..d31e4ba61aeb7259c772752b3e20f72f5ab6a5ce
--- /dev/null
+++ b/NR/NR.h
@@ -0,0 +1,17 @@
+#ifndef _NR_H_
+#define _NR_H_
+
+#include "Gmsh.h"
+#include "nrutil.h"
+
+/* "public" routines used in Gmsh */
+
+void dsvdcmp (double **a, int m, int n, double w[], double **v);
+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));
+void newt(double x[], int n, int *check,
+ void (*vecfunc)(int, double [], double []));
+
+#endif
diff --git a/NR/brent.cpp b/NR/brent.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8dc82e881318dd5f92334a31cc848c16f91f7309
--- /dev/null
+++ b/NR/brent.cpp
@@ -0,0 +1,77 @@
+#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);
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+double brent(double ax, double bx, double cx, double (*f)(double), double tol,
+ double *xmin)
+{
+ int iter;
+ double a,b,d,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. */
diff --git a/NR/dpythag.cpp b/NR/dpythag.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4585295c43027bfa9f6f326556dd424b25f7e6d4
--- /dev/null
+++ b/NR/dpythag.cpp
@@ -0,0 +1,14 @@
+#include <math.h>
+#define NRANSI
+#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+DSQR(absb/absa));
+ else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+DSQR(absa/absb)));
+}
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/NR/dsvdcmp.cpp b/NR/dsvdcmp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..81317ec7db4c1e6ef3558ec6a3c36f3e7f3a4e45
--- /dev/null
+++ b/NR/dsvdcmp.cpp
@@ -0,0 +1,183 @@
+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+
+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);
+}
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/NR/fdjac.cpp b/NR/fdjac.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7767ff1bc40f7a7d562934d745ef11c72b646a64
--- /dev/null
+++ b/NR/fdjac.cpp
@@ -0,0 +1,29 @@
+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define EPS 1.0e-4
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+void fdjac(int n, double x[], double fvec[], double **df,
+ void (*vecfunc)(int, double [], double []))
+{
+ int i,j;
+ double h,temp,*f;
+
+ f=dvector(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_dvector(f,1,n);
+}
+#undef EPS
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/NR/fmin.cpp b/NR/fmin.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f461aa7a05e4cf06261d3c9fd134bf0d62499ce6
--- /dev/null
+++ b/NR/fmin.cpp
@@ -0,0 +1,20 @@
+#define NRANSI
+#include "nrutil.h"
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+extern int nn;
+extern double *fvec;
+extern void (*nrfuncv)(int n, double v[], double f[]);
+
+double fmin(double x[])
+{
+ int i;
+ double 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. */
diff --git a/NR/lnsrch.cpp b/NR/lnsrch.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e4a5346a5f846b570a7c7f105713933923421b24
--- /dev/null
+++ b/NR/lnsrch.cpp
@@ -0,0 +1,65 @@
+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define ALF 1.0e-4
+#define TOLX 1.0e-7
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],
+ double *f, double stpmax, int *check, double (*func)(double []))
+{
+ int i;
+ double 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. */
diff --git a/NR/lubksb.cpp b/NR/lubksb.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..202c826e2f2b0f8b0f6108b2bf152f6837b3cfa0
--- /dev/null
+++ b/NR/lubksb.cpp
@@ -0,0 +1,23 @@
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+void lubksb(double **a, int n, int *indx, double b[])
+{
+ int i,ii=0,ip,j;
+ double 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. */
diff --git a/NR/ludcmp.cpp b/NR/ludcmp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8f4e0c17fca5bfc5ffb10e6e8df231aa0fb07384
--- /dev/null
+++ b/NR/ludcmp.cpp
@@ -0,0 +1,60 @@
+#include <math.h>
+#define NRANSI
+#include "nrutil.h"
+#define TINY 1.0e-20;
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+void ludcmp(double **a, int n, int *indx, double *d)
+{
+ int i,imax,j,k;
+ double big,dum,sum,temp;
+ double *vv;
+
+ vv=dvector(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_dvector(vv,1,n);
+}
+#undef TINY
+#undef NRANSI
+/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
diff --git a/NR/mnbrak.cpp b/NR/mnbrak.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9aedccc886144c9f1d09ed8a694c65e17cdff822
--- /dev/null
+++ b/NR/mnbrak.cpp
@@ -0,0 +1,67 @@
+#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);
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+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. */
diff --git a/NR/newt.cpp b/NR/newt.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b09482177e7bf5f7d97f80c8b430b3dc561c94a7
--- /dev/null
+++ b/NR/newt.cpp
@@ -0,0 +1,92 @@
+#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
+
+/* This file has been modified for inclusion in Gmsh (float->double) */
+
+int nn;
+double *fvec;
+void (*nrfuncv)(int n, double v[], double f[]);
+#define FREERETURN {free_dvector(fvec,1,n);free_dvector(xold,1,n);\
+ free_dvector(p,1,n);free_dvector(g,1,n);free_dmatrix(fjac,1,n,1,n);\
+ free_ivector(indx,1,n);return;}
+
+void newt(double x[], int n, int *check,
+ void (*vecfunc)(int, double [], double []))
+{
+ void fdjac(int n, double x[], double fvec[], double **df,
+ void (*vecfunc)(int, double [], double []));
+ double fmin(double x[]);
+ void lnsrch(int n, double xold[], double fold, double g[], double p[], double x[],
+ double *f, double stpmax, int *check, double (*func)(double []));
+ void lubksb(double **a, int n, int *indx, double b[]);
+ void ludcmp(double **a, int n, int *indx, double *d);
+ int i,its,j,*indx;
+ double d,den,f,fold,stpmax,sum,temp,test,**fjac,*g,*p,*xold;
+
+ indx=ivector(1,n);
+ fjac=dmatrix(1,n,1,n);
+ g=dvector(1,n);
+ p=dvector(1,n);
+ xold=dvector(1,n);
+ fvec=dvector(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),(double)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. */
diff --git a/NR/nrutil.cpp b/NR/nrutil.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e1066b49c052a38cf948957fff72c875a2cb352e
--- /dev/null
+++ b/NR/nrutil.cpp
@@ -0,0 +1,619 @@
+#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
+
+#include <stdio.h>
+#include <stddef.h>
+#include <stdlib.h>
+#define NR_END 1
+#define FREE_ARG char*
+
+#include "Gmsh.h"
+
+void nrerror(char error_text[])
+/* Numerical Recipes standard error handler */
+{
+ Msg(GERROR, "%s", error_text);
+ /*
+ 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));
+}
+
+#else /* ANSI */
+/* traditional - K&R */
+
+#include <stdio.h>
+#define NR_END 1
+#define FREE_ARG char*
+
+void nrerror(error_text)
+char error_text[];
+/* Numerical Recipes standard error handler */
+{
+ void exit();
+
+ 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(nl,nh)
+long nh,nl;
+/* allocate a float vector with subscript range v[nl..nh] */
+{
+ float *v;
+
+ v=(float *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(float)));
+ if (!v) nrerror("allocation failure in vector()");
+ return v-nl+NR_END;
+}
+
+int *ivector(nl,nh)
+long nh,nl;
+/* allocate an int vector with subscript range v[nl..nh] */
+{
+ int *v;
+
+ v=(int *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(int)));
+ if (!v) nrerror("allocation failure in ivector()");
+ return v-nl+NR_END;
+}
+
+unsigned char *cvector(nl,nh)
+long nh,nl;
+/* allocate an unsigned char vector with subscript range v[nl..nh] */
+{
+ unsigned char *v;
+
+ v=(unsigned char *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
+ if (!v) nrerror("allocation failure in cvector()");
+ return v-nl+NR_END;
+}
+
+unsigned long *lvector(nl,nh)
+long nh,nl;
+/* allocate an unsigned long vector with subscript range v[nl..nh] */
+{
+ unsigned long *v;
+
+ v=(unsigned long *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(long)));
+ if (!v) nrerror("allocation failure in lvector()");
+ return v-nl+NR_END;
+}
+
+double *dvector(nl,nh)
+long nh,nl;
+/* allocate a double vector with subscript range v[nl..nh] */
+{
+ double *v;
+
+ v=(double *)malloc((unsigned int) ((nh-nl+1+NR_END)*sizeof(double)));
+ if (!v) nrerror("allocation failure in dvector()");
+ return v-nl+NR_END;
+}
+
+float **matrix(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* 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((unsigned int)((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((unsigned int)((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(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* 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((unsigned int)((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((unsigned int)((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(nrl,nrh,ncl,nch)
+long nch,ncl,nrh,nrl;
+/* 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((unsigned int)((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((unsigned int)((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(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
+float **a;
+long newcl,newrl,oldch,oldcl,oldrh,oldrl;
+/* 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((unsigned int) ((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(a,nrl,nrh,ncl,nch)
+float *a;
+long nch,ncl,nrh,nrl;
+/* 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((unsigned int) ((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(nrl,nrh,ncl,nch,ndl,ndh)
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* 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((unsigned int)((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((unsigned int)((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((unsigned int)((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(v,nl,nh)
+float *v;
+long nh,nl;
+/* free a float vector allocated with vector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_ivector(v,nl,nh)
+int *v;
+long nh,nl;
+/* free an int vector allocated with ivector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_cvector(v,nl,nh)
+long nh,nl;
+unsigned char *v;
+/* free an unsigned char vector allocated with cvector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_lvector(v,nl,nh)
+long nh,nl;
+unsigned long *v;
+/* free an unsigned long vector allocated with lvector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_dvector(v,nl,nh)
+double *v;
+long nh,nl;
+/* free a double vector allocated with dvector() */
+{
+ free((FREE_ARG) (v+nl-NR_END));
+}
+
+void free_matrix(m,nrl,nrh,ncl,nch)
+float **m;
+long nch,ncl,nrh,nrl;
+/* 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(m,nrl,nrh,ncl,nch)
+double **m;
+long nch,ncl,nrh,nrl;
+/* 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(m,nrl,nrh,ncl,nch)
+int **m;
+long nch,ncl,nrh,nrl;
+/* 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(b,nrl,nrh,ncl,nch)
+float **b;
+long nch,ncl,nrh,nrl;
+/* free a submatrix allocated by submatrix() */
+{
+ free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_convert_matrix(b,nrl,nrh,ncl,nch)
+float **b;
+long nch,ncl,nrh,nrl;
+/* free a matrix allocated by convert_matrix() */
+{
+ free((FREE_ARG) (b+nrl-NR_END));
+}
+
+void free_f3tensor(t,nrl,nrh,ncl,nch,ndl,ndh)
+float ***t;
+long nch,ncl,ndh,ndl,nrh,nrl;
+/* 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));
+}
+
+#endif /* ANSI */
diff --git a/NR/nrutil.h b/NR/nrutil.h
new file mode 100644
index 0000000000000000000000000000000000000000..ddfa7d0b51691cfbe767e1ba75d0731e6fea7a3b
--- /dev/null
+++ b/NR/nrutil.h
@@ -0,0 +1,108 @@
+#ifndef _NR_UTILS_H_
+#define _NR_UTILS_H_
+
+/* This file has been modified for inclusion in Gmsh */
+
+/* Gmsh: */
+#include "Numeric.h"
+
+/* Gmsh:
+static float sqrarg;
+#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
+
+static double dsqrarg;
+#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)
+
+static double dmaxarg1,dmaxarg2;
+#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
+ (dmaxarg1) : (dmaxarg2))
+
+static double dminarg1,dminarg2;
+#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
+ (dminarg1) : (dminarg2))
+
+static float maxarg1,maxarg2;
+#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
+ (maxarg1) : (maxarg2))
+
+static float minarg1,minarg2;
+#define FMIN(a,b) (minarg1=(a),minarg2=(b),(minarg1) < (minarg2) ?\
+ (minarg1) : (minarg2))
+
+static long lmaxarg1,lmaxarg2;
+#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
+ (lmaxarg1) : (lmaxarg2))
+
+static long lminarg1,lminarg2;
+#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
+ (lminarg1) : (lminarg2))
+
+static int imaxarg1,imaxarg2;
+#define IMAX(a,b) (imaxarg1=(a),imaxarg2=(b),(imaxarg1) > (imaxarg2) ?\
+ (imaxarg1) : (imaxarg2))
+
+static int iminarg1,iminarg2;
+#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
+ (iminarg1) : (iminarg2))
+
+#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
+*/
+
+#if defined(__STDC__) || defined(ANSI) || defined(NRANSI) /* ANSI */
+
+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);
+
+#else /* ANSI */
+/* traditional - K&R */
+
+void nrerror();
+float *vector();
+float **matrix();
+float **submatrix();
+float **convert_matrix();
+float ***f3tensor();
+double *dvector();
+double **dmatrix();
+int *ivector();
+int **imatrix();
+unsigned char *cvector();
+unsigned long *lvector();
+void free_vector();
+void free_dvector();
+void free_ivector();
+void free_cvector();
+void free_lvector();
+void free_matrix();
+void free_submatrix();
+void free_convert_matrix();
+void free_dmatrix();
+void free_imatrix();
+void free_f3tensor();
+
+#endif /* ANSI */
+
+#endif /* _NR_UTILS_H_ */
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 5e8d725c92d3193e59b25395ac386ed35c198b59..d7214b146254744f0ae62fb4366c8ca3405115f1 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -1,16 +1,14 @@
-# $Id: Makefile,v 1.6 2003-02-11 08:54:56 geuzaine Exp $
+# $Id: Makefile,v 1.7 2003-02-18 05:50:05 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshNumeric.a
INCLUDE = -I../Common -I../DataStr
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Numeric.cpp\
- Newton.cpp \
- Opti.cpp \
- SVD.cpp \
- NRUtil.cpp
+ gsl_newt.cpp\
+ gsl_brent.cpp
OBJ = ${SRC:.cpp=.o}
@@ -37,12 +35,5 @@ depend:
# DO NOT DELETE THIS LINE
Numeric.o: Numeric.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h Numeric.h
-Newton.o: Newton.cpp NRUtil.h Numeric.h
-Opti.o: Opti.cpp NRUtil.h Numeric.h
-SVD.o: SVD.cpp ../Common/Gmsh.h ../Common/Message.h ../DataStr/Malloc.h \
- ../DataStr/List.h ../DataStr/Tree.h ../DataStr/avl.h ../DataStr/Tools.h \
- NRUtil.h Numeric.h
-NRUtil.o: NRUtil.cpp ../Common/Gmsh.h ../Common/Message.h \
- ../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
- ../DataStr/avl.h ../DataStr/Tools.h
+ ../DataStr/avl.h ../DataStr/Tools.h Numeric.h \
+ /usr/local/include/gsl/gsl_version.h /usr/local/include/gsl/gsl_types.h
diff --git a/Numeric/NRUtil.cpp b/Numeric/NRUtil.cpp
deleted file mode 100644
index a8b2b6cf73d611c17f75acaab772a7ca6e63b59c..0000000000000000000000000000000000000000
--- a/Numeric/NRUtil.cpp
+++ /dev/null
@@ -1,299 +0,0 @@
-// $Id: NRUtil.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
-
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
-
-#include <stdio.h>
-#include <stddef.h>
-#include <stdlib.h>
-#include <malloc.h>
-
-#include "Gmsh.h"
-
-
-#define NR_END 1
-#define FREE_ARG char*
-
-void nrerror(char error_text[])
-/* Numerical Recipes standard error handler */
-{
- Msg(GERROR, "%s", error_text);
- /*
- 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));
-}
-
diff --git a/Numeric/NRUtil.h b/Numeric/NRUtil.h
deleted file mode 100644
index 191fc3528e8f13540eedf4f2cb2c706ec1262d69..0000000000000000000000000000000000000000
--- a/Numeric/NRUtil.h
+++ /dev/null
@@ -1,35 +0,0 @@
-#ifndef _NR_UTILS_H_
-#define _NR_UTILS_H_
-
-#include "Numeric.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_ */
diff --git a/Numeric/Newton.cpp b/Numeric/Newton.cpp
deleted file mode 100644
index 22eaaadc0a560b17e6390380a00cbc763bcf8d13..0000000000000000000000000000000000000000
--- a/Numeric/Newton.cpp
+++ /dev/null
@@ -1,329 +0,0 @@
-// $Id: Newton.cpp,v 1.2 2002-06-03 21:31:11 geuzaine Exp $
-
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
-
-#include <math.h>
-
-#define NRANSI
-#include "NRUtil.h"
-
-#define SQR(a) ((a)*(a))
-
-// global variables
-
-int nn;
-double *fvec;
-void (*nrfuncv) (int n, double v[], double f[]);
-
-// fmin.cpp
-
-double
-fmin (double x[])
-{
- int i;
- double sum;
-
- (*nrfuncv) (nn, x, fvec);
- for (sum = 0.0, i = 1; i <= nn; i++)
- sum += SQR (fvec[i]);
- return 0.5 * sum;
-}
-
-// newt.cpp
-
-#define MAXITS 200
-#define TOLF 1.0e-4 /* 1.0e-4 */
-#define TOLMIN 1.0e-6 /* 1.0e-6 */
-#define TOLX 1.0e-7 /* 1.0e-7 */
-#define STPMX 100.0
-#define FREERETURN {free_dvector(fvec,1,n);free_dvector(xold,1,n);\
- free_dvector(p,1,n);free_dvector(g,1,n);free_dmatrix(fjac,1,n,1,n);\
- free_ivector(indx,1,n);return;}
-
-void
-newt (double x[], int n, int *check,
- void (*vecfunc) (int, double[], double[]))
-{
- void fdjac (int n, double x[], double fvec[], double **df,
- void (*vecfunc) (int, double[], double[]));
- double fmin (double x[]);
- void lnsrch (int n, double xold[], double fold, double g[], double p[],
- double x[], double *f, double stpmax, int *check,
- double (*func) (double[]));
- void lubksb (double **a, int n, int *indx, double b[]);
- void ludcmp (double **a, int n, int *indx, double *d);
- int i, its, j, *indx;
- double d, den, f, fold, stpmax, sum, temp, test, **fjac, *g, *p, *xold;
-
- indx = ivector (1, n);
- fjac = dmatrix (1, n, 1, n);
- g = dvector (1, n);
- p = dvector (1, n);
- xold = dvector (1, n);
- fvec = dvector (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 * MAX (sqrt (sum), (double) n);
- for (its = 1; its <= MAXITS; its++) {
- fdjac (n, x, fvec, fjac, vecfunc);
- //vecjac (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 = MAX (f, 0.5 * n);
- for (i = 1; i <= n; i++) {
- temp = fabs (g[i]) * MAX (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])) / MAX (fabs (x[i]), 1.0);
- if (temp > test)
- test = temp;
- }
- if (test < TOLX)
- FREERETURN;
- }
-
- //Msg(WARNING, "MAXITS exceeded in newt");
- *check=1;
-}
-
-#undef MAXITS
-#undef TOLF
-#undef TOLMIN
-#undef TOLX
-#undef STPMX
-#undef FREERETURN
-
-
-// lnsrch.cpp
-
-#define ALF 1.0e-4 /* 1.0e-4 */
-#define TOLX 1.0e-7 /* 1.0e-7 */
-
-void
-lnsrch (int n, double xold[], double fold, double g[], double p[], double x[],
- double *f, double stpmax, int *check, double (*func) (double[]))
-{
- int i;
- double 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]) / MAX (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 = MAX (tmplam, 0.1 * alam);
- }
-}
-
-#undef ALF
-#undef TOLX
-
-
-
-// fdjac.cpp
-
-#define EPS 1.0e-4 /* 1.0e-4 */
-
-void
-fdjac (int n, double x[], double fvec[], double **df,
- void (*vecfunc) (int, double[], double[]))
-{
- int i, j;
- double h, temp, *f;
-
- f = dvector (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_dvector (f, 1, n);
-}
-
-#undef EPS
-
-
-// ludcmp.cpp
-
-#define TINY 1.0e-20;
-
-void
-ludcmp (double **a, int n, int *indx, double *d)
-{
- int i, imax, j, k;
- double big, dum, sum, temp;
- double *vv;
-
- vv = dvector (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_dvector (vv, 1, n);
-}
-
-#undef TINY
-
-
-// lubksb.cpp
-
-void
-lubksb (double **a, int n, int *indx, double b[])
-{
- int i, ii = 0, ip, j;
- double 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];
- }
-}
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 33262d73f13f6e64738f862e91c5585c03dbde10..23db7ae99198b90a3c7b2b3bd42e408e5fddec6f 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.4 2003-01-23 20:19:23 geuzaine Exp $
+// $Id: Numeric.cpp,v 1.5 2003-02-18 05:50:06 geuzaine Exp $
//
// Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
//
@@ -19,11 +19,37 @@
//
// Please report all bugs and problems to "gmsh@geuz.org".
+// this file should contain only purely numerical routines (that do
+// not depend on any Gmsh structures)
+
#include "Gmsh.h"
#include "Numeric.h"
-// this file should contain only purely numerical routines (that do
-// not depend on any Gmsh structures)
+// Check GSL version. We need at least 1.2, since all versions >=
+// 1.1.1 have a buggy SVD routine (caused an infinite loop--fixed on
+// Sun Jun 16 11:45:29 2002 in GSL's cvs tree). We check this at run
+// time since Gmsh could be distributed with the gsl could be
+// dynamically linked.
+
+#if defined(HAVE_GSL)
+#include <gsl/gsl_version.h>
+int check_gsl(){
+ int major, minor;
+ sscanf(gsl_version, "%d.%d", &major, &minor);
+ if(major < 1 || (major == 1 && minor < 2)){
+ Msg(FATAL1, "Your GSL version (%d.%d.X) has a bug in the singular value", major, minor);
+ Msg(FATAL3, "decomposition code. Please upgrade to version 1.2 or above.");
+ return 0;
+ }
+ return 1;
+}
+#else
+int check_gsl(){
+ return 1;
+}
+#endif
+
+// How ? GSL_VERSION is a string "maj.min" or "maj.min.patch"
double myatan2 (double a, double b){
if (a == 0.0 && b == 0)
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index c5bb896a561bd38ae00510c485335f3ae9b24361..25e2136229415b0e07364067834410236b704cd5 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -31,6 +31,7 @@
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)<(b))?(b):(a))
#define SQR(a) ((a)*(a))
+#define SIGN(a,b)((b) >= 0.0 ? fabs(a) : -fabs(a))
#define IMIN MIN
#define LMIN MIN
@@ -43,7 +44,7 @@
#define DMAX MAX
#define DSQR SQR
-#define FSQR SQU
+#define FSQR SQR
#define THRESHOLD(a,b,c) (((a)>(c))?(c):((a)<(b)?(b):(a)))
@@ -53,6 +54,8 @@
#define Succ(x) ((x)->next)
#define square(x) ((x)*(x))
+int check_gsl();
+
double myatan2 (double a, double b);
double myasin (double a);
void prodve (double a[3], double b[3], double c[3]);
@@ -70,4 +73,5 @@ double InterpolateIso(double *X, double *Y, double *Z,
double *XI, double *YI ,double *ZI);
void gradSimplex (double *x, double *y, double *z, double *v, double *grad);
+
#endif
diff --git a/Numeric/Opti.cpp b/Numeric/Opti.cpp
deleted file mode 100644
index a4c9a268e1d3adf2a7dc6276e613a2649158b41a..0000000000000000000000000000000000000000
--- a/Numeric/Opti.cpp
+++ /dev/null
@@ -1,181 +0,0 @@
-// $Id: Opti.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
-
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
-
-#include <math.h>
-#define NRANSI
-#include "NRUtil.h"
-
-// mnbrak
-
-#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
-
-
-// brent
-
-#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
diff --git a/Numeric/SVD.cpp b/Numeric/SVD.cpp
deleted file mode 100644
index 7fa82f79f7e3b509a3f333c8c3a5a5e50bc7d6c3..0000000000000000000000000000000000000000
--- a/Numeric/SVD.cpp
+++ /dev/null
@@ -1,217 +0,0 @@
-// $Id: SVD.cpp,v 1.1 2002-05-18 21:34:29 geuzaine Exp $
-
-/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
-
-#include "Gmsh.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 == 1000){
- Msg(GERROR, "SVD decomposition: no convergence in 1000 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);
-}
-
-void dsvbksb(double **u, double w[], double **v, int m, int n, double b[], double x[])
-{
- int jj,j,i;
- double s,*tmp;
-
- tmp=dvector(1,n);
- for (j=1;j<=n;j++) {
- s=0.0;
- if (w[j]) {
- for (i=1;i<=m;i++) s += u[i][j]*b[i];
- s /= w[j];
- }
- tmp[j]=s;
- }
- for (j=1;j<=n;j++) {
- s=0.0;
- for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
- x[j]=s;
- }
- free_dvector(tmp,1,n);
-}
diff --git a/Numeric/gsl_brent.cpp b/Numeric/gsl_brent.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cc06c81c869a763463537eec90fd57ad2d7b5616
--- /dev/null
+++ b/Numeric/gsl_brent.cpp
@@ -0,0 +1,138 @@
+// $Id: gsl_brent.cpp,v 1.1 2003-02-18 05:50:07 geuzaine Exp $
+//
+// Copyright (C) 1997 - 2003 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".
+
+// This implements a Newton method using the GSL.
+//
+// Author: Nicolas Tardiey <nicolas.tardieu@der.edf.fr>
+
+#if defined(HAVE_GSL)
+
+#include "Gmsh.h"
+#include "Numeric.h"
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_min.h>
+
+static double (*nrfunc)(double);
+
+double fn1 (double x, void * params){
+ return nrfunc(x);
+}
+
+#define MAXITER 100
+
+double brent(double ax, double xx, double bx,
+ double (*f)(double), double tol, double *xmin){
+ int status;
+ int iter = 0;
+ double m = xx, a = ax, b = bx;
+ const gsl_min_fminimizer_type *T;
+ gsl_min_fminimizer *s;
+ gsl_function F;
+
+ nrfunc = f;
+
+ F.function = &fn1;
+ F.params = 0;
+
+ T = gsl_min_fminimizer_brent;
+ s = gsl_min_fminimizer_alloc(T);
+ gsl_min_fminimizer_set(s, &F, m, a, b);
+
+ do{
+ iter++;
+ status = gsl_min_fminimizer_iterate(s);
+
+ m = gsl_min_fminimizer_x_minimum(s);
+ a = gsl_min_fminimizer_x_lower(s);
+ b = gsl_min_fminimizer_x_upper(s);
+
+ status = gsl_min_test_interval(a, b, 0.001, 0.0);
+
+ //if(status == GSL_SUCCESS) printf ("Converged\n");
+ }
+ while(status == GSL_CONTINUE && iter < MAXITER);
+
+ *xmin = m;
+ return fn1(m, NULL);
+}
+
+#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);
+ }
+}
+
+#endif
diff --git a/Numeric/gsl_brent.h b/Numeric/gsl_brent.h
new file mode 100644
index 0000000000000000000000000000000000000000..caa02ad7eff0540b76703f8d36dba73c948f6d9f
--- /dev/null
+++ b/Numeric/gsl_brent.h
@@ -0,0 +1,28 @@
+#ifndef _GSL_BRENT_H_
+#define _GSL_BRENT_H_
+
+// Copyright (C) 1997 - 2003 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".
+
+double brent(double ax, double xx, double bx,
+ double (*f)(double), double tol, double *xmin);
+void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
+ double (*func)(double));
+
+#endif
diff --git a/Numeric/gsl_newt.cpp b/Numeric/gsl_newt.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..966a4febfa132b3297dc6e9f4992e5a52040b1a1
--- /dev/null
+++ b/Numeric/gsl_newt.cpp
@@ -0,0 +1,103 @@
+// $Id: gsl_newt.cpp,v 1.1 2003-02-18 05:50:07 geuzaine Exp $
+//
+// Copyright (C) 1997 - 2003 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".
+
+// This implements a Newton method using the GSL.
+//
+// Author: Nicolas Tardiey <nicolas.tardieu@der.edf.fr>
+
+#if defined(HAVE_GSL)
+
+#include "Gmsh.h"
+#include "Numeric.h"
+
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_multiroots.h>
+#include <gsl/gsl_linalg.h>
+
+#define MAX_DIM_NEWT 10
+#define MAXITER 10000
+#define PREC 1.e-8
+
+int gsl_dim;
+double gsl_u[MAX_DIM_NEWT], gsl_v[MAX_DIM_NEWT];
+static void (*nrfunc)(int n, double x[],double y[]);
+struct gsl_dummy{;};
+
+void convert_vector_from_gsl(const gsl_vector *gv, double * v){
+ int i, m;
+ m=gv->size;
+ for (i=0;i<m;i++){
+ v[i]=gsl_vector_get(gv,i);
+ }
+}
+
+void convert_vector_to_gsl(double *v, int n, gsl_vector *gv){
+ int i;
+ for (i=0;i<n;i++){
+ gsl_vector_set (gv, i, v[i]);
+ }
+}
+
+int gslfunc(const gsl_vector *xx, void *params, gsl_vector *f){
+ convert_vector_from_gsl(xx,gsl_u);
+ (*nrfunc)(gsl_dim,gsl_u,gsl_v);
+ convert_vector_to_gsl(gsl_v,gsl_dim,f);
+ return GSL_SUCCESS;
+}
+
+void newt(double x[], int n, int *check, void (*func)(int, double [],double [])){
+ const gsl_multiroot_fsolver_type *T;
+ gsl_multiroot_fsolver *s;
+ int status;
+ size_t i, iter = 0;
+ struct gsl_dummy p = {};
+ gsl_multiroot_function f = {&gslfunc, n, &p};
+ gsl_vector *xx = gsl_vector_alloc (n);
+
+ if(n > MAX_DIM_NEWT) Msg(FATAL, "Maximum Newton dimension exceeded\n");
+ gsl_dim = n;
+
+ nrfunc = func;
+ convert_vector_to_gsl(x,n,xx);
+
+ T = gsl_multiroot_fsolver_hybrid;
+ s = gsl_multiroot_fsolver_alloc(T, 2);
+ gsl_multiroot_fsolver_set(s, &f, xx);
+
+ do{
+ iter++;
+ status = gsl_multiroot_fsolver_iterate(s);
+ if(status) break; // solver problem
+ status = gsl_multiroot_test_residual(s->f, n*PREC);
+ }
+ while(status == GSL_CONTINUE && iter < MAXITER);
+
+ if (status == GSL_CONTINUE)
+ *check=1;
+ else
+ *check=0;
+
+ gsl_multiroot_fsolver_free(s);
+ gsl_vector_free(xx);
+}
+
+#endif
diff --git a/Numeric/gsl_newt.h b/Numeric/gsl_newt.h
new file mode 100644
index 0000000000000000000000000000000000000000..7539a119418a3fdd4a712fe93816ed28796e2817
--- /dev/null
+++ b/Numeric/gsl_newt.h
@@ -0,0 +1,25 @@
+#ifndef _GSL_NEWT_H_
+#define _GSL_NEWT_H_
+
+// Copyright (C) 1997 - 2003 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".
+
+void newt(double x[], int n, int *check, void (*func)(int, double [],double []));
+
+#endif
diff --git a/Parallel/Makefile b/Parallel/Makefile
index ce471802ea10c4a17a233b1d78d36968156af2bb..ebe509ff534ec696e91974d75b9f40f63ffe0260 100644
--- a/Parallel/Makefile
+++ b/Parallel/Makefile
@@ -1,10 +1,10 @@
-# $Id: Makefile,v 1.12 2003-02-11 08:54:57 geuzaine Exp $
+# $Id: Makefile,v 1.13 2003-02-18 05:50:07 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshParallel.a
INCLUDE =
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = ParUtil.cpp
diff --git a/Parser/Makefile b/Parser/Makefile
index 972fd7651fc76dd52215eee5453c311e81cf5f9a..83a359a6fef6a83edb6651db2869455414368461 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.43 2003-02-11 08:54:57 geuzaine Exp $
+# $Id: Makefile,v 1.44 2003-02-18 05:50:09 geuzaine Exp $
include ../variables
@@ -8,7 +8,7 @@ LEX = flex
LIB = ../lib/libGmshParser.a
INCLUDE = -I../Common -I../DataStr -I../Geo -I../Graphics\
-I../Mesh -I../Numeric -I../Fltk -I../Plugin -I../Parallel
-CFLAGS = ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE} ${GUI_INCLUDE}
+CFLAGS = ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Gmsh.yy.cpp\
Gmsh.tab.cpp\
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 38ee61054d4365885a54624c42f481e84b950a3d..e31edce6d0e7f36a02461e3066a0eb49fd33a064 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -1,11 +1,11 @@
-# $Id: Makefile,v 1.30 2003-02-11 08:54:57 geuzaine Exp $
+# $Id: Makefile,v 1.31 2003-02-18 05:50:09 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshPlugin.a
INCLUDE = -I../Common -I../Graphics -I../DataStr -I../Geo\
-I../Mesh -I../Numeric -I../Triangle
-CFLAGS = -DMPICH_SKIP_MPICXX ${OPT_FLAGS} ${OS_FLAGS} ${VERSION_FLAGS} ${INCLUDE} ${GUI_INCLUDE}
+CFLAGS = -DMPICH_SKIP_MPICXX ${OPTIM} ${FLAGS} ${INCLUDE}
SRC = Plugin.cpp\
LevelsetPlugin.cpp\
diff --git a/Triangle/Makefile b/Triangle/Makefile
index 16bfecef573105ed91f8dbd34b58a096dccc2b0a..f1c2d09771914f784a467b7116940b841288cb7d 100644
--- a/Triangle/Makefile
+++ b/Triangle/Makefile
@@ -1,11 +1,12 @@
-# $Id: Makefile,v 1.12 2003-02-16 19:17:36 geuzaine Exp $
+# $Id: Makefile,v 1.13 2003-02-18 05:50:10 geuzaine Exp $
include ../variables
LIB = ../lib/libGmshTriangle.a
-# Do not optimize: it crashes badly on Linux
-CFLAGS = -O0 ${OS_FLAGS} -DTRILIBRARY
+# Do not optimize triangle: it crashes badly on Linux
+OPTIM = -O0
+CFLAGS = ${OPTIM} ${FLAGS} -DTRILIBRARY
SRC = triangle.c
diff --git a/configure b/configure
index d72dac04b31d0e1ea56f15431d7ba2c6d866c1d6..7f5a791e7589398e012e1b731fc6e1b5fe038e3b 100755
--- a/configure
+++ b/configure
@@ -833,6 +833,7 @@ Optional Features:
--enable-gui build the graphical user interface (default=yes)
--enable-parallel enable parallel version (default=no)
--enable-triangle compile Triangle if available (default=yes)
+ --enable-gsl use GSL as numerical toolkit (default=yes)
Optional Packages:
--with-PACKAGE[=ARG] use PACKAGE [ARG=yes]
@@ -1240,6 +1241,11 @@ fi;
if test "${enable_triangle+set}" = set; then
enableval="$enable_triangle"
+fi;
+# Check whether --enable-gsl or --disable-gsl was given.
+if test "${enable_gsl+set}" = set; then
+ enableval="$enable_gsl"
+
fi;
UNAME=`uname`
@@ -1248,7 +1254,6 @@ if test "x${UNAME}" = "xIRIX64"; then
UNAME="IRIX"
fi
-
ac_ext=c
ac_cpp='$CPP $CPPFLAGS'
ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5'
@@ -2332,6 +2337,9 @@ ac_link='$CC -o conftest$ac_exeext $CFLAGS $CPPFLAGS $LDFLAGS conftest.$ac_ext $
ac_compiler_gnu=$ac_cv_c_compiler_gnu
+FLAGS=""
+OPTIM="${CXXFLAGS}"
+
ac_ext=c
ac_cpp='$CPP $CPPFLAGS'
ac_compile='$CC -c $CFLAGS $CPPFLAGS conftest.$ac_ext >&5'
@@ -2654,6 +2662,67 @@ fi
+echo "$as_me:$LINENO: checking for main in -lm" >&5
+echo $ECHO_N "checking for main in -lm... $ECHO_C" >&6
+if test "${ac_cv_lib_m_main+set}" = set; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+ ac_check_lib_save_LIBS=$LIBS
+LIBS="-lm $LIBS"
+cat >conftest.$ac_ext <<_ACEOF
+#line $LINENO "configure"
+#include "confdefs.h"
+
+
+#ifdef F77_DUMMY_MAIN
+# ifdef __cplusplus
+ extern "C"
+# endif
+ int F77_DUMMY_MAIN() { return 1; }
+#endif
+int
+main ()
+{
+main ();
+ ;
+ return 0;
+}
+_ACEOF
+rm -f conftest.$ac_objext conftest$ac_exeext
+if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
+ (eval $ac_link) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); } &&
+ { ac_try='test -s conftest$ac_exeext'
+ { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+ (eval $ac_try) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); }; }; then
+ ac_cv_lib_m_main=yes
+else
+ echo "$as_me: failed program was:" >&5
+cat conftest.$ac_ext >&5
+ac_cv_lib_m_main=no
+fi
+rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
+LIBS=$ac_check_lib_save_LIBS
+fi
+echo "$as_me:$LINENO: result: $ac_cv_lib_m_main" >&5
+echo "${ECHO_T}$ac_cv_lib_m_main" >&6
+if test $ac_cv_lib_m_main = yes; then
+ cat >>confdefs.h <<_ACEOF
+#define HAVE_LIBM 1
+_ACEOF
+
+ LIBS="-lm $LIBS"
+
+fi
+
+
+
+
if test "x${AR}" = "x:"; then
{ { echo "$as_me:$LINENO: error: Could not find the library archiver" >&5
echo "$as_me: error: Could not find the library archiver" >&2;}
@@ -2667,7 +2736,7 @@ if test "x$enable_gui" != "xno"; then
GMSH_DIRS="${GMSH_DIRS} Graphics Fltk"
GMSH_LIBS="-Llib -lGmshFltk -lGmshParser -lGmshGraphics -lGmshMesh -lGmshGeo -lGmshNumeric -lGmshCommon -lGmshDataStr -lGmshPlugin -lGmshParallel"
- VERSION_FLAGS="-DHAVE_FLTK"
+ FLAGS="${FLAGS} -DHAVE_FLTK"
if test "x${FLTK_PREFIX}" != "x" ; then
# Extract the first word of "fltk-config", so it can be a program name with args.
@@ -2710,8 +2779,8 @@ else
echo "${ECHO_T}no" >&6
fi
- GUI_LIBS="-L${FLTK_PREFIX}/lib"
- GUI_INCLUDE="-I${FLTK_PREFIX}"
+ GMSH_LIBS="${GMSH_LIBS} -L${FLTK_PREFIX}/lib"
+ FLAGS="${FLAGS} -I${FLTK_PREFIX}"
else
# Extract the first word of "fltk-config", so it can be a program name with args.
set dummy fltk-config; ac_word=$2
@@ -2758,329 +2827,201 @@ fi
echo "$as_me: error: Could not find fltk-config. Try --with-fltk-prefix?" >&2;}
{ (exit 1); exit 1; }; }
fi
- GUI_LIBS="${GUI_LIBS} `$FLTKCONFIG --use-gl --use-images --ldflags`"
- GUI_INCLUDE="${GUI_INCLUDE} `$FLTKCONFIG --use-gl --use-images --cxxflags`"
-
- if test "x${JPEG_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} `$FLTKCONFIG --use-gl --use-images --ldflags`"
+ FLAGS="${FLAGS} `$FLTKCONFIG --use-gl --use-images --cxxflags`"
-echo "$as_me:$LINENO: checking for ANSI C header files" >&5
-echo $ECHO_N "checking for ANSI C header files... $ECHO_C" >&6
-if test "${ac_cv_header_stdc+set}" = set; then
+ if test "x${JPEG_PREFIX}" != "x"; then
+ LDFLAGS="-L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib ${LDFLAGS}"
+ fi
+ echo "$as_me:$LINENO: checking for main in -ljpeg" >&5
+echo $ECHO_N "checking for main in -ljpeg... $ECHO_C" >&6
+if test "${ac_cv_lib_jpeg_main+set}" = set; then
echo $ECHO_N "(cached) $ECHO_C" >&6
else
- cat >conftest.$ac_ext <<_ACEOF
+ ac_check_lib_save_LIBS=$LIBS
+LIBS="-ljpeg $LIBS"
+cat >conftest.$ac_ext <<_ACEOF
#line $LINENO "configure"
#include "confdefs.h"
-#include <stdlib.h>
-#include <stdarg.h>
-#include <string.h>
-#include <float.h>
+
+#ifdef F77_DUMMY_MAIN
+# ifdef __cplusplus
+ extern "C"
+# endif
+ int F77_DUMMY_MAIN() { return 1; }
+#endif
+int
+main ()
+{
+main ();
+ ;
+ return 0;
+}
_ACEOF
-if { (eval echo "$as_me:$LINENO: \"$ac_cpp conftest.$ac_ext\"") >&5
- (eval $ac_cpp conftest.$ac_ext) 2>conftest.er1
+rm -f conftest.$ac_objext conftest$ac_exeext
+if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
+ (eval $ac_link) 2>&5
ac_status=$?
- egrep -v '^ *\+' conftest.er1 >conftest.err
- rm -f conftest.er1
- cat conftest.err >&5
echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } >/dev/null; then
- if test -s conftest.err; then
- ac_cpp_err=$ac_c_preproc_warn_flag
- else
- ac_cpp_err=
- fi
+ (exit $ac_status); } &&
+ { ac_try='test -s conftest$ac_exeext'
+ { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+ (eval $ac_try) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); }; }; then
+ ac_cv_lib_jpeg_main=yes
else
- ac_cpp_err=yes
+ echo "$as_me: failed program was:" >&5
+cat conftest.$ac_ext >&5
+ac_cv_lib_jpeg_main=no
fi
-if test -z "$ac_cpp_err"; then
- ac_cv_header_stdc=yes
+rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
+LIBS=$ac_check_lib_save_LIBS
+fi
+echo "$as_me:$LINENO: result: $ac_cv_lib_jpeg_main" >&5
+echo "${ECHO_T}$ac_cv_lib_jpeg_main" >&6
+if test $ac_cv_lib_jpeg_main = yes; then
+ JPEG="yes"
else
- echo "$as_me: failed program was:" >&5
- cat conftest.$ac_ext >&5
- ac_cv_header_stdc=no
+ JPEG="no"
fi
-rm -f conftest.err conftest.$ac_ext
-if test $ac_cv_header_stdc = yes; then
- # SunOS 4.x string.h does not declare mem*, contrary to ANSI.
- cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-#include <string.h>
+ if test "x${JPEG}" = "xyes"; then
+ FLAGS="${FLAGS} -DHAVE_LIBJPEG"
+ if test "x${JPEG_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} -ljpeg"
+ else
+ GMSH_LIBS="${GMSH_LIBS} -L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib -ljpeg"
+ FLAGS="${FLAGS} -I${JPEG_PREFIX} -I${JPEG_PREFIX}/include"
+ fi
+ fi
-_ACEOF
-if (eval "$ac_cpp conftest.$ac_ext") 2>&5 |
- egrep "memchr" >/dev/null 2>&1; then
- :
else
- ac_cv_header_stdc=no
-fi
-rm -f conftest*
-fi
+ GMSH_DIRS="${GMSH_DIRS} Box"
+ GMSH_LIBS="-Llib -lGmshBox -lGmshParser -lGmshMesh -lGmshGeo -lGmshNumeric -lGmshPlugin -lGmshCommon -lGmshDataStr -lGmshParallel"
-if test $ac_cv_header_stdc = yes; then
- # ISC 2.0.2 stdlib.h does not declare free, contrary to ANSI.
- cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-#include <stdlib.h>
+fi
-_ACEOF
-if (eval "$ac_cpp conftest.$ac_ext") 2>&5 |
- egrep "free" >/dev/null 2>&1; then
- :
+echo "$as_me:$LINENO: checking for ./Triangle/triangle.c" >&5
+echo $ECHO_N "checking for ./Triangle/triangle.c... $ECHO_C" >&6
+if test "${ac_cv_file___Triangle_triangle_c+set}" = set; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
else
- ac_cv_header_stdc=no
+ test "$cross_compiling" = yes &&
+ { { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
+echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
+ { (exit 1); exit 1; }; }
+if test -r "./Triangle/triangle.c"; then
+ ac_cv_file___Triangle_triangle_c=yes
+else
+ ac_cv_file___Triangle_triangle_c=no
+fi
+fi
+echo "$as_me:$LINENO: result: $ac_cv_file___Triangle_triangle_c" >&5
+echo "${ECHO_T}$ac_cv_file___Triangle_triangle_c" >&6
+if test $ac_cv_file___Triangle_triangle_c = yes; then
+ TRIANGLE="yes"
+else
+ TRIANGLE="no"
fi
-rm -f conftest*
+if test "x${TRIANGLE}" = "xyes" -a "x$enable_triangle" != "xno"; then
+ GMSH_DIRS="${GMSH_DIRS} Triangle"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshTriangle"
+ FLAGS="${FLAGS} -DHAVE_TRIANGLE"
+else
+ echo "*******************************************************************"
+ echo "If you want to use Jonathan Shewchuk's Triangle as an alternative"
+ echo "isotropic 2D mesh generator, please download Triangle from the"
+ echo "author's web site at http://www.cs.cmu.edu/~quake/triangle.html,"
+ echo "unpack the archive and copy the two files 'triangle.c' and"
+ echo "'triangle.h' in the ./Triangle subdirectory. Then run ./configure"
+ echo "again."
+ echo ""
+ echo "Please note that by doing so, you agree to Triangle's licencing"
+ echo "requirements stated in ./Triangle/README. (Most notably, you may"
+ echo "then only redistribute Gmsh if no compensation is received.)"
+ echo "*******************************************************************"
fi
-if test $ac_cv_header_stdc = yes; then
- # /bin/cc in Irix-4.0.5 gets non-ANSI ctype macros unless using -ansi.
- if test "$cross_compiling" = yes; then
- :
+if test "x$enable_gsl" = "xyes"; then
+ if test "x${GSL_PREFIX}" != "x"; then
+ LDFLAGS="-L${GSL_PREFIX}/lib ${LDFLAGS}"
+ fi
+
+echo "$as_me:$LINENO: checking for main in -lgslcblas" >&5
+echo $ECHO_N "checking for main in -lgslcblas... $ECHO_C" >&6
+if test "${ac_cv_lib_gslcblas_main+set}" = set; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
else
- cat >conftest.$ac_ext <<_ACEOF
+ ac_check_lib_save_LIBS=$LIBS
+LIBS="-lgslcblas $LIBS"
+cat >conftest.$ac_ext <<_ACEOF
#line $LINENO "configure"
#include "confdefs.h"
-#include <ctype.h>
-#if ((' ' & 0x0FF) == 0x020)
-# define ISLOWER(c) ('a' <= (c) && (c) <= 'z')
-# define TOUPPER(c) (ISLOWER(c) ? 'A' + ((c) - 'a') : (c))
-#else
-# define ISLOWER(c) (('a' <= (c) && (c) <= 'i') \
- || ('j' <= (c) && (c) <= 'r') \
- || ('s' <= (c) && (c) <= 'z'))
-# define TOUPPER(c) (ISLOWER(c) ? ((c) | 0x40) : (c))
-#endif
-#define XOR(e, f) (((e) && !(f)) || (!(e) && (f)))
+
+#ifdef F77_DUMMY_MAIN
+# ifdef __cplusplus
+ extern "C"
+# endif
+ int F77_DUMMY_MAIN() { return 1; }
+#endif
int
main ()
{
- int i;
- for (i = 0; i < 256; i++)
- if (XOR (islower (i), ISLOWER (i))
- || toupper (i) != TOUPPER (i))
- exit(2);
- exit (0);
+main ();
+ ;
+ return 0;
}
_ACEOF
-rm -f conftest$ac_exeext
+rm -f conftest.$ac_objext conftest$ac_exeext
if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
(eval $ac_link) 2>&5
ac_status=$?
echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } && { ac_try='./conftest$ac_exeext'
+ (exit $ac_status); } &&
+ { ac_try='test -s conftest$ac_exeext'
{ (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
(eval $ac_try) 2>&5
ac_status=$?
echo "$as_me:$LINENO: \$? = $ac_status" >&5
(exit $ac_status); }; }; then
- :
+ ac_cv_lib_gslcblas_main=yes
else
- echo "$as_me: program exited with status $ac_status" >&5
-echo "$as_me: failed program was:" >&5
+ echo "$as_me: failed program was:" >&5
cat conftest.$ac_ext >&5
-( exit $ac_status )
-ac_cv_header_stdc=no
-fi
-rm -f core core.* *.core conftest$ac_exeext conftest.$ac_objext conftest.$ac_ext
-fi
+ac_cv_lib_gslcblas_main=no
fi
+rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
+LIBS=$ac_check_lib_save_LIBS
fi
-echo "$as_me:$LINENO: result: $ac_cv_header_stdc" >&5
-echo "${ECHO_T}$ac_cv_header_stdc" >&6
-if test $ac_cv_header_stdc = yes; then
-
-cat >>confdefs.h <<\_ACEOF
-#define STDC_HEADERS 1
+echo "$as_me:$LINENO: result: $ac_cv_lib_gslcblas_main" >&5
+echo "${ECHO_T}$ac_cv_lib_gslcblas_main" >&6
+if test $ac_cv_lib_gslcblas_main = yes; then
+ cat >>confdefs.h <<_ACEOF
+#define HAVE_LIBGSLCBLAS 1
_ACEOF
-fi
-
-# On IRIX 5.3, sys/types and inttypes.h are conflicting.
-
-
-
-
-
-
-
+ LIBS="-lgslcblas $LIBS"
+fi
-for ac_header in sys/types.h sys/stat.h stdlib.h string.h memory.h strings.h \
- inttypes.h stdint.h unistd.h
-do
-as_ac_Header=`echo "ac_cv_header_$ac_header" | $as_tr_sh`
-echo "$as_me:$LINENO: checking for $ac_header" >&5
-echo $ECHO_N "checking for $ac_header... $ECHO_C" >&6
-if eval "test \"\${$as_ac_Header+set}\" = set"; then
- echo $ECHO_N "(cached) $ECHO_C" >&6
-else
- cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-$ac_includes_default
-
-#include <$ac_header>
-_ACEOF
-rm -f conftest.$ac_objext
-if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
- (eval $ac_compile) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } &&
- { ac_try='test -s conftest.$ac_objext'
- { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
- (eval $ac_try) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); }; }; then
- eval "$as_ac_Header=yes"
-else
- echo "$as_me: failed program was:" >&5
-cat conftest.$ac_ext >&5
-eval "$as_ac_Header=no"
-fi
-rm -f conftest.$ac_objext conftest.$ac_ext
-fi
-echo "$as_me:$LINENO: result: `eval echo '${'$as_ac_Header'}'`" >&5
-echo "${ECHO_T}`eval echo '${'$as_ac_Header'}'`" >&6
-if test `eval echo '${'$as_ac_Header'}'` = yes; then
- cat >>confdefs.h <<_ACEOF
-#define `echo "HAVE_$ac_header" | $as_tr_cpp` 1
-_ACEOF
-
-fi
-
-done
-
-
-if test "${ac_cv_header_jpeglib_h+set}" = set; then
- echo "$as_me:$LINENO: checking for jpeglib.h" >&5
-echo $ECHO_N "checking for jpeglib.h... $ECHO_C" >&6
-if test "${ac_cv_header_jpeglib_h+set}" = set; then
- echo $ECHO_N "(cached) $ECHO_C" >&6
-fi
-echo "$as_me:$LINENO: result: $ac_cv_header_jpeglib_h" >&5
-echo "${ECHO_T}$ac_cv_header_jpeglib_h" >&6
-else
- # Is the header compilable?
-echo "$as_me:$LINENO: checking jpeglib.h usability" >&5
-echo $ECHO_N "checking jpeglib.h usability... $ECHO_C" >&6
-cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-$ac_includes_default
-#include <jpeglib.h>
-_ACEOF
-rm -f conftest.$ac_objext
-if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
- (eval $ac_compile) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } &&
- { ac_try='test -s conftest.$ac_objext'
- { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
- (eval $ac_try) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); }; }; then
- ac_header_compiler=yes
-else
- echo "$as_me: failed program was:" >&5
-cat conftest.$ac_ext >&5
-ac_header_compiler=no
-fi
-rm -f conftest.$ac_objext conftest.$ac_ext
-echo "$as_me:$LINENO: result: $ac_header_compiler" >&5
-echo "${ECHO_T}$ac_header_compiler" >&6
-
-# Is the header present?
-echo "$as_me:$LINENO: checking jpeglib.h presence" >&5
-echo $ECHO_N "checking jpeglib.h presence... $ECHO_C" >&6
-cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-#include <jpeglib.h>
-_ACEOF
-if { (eval echo "$as_me:$LINENO: \"$ac_cpp conftest.$ac_ext\"") >&5
- (eval $ac_cpp conftest.$ac_ext) 2>conftest.er1
- ac_status=$?
- egrep -v '^ *\+' conftest.er1 >conftest.err
- rm -f conftest.er1
- cat conftest.err >&5
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } >/dev/null; then
- if test -s conftest.err; then
- ac_cpp_err=$ac_c_preproc_warn_flag
- else
- ac_cpp_err=
- fi
-else
- ac_cpp_err=yes
-fi
-if test -z "$ac_cpp_err"; then
- ac_header_preproc=yes
-else
- echo "$as_me: failed program was:" >&5
- cat conftest.$ac_ext >&5
- ac_header_preproc=no
-fi
-rm -f conftest.err conftest.$ac_ext
-echo "$as_me:$LINENO: result: $ac_header_preproc" >&5
-echo "${ECHO_T}$ac_header_preproc" >&6
-
-# So? What about this header?
-case $ac_header_compiler:$ac_header_preproc in
- yes:no )
- { echo "$as_me:$LINENO: WARNING: jpeglib.h: accepted by the compiler, rejected by the preprocessor!" >&5
-echo "$as_me: WARNING: jpeglib.h: accepted by the compiler, rejected by the preprocessor!" >&2;}
- { echo "$as_me:$LINENO: WARNING: jpeglib.h: proceeding with the preprocessor's result" >&5
-echo "$as_me: WARNING: jpeglib.h: proceeding with the preprocessor's result" >&2;};;
- no:yes )
- { echo "$as_me:$LINENO: WARNING: jpeglib.h: present but cannot be compiled" >&5
-echo "$as_me: WARNING: jpeglib.h: present but cannot be compiled" >&2;}
- { echo "$as_me:$LINENO: WARNING: jpeglib.h: check for missing prerequisite headers?" >&5
-echo "$as_me: WARNING: jpeglib.h: check for missing prerequisite headers?" >&2;}
- { echo "$as_me:$LINENO: WARNING: jpeglib.h: proceeding with the preprocessor's result" >&5
-echo "$as_me: WARNING: jpeglib.h: proceeding with the preprocessor's result" >&2;};;
-esac
-echo "$as_me:$LINENO: checking for jpeglib.h" >&5
-echo $ECHO_N "checking for jpeglib.h... $ECHO_C" >&6
-if test "${ac_cv_header_jpeglib_h+set}" = set; then
- echo $ECHO_N "(cached) $ECHO_C" >&6
-else
- ac_cv_header_jpeglib_h=$ac_header_preproc
-fi
-echo "$as_me:$LINENO: result: $ac_cv_header_jpeglib_h" >&5
-echo "${ECHO_T}$ac_cv_header_jpeglib_h" >&6
-
-fi
-if test $ac_cv_header_jpeglib_h = yes; then
- echo "$as_me:$LINENO: checking for jpeg_destroy_decompress in -ljpeg" >&5
-echo $ECHO_N "checking for jpeg_destroy_decompress in -ljpeg... $ECHO_C" >&6
-if test "${ac_cv_lib_jpeg_jpeg_destroy_decompress+set}" = set; then
+ echo "$as_me:$LINENO: checking for main in -lgsl" >&5
+echo $ECHO_N "checking for main in -lgsl... $ECHO_C" >&6
+if test "${ac_cv_lib_gsl_main+set}" = set; then
echo $ECHO_N "(cached) $ECHO_C" >&6
else
ac_check_lib_save_LIBS=$LIBS
-LIBS="-ljpeg $LIBS"
-
+LIBS="-lgsl $LIBS"
cat >conftest.$ac_ext <<_ACEOF
#line $LINENO "configure"
#include "confdefs.h"
-/* Override any gcc2 internal prototype to avoid an error. */
-#ifdef __cplusplus
-extern "C"
-#endif
-/* We use char because int might match the return type of a gcc2
- builtin and then its argument prototype would still apply. */
-char jpeg_destroy_decompress ();
+
#ifdef F77_DUMMY_MAIN
# ifdef __cplusplus
extern "C"
@@ -3090,7 +3031,7 @@ char jpeg_destroy_decompress ();
int
main ()
{
-jpeg_destroy_decompress ();
+main ();
;
return 0;
}
@@ -3107,88 +3048,92 @@ if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
ac_status=$?
echo "$as_me:$LINENO: \$? = $ac_status" >&5
(exit $ac_status); }; }; then
- ac_cv_lib_jpeg_jpeg_destroy_decompress=yes
+ ac_cv_lib_gsl_main=yes
else
echo "$as_me: failed program was:" >&5
cat conftest.$ac_ext >&5
-ac_cv_lib_jpeg_jpeg_destroy_decompress=no
+ac_cv_lib_gsl_main=no
fi
rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
LIBS=$ac_check_lib_save_LIBS
fi
-echo "$as_me:$LINENO: result: $ac_cv_lib_jpeg_jpeg_destroy_decompress" >&5
-echo "${ECHO_T}$ac_cv_lib_jpeg_jpeg_destroy_decompress" >&6
-if test $ac_cv_lib_jpeg_jpeg_destroy_decompress = yes; then
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_LIBJPEG"
- GUI_LIBS="${GUI_LIBS} -ljpeg"
-fi
-
+echo "$as_me:$LINENO: result: $ac_cv_lib_gsl_main" >&5
+echo "${ECHO_T}$ac_cv_lib_gsl_main" >&6
+if test $ac_cv_lib_gsl_main = yes; then
+ GSL="yes"
+else
+ GSL="no"
fi
-
- else
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_LIBJPEG"
- GUI_INCLUDE="${GUI_INCLUDE} -I${JPEG_PREFIX} -I${JPEG_PREFIX}/include"
- GUI_LIBS="${GUI_LIBS} -L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib -ljpeg"
+ if test "x${GSL}" = "xyes"; then
+ FLAGS="${FLAGS} -DHAVE_GSL"
+ if test "x${GSL_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} -lgsl -lgslcblas"
+ else
+ GMSH_LIBS="${GMSH_LIBS} -L${GSL_PREFIX}/lib -lgsl -lgslcblas"
+ FLAGS="${FLAGS} -I${GSL_PREFIX}/include"
+ fi
fi
-
-else
-
- GMSH_DIRS="${GMSH_DIRS} Box"
- GMSH_LIBS="-Llib -lGmshBox -lGmshParser -lGmshMesh -lGmshGeo -lGmshNumeric -lGmshPlugin -lGmshCommon -lGmshDataStr -lGmshParallel"
- VERSION_FLAGS=""
- GUI_LIBS=""
- GUI_INCLUDE=""
-
fi
-
-echo "$as_me:$LINENO: checking for ./Triangle/triangle.c" >&5
-echo $ECHO_N "checking for ./Triangle/triangle.c... $ECHO_C" >&6
-if test "${ac_cv_file___Triangle_triangle_c+set}" = set; then
+if test "x${GSL}" != "xyes"; then
+ echo "$as_me:$LINENO: checking for ./NR/dsvdcmp.cpp" >&5
+echo $ECHO_N "checking for ./NR/dsvdcmp.cpp... $ECHO_C" >&6
+if test "${ac_cv_file___NR_dsvdcmp_cpp+set}" = set; then
echo $ECHO_N "(cached) $ECHO_C" >&6
else
test "$cross_compiling" = yes &&
{ { echo "$as_me:$LINENO: error: cannot check for file existence when cross compiling" >&5
echo "$as_me: error: cannot check for file existence when cross compiling" >&2;}
{ (exit 1); exit 1; }; }
-if test -r "./Triangle/triangle.c"; then
- ac_cv_file___Triangle_triangle_c=yes
+if test -r "./NR/dsvdcmp.cpp"; then
+ ac_cv_file___NR_dsvdcmp_cpp=yes
else
- ac_cv_file___Triangle_triangle_c=no
+ ac_cv_file___NR_dsvdcmp_cpp=no
fi
fi
-echo "$as_me:$LINENO: result: $ac_cv_file___Triangle_triangle_c" >&5
-echo "${ECHO_T}$ac_cv_file___Triangle_triangle_c" >&6
-if test $ac_cv_file___Triangle_triangle_c = yes; then
- HAVE_TRIANGLE="yes"
+echo "$as_me:$LINENO: result: $ac_cv_file___NR_dsvdcmp_cpp" >&5
+echo "${ECHO_T}$ac_cv_file___NR_dsvdcmp_cpp" >&6
+if test $ac_cv_file___NR_dsvdcmp_cpp = yes; then
+ NR="yes"
else
- HAVE_TRIANGLE="no"
+ NR="no"
fi
-if test "x${HAVE_TRIANGLE}" = "xyes" -a "x$enable_triangle" != "xno"; then
- GMSH_DIRS="${GMSH_DIRS} Triangle"
- GMSH_LIBS="${GMSH_LIBS} -lGmshTriangle"
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_TRIANGLE"
+ if test "x${NR}" = "xyes"; then
+ echo "*******************************************************************"
+ echo "You are building a non-free version of Gmsh, using some code"
+ echo "copyrighted by Numerical Recipes."
+ echo "*******************************************************************"
+ GMSH_DIRS="${GMSH_DIRS} NR"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshNR"
+ else
+ echo "*******************************************************************"
+ echo "This is the free version of Gmsh and configure could not find"
+ echo "the GNU Scientific Library (GSL) on your system:"
+ echo "- if it is installed in a non-standard location, please run"
+ echo " ./configure again with the --with-gsl-prefix option"
+ echo "- if it is not installed on your system, you can download it from"
+ echo " http://sources.redhat.com/gsl/"
+ echo ""
+ echo "IMPORTANT NOTE: You need to install GSL 1.2 or above. All versions"
+ echo "<= 1.1.1 have a bug in the singular value decomposition algorithm"
+ echo "that will cause Gmsh to hang during mesh generation."
+ echo "*******************************************************************"
+ { { echo "$as_me:$LINENO: error: " >&5
+echo "$as_me: error: " >&2;}
+ { (exit 1); exit 1; }; }
+ fi
fi
+GMSH_LIBS="${GMSH_LIBS} -lm"
+
if test "x$enable_parallel" = "xyes"; then
- VERSION_FLAGS="${VERSION_FLAGS} -DPARALLEL"
+ FLAGS="${FLAGS} -DPARALLEL"
fi
-
-
-
-OPT_FLAGS="${CXXFLAGS}"
-OS_FLAGS=""
LINKER="${CXX}"
POSTBUILD=""
-if test "x$enable_gui" != "xno"; then
- GMSH_LIBS="${GMSH_LIBS} ${GUI_LIBS}"
-else
- GMSH_LIBS="${GMSH_LIBS} -lm";
-fi
-
case "$UNAME" in
CYGWIN* | MINGW*)
@@ -3199,30 +3144,30 @@ case "$UNAME" in
;;
Darwin*)
- OS_FLAGS="${OS_FLAGS} -D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
POSTBUILD="/Developer/Tools/Rez -t APPL -o bin/gmsh Fltk/MacRes.r"
;;
AIX*)
- OS_FLAGS="${OS_FLAGS} -D_BSD -D_NO_DLL"
+ FLAGS="${FLAGS} -D_BSD -D_NO_DLL"
;;
IRIX*)
CXX="CC"
CC="cc"
- OS_FLAGS="${OS_FLAGS} -mips3 -n32"
- OPT_FLAGS="-O2 -OPT:Olimit=0 -LANG:std"
+ OPTIM="-O2"
+ FLAGS="${FLAGS} -mips3 -n32 -LANG:std -OPT:Olimit=0"
AR="CC -mips3 -n32 -ar -o"
LINKER="CC -O2 -mips3 -n32"
;;
SunOS*)
- OS_FLAGS="${OS_FLAGS} -D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
GMSH_LIBS="${GMSH_LIBS} -lsocket -lnsl -ldl"
;;
HP-UX*)
- OS_FLAGS="-D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
;;
esac
@@ -3414,7 +3359,7 @@ echo "$as_me:$LINENO: result: $ac_cv_c_bigendian" >&5
echo "${ECHO_T}$ac_cv_c_bigendian" >&6
case $ac_cv_c_bigendian in
yes)
- OS_FLAGS="${OS_FLAGS} -D_BIG_ENDIAN" ;;
+ FLAGS="${FLAGS} -D_BIG_ENDIAN" ;;
no)
;;
*)
@@ -3427,88 +3372,6 @@ esac
-
-
-
-
-
-
-
-if test "x${GSL_PREFIX}" = "x"; then
- LDFLAGS="-lgslcblas"
-else
- LDFLAGS="-L${GSL_PREFIX}/lib -lgslcblas"
-fi
-echo "$as_me:$LINENO: checking for gsl_vector_alloc in -lgsl" >&5
-echo $ECHO_N "checking for gsl_vector_alloc in -lgsl... $ECHO_C" >&6
-if test "${ac_cv_lib_gsl_gsl_vector_alloc+set}" = set; then
- echo $ECHO_N "(cached) $ECHO_C" >&6
-else
- ac_check_lib_save_LIBS=$LIBS
-LIBS="-lgsl $LIBS"
-cat >conftest.$ac_ext <<_ACEOF
-#line $LINENO "configure"
-#include "confdefs.h"
-
-/* Override any gcc2 internal prototype to avoid an error. */
-#ifdef __cplusplus
-extern "C"
-#endif
-/* We use char because int might match the return type of a gcc2
- builtin and then its argument prototype would still apply. */
-char gsl_vector_alloc ();
-#ifdef F77_DUMMY_MAIN
-# ifdef __cplusplus
- extern "C"
-# endif
- int F77_DUMMY_MAIN() { return 1; }
-#endif
-int
-main ()
-{
-gsl_vector_alloc ();
- ;
- return 0;
-}
-_ACEOF
-rm -f conftest.$ac_objext conftest$ac_exeext
-if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
- (eval $ac_link) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); } &&
- { ac_try='test -s conftest$ac_exeext'
- { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
- (eval $ac_try) 2>&5
- ac_status=$?
- echo "$as_me:$LINENO: \$? = $ac_status" >&5
- (exit $ac_status); }; }; then
- ac_cv_lib_gsl_gsl_vector_alloc=yes
-else
- echo "$as_me: failed program was:" >&5
-cat conftest.$ac_ext >&5
-ac_cv_lib_gsl_gsl_vector_alloc=no
-fi
-rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
-LIBS=$ac_check_lib_save_LIBS
-fi
-echo "$as_me:$LINENO: result: $ac_cv_lib_gsl_gsl_vector_alloc" >&5
-echo "${ECHO_T}$ac_cv_lib_gsl_gsl_vector_alloc" >&6
-if test $ac_cv_lib_gsl_gsl_vector_alloc = yes; then
- GSL="yes"
-else
- GSL="no"
-fi
-
-if test "x${GSL_PREFIX}" = "x"; then
- GSL_LIBS="-lgsl -lgslcblas"
-else
- GSL_LIBS="-L${GSL_PREFIX}/lib -lgsl -lgslcblas"
- GSL_INCLUDE="-I${GSL_PREFIX}/include"
-fi
-
-
-
echo "$as_me:$LINENO: checking for ANSI C header files" >&5
echo $ECHO_N "checking for ANSI C header files... $ECHO_C" >&6
if test "${ac_cv_header_stdc+set}" = set; then
@@ -3648,6 +3511,64 @@ _ACEOF
fi
+# On IRIX 5.3, sys/types and inttypes.h are conflicting.
+
+
+
+
+
+
+
+
+
+for ac_header in sys/types.h sys/stat.h stdlib.h string.h memory.h strings.h \
+ inttypes.h stdint.h unistd.h
+do
+as_ac_Header=`echo "ac_cv_header_$ac_header" | $as_tr_sh`
+echo "$as_me:$LINENO: checking for $ac_header" >&5
+echo $ECHO_N "checking for $ac_header... $ECHO_C" >&6
+if eval "test \"\${$as_ac_Header+set}\" = set"; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+ cat >conftest.$ac_ext <<_ACEOF
+#line $LINENO "configure"
+#include "confdefs.h"
+$ac_includes_default
+
+#include <$ac_header>
+_ACEOF
+rm -f conftest.$ac_objext
+if { (eval echo "$as_me:$LINENO: \"$ac_compile\"") >&5
+ (eval $ac_compile) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); } &&
+ { ac_try='test -s conftest.$ac_objext'
+ { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+ (eval $ac_try) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); }; }; then
+ eval "$as_ac_Header=yes"
+else
+ echo "$as_me: failed program was:" >&5
+cat conftest.$ac_ext >&5
+eval "$as_ac_Header=no"
+fi
+rm -f conftest.$ac_objext conftest.$ac_ext
+fi
+echo "$as_me:$LINENO: result: `eval echo '${'$as_ac_Header'}'`" >&5
+echo "${ECHO_T}`eval echo '${'$as_ac_Header'}'`" >&6
+if test `eval echo '${'$as_ac_Header'}'` = yes; then
+ cat >>confdefs.h <<_ACEOF
+#define `echo "HAVE_$ac_header" | $as_tr_cpp` 1
+_ACEOF
+
+fi
+
+done
+
+
for ac_header in sys/time.h sys/resource.h
@@ -3764,6 +3685,14 @@ fi
done
+
+
+
+
+
+
+
+
ac_config_files="$ac_config_files variables"
cat >confcache <<\_ACEOF
# This file is a shell script that caches the results of configure
@@ -4342,7 +4271,6 @@ s,@ECHO_C@,$ECHO_C,;t t
s,@ECHO_N@,$ECHO_N,;t t
s,@ECHO_T@,$ECHO_T,;t t
s,@LIBS@,$LIBS,;t t
-s,@UNAME@,$UNAME,;t t
s,@CC@,$CC,;t t
s,@CFLAGS@,$CFLAGS,;t t
s,@LDFLAGS@,$LDFLAGS,;t t
@@ -4358,16 +4286,13 @@ s,@RANLIB@,$RANLIB,;t t
s,@ac_ct_RANLIB@,$ac_ct_RANLIB,;t t
s,@AR@,$AR,;t t
s,@FLTKCONFIG@,$FLTKCONFIG,;t t
-s,@VERSION_FLAGS@,$VERSION_FLAGS,;t t
-s,@GUI_INCLUDE@,$GUI_INCLUDE,;t t
-s,@OPT_FLAGS@,$OPT_FLAGS,;t t
-s,@OS_FLAGS@,$OS_FLAGS,;t t
+s,@UNAME@,$UNAME,;t t
+s,@FLAGS@,$FLAGS,;t t
+s,@OPTIM@,$OPTIM,;t t
s,@LINKER@,$LINKER,;t t
s,@GMSH_DIRS@,$GMSH_DIRS,;t t
s,@GMSH_LIBS@,$GMSH_LIBS,;t t
s,@POSTBUILD@,$POSTBUILD,;t t
-s,@GSL_LIBS@,$GSL_LIBS,;t t
-s,@GSL_INCLUDE@,$GSL_INCLUDE,;t t
CEOF
_ACEOF
@@ -4594,22 +4519,6 @@ fi
echo "*******************************************************************"
echo "Gmsh is configured for"
echo "- Operating system: $UNAME"
-echo "- Version flags: $VERSION_FLAGS"
-echo "- OS flags: $OS_FLAGS"
+echo "- Compiler flags: $FLAGS"
echo "*******************************************************************"
-if test "x${HAVE_TRIANGLE}" = "xno"; then
- echo "*******************************************************************"
- echo "If you want to use Jonathan Shewchuk's Triangle as an alternative"
- echo "isotropic 2D mesh generator, please download Triangle from the"
- echo "author's web site at http://www.cs.cmu.edu/~quake/triangle.html,"
- echo "unpack the archive and copy the two files 'triangle.c' and"
- echo "'triangle.h' in the ./Triangle subdirectory. Then run ./configure"
- echo "again."
- echo ""
- echo "Please note that by doing so, you agree to Triangle's licencing"
- echo "requirements stated in ./Triangle/README. (Most notably, you may"
- echo "then only redistribute Gmsh if no compensation is received.)"
- echo "*******************************************************************"
-fi
-
diff --git a/configure.in b/configure.in
index a3bcfc73ccfe735825b39a006bc7df125f3a1cc0..695179515c9b526bf1f8fe3ca3822ad189f4e34e 100644
--- a/configure.in
+++ b/configure.in
@@ -1,4 +1,4 @@
-dnl "$Id: configure.in,v 1.23 2003-02-17 16:32:37 geuzaine Exp $"
+dnl "$Id: configure.in,v 1.24 2003-02-18 05:50:04 geuzaine Exp $"
dnl
dnl Copyright (C) 1997 - 2003 C. Geuzaine, J.-F. Remacle
dnl
@@ -21,6 +21,10 @@ dnl Please report all bugs and problems to "gmsh@geuz.org".
dnl Process this file with autoconf to produce the configure script.
+dnl Note: each CHECK_LIB adds the found library to LIBS, which serves
+dnl implicitely when checking for the next one. So one should check
+dnl in reverse order! (start with -lm, ...)
+
dnl Check that this is the gmsh source tree
AC_INIT(Parser/Gmsh.y)
@@ -48,6 +52,10 @@ AC_ARG_ENABLE(parallel,
AC_ARG_ENABLE(triangle,
AC_HELP_STRING([--enable-triangle],
[compile Triangle if available (default=yes)]))
+AC_ARG_ENABLE(gsl,
+ AC_HELP_STRING([--enable-gsl],
+ [use GSL as numerical toolkit (default=yes)]))
+dnl FIXME [use GSL as numerical toolkit (default=yes)]))
dnl Get the operating system and version number
UNAME=`uname`
@@ -55,17 +63,23 @@ UVERSION=`uname -r | sed -e '1,$s/[[^0-9]]//g'`
if test "x${UNAME}" = "xIRIX64"; then
UNAME="IRIX"
fi
-AC_SUBST(UNAME)
dnl Check for default compilers
AC_PROG_CC
AC_PROG_CXX
+dnl Set default flags
+FLAGS=""
+OPTIM="${CXXFLAGS}"
+
dnl Check for various programs
AC_PROG_CPP
AC_PROG_RANLIB
AC_PATH_PROG(AR, ar)
+dnl Check for standard math library
+AC_CHECK_LIB(m,main)
+
dnl See if we need a .exe extension on executables
AC_EXEEXT
@@ -83,74 +97,123 @@ if test "x$enable_gui" != "xno"; then
GMSH_DIRS="${GMSH_DIRS} Graphics Fltk"
GMSH_LIBS="-Llib -lGmshFltk -lGmshParser -lGmshGraphics -lGmshMesh -lGmshGeo -lGmshNumeric -lGmshCommon -lGmshDataStr -lGmshPlugin -lGmshParallel"
- VERSION_FLAGS="-DHAVE_FLTK"
+ FLAGS="${FLAGS} -DHAVE_FLTK"
if test "x${FLTK_PREFIX}" != "x" ; then
AC_PATH_PROG(FLTKCONFIG,fltk-config,"",${FLTK_PREFIX})
dnl Find the libs/includes even if fltk is _not_ properly installed (ugly hack!)
- GUI_LIBS="-L${FLTK_PREFIX}/lib"
- GUI_INCLUDE="-I${FLTK_PREFIX}"
+ GMSH_LIBS="${GMSH_LIBS} -L${FLTK_PREFIX}/lib"
+ FLAGS="${FLAGS} -I${FLTK_PREFIX}"
else
AC_PATH_PROG(FLTKCONFIG,fltk-config)
fi
if test "x$FLTKCONFIG" = "x"; then
AC_MSG_ERROR(Could not find fltk-config. Try --with-fltk-prefix?)
fi
- GUI_LIBS="${GUI_LIBS} `$FLTKCONFIG --use-gl --use-images --ldflags`"
- GUI_INCLUDE="${GUI_INCLUDE} `$FLTKCONFIG --use-gl --use-images --cxxflags`"
-
- dnl Check if libjpeg is available (in case fltk didn't find it) to enable/disable gl2jpg
- if test "x${JPEG_PREFIX}" = "x"; then
- AC_CHECK_HEADER(jpeglib.h,
- AC_CHECK_LIB(jpeg, jpeg_destroy_decompress,
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_LIBJPEG"
- GUI_LIBS="${GUI_LIBS} -ljpeg"))
- else
- dnl Should work even if the jpeg libs/includes are _not_ properly installed (ugly hack!)
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_LIBJPEG"
- GUI_INCLUDE="${GUI_INCLUDE} -I${JPEG_PREFIX} -I${JPEG_PREFIX}/include"
- GUI_LIBS="${GUI_LIBS} -L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib -ljpeg"
+ GMSH_LIBS="${GMSH_LIBS} `$FLTKCONFIG --use-gl --use-images --ldflags`"
+ FLAGS="${FLAGS} `$FLTKCONFIG --use-gl --use-images --cxxflags`"
+
+ dnl Check if libjpeg is available to enable/disable gl2jpg
+ if test "x${JPEG_PREFIX}" != "x"; then
+ LDFLAGS="-L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib ${LDFLAGS}"
+ fi
+ AC_CHECK_LIB(jpeg,main,JPEG="yes",JPEG="no")
+ if test "x${JPEG}" = "xyes"; then
+ FLAGS="${FLAGS} -DHAVE_LIBJPEG"
+ if test "x${JPEG_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} -ljpeg"
+ else
+ dnl Find the libs/includes even if libjpeg is _not_ properly installed (ugly hack!)
+ GMSH_LIBS="${GMSH_LIBS} -L${JPEG_PREFIX} -L${JPEG_PREFIX}/lib -ljpeg"
+ FLAGS="${FLAGS} -I${JPEG_PREFIX} -I${JPEG_PREFIX}/include"
+ fi
fi
else
GMSH_DIRS="${GMSH_DIRS} Box"
GMSH_LIBS="-Llib -lGmshBox -lGmshParser -lGmshMesh -lGmshGeo -lGmshNumeric -lGmshPlugin -lGmshCommon -lGmshDataStr -lGmshParallel"
- VERSION_FLAGS=""
- GUI_LIBS=""
- GUI_INCLUDE=""
fi
dnl Check if Triangle is installed
-AC_CHECK_FILE(./Triangle/triangle.c, HAVE_TRIANGLE="yes", HAVE_TRIANGLE="no")
-if test "x${HAVE_TRIANGLE}" = "xyes" -a "x$enable_triangle" != "xno"; then
+AC_CHECK_FILE(./Triangle/triangle.c, TRIANGLE="yes", TRIANGLE="no")
+if test "x${TRIANGLE}" = "xyes" -a "x$enable_triangle" != "xno"; then
GMSH_DIRS="${GMSH_DIRS} Triangle"
GMSH_LIBS="${GMSH_LIBS} -lGmshTriangle"
- VERSION_FLAGS="${VERSION_FLAGS} -DHAVE_TRIANGLE"
+ FLAGS="${FLAGS} -DHAVE_TRIANGLE"
+else
+ echo "*******************************************************************"
+ echo "If you want to use Jonathan Shewchuk's Triangle as an alternative"
+ echo "isotropic 2D mesh generator, please download Triangle from the"
+ echo "author's web site at http://www.cs.cmu.edu/~quake/triangle.html,"
+ echo "unpack the archive and copy the two files 'triangle.c' and"
+ echo "'triangle.h' in the ./Triangle subdirectory. Then run ./configure"
+ echo "again."
+ echo ""
+ echo "Please note that by doing so, you agree to Triangle's licencing"
+ echo "requirements stated in ./Triangle/README. (Most notably, you may"
+ echo "then only redistribute Gmsh if no compensation is received.)"
+ echo "*******************************************************************"
fi
+dnl Check for GSL
+dnl FIXME: if test "x$enable_gsl" != "xno"; then
+if test "x$enable_gsl" = "xyes"; then
+ if test "x${GSL_PREFIX}" != "x"; then
+ LDFLAGS="-L${GSL_PREFIX}/lib ${LDFLAGS}"
+ fi
+ AC_CHECK_LIB(gslcblas,main)
+ AC_CHECK_LIB(gsl,main,GSL="yes",GSL="no")
+ if test "x${GSL}" = "xyes"; then
+ FLAGS="${FLAGS} -DHAVE_GSL"
+ if test "x${GSL_PREFIX}" = "x"; then
+ GMSH_LIBS="${GMSH_LIBS} -lgsl -lgslcblas"
+ else
+ GMSH_LIBS="${GMSH_LIBS} -L${GSL_PREFIX}/lib -lgsl -lgslcblas"
+ FLAGS="${FLAGS} -I${GSL_PREFIX}/include"
+ fi
+ fi
+fi
+if test "x${GSL}" != "xyes"; then
+ dnl Check if non-free numerical recipes routines are in the tree
+ AC_CHECK_FILE(./NR/dsvdcmp.cpp,NR="yes",NR="no")
+ if test "x${NR}" = "xyes"; then
+ echo "*******************************************************************"
+ echo "You are building a non-free version of Gmsh, using some code"
+ echo "copyrighted by Numerical Recipes."
+ echo "*******************************************************************"
+ GMSH_DIRS="${GMSH_DIRS} NR"
+ GMSH_LIBS="${GMSH_LIBS} -lGmshNR"
+ else
+ echo "*******************************************************************"
+ echo "This is the free version of Gmsh and configure could not find"
+ echo "the GNU Scientific Library (GSL) on your system:"
+ echo "- if it is installed in a non-standard location, please run"
+ echo " ./configure again with the --with-gsl-prefix option"
+ echo "- if it is not installed on your system, you can download it from"
+ echo " http://sources.redhat.com/gsl/"
+ echo ""
+ echo "IMPORTANT NOTE: You need to install GSL 1.2 or above. All versions"
+ echo "<= 1.1.1 have a bug in the singular value decomposition algorithm"
+ echo "that will cause Gmsh to hang during mesh generation."
+ echo "*******************************************************************"
+ AC_MSG_ERROR()
+ fi
+fi
+
+dnl Finish link line
+GMSH_LIBS="${GMSH_LIBS} -lm"
+
dnl Check if we should build the parallel version
if test "x$enable_parallel" = "xyes"; then
- VERSION_FLAGS="${VERSION_FLAGS} -DPARALLEL"
+ FLAGS="${FLAGS} -DPARALLEL"
fi
-AC_SUBST(VERSION_FLAGS)
-AC_SUBST(GUI_INCLUDE)
-
-dnl Set default flags, linker and post build action
-OPT_FLAGS="${CXXFLAGS}"
-OS_FLAGS=""
+dnl Set default linker and post build action
LINKER="${CXX}"
POSTBUILD=""
-dnl Set libs to link against
-if test "x$enable_gui" != "xno"; then
- GMSH_LIBS="${GMSH_LIBS} ${GUI_LIBS}"
-else
- GMSH_LIBS="${GMSH_LIBS} -lm";
-fi
-
dnl Modify defaults according to OS
case "$UNAME" in
@@ -162,98 +225,56 @@ case "$UNAME" in
;;
Darwin*)
- OS_FLAGS="${OS_FLAGS} -D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
POSTBUILD="/Developer/Tools/Rez -t APPL -o bin/gmsh Fltk/MacRes.r"
;;
AIX*)
- OS_FLAGS="${OS_FLAGS} -D_BSD -D_NO_DLL"
+ FLAGS="${FLAGS} -D_BSD -D_NO_DLL"
;;
IRIX*)
CXX="CC"
CC="cc"
- OS_FLAGS="${OS_FLAGS} -mips3 -n32"
- OPT_FLAGS="-O2 -OPT:Olimit=0 -LANG:std"
+ OPTIM="-O2"
+ FLAGS="${FLAGS} -mips3 -n32 -LANG:std -OPT:Olimit=0"
AR="CC -mips3 -n32 -ar -o"
LINKER="CC -O2 -mips3 -n32"
;;
SunOS*)
- OS_FLAGS="${OS_FLAGS} -D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
GMSH_LIBS="${GMSH_LIBS} -lsocket -lnsl -ldl"
;;
HP-UX*)
- OS_FLAGS="-D_NO_DLL"
+ FLAGS="${FLAGS} -D_NO_DLL"
;;
esac
dnl Is the machine big or littke endian?
-AC_C_BIGENDIAN(OS_FLAGS="${OS_FLAGS} -D_BIG_ENDIAN")
-
-AC_SUBST(OPT_FLAGS)
-AC_SUBST(OS_FLAGS)
-AC_SUBST(LINKER)
-AC_SUBST(GMSH_DIRS)
-AC_SUBST(GMSH_LIBS)
-AC_SUBST(POSTBUILD)
-AC_SUBST(AR)
-
-dnl Check for GSL
-if test "x${GSL_PREFIX}" = "x"; then
- LDFLAGS="-lgslcblas"
-else
- LDFLAGS="-L${GSL_PREFIX}/lib -lgslcblas"
-fi
-AC_CHECK_LIB(gsl,gsl_vector_alloc,GSL="yes",GSL="no")
-if test "x${GSL_PREFIX}" = "x"; then
- GSL_LIBS="-lgsl -lgslcblas"
-else
- GSL_LIBS="-L${GSL_PREFIX}/lib -lgsl -lgslcblas"
- GSL_INCLUDE="-I${GSL_PREFIX}/include"
-fi
-AC_SUBST(GSL_LIBS)
-AC_SUBST(GSL_INCLUDE)
+AC_C_BIGENDIAN(FLAGS="${FLAGS} -D_BIG_ENDIAN")
dnl Check for header files
AC_HEADER_STDC
AC_CHECK_HEADERS(sys/time.h sys/resource.h)
dnl Write output
+AC_SUBST(UNAME)
+AC_SUBST(FLAGS)
+AC_SUBST(OPTIM)
+AC_SUBST(LINKER)
+AC_SUBST(GMSH_DIRS)
+AC_SUBST(GMSH_LIBS)
+AC_SUBST(POSTBUILD)
+AC_SUBST(AR)
AC_OUTPUT(variables)
dnl Print some information
echo "*******************************************************************"
echo "Gmsh is configured for"
echo "- Operating system: $UNAME"
-echo "- Version flags: $VERSION_FLAGS"
-echo "- OS flags: $OS_FLAGS"
+echo "- Compiler flags: $FLAGS"
echo "*******************************************************************"
-if test "x${HAVE_TRIANGLE}" = "xno"; then
- echo "*******************************************************************"
- echo "If you want to use Jonathan Shewchuk's Triangle as an alternative"
- echo "isotropic 2D mesh generator, please download Triangle from the"
- echo "author's web site at http://www.cs.cmu.edu/~quake/triangle.html,"
- echo "unpack the archive and copy the two files 'triangle.c' and"
- echo "'triangle.h' in the ./Triangle subdirectory. Then run ./configure"
- echo "again."
- echo ""
- echo "Please note that by doing so, you agree to Triangle's licencing"
- echo "requirements stated in ./Triangle/README. (Most notably, you may"
- echo "then only redistribute Gmsh if no compensation is received.)"
- echo "*******************************************************************"
-fi
-
-dnl Don't print the warning at the moment, since we don't use the GSL yet
-dnl if test "x${GSL}" = "xno"; then
-dnl echo "*******************************************************************"
-dnl echo "Configure could not find the GNU Scientific Library (GSL):"
-dnl echo "- if it is installed in a non-standard location, please run"
-dnl echo " ./configure again with the --with-gsl-prefix option"
-dnl echo "- if it is not installed on your system, you can download it from"
-dnl echo " http://sources.redhat.com/gsl/"
-dnl echo "*******************************************************************"
-dnl fi
diff --git a/utils/Makefile b/utils/Makefile
index dad9461d568a7032c2f23de889949ab271d979c4..ae6c1e1c6c6ec11824e61aaf9ded90a08e7c7371 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -1,13 +1,13 @@
-# $Id: Makefile,v 1.11 2003-02-11 22:50:41 geuzaine Exp $
+# $Id: Makefile,v 1.12 2003-02-18 05:50:10 geuzaine Exp $
include ../variables
dxf2geo: dxf2geo.c message.c
- ${CXX} ${OPT_FLAGS} -o ../bin/dxf2geo -I../DataStr\
+ ${CXX} ${OPTIM} -o ../bin/dxf2geo -I../DataStr\
dxf2geo.c message.c ../lib/libGmshDataStr.a -lm
mysolver: mysolver.c
- ${CXX} ${OPT_FLAGS} -o mysolver.exe mysolver.c GmshClient.c
+ ${CXX} ${OPTIM} -o mysolver.exe mysolver.c GmshClient.c
mysolvertgz:
gtar zcvf mysolver.tgz mysolver.c mysolver.opt GmshClient.c GmshClient.h
diff --git a/variables.in b/variables.in
index d2ae9b49b2973806f181ad386b73cd771fd77e1b..9bf9eebca47f9c84a92309c38b1f1fb5de6f1300 100644
--- a/variables.in
+++ b/variables.in
@@ -7,19 +7,13 @@ CXX=@CXX@
LINKER=@LINKER@
# Compiler flags
-OPT_FLAGS=@OPT_FLAGS@
-OS_FLAGS=@OS_FLAGS@
-VERSION_FLAGS=@VERSION_FLAGS@
-GUI_INCLUDE=@GUI_INCLUDE@
+FLAGS=@FLAGS@
+OPTIM=@OPTIM@
# Gmsh subdirectories and libraries
GMSH_DIRS=@GMSH_DIRS@
GMSH_LIBS=@GMSH_LIBS@
-# GSL include and library path
-GSL_INCLUDE=@GSL_INCLUDE@
-GSL_LIBS=@GSL_LIBS@
-
# How you create a static library on this machine
AR=@AR@
RANLIB=@RANLIB@