diff --git a/Numeric/NumericEmbedded.cpp b/Numeric/NumericEmbedded.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ffe3cdd112af8319bdcf1f57321d83bf546ceca
--- /dev/null
+++ b/Numeric/NumericEmbedded.cpp
@@ -0,0 +1,524 @@
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..f6b0368ae12bba405f3bdeefe036415262ab29f7
--- /dev/null
+++ b/Numeric/NumericEmbedded.h
@@ -0,0 +1,76 @@
+// 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/Numeric/gmshAssembler.cpp b/Numeric/gmshAssembler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f28f2ec96b450c1a6bbffe4bdf8bd3d3fc7ff305
--- /dev/null
+++ b/Numeric/gmshAssembler.cpp
@@ -0,0 +1,68 @@
+// 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>.
+
+#include "MVertex.h"
+#include "gmshAssembler.h"
+#include "gmshLinearSystem.h"
+
+void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
+ MVertex *vC, int iCompC, int iFieldC,
+ double val)
+{
+ if (!lsys->isAllocated()){
+ lsys->allocate(numbering.size());
+ }
+
+ std::map<gmshDofKey, int>::iterator itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itR != numbering.end()){
+ std::map<gmshDofKey, int>::iterator itC = numbering.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itC != numbering.end()){
+ lsys->addToMatrix(itR->second, itC->second, val);
+ }
+ else {
+ std::map<gmshDofKey, double>::iterator itF = fixed.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itF != fixed.end()){
+ lsys->addToRightHandSide(itR->second, -val*itF->second);
+ }
+ else{
+ std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator itConstrC =
+ constraints.find(gmshDofKey(vC, iCompC, iFieldC));
+ if (itConstrC != constraints.end()){
+ for (unsigned int i = 0; i < itConstrC->second.size(); i++){
+ gmshDofKey &dofKeyConstrC = itConstrC->second[i].first;
+ double valConstrC = itConstrC->second[i].second;
+ assemble(vR, iCompR, iFieldR,
+ dofKeyConstrC.v,dofKeyConstrC.comp, dofKeyConstrC.field,
+ val * valConstrC);
+ }
+ }
+ }
+ }
+ }
+ else{
+ std::map<gmshDofKey, std::vector<std::pair<gmshDofKey, double> > >::iterator itConstrR =
+ constraints.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itConstrR != constraints.end()){
+ for (unsigned int i = 0; i < itConstrR->second.size(); i++){
+ gmshDofKey &dofKeyConstrR = itConstrR->second[i].first;
+ double valConstrR = itConstrR->second[i].second;
+ assemble(dofKeyConstrR.v,dofKeyConstrR.comp, dofKeyConstrR.field,
+ vC, iCompC, iFieldC,
+ val * valConstrR);
+ }
+ }
+ }
+}
+
+void gmshAssembler::assemble(MVertex *vR, int iCompR, int iFieldR,
+ double val)
+{
+ if (!lsys->isAllocated())lsys->allocate(numbering.size());
+ std::map<gmshDofKey, int>::iterator itR = numbering.find(gmshDofKey(vR, iCompR, iFieldR));
+ if (itR != numbering.end()){
+ lsys->addToRightHandSide(itR->second, val);
+ }
+}
+
diff --git a/Numeric/gmshTermOfFormulation.cpp b/Numeric/gmshTermOfFormulation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..088433923f1e1c4af6a6aa6dd0ea8efecd0be15e
--- /dev/null
+++ b/Numeric/gmshTermOfFormulation.cpp
@@ -0,0 +1,111 @@
+// 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>.
+
+#include "Gmsh.h"
+#include "GModel.h"
+#include "MElement.h"
+#include "GmshMatrix.h"
+#include "gmshTermOfFormulation.h"
+#include "gmshFunction.h"
+#include "gmshLinearSystem.h"
+#include "gmshAssembler.h"
+
+gmshNodalFemTerm::~gmshNodalFemTerm ()
+{
+}
+
+void gmshNodalFemTerm::addDirichlet(int physical,
+ int dim,
+ int comp,
+ int field,
+ const gmshFunction & e,
+ gmshAssembler &lsys)
+{
+}
+
+void gmshNodalFemTerm::addNeumann(int physical,
+ int dim,
+ int comp,
+ int field,
+ const gmshFunction & fct,
+ gmshAssembler &lsys)
+{
+}
+
+void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys) const
+{
+ if (_gm->getNumRegions()){
+ for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
+ addToMatrix(lsys, *it);
+ }
+ }
+ else if(_gm->getNumFaces()){
+ for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
+ addToMatrix(lsys, *it);
+ }
+ }
+}
+
+void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,const std::vector<MElement*> &v) const
+{
+ for (unsigned int i = 0; i < v.size(); i++)
+ addToMatrix(lsys, v[i]);
+}
+
+void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys,
+ Double_Matrix &localMatrix,
+ MElement *e) const
+{
+ const int nbR = sizeOfR(e);
+ const int nbC = sizeOfC(e);
+ for (int j = 0; j < nbR; j++){
+ MVertex *vR;int iCompR,iFieldR;
+ getLocalDofR (e, j, &vR, &iCompR, &iFieldR);
+ for (int k = 0; k < nbC; k++){
+ MVertex *vC;
+ int iCompC, iFieldC;
+ getLocalDofC(e, k, &vC, &iCompC, &iFieldC);
+ lsys.assemble(vR, iCompR, iFieldR,
+ vC, iCompC, iFieldC,
+ localMatrix(j, k));
+ }
+ }
+}
+
+void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, MElement *e) const
+{
+ const int nbR = sizeOfR(e);
+ const int nbC = sizeOfC(e);
+ Double_Matrix localMatrix (nbR, nbC);
+ elementMatrix(e, localMatrix);
+ addToMatrix(lsys, localMatrix, e);
+}
+
+void gmshNodalFemTerm::addToMatrix(gmshAssembler &lsys, GEntity *ge) const
+{
+ for(unsigned int i = 0; i < ge->getNumMeshElements(); i++){
+ MElement *e = ge->getMeshElement(i);
+ addToMatrix(lsys, e);
+ }
+}
+
+void gmshNodalFemTerm::addToRightHandSide(gmshAssembler &lsys) const
+{
+ if (_gm->getNumRegions()){
+ for(GModel::riter it = _gm->firstRegion(); it != _gm->lastRegion(); ++it){
+ addToRightHandSide(lsys, *it);
+ }
+ }
+ else if(_gm->getNumFaces()){
+ for(GModel::fiter it = _gm->firstFace(); it != _gm->lastFace(); ++it){
+ addToRightHandSide(lsys, *it);
+ }
+ }
+}
+
+void gmshNodalFemTerm::addToRightHandSide(gmshAssembler &lsys, GEntity *ge) const
+{
+ throw;
+}
diff --git a/contrib/Chaco/misc/timing.c b/contrib/Chaco/misc/timing.c
new file mode 100644
index 0000000000000000000000000000000000000000..70e9490e1a593a6e834c5d14353bad9247b7a6ad
--- /dev/null
+++ b/contrib/Chaco/misc/timing.c
@@ -0,0 +1,237 @@
+/* This software was developed by Bruce Hendrickson and Robert Leland *
+ * at Sandia National Laboratories under US Department of Energy *
+ * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
+
+#include <stdio.h>
+
+/* Timing parameters. */
+
+double start_time = -1;
+double total_time = 0;
+double input_time = 0;
+double partition_time = 0;
+double sequence_time = 0;
+double kernel_time = 0;
+double reformat_time = 0;
+double check_input_time = 0;
+double count_time = 0;
+double print_assign_time = 0;
+
+double coarsen_time = 0;
+double match_time = 0;
+double make_cgraph_time = 0;
+
+double lanczos_time = 0;
+double splarax_time = 0;
+double orthog_time = 0;
+double ql_time = 0;
+double tevec_time = 0;
+double ritz_time = 0;
+double evec_time = 0;
+double blas_time = 0;
+double init_time = 0;
+double scan_time = 0;
+double debug_time = 0;
+double pause_time = 0;
+
+double rqi_symmlq_time = 0;
+double refine_time = 0;
+
+double kl_total_time = 0;
+double kl_init_time = 0;
+double nway_kl_time = 0;
+double kl_bucket_time = 0;
+
+double inertial_time = 0;
+double inertial_axis_time = 0;
+double median_time = 0;
+
+double sim_time = 0;
+
+
+void time_out(outfile)
+FILE *outfile; /* file to print output to */
+{
+ FILE *tempfile; /* file name for two passes */
+ extern int ECHO; /* parameter for different output styles */
+ extern int OUTPUT_TIME; /* how much timing output should I print? */
+ double time_tol; /* tolerance for ignoring time */
+ double other_time; /* time not accounted for */
+ int i; /* loop counter */
+
+ time_tol = 0.005;
+
+ for (i = 0; i < 2; i++) {
+ if (i == 1) { /* Print to file? */
+ if (ECHO < 0)
+ tempfile = outfile;
+ else
+ break;
+ }
+ else { /* Print to stdout. */
+ tempfile = stdout;
+ }
+
+
+ if (OUTPUT_TIME > 0) {
+ if (total_time != 0) {
+ fprintf(tempfile, "\nTotal time: %g sec.\n", total_time);
+ if (input_time != 0)
+ fprintf(tempfile, " input %g\n", input_time);
+ if (reformat_time != 0)
+ fprintf(tempfile, " reformatting %g\n", reformat_time);
+ if (check_input_time != 0)
+ fprintf(tempfile, " checking input %g\n", check_input_time);
+ if (partition_time != 0)
+ fprintf(tempfile, " partitioning %g\n", partition_time);
+ if (sequence_time != 0)
+ fprintf(tempfile, " sequencing %g\n", sequence_time);
+ if (kernel_time != 0)
+ fprintf(tempfile, " kernel benchmarking %g\n", kernel_time);
+ if (count_time != 0)
+ fprintf(tempfile, " evaluation %g\n", count_time);
+ if (print_assign_time != 0)
+ fprintf(tempfile, " printing assignment file %g\n", print_assign_time);
+
+ if (sim_time != 0)
+ fprintf(tempfile, " simulating %g\n", sim_time);
+ other_time = total_time - input_time - reformat_time -
+ check_input_time - partition_time - count_time -
+ print_assign_time - sim_time - sequence_time - kernel_time;
+ if (other_time > time_tol)
+ fprintf(tempfile, " other %g\n", other_time);
+ }
+ }
+
+ if (OUTPUT_TIME > 1) {
+ if (inertial_time != 0) {
+ if (inertial_time != 0)
+ fprintf(tempfile, "\nInertial time: %g sec.\n", inertial_time);
+ if (inertial_axis_time != 0)
+ fprintf(tempfile, " inertial axis %g\n", inertial_axis_time);
+ if (median_time != 0)
+ fprintf(tempfile, " median finding %g\n", median_time);
+ other_time = inertial_time - inertial_axis_time - median_time;
+ if (other_time > time_tol)
+ fprintf(tempfile, " other %g\n", other_time);
+ }
+
+ if (kl_total_time != 0) {
+ if (kl_total_time != 0)
+ fprintf(tempfile, "\nKL time: %g sec.\n", kl_total_time);
+ if (kl_init_time != 0)
+ fprintf(tempfile, " initialization %g\n", kl_init_time);
+ if (nway_kl_time != 0)
+ fprintf(tempfile, " nway refinement %g\n", nway_kl_time);
+ if (kl_bucket_time != 0)
+ fprintf(tempfile, " bucket sorting %g\n", kl_bucket_time);
+ other_time = kl_total_time - kl_init_time - nway_kl_time;
+ if (other_time > time_tol)
+ fprintf(tempfile, " other %g\n", other_time);
+ }
+
+ if (coarsen_time != 0 && rqi_symmlq_time == 0) {
+ fprintf(tempfile, "\nCoarsening %g sec.\n", coarsen_time);
+ if (match_time != 0)
+ fprintf(tempfile, " maxmatch %g\n", match_time);
+ if (make_cgraph_time != 0)
+ fprintf(tempfile, " makecgraph %g\n", make_cgraph_time);
+ }
+
+ if (lanczos_time != 0) {
+ fprintf(tempfile, "\nLanczos time: %g sec.\n", lanczos_time);
+ if (splarax_time != 0)
+ fprintf(tempfile, " matvec/solve %g\n", splarax_time);
+ if (blas_time != 0)
+ fprintf(tempfile, " vector ops %g\n", blas_time);
+ if (evec_time != 0)
+ fprintf(tempfile, " assemble evec %g\n", evec_time);
+ if (init_time != 0)
+ fprintf(tempfile, " malloc/init/free %g\n", init_time);
+ if (orthog_time != 0)
+ fprintf(tempfile, " maintain orthog %g\n", orthog_time);
+ if (scan_time != 0)
+ fprintf(tempfile, " scan %g\n", scan_time);
+ if (debug_time != 0)
+ fprintf(tempfile, " debug/warning/check %g\n", debug_time);
+ if (ql_time != 0)
+ fprintf(tempfile, " ql/bisection %g\n", ql_time);
+ if (tevec_time != 0)
+ fprintf(tempfile, " tevec %g\n", tevec_time);
+ if (ritz_time != 0)
+ fprintf(tempfile, " compute ritz %g\n", ritz_time);
+ if (pause_time != 0)
+ fprintf(tempfile, " pause %g\n", pause_time);
+ other_time = lanczos_time - splarax_time - orthog_time
+ - ql_time - tevec_time - ritz_time - evec_time
+ - blas_time - init_time - scan_time - debug_time - pause_time;
+ if (other_time > time_tol && other_time != lanczos_time) {
+ fprintf(tempfile, " other %g\n", other_time);
+ }
+ }
+
+ if (rqi_symmlq_time != 0) {
+ fprintf(tempfile, "\nRQI/Symmlq time: %g sec.\n", rqi_symmlq_time);
+ if (coarsen_time != 0)
+ fprintf(tempfile, " coarsening %g\n", coarsen_time);
+ if (match_time != 0)
+ fprintf(tempfile, " maxmatch %g\n", match_time);
+ if (make_cgraph_time != 0)
+ fprintf(tempfile, " makecgraph %g\n", make_cgraph_time);
+ if (refine_time != 0)
+ fprintf(tempfile, " refinement %g\n", refine_time);
+ if (lanczos_time != 0)
+ fprintf(tempfile, " lanczos %g\n", lanczos_time);
+ other_time = rqi_symmlq_time - coarsen_time - refine_time - lanczos_time;
+ if (other_time > time_tol)
+ fprintf(tempfile, " other %g\n", other_time);
+ }
+ }
+ }
+}
+
+
+void clear_timing()
+{
+ start_time = -1;
+ total_time = 0;
+ input_time = 0;
+ partition_time = 0;
+ sequence_time = 0;
+ kernel_time = 0;
+ reformat_time = 0;
+ check_input_time = 0;
+ count_time = 0;
+ print_assign_time = 0;
+
+ coarsen_time = 0;
+ match_time = 0;
+ make_cgraph_time = 0;
+
+ lanczos_time = 0;
+ splarax_time = 0;
+ orthog_time = 0;
+ ql_time = 0;
+ tevec_time = 0;
+ ritz_time = 0;
+ evec_time = 0;
+ blas_time = 0;
+ init_time = 0;
+ scan_time = 0;
+ debug_time = 0;
+ pause_time = 0;
+
+ rqi_symmlq_time = 0;
+ refine_time = 0;
+
+ kl_total_time = 0;
+ kl_init_time = 0;
+ nway_kl_time = 0;
+ kl_bucket_time = 0;
+
+ inertial_time = 0;
+ inertial_axis_time = 0;
+ median_time = 0;
+
+ sim_time = 0;
+}
diff --git a/examples/cube.geo b/examples/cube.geo
deleted file mode 100644
index decdd4a696bb42dae6ae5bda42e91f3bb5b099e1..0000000000000000000000000000000000000000
--- a/examples/cube.geo
+++ /dev/null
@@ -1,83 +0,0 @@
-lcar1 = .0275;
-
-length = 0.25;
-height = 1.0;
-depth = 0.25;
-
-Point(newp) = {length/2,height/2,depth,lcar1}; /* Point 1 */
-Point(newp) = {length/2,height/2,0,lcar1}; /* Point 2 */
-Point(newp) = {-length/2,height/2,depth,lcar1}; /* Point 3 */
-Point(newp) = {-length/2,-height/2,depth,lcar1}; /* Point 4 */
-Point(newp) = {length/2,-height/2,depth,lcar1}; /* Point 5 */
-Point(newp) = {length/2,-height/2,0,lcar1}; /* Point 6 */
-Point(newp) = {-length/2,height/2,0,lcar1}; /* Point 7 */
-Point(newp) = {-length/2,-height/2,0,lcar1}; /* Point 8 */
-Line(1) = {3,1};
-Line(2) = {3,7};
-Line(3) = {7,2};
-Line(4) = {2,1};
-Line(5) = {1,5};
-Line(6) = {5,4};
-Line(7) = {4,8};
-Line(8) = {8,6};
-Line(9) = {6,5};
-Line(10) = {6,2};
-Line(11) = {3,4};
-Line(12) = {8,7};
-Line Loop(13) = {-6,-5,-1,11};
-Plane Surface(14) = {13};
-Line Loop(15) = {4,5,-9,10};
-Plane Surface(16) = {15};
-Line Loop(17) = {-3,-12,8,10};
-Plane Surface(18) = {17};
-Line Loop(19) = {7,12,-2,11};
-Plane Surface(20) = {19};
-Line Loop(21) = {-4,-3,-2,1};
-Plane Surface(22) = {21};
-Line Loop(23) = {8,9,6,7};
-Plane Surface(24) = {23};
-
-Surface Loop(25) = {14,24,-18,22,16,-20};
-Complex Volume(26) = {25};
-
-n = 21;
-
-/*
-sizex = 4;
-sizey = 6;
-sizez = 4;
-*/
-sizex = n*length;
-sizey = n*height;
-sizez = n*depth;
-
-sizex = 9;
-sizey = 33;
-sizez = 9;
-
-
-Transfinite Line {2,4,9,7} = sizez;
-Transfinite Line {11,5,10,12} = sizey;
-Transfinite Line {3,1,6,8} = sizex;
-
-Transfinite Surface {14} = {5,4,3,1};
-Transfinite Surface {16} = {6,2,1,5};
-Transfinite Surface {18} = {6,2,7,8};
-Transfinite Surface {20} = {3,7,8,4};
-Transfinite Surface {22} = {3,7,2,1};
-Transfinite Surface {24} = {4,8,6,5};
-
-Recombine Surface {14,16,18,20,22,24};
-
-Transfinite Volume {26} = {3,7,2,1,4,8,6,5};
-
-Physical Surface(200) = {14,16,18,20,24};
-
-Physical Volume(100) = {26};
-/*
-Mesh.use_cut_plane = 1;
-Mesh.cut_planea = 0;
-Mesh.cut_planeb = 0;
-Mesh.cut_planec = 1;
-Mesh.cut_planed = -0.125;
-*/
\ No newline at end of file