fix bugs in Gmsh_Vector and Gmsh_Matrix (scale & mult were wrong)

Numeric/GmshMatrix.h

@@ -21,7 +21,7 @@ class Gmsh_Vector
   Gmsh_Vector(int R) : r(R)
   {
     data = new SCALAR[r];
-    scale(0);
+    scale(0.);
   }
   Gmsh_Vector(const Gmsh_Vector<SCALAR> &other) : r(other.r)
   {
@@ -52,6 +52,9 @@ class Gmsh_Vector
   }
   inline void scale(const SCALAR s)
   {
+    if (s == 0.) 
+      for (int i = 0; i < r; ++i) data[i] = 0.;
+    else 
       for (int i = 0; i < r; ++i) data[i] *= s;
   }
 };
@@ -142,7 +145,7 @@ class Gmsh_Matrix
   Gmsh_Matrix(int R, int C) : r(R), c(C)
   {
     data = new SCALAR[r * c];
-    scale(0);
+    scale(0.);
   }
   Gmsh_Matrix(const Gmsh_Matrix<SCALAR> &other) : r(other.r), c(other.c)
   {
@@ -172,20 +175,22 @@ class Gmsh_Matrix
   {
     return data[i + r * j];
   }
-  inline void mult(const Gmsh_Matrix<SCALAR> &x, const Gmsh_Matrix<SCALAR> &b)
+  inline void mult(const Gmsh_Matrix<SCALAR> &x, Gmsh_Matrix<SCALAR> &b)
   {
-    for(int i = 0; i < b.size1(); i++)
-      for(int j = 0; j < b.size2(); j++)
-	for(int k = 0; k < size2(); k++)
+    b.scale(0.);
+    for(int i = 0; i < r; i++)
+      for(int j = 0; j < x.size2(); j++)
+	for(int k = 0; k < c; k++)
 	  b.data[i + r * j] += (*this)(i, k) * x(k, j);
   }
   inline void mult(const Gmsh_Vector<SCALAR> &x, Gmsh_Vector<SCALAR> &b)
   {
-    for(int i = 0; i < b.size(); i++)
-      for(int j = 0; j < b.size(); j++)
+    b.scale(0.);
+    for(int i = 0; i < r; i++)
+      for(int j = 0; j < c; j++)
 	b.data[i] += (*this)(i, j) * x(j);
   }
-  inline void blas_dgemm(const Gmsh_Matrix<SCALAR> & x, const Gmsh_Matrix<SCALAR> & b, 
+  inline void blas_dgemm(const Gmsh_Matrix<SCALAR> &x, Gmsh_Matrix<SCALAR> &b, 
 			 const SCALAR c_a = 1.0, const SCALAR c_b = 1.0)
   {
     Gmsh_Matrix<SCALAR> temp (x.size1(), b.size2());
@@ -256,6 +261,9 @@ class Gmsh_Matrix
   }  
   inline void scale(const double s)
   {
+    if(s == 0.)
+      for (int i = 0; i < r * c; ++i) data[i] = 0.;
+    else
       for (int i = 0; i < r * c; ++i) data[i] *= s;
   }
   inline void add(const double &a) 
@@ -308,7 +316,7 @@ class GSL_Vector
   }
   inline void scale(const double &y)
   {
-    if (y == 0.0) gsl_vector_set_zero(data);
+    if (y == 0.) gsl_vector_set_zero(data);
     else gsl_vector_scale(data, y);
   }
 };
@@ -352,7 +360,7 @@ class GSL_Matrix
   {
     return *gsl_matrix_ptr(data, i, j);
   }
-  inline void mult(const GSL_Matrix &x, const GSL_Matrix &b)
+  inline void mult(const GSL_Matrix &x, GSL_Matrix &b)
   {
     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
   }
@@ -426,7 +434,7 @@ class GSL_Matrix
   {
     gsl_blas_dgemv(CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
   }
-  inline void blas_dgemm(const GSL_Matrix & x, const GSL_Matrix& b, 
+  inline void blas_dgemm(const GSL_Matrix &x, GSL_Matrix &b, 
 			 const double c_a = 1.0, const double c_b = 1.0)
   {      
     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, c_a, x.data, b.data, c_b, data);
@@ -441,7 +449,7 @@ class GSL_Matrix
   }
   inline void scale(const double &m) 
   {
-    if (m == 0.0) gsl_matrix_set_zero(data);
+    if (m == 0.) gsl_matrix_set_zero(data);
     else gsl_matrix_scale(data, m);
   }
   inline void add(const double &a) 

Numeric/gmshTermOfFormulation.h

@@ -14,9 +14,9 @@ class gmshLinearSystem;
 class gmshFunction;
 
 #include <math.h>
-#include <GmshMatrix.h>
 #include <map>
 #include <vector>
+#include "GmshMatrix.h"
 
 class gmshAssembler;
 

Post/adaptiveData.cpp

@@ -841,7 +841,6 @@ template <class T>
 int adaptiveElements<T>::_zoomElement(int ielem, int level, GMSH_Post_Plugin *plug)
 {
   const int N = _coeffs->size1();
-  
   Double_Vector val(N),  res(adaptivePoint::all.size());
   Double_Vector valx(N), resx(adaptivePoint::all.size());
   Double_Vector valy(N), resy(adaptivePoint::all.size());
@@ -934,18 +933,18 @@ int adaptiveElements<T>::_zoomElement(int ielem, int level, GMSH_Post_Plugin *pl
 template <class T> 
 void adaptiveElements<T>::_changeResolution(int level, GMSH_Post_Plugin *plug, int *done)
 {
-  const int N = _coeffs->size1();
-  const int nbelm = _posX->size1();
-
-  double sf[100];
   T::create(level, _coeffs, _eexps);
 
+  if(!adaptivePoint::all.size()) return;
+
+  const int N = _coeffs->size1();
   if(_interpolate) delete _interpolate;
   _interpolate = new Double_Matrix(adaptivePoint::all.size(), N);
 
   if(_geometry) delete _geometry;
   _geometry = new Double_Matrix(adaptivePoint::all.size(), _posX->size2());
 
+  double sf[100];
   int kk = 0;
   for(std::set<adaptivePoint>::iterator it = adaptivePoint::all.begin(); 
       it != adaptivePoint::all.end(); ++it) {
@@ -961,6 +960,7 @@ void adaptiveElements<T>::_changeResolution(int level, GMSH_Post_Plugin *plug, i
     kk++;
   }
 
+  const int nbelm = _posX->size1();
   for(int i = 0; i < nbelm; ++i)
     done[i] = _zoomElement(i, level, plug);
 }
@@ -999,6 +999,7 @@ void adaptiveElements<T>::initData(PViewData *data, int step)
 
   int numNodes = getNumNodes();
   int numVal = _coeffs->size1() * numComp;
+  if(!numNodes || !numVal) return;
 
   _minVal = VAL_INF;
   _maxVal = -VAL_INF;
@@ -1083,8 +1084,6 @@ void adaptiveElements<T>::changeResolution(int level, double tol, GMSH_Post_Plug
   List_Reset(_listEle);
   *_numEle = 0;
   
-  if(!_posX->size1()) return;
-
   std::vector<int> done(_posX->size1(), 0);
 
   // We first do the adaptive stuff at level 2 and will only process