Skip to content
Snippets Groups Projects
Commit 913774dc authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

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

parent ffb6c8f9
No related branches found
No related tags found
No related merge requests found
// $Id: Adapt.cpp,v 1.7 2001-04-26 17:58:00 remacle Exp $
#include "Gmsh.h"
#include "Adapt.h"
#include "nrutil.h"
#define TOL 1.e-08
#define MAXDEG 999
static int NN ;
static double MINH , *ERR , *HH , *PP , E0, DIM ;
/* ------------------------------------------------------------------------ */
/* f XXX */
/* ------------------------------------------------------------------------ */
/* h-type version 1 : minimize the number of elements while keeping a
given global error */
double fH1 (double l){
int i;
double val1 = 0.0, val2 = 0.0;
for(i = 1 ; i <= NN ; i++){
val1 += pow(2.*l*DSQR(ERR[i])*PP[i]/DIM, DIM/(2.*PP[i]+DIM));
val2 += DSQR(ERR[i]) * pow(2.*l*DSQR(ERR[i])*PP[i]/DIM, -2.*PP[i]/(2.*PP[i]+DIM));
}
return( -(val1 + l * (val2 - DSQR(E0))) );
}
/* h-type version 2 : minimize the error while keeping a given number
of elements */
double fH2 (double l){
int i;
double val1 = 0.0, val2 = 0.0, qi;
for(i = 1 ; i <= NN ; i++){
qi = pow(DIM*l/(2.*PP[i]*DSQR(ERR[i])), -DIM/(DIM+2.*PP[i]));
val1 += DSQR(ERR[i]) * pow(qi, -2.*PP[i]/DIM);
val2 += qi;
}
return( -(val1 + l * (val2 - E0)) );
}
/* p-type : minimize error by modifying the interpolation order vector */
double fP1 (double l){
int i;
double val1 = 0.0, val2 = 0.0, qi, e;
for(i = 1 ; i <= NN ; i++){
e = ERR[i];
if(e == 0.0) e=1.e-12;
qi = - log(2.*l*log(HH[i]/MINH)*DSQR(e)) / log(HH[i]/MINH);
val1 -= 0.5 * qi;
val2 += pow(HH[i]/MINH, qi) * DSQR(e);
}
return( -(val1 + l * (val2 - DSQR(E0))) );
}
/* ------------------------------------------------------------------------ */
/* A d a p t */
/* ------------------------------------------------------------------------ */
double min1d (int method, 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));
switch(method){
case ADAPT_H1:
case ADAPT_P1: ax=1.e-12; xx=1.e2; break;
default: ax=1.e-15; xx=1.e-12; break;
}
mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,funct);
return( brent(ax,xx,bx,funct,TOL,xmin) );
}
/* Adapt return the constraint (N0 ou e0) for the optimzed problem */
double AdaptMesh (int N, /* Number of elements */
int method, /* ADAPT_H1, ADAPT_H2, ADAPT_P1 or ADAPT_P2 */
int dim, /* 2 or 3 */
double *err, /* elementary errors */
double *h, /* elementary mesh sizes */
double *p, /* elementary exponents */
double e0 /* prescribed error or number of elements */){
int i;
double contr=0.0, pivrai, lambda, minf, qi, ri, pi, obj, obj2, minri, maxri;
double errmin, errmax;
h[N+1] = 1.0;
p[N+1] = 1.0;
NN = N;
ERR = err;
HH = h;
PP = p;
E0 = e0;
DIM = (double)dim;
for(i = 1 ; i <= N ; i++){
if(i == 1)
errmin = errmax = err[i];
else{
errmin = DMIN(errmin, err[i]);
errmax = DMAX(errmax, err[i]);
}
}
switch (method) {
case ADAPT_H1 :
minf = min1d (method, fH1, &lambda);
obj = 0.0;
for(i = 1 ; i <= N ; i++){
qi = pow(2.*lambda*DSQR(err[i])*p[i]/DIM, DIM/(2.*p[i]+DIM));
ri = pow(qi,1./DIM);
if(i==1) minri = maxri = ri;
if(err[i] == 0.0) ri = .5;
minri = DMIN(minri, ri);
maxri = DMAX(maxri, ri);
obj += DSQR(err[i]) * pow(ri, -2.*p[i]) ;
h[i] = sqrt(2.) * h[i]/ri;
p[i] = ri;
}
contr = fabs(minf);
printf("H-Refinement 1, Error %g=>%g, Objective %g, Reduction Factor %g->%g",
e0, sqrt(obj), -minf, minri, maxri);
break;
case ADAPT_H2 :
minf = min1d (method, fH2, &lambda);
obj = 0.0;
for(i = 1 ; i <= N ; i++){
qi = pow((DIM*lambda)/(2.*DSQR(err[i])*p[i]), -DIM/(DIM+2.*p[i]));
ri = pow(qi, 1./DIM);
if(i == 1) minri = maxri = ri;
minri = DMIN(minri, ri);
maxri = DMAX(maxri, ri);
obj += pow(ri,DIM) ;
h[i] = h[i]/ri;
p[i] = p[i];
}
contr = sqrt(fabs(minf));
printf( "H-Refinement 2, Elements %g=>%g, Objective %g, Reduction Factor %g->%g",
e0, obj, 100. * sqrt(fabs(minf)), minri, maxri);
break;
case ADAPT_P1 :
MINH = h[1];
for(i = 1 ; i <= N ; i++) MINH =DMIN(h[i], MINH);
MINH /= 2.;
minf = min1d (method, fP1, &lambda);
obj = obj2 = 0.0;
for(i = 1 ; i <= N ; i++){
qi = -log(2.*lambda*DSQR(err[i])*log(h[i]/MINH)) / log(h[i]/MINH);
pi = p[i] - .5 * qi;
pivrai = DMIN(DMAX(1., (double)(int)(pi+.99)), MAXDEG);
obj2 += pow(h[i]/MINH, 2.*(p[i]-pivrai)) * DSQR(err[i]);
obj += DSQR(err[i]) * pow(h[i]/MINH, qi);
h[i] = h[i];
p[i] = pi;
}
contr = fabs(minf);
printf("P-Refinement, Error %g=%g=>%g, Objective %g",
e0, sqrt(obj), sqrt(obj2), minf);
break;
case ADAPT_P2 :
minf = min1d (method, fH1, &lambda);
break;
default :
printf("Unknown Adaption Method");
}
return(contr) ;
}
#ifndef _ADAPT_H_
#define _ADAPT_H_
#define ADAPT_P1 1
#define ADAPT_P2 2
#define ADAPT_H1 3
#define ADAPT_H2 4
double AdaptMesh (int N, /* Number of elements */
int method, /* ADAPT_H1, ADAPT_H2, ADAPT_P1 or ADAPT_P2 */
int dim, /* 2 or 3 */
double *err, /* elementary errors */
double *h, /* elementary mesh sizes */
double *p, /* elementary exponents */
double e0 /* prescribed error or number of elements */);
#endif
# $Id: Makefile,v 1.21 2001-08-11 23:32:15 geuzaine Exp $
#
# Makefile for "libAdapt.a"
#
.IGNORE:
CC = c++
AR = ar ruvs
RM = rm
RANLIB = ranlib
LIB = ../lib/libAdapt.a
INCLUDE = -I../Common -I../DataStr
C_FLAGS = -g -Wall
OS_FLAGS =
VERSION_FLAGS =
RMFLAGS = -f
CFLAGS = $(C_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE)
SRC = Adapt.cpp \
mnbrak.cpp \
brent.cpp \
nrutil.cpp \
dsvdcmp.cpp \
newt.cpp \
fmin.cpp \
fdjac.cpp \
lnsrch.cpp \
lubksb.cpp \
ludcmp.cpp
OBJ = $(SRC:.cpp=.o)
.SUFFIXES: .o .cpp
$(LIB): $(OBJ)
$(AR) $(LIB) $(OBJ)
$(RANLIB) $(LIB)
.cpp.o:
$(CC) $(CFLAGS) -c $<
clean:
$(RM) $(RMFLAGS) *.o
lint:
$(LINT) $(CFLAGS) $(SRC)
depend:
(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
$(CC) -MM $(CFLAGS) ${SRC} \
) >Makefile.new
cp Makefile Makefile.bak
cp Makefile.new Makefile
$(RM) $(RMFLAGS) Makefile.new
# DO NOT DELETE THIS LINE
Adapt.o: Adapt.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h Adapt.h nrutil.h \
../Common/Numeric.h
mnbrak.o: mnbrak.cpp nrutil.h ../Common/Numeric.h
brent.o: brent.cpp nrutil.h ../Common/Numeric.h
nrutil.o: nrutil.cpp
dsvdcmp.o: dsvdcmp.cpp nrutil.h ../Common/Numeric.h
newt.o: newt.cpp nrutil.h ../Common/Numeric.h
fmin.o: fmin.cpp nrutil.h ../Common/Numeric.h
fdjac.o: fdjac.cpp nrutil.h ../Common/Numeric.h
lnsrch.o: lnsrch.cpp nrutil.h ../Common/Numeric.h
lubksb.o: lubksb.cpp
ludcmp.o: ludcmp.cpp nrutil.h ../Common/Numeric.h
// $Id: brent.cpp,v 1.4 2001-01-08 08:05:39 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 100
#define CGOLD 0.3819660
#define ZEPS 1.0e-10
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
double
brent (double ax, double bx, double cx, double (*f) (double), double tol,
double *xmin)
{
int iter;
double a, b, d=0.0, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w,
x, xm;
double e = 0.0;
a = (ax < cx ? ax : cx);
b = (ax > cx ? ax : cx);
x = w = v = bx;
fw = fv = fx = (*f) (x);
for (iter = 1; iter <= ITMAX; iter++)
{
xm = 0.5 * (a + b);
tol2 = 2.0 * (tol1 = tol * fabs (x) + ZEPS);
if (fabs (x - xm) <= (tol2 - 0.5 * (b - a)))
{
*xmin = x;
return fx;
}
if (fabs (e) > tol1)
{
r = (x - w) * (fx - fv);
q = (x - v) * (fx - fw);
p = (x - v) * q - (x - w) * r;
q = 2.0 * (q - r);
if (q > 0.0)
p = -p;
q = fabs (q);
etemp = e;
e = d;
if (fabs (p) >= fabs (0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x))
d = CGOLD * (e = (x >= xm ? a - x : b - x));
else
{
d = p / q;
u = x + d;
if (u - a < tol2 || b - u < tol2)
d = SIGN (tol1, xm - x);
}
}
else
{
d = CGOLD * (e = (x >= xm ? a - x : b - x));
}
u = (fabs (d) >= tol1 ? x + d : x + SIGN (tol1, d));
fu = (*f) (u);
if (fu <= fx)
{
if (u >= x)
a = x;
else
b = x;
SHFT (v, w, x, u)
SHFT (fv, fw, fx, fu)
}
else
{
if (u < x)
a = u;
else
b = u;
if (fu <= fw || w == x)
{
v = w;
w = u;
fv = fw;
fw = fu;
}
else if (fu <= fv || v == x || v == w)
{
v = u;
fv = fu;
}
}
}
nrerror ("Too many iterations in brent");
*xmin = x;
return fx;
}
#undef ITMAX
#undef CGOLD
#undef ZEPS
#undef SHFT
#undef NRANSI
/* (C) Copr. 1986-92 Numerical Recipes Software J!0. */
// $Id: dsvdcmp.cpp,v 1.4 2001-01-08 08:05:39 geuzaine Exp $
#include <math.h>
#include "nrutil.h"
double dpythag(double a, double b)
{
double absa,absb;
absa=fabs(a);
absb=fabs(b);
if (absa > absb) return absa*sqrt(1.0+SQR(absb/absa));
else return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+SQR(absa/absb)));
}
void dsvdcmp(double **a, int m, int n, double w[], double **v)
{
double dpythag(double a, double b);
int flag,i,its,j,jj,k,l,nm;
double anorm,c,f,g,h,s,scale,x,y,z,*rv1;
rv1=dvector(1,n);
g=scale=anorm=0.0;
for (i=1;i<=n;i++) {
l=i+1;
rv1[i]=scale*g;
g=s=scale=0.0;
if (i <= m) {
for (k=i;k<=m;k++) scale += fabs(a[k][i]);
if (scale) {
for (k=i;k<=m;k++) {
a[k][i] /= scale;
s += a[k][i]*a[k][i];
}
f=a[i][i];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][i]=f-g;
for (j=l;j<=n;j++) {
for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j];
f=s/h;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (k=i;k<=m;k++) a[k][i] *= scale;
}
}
w[i]=scale *g;
g=s=scale=0.0;
if (i <= m && i != n) {
for (k=l;k<=n;k++) scale += fabs(a[i][k]);
if (scale) {
for (k=l;k<=n;k++) {
a[i][k] /= scale;
s += a[i][k]*a[i][k];
}
f=a[i][l];
g = -SIGN(sqrt(s),f);
h=f*g-s;
a[i][l]=f-g;
for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
for (j=l;j<=m;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k];
for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
}
for (k=l;k<=n;k++) a[i][k] *= scale;
}
}
anorm=DMAX(anorm,(fabs(w[i])+fabs(rv1[i])));
}
for (i=n;i>=1;i--) {
if (i < n) {
if (g) {
for (j=l;j<=n;j++) v[j][i]=(a[i][j]/a[i][l])/g;
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j];
for (k=l;k<=n;k++) v[k][j] += s*v[k][i];
}
}
for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
g=rv1[i];
l=i;
}
for (i=IMIN(m,n);i>=1;i--) {
l=i+1;
g=w[i];
for (j=l;j<=n;j++) a[i][j]=0.0;
if (g) {
g=1.0/g;
for (j=l;j<=n;j++) {
for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j];
f=(s/a[i][i])*g;
for (k=i;k<=m;k++) a[k][j] += f*a[k][i];
}
for (j=i;j<=m;j++) a[j][i] *= g;
} else for (j=i;j<=m;j++) a[j][i]=0.0;
++a[i][i];
}
for (k=n;k>=1;k--) {
for (its=1;its<=30;its++) {
flag=1;
for (l=k;l>=1;l--) {
nm=l-1;
if ((double)(fabs(rv1[l])+anorm) == anorm) {
flag=0;
break;
}
if ((double)(fabs(w[nm])+anorm) == anorm) break;
}
if (flag) {
c=0.0;
s=1.0;
for (i=l;i<=k;i++) {
f=s*rv1[i];
rv1[i]=c*rv1[i];
if ((double)(fabs(f)+anorm) == anorm) break;
g=w[i];
h=dpythag(f,g);
w[i]=h;
h=1.0/h;
c=g*h;
s = -f*h;
for (j=1;j<=m;j++) {
y=a[j][nm];
z=a[j][i];
a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
z=w[k];
if (l == k) {
if (z < 0.0) {
w[k] = -z;
for (j=1;j<=n;j++) v[j][k] = -v[j][k];
}
break;
}
if (its == 30) nrerror("no convergence in 30 dsvdcmp iterations");
x=w[l];
nm=k-1;
y=w[nm];
g=rv1[nm];
h=rv1[k];
f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=dpythag(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
c=s=1.0;
for (j=l;j<=nm;j++) {
i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g=c*g;
z=dpythag(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g = g*c-x*s;
h=y*s;
y *= c;
for (jj=1;jj<=n;jj++) {
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=dpythag(f,h);
w[j]=z;
if (z) {
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for (jj=1;jj<=m;jj++) {
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
free_dvector(rv1,1,n);
}
/* cf. Numerical Recipes in C, p. 62 */
#define PREC 1.e-16
void invert_singular_matrix(double **M, int n, double **I){
double **V, **T, *W;
int i, j, k;
V = dmatrix(1,n,1,n);
T = dmatrix(1,n,1,n);
W = dvector(1,n);
dsvdcmp(M, n, n, W, V);
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
I[i][j] = 0.0 ;
T[i][j] = 0.0 ;
}
}
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
if(fabs(W[i]) > PREC){
T[i][j] += M[j][i] / W[i] ;
}
/*
else{
T[i][j] += 0.0 ;
}
*/
}
}
for(i=1 ; i<=n ; i++){
for(j=1 ; j<=n ; j++){
for(k=1 ; k<=n ; k++){
I[i][j] += V[i][k] * T[k][j] ;
}
}
}
free_dmatrix(V,1,n,1,n);
free_dmatrix(T,1,n,1,n);
free_dvector(W,1,n);
}
#undef PREC
// $Id: fdjac.cpp,v 1.4 2001-01-08 08:05:39 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define EPS 1.0e-4
void
fdjac (int n, float x[], float fvec[], float **df,
void (*vecfunc) (int, float[], float[]))
{
int i, j;
float h, temp, *f;
f = vector (1, n);
for (j = 1; j <= n; j++)
{
temp = x[j];
h = EPS * fabs (temp);
if (h == 0.0)
h = EPS;
x[j] = temp + h;
h = x[j] - temp;
(*vecfunc) (n, x, f);
x[j] = temp;
for (i = 1; i <= n; i++)
df[i][j] = (f[i] - fvec[i]) / h;
}
free_vector (f, 1, n);
}
#undef EPS
#undef NRANSI
// $Id: fmin.cpp,v 1.3 2001-01-08 08:05:40 geuzaine Exp $
#define NRANSI
#include "nrutil.h"
extern int nn;
extern float *fvec;
extern void (*nrfuncv) (int n, float v[], float f[]);
float
fmin (float x[])
{
int i;
float sum;
(*nrfuncv) (nn, x, fvec);
for (sum = 0.0, i = 1; i <= nn; i++)
sum += SQR (fvec[i]);
return 0.5 * sum;
}
#undef NRANSI
// $Id: lnsrch.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define ALF 1.0e-4
#define TOLX 1.0e-7
void
lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
float *f, float stpmax, int *check, float (*func) (float[]))
{
int i;
float a, alam, alam2, alamin, b, disc, f2, fold2, rhs1, rhs2, slope,
sum, temp, test, tmplam;
*check = 0;
for (sum = 0.0, i = 1; i <= n; i++)
sum += p[i] * p[i];
sum = sqrt (sum);
if (sum > stpmax)
for (i = 1; i <= n; i++)
p[i] *= stpmax / sum;
for (slope = 0.0, i = 1; i <= n; i++)
slope += g[i] * p[i];
test = 0.0;
for (i = 1; i <= n; i++)
{
temp = fabs (p[i]) / FMAX (fabs (xold[i]), 1.0);
if (temp > test)
test = temp;
}
alamin = TOLX / test;
alam = 1.0;
for (;;)
{
for (i = 1; i <= n; i++)
x[i] = xold[i] + alam * p[i];
*f = (*func) (x);
if (alam < alamin)
{
for (i = 1; i <= n; i++)
x[i] = xold[i];
*check = 1;
return;
}
else if (*f <= fold + ALF * alam * slope)
return;
else
{
if (alam == 1.0)
tmplam = -slope / (2.0 * (*f - fold - slope));
else
{
rhs1 = *f - fold - alam * slope;
rhs2 = f2 - fold2 - alam2 * slope;
a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
if (a == 0.0)
tmplam = -slope / (2.0 * b);
else
{
disc = b * b - 3.0 * a * slope;
if (disc < 0.0)
nrerror ("Roundoff problem in lnsrch.");
else
tmplam = (-b + sqrt (disc)) / (3.0 * a);
}
if (tmplam > 0.5 * alam)
tmplam = 0.5 * alam;
}
}
alam2 = alam;
f2 = *f;
fold2 = fold;
alam = FMAX (tmplam, 0.1 * alam);
}
}
#undef ALF
#undef TOLX
#undef NRANSI
// $Id: lubksb.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
void
lubksb (float **a, int n, int *indx, float b[])
{
int i, ii = 0, ip, j;
float sum;
for (i = 1; i <= n; i++)
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if (ii)
for (j = ii; j <= i - 1; j++)
sum -= a[i][j] * b[j];
else if (sum)
ii = i;
b[i] = sum;
}
for (i = n; i >= 1; i--)
{
sum = b[i];
for (j = i + 1; j <= n; j++)
sum -= a[i][j] * b[j];
b[i] = sum / a[i][i];
}
}
// $Id: ludcmp.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define TINY 1.0e-20;
void
ludcmp (float **a, int n, int *indx, float *d)
{
int i, imax, j, k;
float big, dum, sum, temp;
float *vv;
vv = vector (1, n);
*d = 1.0;
for (i = 1; i <= n; i++)
{
big = 0.0;
for (j = 1; j <= n; j++)
if ((temp = fabs (a[i][j])) > big)
big = temp;
if (big == 0.0)
nrerror ("Singular matrix in routine ludcmp");
vv[i] = 1.0 / big;
}
for (j = 1; j <= n; j++)
{
for (i = 1; i < j; i++)
{
sum = a[i][j];
for (k = 1; k < i; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
}
big = 0.0;
for (i = j; i <= n; i++)
{
sum = a[i][j];
for (k = 1; k < j; k++)
sum -= a[i][k] * a[k][j];
a[i][j] = sum;
if ((dum = vv[i] * fabs (sum)) >= big)
{
big = dum;
imax = i;
}
}
if (j != imax)
{
for (k = 1; k <= n; k++)
{
dum = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j] = imax;
if (a[j][j] == 0.0)
a[j][j] = TINY;
if (j != n)
{
dum = 1.0 / (a[j][j]);
for (i = j + 1; i <= n; i++)
a[i][j] *= dum;
}
}
free_vector (vv, 1, n);
}
#undef TINY
#undef NRANSI
// $Id: mnbrak.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define GOLD 1.618034
#define GLIMIT 100.0
#define TINY 1.0e-20
#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
void
mnbrak (double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
double (*func) (double))
{
double ulim, u, r, q, fu, dum;
*fa = (*func) (*ax);
*fb = (*func) (*bx);
if (*fb > *fa)
{
SHFT (dum, *ax, *bx, dum)
SHFT (dum, *fb, *fa, dum)
}
*cx = (*bx) + GOLD * (*bx - *ax);
*fc = (*func) (*cx);
while (*fb > *fc)
{
r = (*bx - *ax) * (*fb - *fc);
q = (*bx - *cx) * (*fb - *fa);
u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
(2.0 * SIGN (FMAX (fabs (q - r), TINY), q - r));
ulim = (*bx) + GLIMIT * (*cx - *bx);
if ((*bx - u) * (u - *cx) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
*ax = (*bx);
*bx = u;
*fa = (*fb);
*fb = fu;
return;
}
else if (fu > *fb)
{
*cx = u;
*fc = fu;
return;
}
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
else if ((*cx - u) * (u - ulim) > 0.0)
{
fu = (*func) (u);
if (fu < *fc)
{
SHFT (*bx, *cx, u, *cx + GOLD * (*cx - *bx))
SHFT (*fb, *fc, fu, (*func) (u))
}
}
else if ((u - ulim) * (ulim - *cx) >= 0.0)
{
u = ulim;
fu = (*func) (u);
}
else
{
u = (*cx) + GOLD * (*cx - *bx);
fu = (*func) (u);
}
SHFT (*ax, *bx, *cx, u)
SHFT (*fa, *fb, *fc, fu)
}
}
#undef GOLD
#undef GLIMIT
#undef TINY
#undef SHFT
#undef NRANSI
// $Id: newt.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
#include <math.h>
#define NRANSI
#include "nrutil.h"
#define MAXITS 200
#define TOLF 1.0e-4
#define TOLMIN 1.0e-6
#define TOLX 1.0e-7
#define STPMX 100.0
int nn;
float *fvec;
void (*nrfuncv) (int n, float v[], float f[]);
#define FREERETURN {free_vector(fvec,1,n);free_vector(xold,1,n);\
free_vector(p,1,n);free_vector(g,1,n);free_matrix(fjac,1,n,1,n);\
free_ivector(indx,1,n);return;}
void
newt (float x[], int n, int *check,
void (*vecfunc) (int, float[], float[]))
{
void fdjac (int n, float x[], float fvec[], float **df,
void (*vecfunc) (int, float[], float[]));
float fmin (float x[]);
void lnsrch (int n, float xold[], float fold, float g[], float p[], float x[],
float *f, float stpmax, int *check, float (*func) (float[]));
void lubksb (float **a, int n, int *indx, float b[]);
void ludcmp (float **a, int n, int *indx, float *d);
int i, its, j, *indx;
float d, den, f, fold, stpmax, sum, temp, test, **fjac, *g, *p, *xold;
indx = ivector (1, n);
fjac = matrix (1, n, 1, n);
g = vector (1, n);
p = vector (1, n);
xold = vector (1, n);
fvec = vector (1, n);
nn = n;
nrfuncv = vecfunc;
f = fmin (x);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < 0.01 * TOLF)
FREERETURN
for (sum = 0.0, i = 1; i <= n; i++)
sum += SQR (x[i]);
stpmax = STPMX * FMAX (sqrt (sum), (float) n);
for (its = 1; its <= MAXITS; its++)
{
fdjac (n, x, fvec, fjac, vecfunc);
for (i = 1; i <= n; i++)
{
for (sum = 0.0, j = 1; j <= n; j++)
sum += fjac[j][i] * fvec[j];
g[i] = sum;
}
for (i = 1; i <= n; i++)
xold[i] = x[i];
fold = f;
for (i = 1; i <= n; i++)
p[i] = -fvec[i];
ludcmp (fjac, n, indx, &d);
lubksb (fjac, n, indx, p);
lnsrch (n, xold, fold, g, p, x, &f, stpmax, check, fmin);
test = 0.0;
for (i = 1; i <= n; i++)
if (fabs (fvec[i]) > test)
test = fabs (fvec[i]);
if (test < TOLF)
{
*check = 0;
FREERETURN
}
if (*check)
{
test = 0.0;
den = FMAX (f, 0.5 * n);
for (i = 1; i <= n; i++)
{
temp = fabs (g[i]) * FMAX (fabs (x[i]), 1.0) / den;
if (temp > test)
test = temp;
}
*check = (test < TOLMIN ? 1 : 0);
FREERETURN
}
test = 0.0;
for (i = 1; i <= n; i++)
{
temp = (fabs (x[i] - xold[i])) / FMAX (fabs (x[i]), 1.0);
if (temp > test)
test = temp;
}
if (test < TOLX)
FREERETURN
}
nrerror ("MAXITS exceeded in newt");
}
#undef MAXITS
#undef TOLF
#undef TOLMIN
#undef TOLX
#undef STPMX
#undef FREERETURN
#undef NRANSI
// $Id: nrutil.cpp,v 1.4 2001-01-08 08:05:40 geuzaine Exp $
#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <malloc.h>
#define NR_END 1
#define FREE_ARG char*
void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;
v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}
int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;
v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}
unsigned char *cvector(long nl, long nh)
/* allocate an unsigned char vector with subscript range v[nl..nh] */
{
unsigned char *v;
v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char)));
if (!v) nrerror("allocation failure in cvector()");
return v-nl+NR_END;
}
unsigned long *lvector(long nl, long nh)
/* allocate an unsigned long vector with subscript range v[nl..nh] */
{
unsigned long *v;
v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long)));
if (!v) nrerror("allocation failure in lvector()");
return v-nl+NR_END;
}
double *dvector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
double *v;
v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
if (!v) nrerror("allocation failure in dvector()");
return v-nl+NR_END;
}
float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
double **dmatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
double **m;
/* allocate pointers to rows */
m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
int **imatrix(long nrl, long nrh, long ncl, long nch)
/* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
int **m;
/* allocate pointers to rows */
m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;
/* allocate rows and set pointers to them */
m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int)));
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
float **m;
/* allocate array of pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in submatrix()");
m += NR_END;
m -= newrl;
/* set pointers to rows */
for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix
declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1
and ncol=nch-ncl+1. The routine should be called with the address
&a[0][0] as the first argument. */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;
/* allocate pointers to rows */
m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure in convert_matrix()");
m += NR_END;
m -= nrl;
/* set pointers to rows */
m[nrl]=a-ncl;
for(i=1,j=nrl+1;i<nrow;i++,j++) m[j]=m[j-1]+ncol;
/* return pointer to array of pointers to rows */
return m;
}
float ***f3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
/* allocate a float 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] */
{
long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
float ***t;
/* allocate pointers to pointers to rows */
t=(float ***) malloc((size_t)((nrow+NR_END)*sizeof(float**)));
if (!t) nrerror("allocation failure 1 in f3tensor()");
t += NR_END;
t -= nrl;
/* allocate pointers to rows and set pointers to them */
t[nrl]=(float **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float*)));
if (!t[nrl]) nrerror("allocation failure 2 in f3tensor()");
t[nrl] += NR_END;
t[nrl] -= ncl;
/* allocate rows and set pointers to them */
t[nrl][ncl]=(float *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(float)));
if (!t[nrl][ncl]) nrerror("allocation failure 3 in f3tensor()");
t[nrl][ncl] += NR_END;
t[nrl][ncl] -= ndl;
for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
for(i=nrl+1;i<=nrh;i++) {
t[i]=t[i-1]+ncol;
t[i][ncl]=t[i-1][ncl]+ncol*ndep;
for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
}
/* return pointer to array of pointers to rows */
return t;
}
void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_cvector(unsigned char *v, long nl, long nh)
/* free an unsigned char vector allocated with cvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_lvector(unsigned long *v, long nl, long nh)
/* free an unsigned long vector allocated with lvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_dvector(double *v, long nl, long nh)
/* free a double vector allocated with dvector() */
{
free((FREE_ARG) (v+nl-NR_END));
}
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_dmatrix(double **m, long nrl, long nrh, long ncl, long nch)
/* free a double matrix allocated by dmatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_imatrix(int **m, long nrl, long nrh, long ncl, long nch)
/* free an int matrix allocated by imatrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_submatrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a submatrix allocated by submatrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_convert_matrix(float **b, long nrl, long nrh, long ncl, long nch)
/* free a matrix allocated by convert_matrix() */
{
free((FREE_ARG) (b+nrl-NR_END));
}
void free_f3tensor(float ***t, long nrl, long nrh, long ncl, long nch,
long ndl, long ndh)
/* free a float f3tensor allocated by f3tensor() */
{
free((FREE_ARG) (t[nrl][ncl]+ndl-NR_END));
free((FREE_ARG) (t[nrl]+ncl-NR_END));
free((FREE_ARG) (t+nrl-NR_END));
}
#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_ */
// $Id: Main.cpp,v 1.10 2001-08-11 23:28:31 geuzaine Exp $
#include <signal.h>
#include "ParUtil.h"
#include <signal.h>
#if !defined(WIN32) || defined(__CYGWIN__)
#include <sys/resource.h>
#endif
#include "PluginManager.h"
#include "Gmsh.h"
#include "GmshVersion.h"
#include "Numeric.h"
#include "Geo.h"
#include "Mesh.h"
#include "Views.h"
#include "Parser.h"
#include "Context.h"
#include "Options.h"
#include "OpenFile.h"
#include "GetOptions.h"
#include "MinMax.h"
#include "Static.h"
/* dummy defs for link purposes */
void AddViewInUI(int, char *, int){}
void draw_polygon_2d (double, double, double, int, double *, double *, double *){}
void set_r(int, double){}
void Draw(void){}
void DrawUI(void){}
void Replot(void){}
void CreateOutputFile(char *, int){}
/* ------------------------------------------------------------------------ */
/* I n f o */
/* ------------------------------------------------------------------------ */
void Info (int level, char *arg0){
switch(level){
case 0 :
if(ParUtil::Instance()->master())
{
fprintf(stderr, "%s\n", gmsh_progname);
fprintf(stderr, "%s\n", gmsh_copyright);
Print_Usage(arg0);
}
ParUtil::Instance()->Exit();
case 1:
if(ParUtil::Instance()->master())
fprintf(stderr, "%.2f\n", GMSH_VERSION);
ParUtil::Instance()->Exit();
case 2:
if(ParUtil::Instance()->master())
{
fprintf(stderr, "%s%.2f\n", gmsh_version, GMSH_VERSION);
fprintf(stderr, "%s\n", gmsh_os);
fprintf(stderr, "%s\n", gmsh_date);
fprintf(stderr, "%s\n", gmsh_host);
fprintf(stderr, "%s\n", gmsh_packager);
fprintf(stderr, "%s\n", gmsh_url);
fprintf(stderr, "%s\n", gmsh_email);
}
ParUtil::Instance()->Exit();
default :
break;
}
}
/* ------------------------------------------------------------------------ */
/* m a i n */
/* ------------------------------------------------------------------------ */
int main(int argc, char *argv[]){
int i, nbf;
ParUtil::Instance()->init(argc,argv);
if(argc < 2) Info(0,argv[0]);
Init_Options(0);
Get_Options(argc, argv, &nbf);
M.Vertices = NULL ;
M.VertexEdges = NULL ;
M.Simplexes = NULL ;
M.Points = NULL ;
M.Curves = NULL ;
M.SurfaceLoops = NULL ;
M.EdgeLoops = NULL ;
M.Surfaces = NULL ;
M.Volumes = NULL ;
M.PhysicalGroups = NULL ;
M.Metric = NULL ;
signal(SIGINT, Signal);
signal(SIGSEGV, Signal);
signal(SIGFPE, Signal);
if(CTX.default_plugins)
GMSH_PluginManager::Instance()->RegisterDefaultPlugins();
OpenProblem(CTX.filename);
if(yyerrorstate)
ParUtil::Instance()->Abort();
else{
for(i=1;i<nbf;i++) MergeProblem(TheFileNameTab[i]);
if(TheBgmFileName){
MergeProblem(TheBgmFileName);
if(List_Nbr(Post_ViewList))
BGMWithView((Post_View*)List_Pointer(Post_ViewList, List_Nbr(Post_ViewList)-1));
else
fprintf(stderr, ERROR_STR "Invalid background mesh (no view)\n"); exit(1);
}
if(CTX.batch > 0){
mai3d(THEM, CTX.batch);
Print_Mesh(THEM,CTX.output_filename,CTX.mesh.format);
}
else
Print_Geo(THEM, CTX.output_filename);
if(CTX.mesh.histogram)
Print_Histogram(THEM->Histogram[0]);
ParUtil::Instance()->Barrier(__LINE__,__FILE__);
ParUtil::Instance()->Exit();
return 1;
}
ParUtil::Instance()->Barrier(__LINE__,__FILE__);
ParUtil::Instance()->Exit();
return 1;
}
/* ------------------------------------------------------------------------ */
/* S i g n a l */
/* ------------------------------------------------------------------------ */
void Signal (int sig_num){
switch (sig_num){
case SIGSEGV : Msg(FATAL, "Segmentation violation (invalid memory reference)"); break;
case SIGFPE : Msg(FATAL, "Floating point exception (division by zero?)"); break;
case SIGINT : Msg(FATAL, "Interrupt (generated from terminal special char)"); break;
default : Msg(FATAL, "Unknown signal"); break;
}
}
/* ------------------------------------------------------------------------ */
/* M s g */
/* ------------------------------------------------------------------------ */
void Msg(int level, char *fmt, ...){
va_list args;
int abort=0;
int nb, nbvis;
va_start (args, fmt);
switch(level){
case DIRECT :
if(ParUtil::Instance()->master())
vfprintf(stdout, fmt, args); fprintf(stdout, "\n");
break;
case FATAL :
case FATAL1 :
case FATAL2 :
case FATAL3 :
fprintf(stderr,"On processor %d : ",ParUtil::Instance()->rank());
fprintf(stderr, FATAL_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
abort = 1 ;
break ;
case GERROR :
case GERROR1 :
case GERROR2 :
case GERROR3 :
fprintf(stderr,"On processor %d : ",ParUtil::Instance()->rank());
fprintf(stderr, ERROR_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
abort = 1 ;
break ;
case WARNING :
case WARNING1 :
case WARNING2 :
case WARNING3 :
fprintf(stderr,"On processor %d : ",ParUtil::Instance()->rank());
fprintf(stderr, WARNING_STR);
vfprintf(stderr, fmt,args); fprintf(stderr, "\n");
break;
case PARSER_ERROR :
if(ParUtil::Instance()->master())
{
fprintf(stderr, PARSER_ERROR_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
}
break ;
case PARSER_INFO :
if(CTX.verbosity == 5 && ParUtil::Instance()->master()){
fprintf(stderr, PARSER_INFO_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
}
break ;
case DEBUG :
case DEBUG1 :
case DEBUG2 :
case DEBUG3 :
if(ParUtil::Instance()->master() && CTX.verbosity > 2){
fprintf(stderr, DEBUG_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
}
break;
default :
if(ParUtil::Instance()->master() && CTX.verbosity > 0){
fprintf(stderr, INFO_STR);
vfprintf(stderr, fmt, args); fprintf(stderr, "\n");
}
break;
}
va_end (args);
if(abort) exit(1);
}
/* ------------------------------------------------------------------------ */
/* C p u */
/* ------------------------------------------------------------------------ */
void GetResources(long *s, long *us, long *mem){
#if !defined(WIN32) || defined(__CYGWIN__)
static struct rusage r;
getrusage(RUSAGE_SELF,&r);
*s = (long)r.ru_utime.tv_sec ;
*us = (long)r.ru_utime.tv_usec ;
*mem = (long)r.ru_maxrss ;
#else
*s = *us = *mem = 0;
#endif
}
double Cpu(void){
long s, us, mem;
GetResources(&s, &us, &mem);
return (double)s + (double)us/1.e6 ;
}
/* ------------------------------------------------------------------------ */
/* P r o g r e s s */
/* ------------------------------------------------------------------------ */
void Progress(int i){
}
void AddALineInTheEditGeometryForm (char* line){
};
# $Id: Makefile,v 1.10 2001-08-12 13:08:20 geuzaine Exp $
#
# Makefile for "libBox.a"
#
.IGNORE:
CC = c++
AR = ar ruvs
RM = rm
RANLIB = ranlib
LIB = ../lib/libBox.a
INCLUDE = -I../Common -I../DataStr -I../Geo\
-I../Graphics -I../Mesh -I../Parser -I../Motif -I../Fltk -I../Plugin -I../Parallel
C_FLAGS = -g
OS_FLAGS =
VERSION_FLAGS =
RMFLAGS = -f
CFLAGS = $(C_FLAGS) $(OS_FLAGS) $(VERSION_FLAGS) $(INCLUDE)
SRC = Main.cpp
OBJ = $(SRC:.cpp=.o)
.SUFFIXES: .o .cpp
$(LIB): $(OBJ)
$(AR) $(LIB) $(OBJ)
$(RANLIB) $(LIB)
.cpp.o:
$(CC) $(CFLAGS) -c $<
clean:
$(RM) $(RMFLAGS) *.o
lint:
$(LINT) $(CFLAGS) $(SRC)
depend:
(sed '/^# DO NOT DELETE THIS LINE/q' Makefile && \
$(CC) -MM $(CFLAGS) ${SRC} \
) >Makefile.new
cp Makefile Makefile.bak
cp Makefile.new Makefile
$(RM) $(RMFLAGS) Makefile.new
# DO NOT DELETE THIS LINE
Main.o: Main.cpp ../Common/Gmsh.h ../Common/Message.h \
../DataStr/Malloc.h ../DataStr/List.h ../DataStr/Tree.h \
../DataStr/avl.h ../DataStr/Tools.h ../Common/GmshVersion.h \
../Common/Numeric.h ../Geo/Geo.h ../Mesh/Mesh.h ../Mesh/Vertex.h \
../Mesh/Simplex.h ../Mesh/Edge.h ../Geo/ExtrudeParams.h \
../Mesh/Metric.h ../Common/Views.h ../Common/ColorTable.h \
../Parser/Parser.h ../Common/Context.h ../Parser/OpenFile.h \
../Common/GetOptions.h ../Geo/MinMax.h ../Common/Static.h
#ifndef _BITMAPS_H_
#define _BITMAPS_H_
// 'Gmsh' (Unix) icon
#define g1_width 66
#define g1_height 29
static char g1_bits[] = {
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x03,0xfc,0x00,0x00,0x00,
0x00,0x00,0x00,0x80,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,
0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0xe0,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x70,0x00,0xfc,0x00,0x00,0x00,
0x00,0x00,0x10,0x30,0x00,0xfc,0x00,0x08,0x00,0x00,0x00,0x18,0x18,0x00,0xfc,
0x00,0x0c,0x00,0x04,0x00,0x3c,0x08,0x06,0xfc,0xc0,0x0f,0x04,0x8f,0x01,0x3e,
0x04,0x0f,0xfc,0x30,0x0e,0x8e,0xcf,0x03,0x7b,0x84,0x0f,0xfc,0x18,0x0e,0xcf,
0xe7,0x83,0x70,0x42,0x0f,0xfc,0x1c,0x8e,0x2f,0x97,0x43,0x70,0x22,0x0e,0xfd,
0x1e,0xcd,0x1e,0x8f,0x23,0x21,0x17,0x8e,0xfc,0xfe,0x68,0x8e,0x87,0x97,0x31,
0x0f,0x4e,0xfc,0x7e,0x38,0x87,0x83,0x8f,0x1b,0x07,0x2e,0xfc,0x1c,0x1c,0x83,
0x01,0x03,0x07,0x02,0x1e,0xfc,0x00,0x0e,0x00,0x00,0x00,0x00,0x00,0x0c,0xfc,
0x00,0x07,0x00,0x00,0x00,0x04,0x01,0x00,0xfc,0x80,0x03,0x00,0x00,0x00,0x8c,
0x01,0x00,0xfc,0xc0,0x03,0x00,0x00,0x00,0x8c,0x99,0x26,0xfd,0xe0,0x01,0x00,
0x00,0x00,0x54,0xa5,0x29,0xfd,0xf0,0x01,0x00,0x00,0x00,0x54,0xbd,0x28,0xfd,
0xf0,0x00,0x00,0x00,0x00,0x24,0x85,0x28,0xfd,0x78,0x00,0x00,0x00,0x00,0x24,
0xa5,0x28,0xfd,0x30,0x00,0x00,0x00,0x00,0x24,0x99,0xc8,0xfd,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0xfc};
// 'Gmsh Menu' (Unix) icon
#define g2_width 66
#define g2_height 29
static char g2_bits[] = {
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x03,0xfc,0x00,0x00,0x00,
0x00,0x00,0x00,0x80,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,
0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0xe0,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x70,0x00,0xfc,0x00,0x00,0x00,
0x00,0x00,0x10,0x30,0x00,0xfc,0x00,0x08,0x00,0x00,0x00,0x18,0x18,0x00,0xfc,
0x00,0x0c,0x00,0x04,0x00,0x3c,0x08,0x06,0xfc,0xc0,0x0f,0x04,0x8f,0x01,0x3e,
0x04,0x0f,0xfc,0x30,0x0e,0x8e,0xcf,0x03,0x7b,0x84,0x0f,0xfc,0x18,0x0e,0xcf,
0xe7,0x83,0x70,0x42,0x0f,0xfc,0x1c,0x8e,0x2f,0x97,0x43,0x70,0x22,0x0e,0xfd,
0x1e,0xcd,0x1e,0x8f,0x23,0x21,0x17,0x8e,0xfc,0xfe,0x68,0x8e,0x87,0x97,0x31,
0x0f,0x4e,0xfc,0x7e,0x38,0x87,0x83,0x8f,0x1b,0x07,0x2e,0xfc,0x1c,0x1c,0x83,
0x01,0x03,0x07,0x02,0x1e,0xfc,0x00,0x0e,0x00,0x00,0x00,0x00,0x00,0x0c,0xfc,
0x00,0x07,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x80,0x03,0x00,0x00,0x00,0x00,
0x00,0x00,0xfc,0xc0,0x03,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0xe0,0x01,0x00,
0x00,0x00,0x00,0x00,0x00,0xfc,0xf0,0x01,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,
0xf0,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x78,0x00,0x00,0x00,0x00,0x00,
0x00,0x00,0xfc,0x30,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0xfc};
// 'Gmsh command' (Unix Motif) icon
#ifdef _XMOTIF
#define g3_width 66
#define g3_height 29
static char g3_bits[] = {
0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0x00,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x00,0x03,0xfc,0x00,0x00,0x00,
0x00,0x00,0x00,0x80,0x03,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,
0x00,0x00,0x00,0x00,0x00,0x00,0xc0,0x01,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,
0xe0,0x00,0xfc,0x00,0x00,0x00,0x00,0x00,0x00,0x70,0x00,0xfc,0x00,0x00,0x00,
0x00,0x00,0x10,0x30,0x00,0xfc,0x00,0x08,0x00,0x00,0x00,0x18,0x18,0x00,0xfc,
0x00,0x0c,0x00,0x04,0x00,0x3c,0x08,0x06,0xfc,0xc0,0x0f,0x04,0x8f,0x01,0x3e,
0x04,0x0f,0xfc,0x30,0x0e,0x8e,0xcf,0x03,0x7b,0x84,0x0f,0xfc,0x18,0x0e,0xcf,
0xe7,0x83,0x70,0x42,0x0f,0xfc,0x1c,0x8e,0x2f,0x97,0x43,0x70,0x22,0x0e,0xfd,
0x1e,0xcd,0x1e,0x8f,0x23,0x21,0x17,0x8e,0xfc,0xfe,0x68,0x8e,0x87,0x97,0x31,
0x0f,0x4e,0xfc,0x7e,0x38,0x87,0x83,0x8f,0x1b,0x07,0x2e,0xfc,0x1c,0x1c,0x83,
0x01,0x03,0x07,0x02,0x1e,0xfc,0x00,0x0e,0x00,0x00,0x00,0x00,0x00,0x0c,0xfc,
0x00,0x07,0xe0,0x01,0x00,0x00,0x00,0x00,0xfd,0x80,0x03,0x10,0x02,0x00,0x00,
0x00,0x00,0xfd,0xc0,0x03,0x10,0x70,0x6e,0x6e,0x4e,0x63,0xfd,0xe0,0x01,0x10,
0x88,0x92,0x92,0xd0,0x94,0xfd,0xf0,0x01,0x10,0x88,0x92,0x92,0x5c,0x14,0xfd,
0xf0,0x00,0x10,0x88,0x92,0x92,0x52,0x14,0xfd,0x78,0x00,0x10,0x8a,0x92,0x92,
0x52,0x94,0xfd,0x30,0x00,0xe0,0x71,0x92,0x92,0x6c,0x64,0xfd,0x00,0x00,0x00,
0x00,0x00,0x00,0x00,0x00,0xfc};
#endif
// 'About Gmsh' bitmap
#define about_width 49
#define about_height 111
static char about_bits[] = {
0x00,0x00,0x00,0x00,0x00,0x00,0xfe,0x00,0x00,0x00,0x00,0x00,0x00,0xfe,0x00,
0x00,0x00,0x08,0x00,0x00,0xfe,0x00,0x00,0x00,0x10,0x00,0x00,0xfe,0x00,0x00,
0x00,0x20,0x00,0x00,0xfe,0x00,0x00,0x00,0x40,0x00,0x00,0xfe,0x00,0x00,0x00,
0xc0,0x00,0x00,0xfe,0x00,0x00,0x00,0xc0,0x01,0x00,0xfe,0x00,0x00,0x00,0xc0,
0x01,0x00,0xfe,0x00,0x00,0x00,0xff,0x03,0x00,0xfe,0x00,0x00,0xf0,0xff,0x03,
0x00,0xfe,0x00,0x00,0xf0,0xff,0x03,0x00,0xfe,0x70,0x00,0xf0,0xff,0x03,0x00,
0xfe,0xfc,0x03,0xe0,0xff,0x00,0x00,0xfe,0xf8,0x0f,0xc0,0x01,0x00,0x00,0xfe,
0xf0,0x3f,0xc0,0x00,0x00,0x00,0xfe,0xe0,0x7f,0x00,0x01,0x00,0x00,0xfe,0x80,
0xff,0x00,0x01,0x00,0x00,0xfe,0x00,0xfe,0x01,0x02,0x00,0x00,0xfe,0x00,0xf8,
0x03,0x0c,0x00,0x00,0xfe,0x00,0xf0,0x07,0x08,0x00,0x00,0xfe,0x00,0xc0,0x0f,
0x10,0x00,0x00,0xfe,0x00,0x00,0x1f,0x70,0x00,0x00,0xfe,0x00,0x00,0x7c,0xf8,
0x00,0x00,0xfe,0x00,0x00,0xf0,0xff,0x00,0x00,0xfe,0x00,0x00,0xc0,0xff,0x01,
0x00,0xfe,0x00,0x00,0x00,0xfe,0x01,0x00,0xfe,0x00,0x00,0x00,0xf8,0x00,0x00,
0xfe,0x00,0x00,0x00,0x00,0x00,0x00,0xfe,0x00,0x00,0x00,0x00,0x00,0x00,0xfe,
0x00,0x00,0x00,0x00,0x00,0x00,0xfe,0x00,0x00,0x00,0x00,0x00,0x00,0xfe,0x00,
0x00,0xe0,0x03,0x00,0x00,0xfe,0x00,0x00,0xf0,0x0f,0x00,0x00,0xfe,0x00,0x80,
0xf8,0x1f,0x00,0x00,0xfe,0x00,0x00,0xff,0x21,0x00,0x00,0xfe,0x00,0x00,0xff,
0x40,0x00,0x00,0xfe,0x00,0x00,0x7e,0x80,0x00,0x00,0xfe,0x00,0x00,0x3c,0x80,
0x00,0x00,0xfe,0x00,0x00,0x38,0xc0,0x01,0x00,0xfe,0x00,0x00,0x20,0xf0,0x01,
0x00,0xfe,0x00,0x00,0x40,0xf8,0x01,0x00,0xfe,0x00,0x00,0x80,0xf8,0x01,0x00,
0xfe,0x00,0x00,0x00,0xf1,0x00,0x00,0xfe,0x00,0x00,0x00,0x04,0x00,0x00,0xfe,
0x00,0x00,0x00,0x08,0x00,0x00,0xfe,0x00,0x00,0x00,0x08,0x00,0x00,0xfe,0x00,
0x00,0x00,0x10,0x00,0x00,0xfe,0x00,0x00,0x00,0x30,0x00,0x00,0xfe,0x00,0x00,
0x00,0x70,0x00,0x00,0xfe,0x00,0x00,0x00,0x70,0x00,0x00,0xfe,0x00,0x00,0x00,
0xf8,0x00,0x00,0xfe,0x00,0x00,0xc0,0xff,0x01,0x00,0xfe,0x00,0x00,0xf0,0xff,
0x01,0x00,0xfe,0x00,0x00,0xf0,0xff,0x01,0x00,0xfe,0x00,0x00,0xf0,0xff,0x00,
0x00,0xfe,0x00,0x00,0xe0,0x1f,0x00,0x00,0xfe,0x00,0x00,0xc0,0x01,0x00,0x00,
0xfe,0x00,0x00,0xc0,0x00,0x00,0x00,0xfe,0x00,0x00,0x80,0x00,0x00,0x00,0xfe,
0x00,0x00,0x00,0x01,0x00,0x00,0xfe,0x00,0x00,0x00,0x02,0x00,0x00,0xfe,0x00,
0x00,0x00,0x04,0x00,0x00,0xfe,0x00,0x00,0x38,0x08,0x00,0x00,0xfe,0x00,0x00,
0xf8,0x18,0x00,0x00,0xfe,0x00,0x00,0xf8,0x3f,0x00,0x00,0xfe,0x00,0x00,0xf8,
0x7f,0x00,0x00,0xfe,0x00,0x00,0xf0,0xff,0x00,0x00,0xfe,0x00,0x00,0xf0,0xff,
0x00,0x00,0xfe,0x00,0x00,0x60,0xfc,0x01,0x00,0xfe,0x00,0x00,0x40,0xf0,0x01,
0x00,0xfe,0x00,0x00,0x80,0x00,0x00,0x00,0xfe,0x00,0x00,0x80,0x00,0x00,0x00,
0xfe,0x00,0x00,0x00,0x01,0x00,0x00,0xfe,0x00,0x00,0x00,0x02,0x00,0x00,0xfe,
0x00,0x00,0x0c,0x06,0x00,0x00,0xfe,0x00,0x00,0x3c,0x04,0x00,0x00,0xfe,0x00,
0x00,0xf8,0x08,0x00,0x00,0xfe,0x00,0x00,0xf8,0x1f,0x00,0x00,0xfe,0x00,0x00,
0xf0,0x3f,0x00,0x00,0xfe,0x00,0x00,0xe0,0x7f,0x00,0x00,0xfe,0x00,0x00,0xc0,
0xff,0x00,0x00,0xfe,0x00,0x00,0x00,0xff,0x00,0x00,0xfe,0x00,0x00,0x00,0xfc,
0x00,0x00,0xfe,0x00,0x00,0x00,0x0c,0x00,0x00,0xfe,0x00,0x00,0x00,0x10,0x00,
0x00,0xfe,0x00,0x00,0x00,0x10,0x00,0x00,0xfe,0x00,0x00,0x00,0x20,0x00,0x00,
0xfe,0x00,0x00,0x00,0x60,0x00,0x00,0xfe,0x00,0x00,0x00,0xcf,0x00,0x00,0xfe,
0x00,0x00,0xde,0xbf,0x01,0x00,0xfe,0x00,0x00,0xfe,0xff,0x03,0x00,0xfe,0x00,
0x00,0xfe,0x1f,0x0f,0x00,0xfe,0x00,0x00,0xf8,0x07,0x1e,0x00,0xfe,0x00,0x00,
0xf0,0x03,0x7e,0x00,0xfe,0x00,0x00,0xd0,0x07,0xfc,0x00,0xfe,0x00,0x00,0x10,
0x08,0xf8,0x03,0xfe,0x00,0x00,0x10,0x18,0xf0,0x07,0xfe,0x00,0x00,0x10,0x30,
0xe0,0x0f,0xfe,0x00,0x00,0x10,0x30,0xc0,0x1f,0xfe,0x00,0x00,0x10,0x70,0x80,
0x3f,0xfe,0x00,0x00,0x20,0xf8,0x00,0x7f,0xfe,0x00,0x00,0x60,0xfc,0x00,0xfe,
0xfe,0x00,0x00,0xc0,0xff,0x01,0xfc,0xfe,0x00,0x00,0xc0,0xff,0x01,0x78,0xfe,
0x00,0x00,0x80,0xff,0x01,0x00,0xfe,0x00,0x00,0x00,0xff,0x01,0x00,0xfe,0x00,
0x00,0x00,0xfc,0x01,0x00,0xfe,0x00,0x00,0x00,0xf8,0x00,0x00,0xfe,0x00,0x00,
0x00,0x00,0x00,0x00,0xfe,0x00,0x00,0x00,0x00,0x00,0x00,0xfe};
// 'Abort' bitmap
#define abort_width 13
#define abort_height 13
static char abort_bits[] = {
0x00,0xe0,0x40,0xe0,0x40,0xe0,0x50,0xe1,0x48,0xe2,0x44,0xe4,0x44,0xe4,0x44,
0xe4,0x04,0xe4,0x04,0xe4,0x08,0xe2,0xf0,0xe1,0x00,0xe0};
// 'Play button' bitmap
#define start_width 9
#define start_height 13
static char start_bits[] = {
0x00,0xfe,0x06,0xfe,0x0a,0xfe,0x12,0xfe,0x22,0xfe,0x42,0xfe,0x82,0xfe,0x42,
0xfe,0x22,0xfe,0x12,0xfe,0x0a,0xfe,0x06,0xfe,0x00,0xfe};
// 'Pause button' bitmap
#define stop_width 9
#define stop_height 13
static char stop_bits[] = {
0x00,0xfe,0xee,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,
0xfe,0xaa,0xfe,0xaa,0xfe,0xaa,0xfe,0xee,0xfe,0x00,0xfe};
#endif
// $Id: ColorTable.cpp,v 1.3 2001-02-12 17:38:02 geuzaine Exp $
#include "Gmsh.h"
#include "ColorTable.h"
#include "Context.h"
extern Context_T CTX ;
void ColorTable_InitParam(int number, ColorTable *ct,
int rgb_flag, int alpha_flag){
ct->ipar[COLORTABLE_NUMBER] = number;
if(rgb_flag) {
ct->ipar[COLORTABLE_INVERT] = 0;
ct->ipar[COLORTABLE_SWAP] = 0;
ct->ipar[COLORTABLE_ROTATE] = 0;
ct->fpar[COLORTABLE_CURVE] = 0.0;
ct->fpar[COLORTABLE_BIAS] = 0.0;
ct->fpar[COLORTABLE_BETA] = 0.0;
}
if(alpha_flag) {
ct->fpar[COLORTABLE_ALPHAPOW] = 2.;
ct->fpar[COLORTABLE_ALPHAVAL] = 255.;
}
}
void ColorTable_Recompute(ColorTable *ct, int rgb_flag, int alpha_flag){
float curve, bias;
double gamma;
int i,r,g,b,a,rotate;
float s,t;
ct->ipar[COLORTABLE_CHANGED] = 1 ;
bias = ct->fpar[COLORTABLE_BIAS];
curve = ct->fpar[COLORTABLE_CURVE];
rotate = ct->ipar[COLORTABLE_ROTATE];
for (i=0 ; i<ct->size ; i++) {
if(i+rotate<0)
s = (float) (i+rotate+ct->size) / (float) (ct->size-1);
else if(i+rotate>ct->size-1)
s = (float) (i+rotate-ct->size) / (float) (ct->size-1);
else
s = (float) (i+rotate) / (float) (ct->size-1);
if(ct->ipar[COLORTABLE_SWAP]) s = 1.0 - s;
if (rgb_flag) {
switch(ct->ipar[COLORTABLE_NUMBER]){
case 1 : /* vis5d */
t = (curve+1.4) * (s - (1.+bias)/2.);
r = (int)( 128.0 + 127.0 * atan( 7.0*t ) / 1.57 );
g = (int)( 128.0 + 127.0 * (2 * exp(-7*t*t) - 1) );
b = (int)( 128.0 + 127.0 * atan( -7.0*t ) / 1.57 );
break;
case 2: /* samcef */
if (s-bias<=0.00){
r = 0 ; g = 0 ; b = 255 ;
}
else if(s-bias<=0.40){
r = 0 ; g = (int)((s-bias)*637.5) ; b = (int)(255.-(s-bias)*637.5) ;
}
else if(s-bias<=0.60){
r = (int)(1275.*(s-bias-0.4)) ; g = 255 ; b = 0 ;
}
else if(s-bias<=1.00){
r = 255 ; g = (int)(255.-637.5*(s-bias-0.6)) ; b = 0 ;
}
else {
r = 255 ; g = 0 ; b = 0 ;
}
break;
case 3: /* rainbow */
if (s-bias<=0.00) {
r = 0 ; g = 0 ; b = 255 ;
}
else if(s-bias<=0.25+curve){
r = 0 ; g = (int)((s-bias)*(255./(0.25+curve))) ; b = 255 ;
}
else if(s-bias<=0.50) {
r = 0 ; g = 255 ; b = (int)(255.-(255./(0.25-curve))*(s-bias-0.25-curve));
}
else if(s-bias<=0.75-curve){
r = (int)((s-bias-0.5)*(255./(0.25-curve))); g = 255 ; b = 0 ;
}
else if(s-bias<=1.00) {
r = 255; g = (int)(255.-(255./(0.25+curve))*(s-bias-0.75+curve)) ; b = 0 ;
}
else {
r = 255 ; g = 0 ; b = 0 ;
}
break;
case 4: /* blue-yellow-white */
#define myfct(a,b,c,d) ((a)+\
(b)*(s-bias)+\
(c)*(s-bias)*(s-bias)+\
(d)*(s-bias)*(s-bias)*(s-bias))
#define clamp(x) x = (x)<0?0:((x)>255?255:(x))
r = (int)(255. * myfct(-0.0506169, 2.81633, -1.87033, 0.0524573) );
g = (int)(255. * myfct(0.0485868, -1.26109, 6.3074, -4.12498) );
b = (int)(255. * myfct(0.364662, 1.50814, -7.36756, 6.51847 ) );
clamp(r); clamp(g); clamp(b);
#undef myfct
#undef clamp
break;
case 5: /* grayscale */
if (s-bias<=0.00){ r = g = b = 0 ; }
else if (s-bias<=1.00){ r = g = b = (int)(255*(s-bias)); }
else { r = g = b = 255 ; }
break;
case 6: /* monochrome */
r = g = b = 0 ;
break;
default: /* grayscale without white */
if (s-bias<=0.00){ r = g = b = 0 ; }
else if (s-bias<=1.00){ r = g = b = (int)(220*(s-bias)); }
else { r = g = b = 220 ; }
break;
}
if(ct->fpar[COLORTABLE_BETA]){
if(ct->fpar[COLORTABLE_BETA] > 0.0)
gamma = 1. - ct->fpar[COLORTABLE_BETA];
else
gamma = 1./(1.001 + ct->fpar[COLORTABLE_BETA]);
r = (int)( 255. * pow((double)r/255.,gamma) );
g = (int)( 255. * pow((double)g/255.,gamma) );
b = (int)( 255. * pow((double)b/255.,gamma) );
}
if(ct->ipar[COLORTABLE_INVERT]){
r = 255-r ;
g = 255-g ;
b = 255-b ;
}
}
else {
r = UNPACK_RED( ct->table[i] );
g = UNPACK_GREEN( ct->table[i] );
b = UNPACK_BLUE( ct->table[i] );
}
if (alpha_flag) {
if (ct->fpar[COLORTABLE_ALPHAVAL]<0) {
a = (int)( 255.0 * pow( s, ct->fpar[COLORTABLE_ALPHAPOW] ) );
}
else {
a = (int)( ct->fpar[COLORTABLE_ALPHAVAL] );
}
}
else {
a = UNPACK_ALPHA( ct->table[i] );
}
ct->table[i] = PACK_COLOR( r, g, b, a );
}
}
static ColorTable clip;
void ColorTable_Copy(ColorTable *ct){
clip.size = ct->size;
memcpy(clip.table, ct->table, ct->size * sizeof(unsigned int));
memcpy(clip.ipar, ct->ipar, COLORTABLE_NBMAX_PARAM * sizeof(int));
memcpy(clip.fpar, ct->fpar, COLORTABLE_NBMAX_PARAM * sizeof(float));
}
void ColorTable_Paste(ColorTable *ct){
ct->size = clip.size;
memcpy(ct->table, clip.table, clip.size * sizeof(unsigned int));
memcpy(ct->ipar, clip.ipar, COLORTABLE_NBMAX_PARAM * sizeof(int));
memcpy(ct->fpar, clip.fpar, COLORTABLE_NBMAX_PARAM * sizeof(float));
}
void ColorTable_Print(ColorTable *ct, FILE *fp){
int i, r, g, b, a;
char tmp1[1024],tmp2[1024];
strcpy(tmp1, "");
for (i=0;i<ct->size;i++) {
r = UNPACK_RED( ct->table[i] );
g = UNPACK_GREEN( ct->table[i] );
b = UNPACK_BLUE( ct->table[i] );
a = UNPACK_ALPHA( ct->table[i] );
if(i && !(i%4)){
if(fp) fprintf(fp, "%s\n", tmp1); else Msg(DIRECT, tmp1);
strcpy(tmp1, "");
}
sprintf(tmp2, "{%d, %d, %d, %d}", r, g, b, a);
strcat(tmp1, tmp2);
if(i!=ct->size-1) strcat(tmp1, ", ");
}
if(fp) fprintf(fp, "%s\n", tmp1); else Msg(DIRECT, tmp1);
}
#ifndef _COLORTABLE_H_
#define _COLORTABLE_H_
#define COLORTABLE_NBMAX_PARAM 10
#define COLORTABLE_NBMAX_COLOR 255
typedef struct{
unsigned int table[COLORTABLE_NBMAX_COLOR];
int size;
int ipar[COLORTABLE_NBMAX_PARAM];
float fpar[COLORTABLE_NBMAX_PARAM];
}ColorTable;
/* COLORTABLE_MODE */
#define COLORTABLE_RGB 1
#define COLORTABLE_HSV 2
/* integrer parameters indices */
#define COLORTABLE_NUMBER 0 /* predefined curve index */
#define COLORTABLE_CHANGED 1 /* did the colortable change ? */
#define COLORTABLE_INVERT 2 /* invert (rbg<->255-rgb) */
#define COLORTABLE_SWAP 3 /* swap (min<->max) */
#define COLORTABLE_ROTATE 4 /* rotate */
#define COLORTABLE_MODE 5 /* mode (rgb, hsv) */
/* float parameters indices */
#define COLORTABLE_CURVE 0 /* curvature */
#define COLORTABLE_BIAS 1 /* offset */
#define COLORTABLE_ALPHAPOW 2 /* alpha channel power */
#define COLORTABLE_ALPHAVAL 3 /* alpha channel value */
#define COLORTABLE_BETA 4 /* beta coeff for brighten */
void ColorTable_InitParam (int number, ColorTable * ct, int rgb_flag, int alpha_flag);
void ColorTable_Recompute (ColorTable * ct, int rgb_flag, int alpha_flag);
void ColorTable_Copy(ColorTable *ct);
void ColorTable_Paste(ColorTable *ct);
void ColorTable_Print(ColorTable *ct, FILE *fp) ;
#endif
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment