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

Numeric/NumericEmbedded.cpp

+// 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;
+}

Numeric/NumericEmbedded.h

+// 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

Numeric/gmshAssembler.cpp

+// 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);
+  }
+}
+

Numeric/gmshTermOfFormulation.cpp

+// 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;
+}

contrib/Chaco/misc/timing.c

+/* 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;
+}

examples/cube.geo

-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