diff --git a/Fltk/aboutWindow.cpp b/Fltk/aboutWindow.cpp index 0ccda22976f90be434023c7a586adeb0b837e722..94abcdb16cc120dfb8c0d68095824088cc2579b7 100644 --- a/Fltk/aboutWindow.cpp +++ b/Fltk/aboutWindow.cpp @@ -76,7 +76,7 @@ aboutWindow::aboutWindow() std::vector<std::string> lines = SplitWhiteSpace(Get_GmshBuildOptions(), 30); for(unsigned int i = 0; i < lines.size(); i++){ if(!i) - sprintf(buffer, "@c@.Build options: %s", lines[i].c_str()); + sprintf(buffer, "@c@.Build options:%s", lines[i].c_str()); else sprintf(buffer, "@c@.%s", lines[i].c_str()); o->add(buffer); diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index 08549a695a5dd96b9ef71d770c08990206f6970d..f0454bf757c619827fcd292588a03d9669730824 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -12,11 +12,11 @@ #include "MElement.h" #include "VertexArray.h" #include "GmshMatrix.h" +#include "Numeric.h" #if defined(HAVE_GMSH_EMBEDDED) #include "GmshEmbedded.h" #else -#include "Numeric.h" #include "GaussLegendre1D.h" #include "Context.h" #endif diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp index c0aa06b51c4e9ba0815c7af406329b3069a1e4fd..d1c673727979984b3ecab2feeed9025c79ec924e 100644 --- a/Geo/MElement.cpp +++ b/Geo/MElement.cpp @@ -11,11 +11,11 @@ #include "GEntity.h" #include "GFace.h" #include "StringUtils.h" +#include "Numeric.h" #if defined(HAVE_GMSH_EMBEDDED) #include "GmshEmbedded.h" #else -#include "Numeric.h" #include "GaussLegendre1D.h" #include "Context.h" #include "qualityMeasures.h" diff --git a/Geo/MFace.cpp b/Geo/MFace.cpp index 07c852b572d3273c8a5c4a9dc8bc36c454af6387..4ef4dbfa13d3ac13706c1d36a9e07f5c223a07c6 100644 --- a/Geo/MFace.cpp +++ b/Geo/MFace.cpp @@ -5,11 +5,11 @@ #include "GmshConfig.h" #include "MFace.h" +#include "Numeric.h" #if defined(HAVE_GMSH_EMBEDDED) #include "GmshEmbedded.h" #else -#include "Numeric.h" #include "Context.h" #endif diff --git a/Numeric/NumericEmbedded.cpp b/Numeric/NumericEmbedded.cpp deleted file mode 100644 index 5ffe3cdd112af8319bdcf1f57321d83bf546ceca..0000000000000000000000000000000000000000 --- a/Numeric/NumericEmbedded.cpp +++ /dev/null @@ -1,524 +0,0 @@ -// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle -// -// See the LICENSE.txt file for license information. 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 or external functions, expect -// Msg) - -#include "NumericEmbedded.h" -#include "GmshMessage.h" - -#define SQU(a) ((a)*(a)) - -double myatan2(double a, double b) -{ - if(a == 0.0 && b == 0) - return 0.0; - return atan2(a, b); -} - -double myasin(double a) -{ - if(a <= -1.) - return -M_PI / 2.; - else if(a >= 1.) - return M_PI / 2.; - else - return asin(a); -} - -double myacos(double a) -{ - if(a <= -1.) - return M_PI; - else if(a >= 1.) - return 0.; - else - return acos(a); -} - -void matvec(double mat[3][3], double vec[3], double res[3]) -{ - res[0] = mat[0][0] * vec[0] + mat[0][1] * vec[1] + mat[0][2] * vec[2]; - res[1] = mat[1][0] * vec[0] + mat[1][1] * vec[1] + mat[1][2] * vec[2]; - res[2] = mat[2][0] * vec[0] + mat[2][1] * vec[1] + mat[2][2] * vec[2]; -} - -void normal3points(double x0, double y0, double z0, - double x1, double y1, double z1, - double x2, double y2, double z2, - double n[3]) -{ - double t1[3] = {x1 - x0, y1 - y0, z1 - z0}; - double t2[3] = {x2 - x0, y2 - y0, z2 - z0}; - prodve(t1, t2, n); - norme(n); -} - -int sys2x2(double mat[2][2], double b[2], double res[2]) -{ - double det, ud, norm; - int i; - - norm = SQU(mat[0][0]) + SQU(mat[1][1]) + SQU(mat[0][1]) + SQU(mat[1][0]); - det = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; - - // TOLERANCE ! WARNING WARNING - if(norm == 0.0 || fabs(det) / norm < 1.e-12) { - if(norm) - Msg::Debug("Assuming 2x2 matrix is singular (det/norm == %.16g)", - fabs(det) / norm); - res[0] = res[1] = 0.0; - return 0; - } - ud = 1. / det; - - res[0] = b[0] * mat[1][1] - mat[0][1] * b[1]; - res[1] = mat[0][0] * b[1] - mat[1][0] * b[0]; - - for(i = 0; i < 2; i++) - res[i] *= ud; - - return (1); -} - -double det3x3(double mat[3][3]) -{ - return (mat[0][0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - - mat[0][1] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) + - mat[0][2] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0])); -} - -double trace3x3(double mat[3][3]) -{ - return mat[0][0] + mat[1][1] + mat[2][2]; -} - -double trace2 (double mat[3][3]) -{ - double a00 = mat[0][0] * mat[0][0] + mat[1][0] * mat[0][1] + mat[2][0] * mat[0][2]; - double a11 = mat[1][0] * mat[0][1] + mat[1][1] * mat[1][1] + mat[1][2] * mat[2][1]; - double a22 = mat[2][0] * mat[0][2] + mat[2][1] * mat[1][2] + mat[2][2] * mat[2][2]; - - return a00 + a11 + a22; -} - -int sys3x3(double mat[3][3], double b[3], double res[3], double *det) -{ - double ud; - int i; - - *det = det3x3(mat); - - if(*det == 0.0) { - res[0] = res[1] = res[2] = 0.0; - return (0); - } - - ud = 1. / (*det); - - res[0] = b[0] * (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) - - mat[0][1] * (b[1] * mat[2][2] - mat[1][2] * b[2]) + - mat[0][2] * (b[1] * mat[2][1] - mat[1][1] * b[2]); - - res[1] = mat[0][0] * (b[1] * mat[2][2] - mat[1][2] * b[2]) - - b[0] * (mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) + - mat[0][2] * (mat[1][0] * b[2] - b[1] * mat[2][0]); - - res[2] = mat[0][0] * (mat[1][1] * b[2] - b[1] * mat[2][1]) - - mat[0][1] * (mat[1][0] * b[2] - b[1] * mat[2][0]) + - b[0] * (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]); - - for(i = 0; i < 3; i++) - res[i] *= ud; - return (1); -} - -int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det) -{ - int out; - double norm; - - out = sys3x3(mat, b, res, det); - norm = - SQU(mat[0][0]) + SQU(mat[0][1]) + SQU(mat[0][2]) + - SQU(mat[1][0]) + SQU(mat[1][1]) + SQU(mat[1][2]) + - SQU(mat[2][0]) + SQU(mat[2][1]) + SQU(mat[2][2]); - - // TOLERANCE ! WARNING WARNING - if(norm == 0.0 || fabs(*det) / norm < 1.e-12) { - if(norm) - Msg::Debug("Assuming 3x3 matrix is singular (det/norm == %.16g)", - fabs(*det) / norm); - res[0] = res[1] = res[2] = 0.0; - return 0; - } - - return out; -} - -double det2x2(double mat[2][2]) -{ - return mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]; -} - -double det2x3(double mat[2][3]) -{ - double v1[3] = {mat[0][0], mat[0][1], mat[0][2]}; - double v2[3] = {mat[1][0], mat[1][1], mat[1][2]}; - double n[3]; - - prodve(v1, v2, n); - return norm3(n); -} - -double inv2x2(double mat[2][2], double inv[2][2]) -{ - const double det = det2x2(mat); - if(det){ - double ud = 1. / det; - inv[0][0] = mat[1][1] * ud; - inv[1][0] = -mat[1][0] * ud; - inv[0][1] = -mat[0][1] * ud; - inv[1][1] = mat[0][0] * ud; - } - else{ - Msg::Error("Singular matrix"); - for(int i = 0; i < 2; i++) - for(int j = 0; j < 2; j++) - inv[i][j] = 0.; - } - return det; -} - -double inv3x3(double mat[3][3], double inv[3][3]) -{ - double det = det3x3(mat); - if(det){ - double ud = 1. / det; - inv[0][0] = (mat[1][1] * mat[2][2] - mat[1][2] * mat[2][1]) * ud; - inv[1][0] = -(mat[1][0] * mat[2][2] - mat[1][2] * mat[2][0]) * ud; - inv[2][0] = (mat[1][0] * mat[2][1] - mat[1][1] * mat[2][0]) * ud; - inv[0][1] = -(mat[0][1] * mat[2][2] - mat[0][2] * mat[2][1]) * ud; - inv[1][1] = (mat[0][0] * mat[2][2] - mat[0][2] * mat[2][0]) * ud; - inv[2][1] = -(mat[0][0] * mat[2][1] - mat[0][1] * mat[2][0]) * ud; - inv[0][2] = (mat[0][1] * mat[1][2] - mat[0][2] * mat[1][1]) * ud; - inv[1][2] = -(mat[0][0] * mat[1][2] - mat[0][2] * mat[1][0]) * ud; - inv[2][2] = (mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) * ud; - } - else{ - Msg::Error("Singular matrix"); - for(int i = 0; i < 3; i++) - for(int j = 0; j < 3; j++) - inv[i][j] = 0.; - } - return det; -} - -double angle_02pi(double A3) -{ - double DP = 2 * M_PI; - while(A3 > DP || A3 < 0.) { - if(A3 > 0) - A3 -= DP; - else - A3 += DP; - } - return A3; -} - -double angle_plan(double V[3], double P1[3], double P2[3], double n[3]) -{ - double PA[3], PB[3], angplan; - double cosc, sinc, c[3]; - - PA[0] = P1[0] - V[0]; - PA[1] = P1[1] - V[1]; - PA[2] = P1[2] - V[2]; - - PB[0] = P2[0] - V[0]; - PB[1] = P2[1] - V[1]; - PB[2] = P2[2] - V[2]; - - norme(PA); - norme(PB); - - prodve(PA, PB, c); - - prosca(PA, PB, &cosc); - prosca(c, n, &sinc); - angplan = myatan2(sinc, cosc); - - return angplan; -} - -double triangle_area(double p0[3], double p1[3], double p2[3]) -{ - double a[3], b[3], c[3]; - - a[0] = p2[0] - p1[0]; - a[1] = p2[1] - p1[1]; - a[2] = p2[2] - p1[2]; - - b[0] = p0[0] - p1[0]; - b[1] = p0[1] - p1[1]; - b[2] = p0[2] - p1[2]; - - prodve(a, b, c); - return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]); -} - -char float2char(float f) -{ - // float normalized in [-1, 1], char in [-127, 127] - f *= 127.; - if(f > 127.) return 127; - else if(f < -127.) return -127; - else return (char)f; -} - -float char2float(char c) -{ - float f = c; - f /= 127.; - if(f > 1.) return 1.; - else if(f < -1) return -1.; - else return f; -} - -double InterpolateIso(double *X, double *Y, double *Z, - double *Val, double V, int I1, int I2, - double *XI, double *YI, double *ZI) -{ - if(Val[I1] == Val[I2]) { - *XI = X[I1]; - *YI = Y[I1]; - *ZI = Z[I1]; - return 0; - } - else { - double coef = (V - Val[I1]) / (Val[I2] - Val[I1]); - *XI = coef * (X[I2] - X[I1]) + X[I1]; - *YI = coef * (Y[I2] - Y[I1]) + Y[I1]; - *ZI = coef * (Z[I2] - Z[I1]) + Z[I1]; - return coef; - } -} - -void gradSimplex(double *x, double *y, double *z, double *v, double *grad) -{ - // p = p1 * (1-u-v-w) + p2 u + p3 v + p4 w - - double mat[3][3]; - double det, b[3]; - mat[0][0] = x[1] - x[0]; - mat[1][0] = x[2] - x[0]; - mat[2][0] = x[3] - x[0]; - mat[0][1] = y[1] - y[0]; - mat[1][1] = y[2] - y[0]; - mat[2][1] = y[3] - y[0]; - mat[0][2] = z[1] - z[0]; - mat[1][2] = z[2] - z[0]; - mat[2][2] = z[3] - z[0]; - b[0] = v[1] - v[0]; - b[1] = v[2] - v[0]; - b[2] = v[3] - v[0]; - sys3x3(mat, b, grad, &det); -} - -double ComputeVonMises(double *V) -{ - double tr = (V[0] + V[4] + V[8]) / 3.; - double v11 = V[0] - tr, v12 = V[1], v13 = V[2]; - double v21 = V[3], v22 = V[4] - tr, v23 = V[5]; - double v31 = V[6], v32 = V[7], v33 = V[8] - tr; - return sqrt(1.5 * (v11 * v11 + v12 * v12 + v13 * v13 + - v21 * v21 + v22 * v22 + v23 * v23 + - v31 * v31 + v32 * v32 + v33 * v33)); -} - -double ComputeScalarRep(int numComp, double *val) -{ - if(numComp == 1) - return val[0]; - else if(numComp == 3) - return sqrt(val[0] * val[0] + val[1] * val[1] + val[2] * val[2]); - else if(numComp == 9) - return ComputeVonMises(val); - return 0.; -} - -void eigenvalue(double mat[3][3], double v[3]) -{ - // characteristic polynomial of T : find v root of - // v^3 - I1 v^2 + I2 T - I3 = 0 - // I1 : first invariant , trace(T) - // I2 : second invariant , 1/2 (I1^2 -trace(T^2)) - // I3 : third invariant , det T - - double c[4]; - c[3] = 1.0; - c[2] = - trace3x3(mat); - c[1] = 0.5 * (c[2] * c[2] - trace2(mat)); - c[0] = - det3x3(mat); - - //printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]); - //printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]); - //printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]); - //printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]); - - double imag[3]; - FindCubicRoots(c, v, imag); - eigsort(v); -} - -void FindCubicRoots(const double coef[4], double real[3], double imag[3]) -{ - double a = coef[3]; - double b = coef[2]; - double c = coef[1]; - double d = coef[0]; - - if(!a || !d){ - Msg::Error("Degenerate cubic: use a second degree solver!"); - return; - } - - b /= a; - c /= a; - d /= a; - - double q = (3.0*c - (b*b))/9.0; - double r = -(27.0*d) + b*(9.0*c - 2.0*(b*b)); - r /= 54.0; - - double discrim = q*q*q + r*r; - imag[0] = 0.0; // The first root is always real. - double term1 = (b/3.0); - - if (discrim > 0) { // one root is real, two are complex - double s = r + sqrt(discrim); - s = ((s < 0) ? -pow(-s, (1.0/3.0)) : pow(s, (1.0/3.0))); - double t = r - sqrt(discrim); - t = ((t < 0) ? -pow(-t, (1.0/3.0)) : pow(t, (1.0/3.0))); - real[0] = -term1 + s + t; - term1 += (s + t)/2.0; - real[1] = real[2] = -term1; - term1 = sqrt(3.0)*(-t + s)/2; - imag[1] = term1; - imag[2] = -term1; - return; - } - - // The remaining options are all real - imag[1] = imag[2] = 0.0; - - double r13; - if (discrim == 0){ // All roots real, at least two are equal. - r13 = ((r < 0) ? -pow(-r,(1.0/3.0)) : pow(r,(1.0/3.0))); - real[0] = -term1 + 2.0*r13; - real[1] = real[2] = -(r13 + term1); - return; - } - - // Only option left is that all roots are real and unequal (to get - // here, q < 0) - q = -q; - double dum1 = q*q*q; - dum1 = acos(r/sqrt(dum1)); - r13 = 2.0*sqrt(q); - real[0] = -term1 + r13*cos(dum1/3.0); - real[1] = -term1 + r13*cos((dum1 + 2.0*M_PI)/3.0); - real[2] = -term1 + r13*cos((dum1 + 4.0*M_PI)/3.0); -} - -void eigsort(double d[3]) -{ - int k, j, i; - double p; - - for (i=0; i<3; i++) { - p=d[k=i]; - for (j=i+1;j<3;j++) - if (d[j] >= p) p=d[k=j]; - if (k != i) { - d[k]=d[i]; - d[i]=p; - } - } -} - -void invert_singular_matrix3x3(double MM[3][3], double II[3][3]) -{ - int i, j, k, n = 3; - double TT[3][3]; - - for(i = 1; i <= n; i++) { - for(j = 1; j <= n; j++) { - II[i - 1][j - 1] = 0.0; - TT[i - 1][j - 1] = 0.0; - } - } - - Double_Matrix M(3, 3), V(3, 3); - Double_Vector W(3); - for(i = 1; i <= n; i++){ - for(j = 1; j <= n; j++){ - M(i - 1, j - 1) = MM[i - 1][j - 1]; - } - } - M.svd(V, W); - for(i = 1; i <= n; i++) { - for(j = 1; j <= n; j++) { - double ww = W(i - 1); - if(fabs(ww) > 1.e-16) { // singular value precision - TT[i - 1][j - 1] += M(j - 1, i - 1) / ww; - } - } - } - for(i = 1; i <= n; i++) { - for(j = 1; j <= n; j++) { - for(k = 1; k <= n; k++) { - II[i - 1][j - 1] += V(i - 1, k - 1) * TT[k - 1][j - 1]; - } - } - } -} - -bool newton_fd(void (*func)(Double_Vector &, Double_Vector &, void *), - Double_Vector &x, void *data, double relax, double tolx) -{ - const int MAXIT = 50; - const double EPS = 1.e-4; - const int N = x.size(); - - Double_Matrix J(N, N); - Double_Vector f(N), feps(N), dx(N); - - for (int iter = 0; iter < MAXIT; iter++){ - func(x, f, data); - - for (int j = 0; j < N; j++){ - double h = EPS * fabs(x(j)); - if(h == 0.) h = EPS; - x(j) += h; - func(x, feps, data); - for (int i = 0; i < N; i++) - J(i, j) = (feps(i) - f(i)) / h; - x(j) -= h; - } - - if (N == 1) - dx(0) = f(0) / J(0, 0); - else - J.lu_solve(f, dx); - - for (int i = 0; i < N; i++) - x(i) -= relax * dx(i); - - if(dx.norm() < tolx) return true; - } - return false; -} diff --git a/Numeric/NumericEmbedded.h b/Numeric/NumericEmbedded.h deleted file mode 100644 index f6b0368ae12bba405f3bdeefe036415262ab29f7..0000000000000000000000000000000000000000 --- a/Numeric/NumericEmbedded.h +++ /dev/null @@ -1,76 +0,0 @@ -// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle -// -// See the LICENSE.txt file for license information. Please report all -// bugs and problems to <gmsh@geuz.org>. - -#ifndef _NUMERIC_EMBEDDED_H_ -#define _NUMERIC_EMBEDDED_H_ - -#include <math.h> -#include "GmshMatrix.h" - -#define myhypot(a,b) (sqrt((a)*(a)+(b)*(b))) -#define sign(x) (((x)>=0)?1:-1) - -double myatan2(double a, double b); -double myasin(double a); -double myacos(double a); -inline void prodve(double a[3], double b[3], double c[3]) -{ - c[2] = a[0] * b[1] - a[1] * b[0]; - c[1] = -a[0] * b[2] + a[2] * b[0]; - c[0] = a[1] * b[2] - a[2] * b[1]; -} -inline void prosca(double a[3], double b[3], double *c) -{ - *c = a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; -} -void matvec(double mat[3][3], double vec[3], double res[3]); -inline double norm3(double a[3]) -{ - return sqrt(a[0] * a[0] + a[1] * a[1] + a[2] * a[2]); -} -inline double norme(double a[3]) -{ - const double mod = norm3(a); - if(mod != 0.0){ - const double one_over_mod = 1./mod; - a[0] *= one_over_mod; - a[1] *= one_over_mod; - a[2] *= one_over_mod; - } - return mod; -} -void normal3points(double x0, double y0, double z0, - double x1, double y1, double z1, - double x2, double y2, double z2, - double n[3]); -int sys2x2(double mat[2][2], double b[2], double res[2]); -int sys3x3(double mat[3][3], double b[3], double res[3], double *det); -int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det); -double det2x2(double mat[2][2]); -double det2x3(double mat[2][3]); -double det3x3(double mat[3][3]); -double trace3x3(double mat[3][3]); -double trace2 (double mat[3][3]); -double inv3x3(double mat[3][3], double inv[3][3]); -double inv2x2(double mat[2][2], double inv[2][2]); -double angle_02pi(double A3); -double angle_plan(double V[3], double P1[3], double P2[3], double n[3]); -double triangle_area(double p0[3], double p1[3], double p2[3]); -char float2char(float f); -float char2float(char c); -void eigenvalue(double mat[3][3], double re[3]); -void FindCubicRoots(const double coeff[4], double re[3], double im[3]); -void eigsort(double d[3]); -double InterpolateIso(double *X, double *Y, double *Z, - double *Val, double V, int I1, int I2, - double *XI, double *YI ,double *ZI); -void gradSimplex(double *x, double *y, double *z, double *v, double *grad); -double ComputeVonMises(double *val); -double ComputeScalarRep(int numComp, double *val); -void invert_singular_matrix3x3(double MM[3][3], double II[3][3]); -bool newton_fd(void (*func)(Double_Vector &, Double_Vector &, void *), - Double_Vector &x, void *data, double relax=1., double tolx=1.e-6); - -#endif diff --git a/utils/embed/GmshEmbedded.h b/utils/embed/GmshEmbedded.h index 8795aabd45573b14a9ca411e50aac75ad1b97d7c..779904bb7522644174aa5313e60774439a6c537c 100644 --- a/utils/embed/GmshEmbedded.h +++ b/utils/embed/GmshEmbedded.h @@ -6,8 +6,6 @@ #ifndef _GMSH_EMBEDDED_H_ #define _GMSH_EMBEDDED_H_ -#include "NumericEmbedded.h" - class CTX{ private: static CTX *_instance; diff --git a/utils/embed/Makefile b/utils/embed/Makefile index 684b115c13ed12c661736773820697f0de0a84cb..59b5a50841054720c960a2ad7d78c2e354854d10 100644 --- a/utils/embed/Makefile +++ b/utils/embed/Makefile @@ -18,7 +18,7 @@ SRC = GModel.cpp\ discreteEdge.cpp discreteFace.cpp discreteRegion.cpp\ MElement.cpp\ MFace.cpp MVertex.cpp\ - NumericEmbedded.cpp\ + Numeric.cpp\ FunctionSpace.cpp\ StringUtils.cpp\ GmshEmbedded.cpp