*** empty log message ***

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();
 }

Common/GmshMatrix.h

 #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

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();