diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index 8e2daee1a296ab8242786b22231eacef87f69afb..698e34ee3a0cee9628ffc9a4fbe3a76fbd7ce422 100644
--- a/Numeric/GmshMatrix.h
+++ b/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,7 +52,10 @@ class Gmsh_Vector
}
inline void scale(const SCALAR s)
{
- for (int i = 0; i < r; ++i) data[i] *= 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.data[i + r *j] += (*this)(i, k) * x(k, j);
+ 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,7 +261,10 @@ class Gmsh_Matrix
}
inline void scale(const double s)
{
- for (int i = 0; i < r * c; ++i) data[i] *= 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,10 +434,10 @@ 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);
+ gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, c_a, x.data, b.data, c_b, data);
}
inline gsl_matrix_view touchSubmatrix(int i0, int ni, int j0, int nj)
{
@@ -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)
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index 83f2eb8a5a4c06a4c1fdd9c69249d7c5ff38c22c..00d63bee781eb46e4f3b0e45b8b51f8a3f033847 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/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;
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index fe1d0a3debb7b8a722fa742159ba4d31e08b91e1..92b10137d4b346ef47229403f23fab96a5c4ae6d 100644
--- a/Post/adaptiveData.cpp
+++ b/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