diff --git a/Common/AdaptiveViews.cpp b/Common/AdaptiveViews.cpp
index 0c21eb1661aad2693268e69a02f11b28198b41f4..68649cb22bcc7ce14711e65338f9449a6e57ebb4 100644
--- a/Common/AdaptiveViews.cpp
+++ b/Common/AdaptiveViews.cpp
@@ -64,6 +64,7 @@ void _triangle::clean ()
}
void _triangle::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps)
{
+ printf("creating the sub-elements\n");
int level = 0;
clean();
_point *p1 = _point::New ( 0,0,0, coeffs, eexps);
@@ -71,6 +72,8 @@ void _triangle::Create (int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexp
_point *p3 = _point::New ( 1,0,0, coeffs, eexps);
_triangle *t = new _triangle(p1,p2,p3);
Recur_Create (t, maxlevel,level,coeffs,eexps) ;
+
+ printf("%d _triangle and %d _point created\n",_triangle::all_triangles.size(),_point::all_points.size());
}
void _triangle::Recur_Create (_triangle *t, int maxlevel, int level , Double_Matrix *coeffs, Double_Matrix *eexps)
{
@@ -93,6 +96,7 @@ void _triangle::Recur_Create (_triangle *t, int maxlevel, int level , Double_Mat
_triangle *t4 = new _triangle (p12,p23,p13);
Recur_Create (t4, maxlevel,level,coeffs,eexps);
t->t[0]=t1;t->t[1]=t2;t->t[2]=t3;t->t[3]=t4;
+
}
void _triangle::Error ( double AVG , double tol )
@@ -173,31 +177,63 @@ void _triangle::Recur_Error ( _triangle *t, double AVG, double tol )
}
}
+static double t0,t1,t2,t3;
+
void Adaptive_Post_View:: zoomElement (Post_View * view ,
- int ielem , int level, GMSH_Post_Plugin *plug)
+ int ielem ,
+ int level,
+ GMSH_Post_Plugin *plug)
{
std::set<_point>::iterator it = _point::all_points.begin();
std::set<_point>::iterator ite = _point::all_points.end();
+
+ clock_t c0 = clock();
+
+ const int N = _coefs->size1();
+ Double_Vector val ( N ), res(_point::all_points.size());
+ Double_Matrix xyz (3,3);
+ Double_Matrix XYZ (_point::all_points.size(),3);
+
+ for ( int k=0;k<3;++k)
+ {
+ xyz(k,0) = (*_STposX) ( ielem , k );
+ xyz(k,1) = (*_STposY) ( ielem , k );
+ xyz(k,2) = (*_STposZ) ( ielem , k );
+ }
+
+ for ( int k=0;k<N;++k)
+ {
+ val(k) = (*_STval )( ielem , k );
+ }
+ _Interpolate->mult(val,res);
+ _Geometry->mult(xyz,XYZ);
+ clock_t c1 = clock();
+
+ int kk=0;
for ( ; it !=ite ; ++it)
{
_point *p = (_point*) &(*it);
- p->val = 0;
- for ( int k=0;k<_coefs->size1();++k)
- {
- p->val += it->shape_functions[k] * (*_STval )( ielem , k );
- }
- p->X = (*_STposX) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposX) ( ielem , 1 ) * p->x + (*_STposX) ( ielem , 2 ) * p->y;
- p->Y = (*_STposY) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposY) ( ielem , 1 ) * p->x + (*_STposY) ( ielem , 2 ) * p->y;
- p->Z = (*_STposZ) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposZ) ( ielem , 1 ) * p->x + (*_STposZ) ( ielem , 2 ) * p->y;
-
- if (level == 0)
- {
- if (min > p->val) min = p->val;
- if (max < p->val) max = p->val;
- }
+ p->val = res(kk);
+ // for ( int k=0;k<N;++k)
+ // {
+ // p->val += p->shape_functions[k] * (*_STval )( ielem , k );
+ // }
+ // p->X = (*_STposX) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposX) ( ielem , 1 ) * p->x + (*_STposX) ( ielem , 2 ) * p->y;
+ // p->Y = (*_STposY) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposY) ( ielem , 1 ) * p->x + (*_STposY) ( ielem , 2 ) * p->y;
+ // p->Z = (*_STposZ) ( ielem , 0 ) * ( 1.-p->x-p->y) + (*_STposZ) ( ielem , 1 ) * p->x + (*_STposZ) ( ielem , 2 ) * p->y;
+
+ p->X = XYZ(kk,0);
+ p->Y = XYZ(kk,1);
+ p->Z = XYZ(kk,2);
+
+ if (min > p->val) min = p->val;
+ if (max < p->val) max = p->val;
+ kk++;
}
+ clock_t c2 = clock();
+
std::list<_triangle*>::iterator itt = _triangle::all_triangles.begin();
std::list<_triangle*>::iterator itte = _triangle::all_triangles.end();
@@ -210,8 +246,9 @@ void Adaptive_Post_View:: zoomElement (Post_View * view ,
if (plug)
plug->assign_specific_visibility ();
else
- _triangle::Error ( max-min, tol );
-
+ _triangle::Error ( max-min, tol );
+
+ clock_t c3 = clock();
itt = _triangle::all_triangles.begin();
for ( ;itt != itte ; itt++)
{
@@ -235,6 +272,13 @@ void Adaptive_Post_View:: zoomElement (Post_View * view ,
view->NbST++;
}
}
+ clock_t c4 = clock();
+
+ t0 += (double)(c1-c0)/CLOCKS_PER_SEC;
+ t1 += (double)(c2-c1)/CLOCKS_PER_SEC;
+ t2 += (double)(c3-c2)/CLOCKS_PER_SEC;
+ t3 += (double)(c4-c3)/CLOCKS_PER_SEC;
+
}
void Adaptive_Post_View:: setAdaptiveResolutionLevel (Post_View * view , int level, GMSH_Post_Plugin *plug)
@@ -245,17 +289,48 @@ void Adaptive_Post_View:: setAdaptiveResolutionLevel (Post_View * view , int lev
_triangle::Create ( level, _coefs, _eexps );
+ std::set<_point>::iterator it = _point::all_points.begin();
+ std::set<_point>::iterator ite = _point::all_points.end();
+
+ const int N = _coefs->size1();
+ if (_Interpolate)
+ delete _Interpolate;
+ if (_Geometry)
+ delete _Geometry;
+ _Interpolate = new Double_Matrix ( _point::all_points.size(), N);
+ _Geometry = new Double_Matrix ( _point::all_points.size(), 3);
+
+ int kk=0;
+ for ( ; it !=ite ; ++it)
+ {
+ _point *p = (_point*) &(*it);
+ for ( int k=0;k<N;++k)
+ {
+ (*_Interpolate)(kk,k) = p->shape_functions[k];
+ }
+ (*_Geometry)(kk,0) = ( 1.-p->x-p->y);
+ (*_Geometry)(kk,1) = p->x;
+ (*_Geometry)(kk,2) = p->y;
+ kk++;
+ }
+
List_Delete(view->ST); view->ST = 0;
view->NbST = 0;
/// for now, that's all we do, 1 TS
view->NbTimeStep=1;
int nbelm = _STposX->size1();
- view->ST = List_Create ( nbelm * (int) pow(4.,level), nbelm *12, sizeof(double));
+ view->ST = List_Create ( nbelm * 4, nbelm , sizeof(double));
+
+
+ t0=t1 = t2 = t3 = 0;
for ( int i=0;i<nbelm;++i)
{
zoomElement ( view , i , level, plug);
}
+
+ printf("finished %g %g %g %g\n",t0,t1,t2,t3);
+
view->Changed = 1;
presentZoomLevel = level;
presentTol = tol;
@@ -322,6 +397,7 @@ Adaptive_Post_View:: Adaptive_Post_View (Post_View *view, List_T *_c , List_T *_
: tol(1.e-3)
{
+ _Interpolate = _Geometry = 0;
_coefs = new Double_Matrix ( List_Nbr (_c) , List_Nbr (_c) );
_eexps = new Double_Matrix ( List_Nbr (_c) , 3 );
@@ -359,5 +435,7 @@ Adaptive_Post_View::~Adaptive_Post_View()
delete _STposY;
delete _STposZ;
delete _STval;
+ if(_Interpolate)delete _Interpolate;
+ if(_Geometry)delete _Geometry;
_triangle::clean();
}
diff --git a/Common/GmshMatrix.h b/Common/GmshMatrix.h
index d20a53fda3e0efd6a35a791bd78edd8ee2c42e9b..0658f54879bb3cae514197d911b9090cf7c57b14 100644
--- a/Common/GmshMatrix.h
+++ b/Common/GmshMatrix.h
@@ -1,7 +1,6 @@
#ifndef _GMSH_BOOSTMATRIX_
#define _GMSH_BOOSTMATRIX_
-#ifndef HAVE_BOOST
template <class SCALAR>
class Gmsh_Matrix
{
@@ -46,17 +45,124 @@ public:
}
};
-typedef Gmsh_Matrix<double> Double_Matrix;
-typedef Gmsh_Matrix<int> Int_Matrix;
+template <class SCALAR>
+class Gmsh_Vector
+{
+private:
+ int r;
+public:
+ inline int size() const {return r;}
+ SCALAR *data;
+ ~Gmsh_Vector() {delete [] data;}
+ Gmsh_Vector(int R)
+ : r(R)
+ {
+ data = new SCALAR [r];
+ scal(0);
+ }
+ Gmsh_Vector(const Gmsh_Vector<SCALAR> &other)
+ : r(other.r)
+ {
+ data = new double [r];
+ copy ( other.data ) ;
+ }
+ inline SCALAR operator () (int i) const
+ {
+ return data[i];
+ }
+ inline SCALAR & operator () (int i)
+ {
+ return data[i];
+ }
+ inline SCALAR operator *(const Gmsh_Vector<SCALAR> & other)
+ {
+ throw;
+ }
+ inline void scal ( const SCALAR s)
+ {
+ for (int i=0;i<r;++i)data[i]*=s;
+ }
+ inline void copy ( const SCALAR **other)
+ {
+ for (int i=0;i<r;++i)data[i]=other.data[i];
+ }
+};
+
+#ifdef HAVE_GSL
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+class GSL_Vector
+{
+private:
+ int r,c;
+public:
+ inline int size() const {return r;}
+ gsl_vector *data;
+ ~GSL_Vector() {gsl_vector_free (data);}
+ GSL_Vector(int R)
+ : r(R)
+ {
+ data = gsl_vector_calloc (r);
+ }
+ GSL_Vector(const GSL_Vector&other)
+ : r(other.r)
+ {
+ data = gsl_vector_calloc (r);
+ gsl_vector_memcpy (data, other.data);
+ }
+ inline double operator () (int i) const
+ {
+ return gsl_vector_get (data,i);
+ }
+ inline double & operator () (int i)
+ {
+ return *gsl_vector_ptr (data,i);
+ }
+};
+class GSL_Matrix
+{
+private:
+ int r,c;
+public:
+ inline int size1() const {return r;}
+ inline int size2() const {return c;}
+ gsl_matrix *data;
+ ~GSL_Matrix() {gsl_matrix_free (data);}
+ GSL_Matrix(int R,int C)
+ : r(R),c(C)
+ {
+ data = gsl_matrix_calloc (r, c);
+ }
+ GSL_Matrix(const GSL_Matrix&other)
+ : r(other.r),c(other.c)
+ {
+ data = gsl_matrix_calloc (r, c);
+ gsl_matrix_memcpy (data, other.data);
+ }
+ inline double operator () (int i, int j) const
+ {
+ return gsl_matrix_get (data,i,j);
+ }
+ inline double & operator () (int i, int j)
+ {
+ return *gsl_matrix_ptr (data,i,j);
+ }
+ inline void mult(const GSL_Matrix & x, const GSL_Matrix & b)
+ {
+ gsl_blas_dgemm (CblasNoTrans,CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
+ }
+ inline void mult (const GSL_Vector & x, GSL_Vector & b )
+ {
+ gsl_blas_dgemv (CblasNoTrans, 1.0, data, x.data, 1.0, b.data);
+ }
+};
+typedef GSL_Matrix Double_Matrix;
+typedef GSL_Vector Double_Vector;
#else
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/blas.hpp>
-#include <boost/numeric/ublas/io.hpp>
-typedef numerics::matrix<double, numerics::column_major> Double_Matrix;
-typedef numerics::vector<double> Double_Vector;
-typedef numerics::matrix<int, numerics::column_major> Int_Matrix;
-typedef numerics::vector<int> Int_Vector;
+typedef Gmsh_Matrix<double> Double_Matrix;
+typedef Gmsh_Vector<double> Double_Vector;
#endif
+typedef Gmsh_Matrix<int> Int_Matrix;
+typedef Gmsh_Vector<int> Int_Vector;
#endif
diff --git a/Common/Views.h b/Common/Views.h
index eacb34079b59be7faf12498f60b0dd3b6df71b68..9c801221891da1b8869971f8c2fcfdbeeb73383e 100644
--- a/Common/Views.h
+++ b/Common/Views.h
@@ -108,6 +108,8 @@ class Adaptive_Post_View
Double_Matrix * _STposY;
Double_Matrix * _STposZ;
Double_Matrix * _STval;
+ Double_Matrix * _Interpolate;
+ Double_Matrix * _Geometry;
public:
Adaptive_Post_View (Post_View *view, List_T *_coeffs, List_T *_eexps);
~Adaptive_Post_View();