diff --git a/Fltk/classificationEditor.cpp b/Fltk/classificationEditor.cpp
index edee9e5e628e0eea0d404bcda5c7bfb34c3ab39c..42f7352717e83e96d6648730a07b0828c2a513fa 100644
--- a/Fltk/classificationEditor.cpp
+++ b/Fltk/classificationEditor.cpp
@@ -331,7 +331,7 @@ static void class_color_cb(Fl_Widget* w, void* data)
     while(it != e->_faces.end()){
       GFace *gf = *it;
       for(unsigned int i = 0; i < gf->triangles.size(); i++){
-            tris.push_back(new MTri3(gf->triangles [i], 0));
+        tris.push_back(new MTri3(gf->triangles[i], 0));
       }
       gf->triangles.clear();
       ++it;
diff --git a/Fltk/partitionDialog.cpp b/Fltk/partitionDialog.cpp
index d49f4e6dd7b7520afdd08f583479141bbb90725c..c268cd6bd5f1a130a90a39acde6f5cb74c653143 100644
--- a/Fltk/partitionDialog.cpp
+++ b/Fltk/partitionDialog.cpp
@@ -30,8 +30,6 @@
 
 extern Context_T CTX;
 
-#define GMSH_WINDOW_BOX FL_FLAT_BOX
-
 #if defined(HAVE_CHACO) || defined(HAVE_METIS)
 
 // Forward declarations of some callbacks
diff --git a/Graphics/drawContext.h b/Graphics/drawContext.h
index 115d8100c7001a05622f9eb688e262acc4d6c019..d00de75966d74f9a5f4c5e98d013f037db3ea766 100644
--- a/Graphics/drawContext.h
+++ b/Graphics/drawContext.h
@@ -16,6 +16,8 @@ class drawTransform {
   drawTransform(){}
   virtual ~drawTransform(){}
   virtual void transform(double &x, double &y, double &z){}
+  virtual void transformOneForm(double &x, double &y, double &z){}
+  virtual void transformTwoForm(double &x, double &y, double &z){}
   virtual void setMatrix(double mat[3][3], double tra[3]){}
 };
 
@@ -77,9 +79,17 @@ class drawContext {
   void setTransform(drawTransform *transform){ _transform = transform; }
   drawTransform *getTransform(){ return _transform; }
   void transform(double &x, double &y, double &z)
-  {
+  { 
     if(_transform) _transform->transform(x, y, z); 
   }
+  void transformOneForm(double &x, double &y, double &z)
+  {
+    if(_transform) _transform->transformOneForm(x, y, z); 
+  }
+  void transformTwoForm(double &x, double &y, double &z)
+  {
+    if(_transform) _transform->transformTwoForm(x, y, z); 
+  }
   void buildRotationMatrix();
   void setQuaternion(double p1x, double p1y, double p2x, double p2y);
   void addQuaternion(double p1x, double p1y, double p2x, double p2y);
diff --git a/Graphics/drawGeom.cpp b/Graphics/drawGeom.cpp
index bc80f4e226a78c920b12018cedc8424855dfb5e1..a461b67c3a0bf69a48695f4091cb3579719e030d 100644
--- a/Graphics/drawGeom.cpp
+++ b/Graphics/drawGeom.cpp
@@ -179,7 +179,7 @@ class drawGEdge {
       glColor4ubv((GLubyte *) & CTX.color.geom.tangents);
       double x = p.x(), y = p.y(), z = p.z();
       _ctx->transform(x, y, z);
-      // FIXME: transform the tangent
+      _ctx->transformOneForm(der[0], der[1], der[2]);
       _ctx->drawVector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
                        CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
                        x, y, z, der[0], der[1], der[2], CTX.geom.light);
@@ -294,7 +294,7 @@ class drawGFace {
       glColor4ubv((GLubyte *) & CTX.color.geom.normals);
       double x = p.x(), y = p.y(), z = p.z();
       _ctx->transform(x, y, z);
-      // FIXME: transform the normal
+      _ctx->transformTwoForm(n[0], n[1], n[2]);
       _ctx->drawVector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
                        CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius,
                        x, y, z, n[0], n[1], n[2], CTX.geom.light);
@@ -364,7 +364,7 @@ class drawGFace {
       glColor4ubv((GLubyte *) & CTX.color.geom.normals);
       double x = p.x(), y = p.y(), z = p.z();
       _ctx->transform(x, y, z);
-      // FIXME: transform the normal
+      _ctx->transformTwoForm(n[0], n[1], n[2]);
       _ctx->drawVector(CTX.vector_type, 0, CTX.arrow_rel_head_radius, 
                        CTX.arrow_rel_stem_length, CTX.arrow_rel_stem_radius, 
                        x, y, z, n[0], n[1], n[2], CTX.geom.light);
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index cebcaf348f58fd09a8fadcf00994d2a4eb910301..ca16a8c848155b5951688aa4aa3148e24eef11b9 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -17,8 +17,6 @@
 #include <functional>
 #include <list>
 #include <math.h>
-#include "GFace.h"
-#include "PView.h"
 #include "GmshMessage.h"
 
 class BDS_Edge;
@@ -100,9 +98,9 @@ public:
   }
   BDS_Vector& operator /= (const double &v)
   {
-    x/=v;
-    y/=v;
-    z/=v;
+    x /= v;
+    y /= v;
+    z /= v;
     return *this;
   }
   BDS_Vector operator / (const double &v)
@@ -135,11 +133,8 @@ public:
     return (x * v.x + y * v.y + z * v.z);
   }
   BDS_Vector(const BDS_Point &p2, const BDS_Point &p1);
-  
   BDS_Vector(const double X=0., const double Y=0., const double Z=0.)
-    : x(X), y(Y), z(Z)
-  {
-  }
+    : x(X), y(Y), z(Z) {}
   static double t;
 };
 
@@ -149,7 +144,7 @@ class BDS_Point
   // second one is dictated by charecteristic lengths at points and is
   // propagated
   double _lcBGM, _lcPTS;
-public:
+ public:
   double X, Y, Z;
   double u, v;
   bool config_modified;
@@ -178,16 +173,14 @@ public:
   void getTriangles(std::list<BDS_Face *> &t) const;
   BDS_Point(int id, double x=0, double y=0, double z=0)
     : _lcBGM(1.e22), _lcPTS(1.e22), X(x), Y(y), Z(z), u(0), v(0),
-      config_modified(true), iD(id), g(0)
-  {         
-  }
+      config_modified(true), iD(id), g(0) {}
 };
 
 class BDS_Edge
 {
   double _length;
   std::vector<BDS_Face*> _faces;
-public:
+ public:
   bool deleted;
   BDS_Point *p1, *p2;
   BDS_GeomEntity *g;
@@ -270,7 +263,7 @@ public:
 
 class BDS_Face
 {
-public:
+ public:
   bool deleted;
   BDS_Edge *e1, *e2, *e3, *e4;
   BDS_GeomEntity *g;
@@ -328,7 +321,7 @@ public:
 
 class GeomLessThan
 {
-public:
+ public:
   bool operator()(const BDS_GeomEntity *ent1, const BDS_GeomEntity *ent2) const
   {
     return *ent1 < *ent2;
@@ -337,7 +330,7 @@ public:
 
 class PointLessThan
 {
-public:
+ public:
   bool operator()(const BDS_Point *ent1, const BDS_Point *ent2) const
   {
     return *ent1 < *ent2;
@@ -346,7 +339,7 @@ public:
 
 class PointLessThanLexicographic
 {
-public:
+ public:
   static double t;
   bool operator()(const BDS_Point *ent1, const BDS_Point *ent2) const
   {
@@ -361,7 +354,7 @@ public:
 
 class EdgeLessThan
 {
-public:
+ public:
   bool operator()(const BDS_Edge *ent1, const BDS_Edge *ent2) const
   {
     return *ent1 < *ent2;
@@ -420,7 +413,7 @@ struct EdgeToRecover
 
 class BDS_Mesh 
 {    
-public:
+ public:
   int MAXPOINTNUMBER, SNAP_SUCCESS, SNAP_FAILURE;
   double Min[3], Max[3], LC;
   double scalingU, scalingV;
diff --git a/Mesh/meshGFaceQuadrilateralize.cpp b/Mesh/meshGFaceQuadrilateralize.cpp
index 16ae863c46ea3c6ba2be01fe04d08c71cf614f9f..7adff0cb6d4e5d9bb4eec8c5f1d33e7feabd0d0d 100644
--- a/Mesh/meshGFaceQuadrilateralize.cpp
+++ b/Mesh/meshGFaceQuadrilateralize.cpp
@@ -6,6 +6,7 @@
 #include "meshGFaceQuadrilateralize.h"
 #include "GmshMessage.h"
 #include "Numeric.h"
+#include "GFace.h"
 #include "meshGFaceDelaunayInsertion.h"
 #include "meshGFaceOptimize.h"
 #include "meshGFaceBDS.h"
diff --git a/Plugin/ExtractEdges.cpp b/Plugin/ExtractEdges.cpp
index 03a8f60fd49e7cd11bb9b074a7af4312c0531107..1b196b02708a49a2ce60341cad6095e0076cfcb1 100644
--- a/Plugin/ExtractEdges.cpp
+++ b/Plugin/ExtractEdges.cpp
@@ -95,7 +95,6 @@ PView *GMSH_ExtractEdgesPlugin::execute(PView *v)
     ++it;
   }
 
-
   data2->setName(data1->getName() + "_ExtractEdges");
   data2->setFileName(data1->getName() + "_ExtractEdges.pos");
   data2->finalize();
diff --git a/Post/adaptiveData.cpp b/Post/adaptiveData.cpp
index d821604fb309ce1e383b7df5355e0bf8d1865b1d..8df5037f3d531009fa2c3c33ce8dabb3d2329af5 100644
--- a/Post/adaptiveData.cpp
+++ b/Post/adaptiveData.cpp
@@ -8,8 +8,15 @@
 #include <set>
 #include "Plugin.h"
 #include "adaptiveData.h"
+#include "OS.h"
+
+std::set<adaptivePoint> adaptiveLine::allPoints;
+std::set<adaptivePoint> adaptiveTriangle::allPoints;
+std::set<adaptivePoint> adaptiveQuadrangle::allPoints;
+std::set<adaptivePoint> adaptiveTetrahedron::allPoints;
+std::set<adaptivePoint> adaptiveHexahedron::allPoints;
+std::set<adaptivePoint> adaptivePrism::allPoints;
 
-std::set<adaptivePoint> adaptivePoint::all;
 std::list<adaptiveLine*> adaptiveLine::all;
 std::list<adaptiveTriangle*> adaptiveTriangle::all;
 std::list<adaptiveQuadrangle*> adaptiveQuadrangle::all;
@@ -37,7 +44,7 @@ static void cleanElement()
   for(typename std::list<T*>::iterator it = T::all.begin(); it != T::all.end(); ++it)
     delete *it;
   T::all.clear();
-  adaptivePoint::all.clear();
+  T::allPoints.clear();
 }
 
 static void computeShapeFunctions(Double_Matrix *coeffs, Double_Matrix *eexps,
@@ -57,34 +64,31 @@ static void computeShapeFunctions(Double_Matrix *coeffs, Double_Matrix *eexps,
   }
 }
 
-adaptivePoint *adaptivePoint::create(double x, double y, double z,
-				     Double_Matrix *coeffs, Double_Matrix *eexps)
+adaptivePoint *adaptivePoint::add(double x, double y, double z,
+                                  std::set<adaptivePoint> &allPoints)
 {
   adaptivePoint p;
   p.x = x;
   p.y = y;
   p.z = z;
-  std::set<adaptivePoint>::iterator it = all.find(p);
-  if(it == all.end()) {
-    all.insert(p);
-    it = all.find(p);
-    computeShapeFunctions(coeffs, eexps, x, y, z, (double*)it->shapeFunctions);
+  std::set<adaptivePoint>::iterator it = allPoints.find(p);
+  if(it == allPoints.end()){
+    allPoints.insert(p);
+    it = allPoints.find(p);
   }
   return (adaptivePoint*)&(*it);
 }
 
-void adaptiveLine::create(int maxlevel, 
-			  Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveLine::create(int maxlevel)
 {
   cleanElement<adaptiveLine>();
-  adaptivePoint *p1 = adaptivePoint::create(-1, 0, 0, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(1, 0, 0, coeffs, eexps);
+  adaptivePoint *p1 = adaptivePoint::add(-1, 0, 0, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(1, 0, 0, allPoints);
   adaptiveLine *t = new adaptiveLine(p1, p2);
-  recurCreate(t, maxlevel, 0, coeffs, eexps);
+  recurCreate(t, maxlevel, 0);
 }
 
-void adaptiveLine::recurCreate(adaptiveLine *e, int maxlevel, int level,
-			       Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveLine::recurCreate(adaptiveLine *e, int maxlevel, int level)
 {
   all.push_back(e);
   if(level++ >= maxlevel) return;
@@ -92,13 +96,13 @@ void adaptiveLine::recurCreate(adaptiveLine *e, int maxlevel, int level,
   // p1    p12    p2
   adaptivePoint *p1 = e->p[0];
   adaptivePoint *p2 = e->p[1];
-  adaptivePoint *p12 = adaptivePoint::create
+  adaptivePoint *p12 = adaptivePoint::add
     ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
-     coeffs, eexps);
+     allPoints);
   adaptiveLine *e1 = new adaptiveLine(p1, p12);
-  recurCreate(e1, maxlevel, level, coeffs, eexps);
+  recurCreate(e1, maxlevel, level);
   adaptiveLine *e2 = new adaptiveLine(p12, p2);
-  recurCreate(e2, maxlevel, level, coeffs, eexps);
+  recurCreate(e2, maxlevel, level);
   e->e[0] = e1;
   e->e[1] = e2;
 }
@@ -149,19 +153,17 @@ void adaptiveLine::recurError(adaptiveLine *e, double AVG, double tol)
   }
 }
 
-void adaptiveTriangle::create(int maxlevel, 
-			      Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveTriangle::create(int maxlevel)
 {
   cleanElement<adaptiveTriangle>();
-  adaptivePoint *p1 = adaptivePoint::create(0, 0, 0, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(0, 1, 0, coeffs, eexps);
-  adaptivePoint *p3 = adaptivePoint::create(1, 0, 0, coeffs, eexps);
+  adaptivePoint *p1 = adaptivePoint::add(0, 0, 0, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(0, 1, 0, allPoints);
+  adaptivePoint *p3 = adaptivePoint::add(1, 0, 0, allPoints);
   adaptiveTriangle *t = new adaptiveTriangle(p1, p2, p3);
-  recurCreate(t, maxlevel, 0, coeffs, eexps);
+  recurCreate(t, maxlevel, 0);
 }
 
-void adaptiveTriangle::recurCreate(adaptiveTriangle *t, int maxlevel, int level,
-				   Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveTriangle::recurCreate(adaptiveTriangle *t, int maxlevel, int level)
 {
   all.push_back(t);
   if(level++ >= maxlevel) return;
@@ -172,23 +174,23 @@ void adaptiveTriangle::recurCreate(adaptiveTriangle *t, int maxlevel, int level,
   adaptivePoint *p1 = t->p[0];
   adaptivePoint *p2 = t->p[1];
   adaptivePoint *p3 = t->p[2];
-  adaptivePoint *p12 = adaptivePoint::create
+  adaptivePoint *p12 = adaptivePoint::add
     ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p13 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p13 = adaptivePoint::add
     ((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5, (p1->z + p3->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p23 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p23 = adaptivePoint::add
     ((p3->x + p2->x) * 0.5, (p3->y + p2->y) * 0.5, (p3->z + p2->z) * 0.5, 
-     coeffs, eexps);
+     allPoints);
   adaptiveTriangle *t1 = new adaptiveTriangle(p1, p12, p13);
-  recurCreate(t1, maxlevel, level, coeffs, eexps);
+  recurCreate(t1, maxlevel, level);
   adaptiveTriangle *t2 = new adaptiveTriangle(p2, p23, p12);
-  recurCreate(t2, maxlevel, level, coeffs, eexps);
+  recurCreate(t2, maxlevel, level);
   adaptiveTriangle *t3 = new adaptiveTriangle(p3, p13, p23);
-  recurCreate(t3, maxlevel, level, coeffs, eexps);
+  recurCreate(t3, maxlevel, level);
   adaptiveTriangle *t4 = new adaptiveTriangle(p12, p23, p13);
-  recurCreate(t4, maxlevel, level, coeffs, eexps);
+  recurCreate(t4, maxlevel, level);
   t->e[0] = t1;
   t->e[1] = t2;
   t->e[2] = t3;
@@ -263,20 +265,18 @@ void adaptiveTriangle::recurError(adaptiveTriangle *t, double AVG, double tol)
   }
 }
 
-void adaptiveQuadrangle::create(int maxlevel, 
-				Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveQuadrangle::create(int maxlevel)
 {
   cleanElement<adaptiveQuadrangle>();
-  adaptivePoint *p1 = adaptivePoint::create(-1, -1, 0, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(1, -1, 0, coeffs, eexps);
-  adaptivePoint *p3 = adaptivePoint::create(1, 1, 0, coeffs, eexps);
-  adaptivePoint *p4 = adaptivePoint::create(-1, 1, 0, coeffs, eexps);
+  adaptivePoint *p1 = adaptivePoint::add(-1, -1, 0, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(1, -1, 0, allPoints);
+  adaptivePoint *p3 = adaptivePoint::add(1, 1, 0, allPoints);
+  adaptivePoint *p4 = adaptivePoint::add(-1, 1, 0, allPoints);
   adaptiveQuadrangle *q = new adaptiveQuadrangle(p1, p2, p3, p4);
-  recurCreate(q, maxlevel, 0, coeffs, eexps);
+  recurCreate(q, maxlevel, 0);
 }
 
-void adaptiveQuadrangle::recurCreate(adaptiveQuadrangle *q, int maxlevel, int level,
-				     Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveQuadrangle::recurCreate(adaptiveQuadrangle *q, int maxlevel, int level)
 {
   all.push_back(q);
   if(level++ >= maxlevel) return;
@@ -288,29 +288,29 @@ void adaptiveQuadrangle::recurCreate(adaptiveQuadrangle *q, int maxlevel, int le
   adaptivePoint *p2 = q->p[1];
   adaptivePoint *p3 = q->p[2];
   adaptivePoint *p4 = q->p[3];
-  adaptivePoint *p12 = adaptivePoint::create
+  adaptivePoint *p12 = adaptivePoint::add
     ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p23 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p23 = adaptivePoint::add
     ((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, (p2->z + p3->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p34 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p34 = adaptivePoint::add
     ((p3->x + p4->x) * 0.5, (p3->y + p4->y) * 0.5, (p3->z + p4->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p14 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p14 = adaptivePoint::add
     ((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5, (p1->z + p4->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *pc = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pc = adaptivePoint::add
     ((p1->x + p2->x + p3->x + p4->x) * 0.25, (p1->y + p2->y + p3->y + p4->y) * 0.25, 
-     (p1->z + p2->z + p3->z + p4->z) * 0.25, coeffs, eexps);
+     (p1->z + p2->z + p3->z + p4->z) * 0.25, allPoints);
   adaptiveQuadrangle *q1 = new adaptiveQuadrangle(p1, p12, pc, p14);
-  recurCreate(q1, maxlevel, level, coeffs, eexps);
+  recurCreate(q1, maxlevel, level);
   adaptiveQuadrangle *q2 = new adaptiveQuadrangle(p2, p23, pc, p12);
-  recurCreate(q2, maxlevel, level, coeffs, eexps);
+  recurCreate(q2, maxlevel, level);
   adaptiveQuadrangle *q3 = new adaptiveQuadrangle(p3, p34, pc, p23);
-  recurCreate(q3, maxlevel, level, coeffs, eexps);
+  recurCreate(q3, maxlevel, level);
   adaptiveQuadrangle *q4 = new adaptiveQuadrangle(p4, p14, pc, p34);
-  recurCreate(q4, maxlevel, level, coeffs, eexps);
+  recurCreate(q4, maxlevel, level);
   q->e[0] = q1;
   q->e[1] = q2;
   q->e[2] = q3;
@@ -385,170 +385,18 @@ void adaptiveQuadrangle::recurError(adaptiveQuadrangle *q, double AVG, double to
   }
 }
 
-void adaptivePrism::create(int maxlevel, 
-				Double_Matrix *coeffs, Double_Matrix *eexps)
-{
-  cleanElement<adaptivePrism>();
-  adaptivePoint *p1 = adaptivePoint::create(0, 0, -1, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(1, 0, -1, coeffs, eexps);
-  adaptivePoint *p3 = adaptivePoint::create(0, 1, -1, coeffs, eexps);
-  adaptivePoint *p4 = adaptivePoint::create(0, 0, 1, coeffs, eexps);
-  adaptivePoint *p5 = adaptivePoint::create(1, 0, 1, coeffs, eexps);
-  adaptivePoint *p6 = adaptivePoint::create(0, 1, 1, coeffs, eexps);
-  adaptivePrism *p = new adaptivePrism(p1, p2, p3, p4, p5, p6);
-  recurCreate(p, maxlevel, 0, coeffs, eexps);
-}
-
-void adaptivePrism::recurCreate(adaptivePrism *p, int maxlevel, int level,
-				     Double_Matrix *coeffs, Double_Matrix *eexps)
-{
-  all.push_back(p);
-  if(level++ >= maxlevel) return;
-
-  // p4   p34    p3
-  // p14  pc     p23
-  // p1   p12    p2
-  adaptivePoint *p1 = p->p[0];
-  adaptivePoint *p2 = p->p[1];
-  adaptivePoint *p3 = p->p[2];
-  adaptivePoint *p4 = p->p[3];
-  adaptivePoint *p5 = p->p[4];
-  adaptivePoint *p6 = p->p[5];
-  adaptivePoint *p14 = adaptivePoint::create
-    ((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5, (p1->z + p4->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p25 = adaptivePoint::create
-    ((p2->x + p5->x) * 0.5, (p2->y + p5->y) * 0.5, (p2->z + p5->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p36 = adaptivePoint::create
-    ((p3->x + p6->x) * 0.5, (p3->y + p6->y) * 0.5, (p3->z + p6->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p12 = adaptivePoint::create
-    ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p23 = adaptivePoint::create
-    ((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, (p2->z + p3->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p31 = adaptivePoint::create
-    ((p3->x + p1->x) * 0.5, (p3->y + p1->y) * 0.5, (p3->z + p1->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p1425 = adaptivePoint::create
-    ((p14->x + p25->x) * 0.5, (p14->y + p25->y) * 0.5, (p14->z + p25->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p2536 = adaptivePoint::create
-    ((p25->x + p36->x) * 0.5, (p25->y + p36->y) * 0.5, (p25->z + p36->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p3614 = adaptivePoint::create
-    ((p36->x + p14->x) * 0.5, (p36->y + p14->y) * 0.5, (p36->z + p14->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p45 = adaptivePoint::create
-    ((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5, (p4->z + p5->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p56 = adaptivePoint::create
-    ((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5, (p5->z + p6->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p64 = adaptivePoint::create
-    ((p6->x + p4->x) * 0.5, (p6->y + p4->y) * 0.5, (p6->z + p4->z) * 0.5, 
-     coeffs, eexps);
-  p->e[0] = new adaptivePrism(p1, p12, p31, p14, p1425, p3614);
-  recurCreate(p->e[0], maxlevel, level, coeffs, eexps);
-  p->e[1] = new adaptivePrism(p2, p23, p12, p25, p2536, p1425);
-  recurCreate(p->e[1], maxlevel, level, coeffs, eexps);
-  p->e[2] = new adaptivePrism(p3, p31, p23, p36, p3614, p2536);
-  recurCreate(p->e[2], maxlevel, level, coeffs, eexps);
-  p->e[3] = new adaptivePrism(p12, p23, p31, p1425, p2536, p3614);
-  recurCreate(p->e[3], maxlevel, level, coeffs, eexps);
-  p->e[4] = new adaptivePrism(p14, p1425, p3614, p4, p45, p64);
-  recurCreate(p->e[4], maxlevel, level, coeffs, eexps);
-  p->e[5] = new adaptivePrism(p25, p2536, p1425, p5, p56, p45);
-  recurCreate(p->e[5], maxlevel, level, coeffs, eexps);
-  p->e[6] = new adaptivePrism(p36, p3614, p2536, p6, p64, p56);
-  recurCreate(p->e[6], maxlevel, level, coeffs, eexps);
-  p->e[7] = new adaptivePrism(p1425, p2536, p3614, p45, p56, p64);
-  recurCreate(p->e[7], maxlevel, level, coeffs, eexps);
-}
-
-void adaptivePrism::error(double AVG, double tol)
-{
-  adaptivePrism *p = *all.begin();
-  recurError(p, AVG, tol);
-}
-
-void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
-{
-  if(!p->e[0])
-    p->visible = true;
-  else {
-    double vr;
-    if(!p->e[0]->e[0]) {
-      double v1 = p->e[0]->V();
-      double v2 = p->e[1]->V();
-      double v3 = p->e[2]->V();
-      double v4 = p->e[3]->V();
-      double v5 = p->e[4]->V();
-      double v6 = p->e[5]->V();
-      double v7 = p->e[6]->V();
-      double v8 = p->e[7]->V();
-      vr = (v1 + v2 + v3 + v4/2 +v5 +v6 +v7 +v8/2) / 7;
-      double v = p->V();
-      if(fabs(v - vr) > AVG * tol){
-        p->visible = false;
-        recurError(p->e[0], AVG, tol);
-        recurError(p->e[1], AVG, tol);
-        recurError(p->e[2], AVG, tol);
-        recurError(p->e[3], AVG, tol);
-        recurError(p->e[4], AVG, tol);
-        recurError(p->e[5], AVG, tol);
-        recurError(p->e[6], AVG, tol);
-        recurError(p->e[7], AVG, tol);
-      }
-      else
-        p->visible = true;
-    }
-    else {
-      bool err=false;
-      double ve[8];
-      for(int i=0; i<8; i++){
-        double v1 = p->e[i]->e[0]->V();
-        double v2 = p->e[i]->e[1]->V();
-        double v3 = p->e[i]->e[2]->V();
-        double v4 = p->e[i]->e[3]->V();
-        double v5 = p->e[i]->e[4]->V();
-        double v6 = p->e[i]->e[5]->V();
-        double v7 = p->e[i]->e[6]->V();
-        double v8 = p->e[i]->e[7]->V();
-        double vr = (v1 + v2 + v3 + v4/2 +v5 +v6 +v7 +v8/2) / 7;
-        ve[i] = p->e[i]->V();
-        err |= (fabs((ve[i] - vr)) > AVG*tol);
-      }
-      double vr = (ve[0] + ve[1] + ve[2] + ve[3] / 2 + 
-                   ve[4] + ve[5] + ve[6] + ve[7] / 2) / 7;
-      err |= (fabs((p->V() - vr))>AVG*tol);
-      if(err) {
-        p->visible = false;
-        for(int i=0;i<8;i++)
-          recurError(p->e[i], AVG, tol);
-      }
-      else
-        p->visible = true;
-    }
-  }
-}
-
-void adaptiveTetrahedron::create(int maxlevel,
-				 Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveTetrahedron::create(int maxlevel)
 {
   cleanElement<adaptiveTetrahedron>();
-  adaptivePoint *p1 = adaptivePoint::create(0, 0, 0, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(0, 1, 0, coeffs, eexps);
-  adaptivePoint *p3 = adaptivePoint::create(1, 0, 0, coeffs, eexps);
-  adaptivePoint *p4 = adaptivePoint::create(0, 0, 1, coeffs, eexps);
+  adaptivePoint *p1 = adaptivePoint::add(0, 0, 0, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(0, 1, 0, allPoints);
+  adaptivePoint *p3 = adaptivePoint::add(1, 0, 0, allPoints);
+  adaptivePoint *p4 = adaptivePoint::add(0, 0, 1, allPoints);
   adaptiveTetrahedron *t = new adaptiveTetrahedron(p1, p2, p3, p4);
-  recurCreate(t, maxlevel, 0, coeffs, eexps);
+  recurCreate(t, maxlevel, 0);
 }
 
-void adaptiveTetrahedron::recurCreate(adaptiveTetrahedron *t, int maxlevel, int level,
-			              Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveTetrahedron::recurCreate(adaptiveTetrahedron *t, int maxlevel, int level)
 {
   all.push_back(t);
   if(level++ >= maxlevel) return;
@@ -557,40 +405,40 @@ void adaptiveTetrahedron::recurCreate(adaptiveTetrahedron *t, int maxlevel, int
   adaptivePoint *p1 = t->p[1];
   adaptivePoint *p2 = t->p[2];
   adaptivePoint *p3 = t->p[3];
-  adaptivePoint *pe0 = adaptivePoint::create
+  adaptivePoint *pe0 = adaptivePoint::add
     ((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5, (p0->z + p1->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *pe1 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pe1 = adaptivePoint::add
     ((p0->x + p2->x) * 0.5, (p0->y + p2->y) * 0.5, (p0->z + p2->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *pe2 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pe2 = adaptivePoint::add
     ((p0->x + p3->x) * 0.5, (p0->y + p3->y) * 0.5, (p0->z + p3->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *pe3 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pe3 = adaptivePoint::add
     ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *pe4 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pe4 = adaptivePoint::add
     ((p1->x + p3->x) * 0.5, (p1->y + p3->y) * 0.5, (p1->z + p3->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *pe5 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pe5 = adaptivePoint::add
     ((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, (p2->z + p3->z) * 0.5,
-     coeffs, eexps);
+     allPoints);
   adaptiveTetrahedron *t1 = new adaptiveTetrahedron(p0, pe0, pe2, pe1);
-  recurCreate(t1, maxlevel, level, coeffs, eexps);
+  recurCreate(t1, maxlevel, level);
   adaptiveTetrahedron *t2 = new adaptiveTetrahedron(p1, pe0, pe3, pe4);
-  recurCreate(t2, maxlevel, level, coeffs, eexps);
+  recurCreate(t2, maxlevel, level);
   adaptiveTetrahedron *t3 = new adaptiveTetrahedron(p2, pe3, pe1, pe5);
-  recurCreate(t3, maxlevel, level, coeffs, eexps);
+  recurCreate(t3, maxlevel, level);
   adaptiveTetrahedron *t4 = new adaptiveTetrahedron(p3, pe2, pe4, pe5);
-  recurCreate(t4, maxlevel, level, coeffs, eexps);
+  recurCreate(t4, maxlevel, level);
   adaptiveTetrahedron *t5 = new adaptiveTetrahedron(pe3, pe5, pe2, pe4);
-  recurCreate(t5, maxlevel, level, coeffs, eexps);
+  recurCreate(t5, maxlevel, level);
   adaptiveTetrahedron *t6 = new adaptiveTetrahedron(pe3, pe2, pe0, pe4);
-  recurCreate(t6, maxlevel, level, coeffs, eexps);
+  recurCreate(t6, maxlevel, level);
   adaptiveTetrahedron *t7 = new adaptiveTetrahedron(pe2, pe5, pe3, pe1);
-  recurCreate(t7, maxlevel, level, coeffs, eexps);
+  recurCreate(t7, maxlevel, level);
   adaptiveTetrahedron *t8 = new adaptiveTetrahedron(pe0, pe2, pe3, pe1);
-  recurCreate(t8, maxlevel, level, coeffs, eexps);
+  recurCreate(t8, maxlevel, level);
   t->e[0] = t1;
   t->e[1] = t2;
   t->e[2] = t3;
@@ -675,24 +523,22 @@ void adaptiveTetrahedron::recurError(adaptiveTetrahedron *t, double AVG, double
   }
 }
 
-void adaptiveHexahedron::create(int maxlevel, 
-				Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveHexahedron::create(int maxlevel)
 {
   cleanElement<adaptiveHexahedron>();
-  adaptivePoint *p1 = adaptivePoint::create(-1, -1, -1, coeffs, eexps);
-  adaptivePoint *p2 = adaptivePoint::create(-1, 1, -1, coeffs, eexps);
-  adaptivePoint *p3 = adaptivePoint::create(1, 1, -1, coeffs, eexps);
-  adaptivePoint *p4 = adaptivePoint::create(1, -1, -1, coeffs, eexps);
-  adaptivePoint *p11 = adaptivePoint::create(-1, -1, 1, coeffs, eexps);
-  adaptivePoint *p21 = adaptivePoint::create(-1, 1, 1, coeffs, eexps);
-  adaptivePoint *p31 = adaptivePoint::create(1, 1, 1, coeffs, eexps);
-  adaptivePoint *p41 = adaptivePoint::create(1, -1, 1, coeffs, eexps);
+  adaptivePoint *p1 = adaptivePoint::add(-1, -1, -1, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(-1, 1, -1, allPoints);
+  adaptivePoint *p3 = adaptivePoint::add(1, 1, -1, allPoints);
+  adaptivePoint *p4 = adaptivePoint::add(1, -1, -1, allPoints);
+  adaptivePoint *p11 = adaptivePoint::add(-1, -1, 1, allPoints);
+  adaptivePoint *p21 = adaptivePoint::add(-1, 1, 1, allPoints);
+  adaptivePoint *p31 = adaptivePoint::add(1, 1, 1, allPoints);
+  adaptivePoint *p41 = adaptivePoint::add(1, -1, 1, allPoints);
   adaptiveHexahedron *h = new adaptiveHexahedron(p1, p2, p3, p4, p11, p21, p31, p41);
-  recurCreate(h, maxlevel, 0, coeffs, eexps);
+  recurCreate(h, maxlevel, 0);
 }
 
-void adaptiveHexahedron::recurCreate(adaptiveHexahedron *h, int maxlevel, int level,
-				     Double_Matrix *coeffs, Double_Matrix *eexps)
+void adaptiveHexahedron::recurCreate(adaptiveHexahedron *h, int maxlevel, int level)
 {
   all.push_back(h);
   if(level++ >= maxlevel) return;
@@ -705,90 +551,90 @@ void adaptiveHexahedron::recurCreate(adaptiveHexahedron *h, int maxlevel, int le
   adaptivePoint *p5 = h->p[5];
   adaptivePoint *p6 = h->p[6];
   adaptivePoint *p7 = h->p[7];
-  adaptivePoint *p01 = adaptivePoint::create
+  adaptivePoint *p01 = adaptivePoint::add
     ((p0->x + p1->x) * 0.5, (p0->y + p1->y) * 0.5, (p0->z + p1->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p12 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p12 = adaptivePoint::add
     ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p23 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p23 = adaptivePoint::add
     ((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, (p2->z + p3->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p03 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p03 = adaptivePoint::add
     ((p3->x + p0->x) * 0.5, (p3->y + p0->y) * 0.5, (p3->z + p0->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p45 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p45 = adaptivePoint::add
     ((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5, (p4->z + p5->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p56 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p56 = adaptivePoint::add
     ((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5, (p5->z + p6->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p67 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p67 = adaptivePoint::add
     ((p6->x + p7->x) * 0.5, (p6->y + p7->y) * 0.5, (p6->z + p7->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p47 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p47 = adaptivePoint::add
     ((p7->x + p4->x) * 0.5, (p7->y + p4->y) * 0.5, (p7->z + p4->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p04 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p04 = adaptivePoint::add
     ((p4->x + p0->x) * 0.5, (p4->y + p0->y) * 0.5, (p4->z + p0->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p15 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p15 = adaptivePoint::add
     ((p5->x + p1->x) * 0.5, (p5->y + p1->y) * 0.5, (p5->z + p1->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p26 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p26 = adaptivePoint::add
     ((p6->x + p2->x) * 0.5, (p6->y + p2->y) * 0.5, (p6->z + p2->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p37 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p37 = adaptivePoint::add
     ((p7->x + p3->x) * 0.5, (p7->y + p3->y) * 0.5, (p7->z + p3->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p0145 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p0145 = adaptivePoint::add
     ((p45->x + p01->x) * 0.5, (p45->y + p01->y) * 0.5,(p45->z + p01->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p1256 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p1256 = adaptivePoint::add
     ((p12->x + p56->x) * 0.5, (p12->y + p56->y) * 0.5, (p12->z + p56->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p2367 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p2367 = adaptivePoint::add
     ((p23->x + p67->x) * 0.5, (p23->y + p67->y) * 0.5, (p23->z + p67->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p0347 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p0347 = adaptivePoint::add
     ((p03->x + p47->x) * 0.5, (p03->y + p47->y) * 0.5, (p03->z + p47->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *p4756 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p4756 = adaptivePoint::add
     ((p47->x + p56->x) * 0.5, (p47->y + p56->y) * 0.5, (p47->z + p56->z) * 0.5, 
-     coeffs, eexps);
-  adaptivePoint *p0312 = adaptivePoint::create
+     allPoints);
+  adaptivePoint *p0312 = adaptivePoint::add
     ((p03->x + p12->x) * 0.5, (p03->y + p12->y) * 0.5, (p03->z + p12->z) * 0.5,
-     coeffs, eexps);
-  adaptivePoint *pc = adaptivePoint::create
+     allPoints);
+  adaptivePoint *pc = adaptivePoint::add
     ((p0->x + p1->x + p2->x + p3->x + p4->x + p5->x + p6->x + p7->x) * 0.125,
      (p0->y + p1->y + p2->y + p3->y + p4->y + p5->y + p6->y + p7->y) * 0.125,
      (p0->z + p1->z + p2->z + p3->z + p4->z + p5->z + p6->z + p7->z) * 0.125,
-     coeffs, eexps);
+     allPoints);
 
   adaptiveHexahedron *h1 = new adaptiveHexahedron
     (p0, p01, p0312, p03, p04, p0145, pc, p0347); // p0
-  recurCreate(h1, maxlevel, level, coeffs, eexps);
+  recurCreate(h1, maxlevel, level);
   adaptiveHexahedron *h2 = new adaptiveHexahedron
     (p01, p0145, p15, p1, p0312, pc, p1256, p12); // p1
-  recurCreate(h2, maxlevel, level, coeffs, eexps);
+  recurCreate(h2, maxlevel, level);
   adaptiveHexahedron *h3 = new adaptiveHexahedron
     (p04, p4, p45, p0145, p0347, p47, p4756, pc); // p4
-  recurCreate(h3, maxlevel, level, coeffs, eexps);
+  recurCreate(h3, maxlevel, level);
   adaptiveHexahedron *h4 = new adaptiveHexahedron
     (p0145, p45, p5, p15, pc, p4756, p56, p1256); // p5
-  recurCreate(h4, maxlevel, level, coeffs, eexps);
+  recurCreate(h4, maxlevel, level);
   adaptiveHexahedron *h5 = new adaptiveHexahedron
     (p0347, p47, p4756, pc, p37, p7, p67, p2367); // p7
-  recurCreate(h5, maxlevel, level, coeffs, eexps);
+  recurCreate(h5, maxlevel, level);
   adaptiveHexahedron *h6 = new adaptiveHexahedron
     (pc, p4756, p56, p1256, p2367, p67, p6, p26); // p6
-  recurCreate(h6, maxlevel, level, coeffs, eexps);
+  recurCreate(h6, maxlevel, level);
   adaptiveHexahedron *h7 = new adaptiveHexahedron
     (p03, p0347, pc, p0312, p3, p37, p2367, p23); // p3
-  recurCreate(h7, maxlevel, level, coeffs, eexps);
+  recurCreate(h7, maxlevel, level);
   adaptiveHexahedron *h8 = new adaptiveHexahedron
     (p0312, pc, p1256, p12, p23, p2367, p26, p2); //p2
-  recurCreate(h8, maxlevel, level, coeffs, eexps);
+  recurCreate(h8, maxlevel, level);
   h->e[0] = h1;
   h->e[1] = h2;
   h->e[2] = h3;
@@ -837,16 +683,164 @@ void adaptiveHexahedron::recurError(adaptiveHexahedron *h, double AVG, double to
   }
 }
 
+void adaptivePrism::create(int maxlevel)
+{
+  cleanElement<adaptivePrism>();
+  adaptivePoint *p1 = adaptivePoint::add(0, 0, -1, allPoints);
+  adaptivePoint *p2 = adaptivePoint::add(1, 0, -1, allPoints);
+  adaptivePoint *p3 = adaptivePoint::add(0, 1, -1, allPoints);
+  adaptivePoint *p4 = adaptivePoint::add(0, 0, 1, allPoints);
+  adaptivePoint *p5 = adaptivePoint::add(1, 0, 1, allPoints);
+  adaptivePoint *p6 = adaptivePoint::add(0, 1, 1, allPoints);
+  adaptivePrism *p = new adaptivePrism(p1, p2, p3, p4, p5, p6);
+  recurCreate(p, maxlevel, 0);
+}
+
+void adaptivePrism::recurCreate(adaptivePrism *p, int maxlevel, int level)
+{
+  all.push_back(p);
+  if(level++ >= maxlevel) return;
+
+  // p4   p34    p3
+  // p14  pc     p23
+  // p1   p12    p2
+  adaptivePoint *p1 = p->p[0];
+  adaptivePoint *p2 = p->p[1];
+  adaptivePoint *p3 = p->p[2];
+  adaptivePoint *p4 = p->p[3];
+  adaptivePoint *p5 = p->p[4];
+  adaptivePoint *p6 = p->p[5];
+  adaptivePoint *p14 = adaptivePoint::add
+    ((p1->x + p4->x) * 0.5, (p1->y + p4->y) * 0.5, (p1->z + p4->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p25 = adaptivePoint::add
+    ((p2->x + p5->x) * 0.5, (p2->y + p5->y) * 0.5, (p2->z + p5->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p36 = adaptivePoint::add
+    ((p3->x + p6->x) * 0.5, (p3->y + p6->y) * 0.5, (p3->z + p6->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p12 = adaptivePoint::add
+    ((p1->x + p2->x) * 0.5, (p1->y + p2->y) * 0.5, (p1->z + p2->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p23 = adaptivePoint::add
+    ((p2->x + p3->x) * 0.5, (p2->y + p3->y) * 0.5, (p2->z + p3->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p31 = adaptivePoint::add
+    ((p3->x + p1->x) * 0.5, (p3->y + p1->y) * 0.5, (p3->z + p1->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p1425 = adaptivePoint::add
+    ((p14->x + p25->x) * 0.5, (p14->y + p25->y) * 0.5, (p14->z + p25->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p2536 = adaptivePoint::add
+    ((p25->x + p36->x) * 0.5, (p25->y + p36->y) * 0.5, (p25->z + p36->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p3614 = adaptivePoint::add
+    ((p36->x + p14->x) * 0.5, (p36->y + p14->y) * 0.5, (p36->z + p14->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p45 = adaptivePoint::add
+    ((p4->x + p5->x) * 0.5, (p4->y + p5->y) * 0.5, (p4->z + p5->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p56 = adaptivePoint::add
+    ((p5->x + p6->x) * 0.5, (p5->y + p6->y) * 0.5, (p5->z + p6->z) * 0.5, 
+     allPoints);
+  adaptivePoint *p64 = adaptivePoint::add
+    ((p6->x + p4->x) * 0.5, (p6->y + p4->y) * 0.5, (p6->z + p4->z) * 0.5, 
+     allPoints);
+  p->e[0] = new adaptivePrism(p1, p12, p31, p14, p1425, p3614);
+  recurCreate(p->e[0], maxlevel, level);
+  p->e[1] = new adaptivePrism(p2, p23, p12, p25, p2536, p1425);
+  recurCreate(p->e[1], maxlevel, level);
+  p->e[2] = new adaptivePrism(p3, p31, p23, p36, p3614, p2536);
+  recurCreate(p->e[2], maxlevel, level);
+  p->e[3] = new adaptivePrism(p12, p23, p31, p1425, p2536, p3614);
+  recurCreate(p->e[3], maxlevel, level);
+  p->e[4] = new adaptivePrism(p14, p1425, p3614, p4, p45, p64);
+  recurCreate(p->e[4], maxlevel, level);
+  p->e[5] = new adaptivePrism(p25, p2536, p1425, p5, p56, p45);
+  recurCreate(p->e[5], maxlevel, level);
+  p->e[6] = new adaptivePrism(p36, p3614, p2536, p6, p64, p56);
+  recurCreate(p->e[6], maxlevel, level);
+  p->e[7] = new adaptivePrism(p1425, p2536, p3614, p45, p56, p64);
+  recurCreate(p->e[7], maxlevel, level);
+}
+
+void adaptivePrism::error(double AVG, double tol)
+{
+  adaptivePrism *p = *all.begin();
+  recurError(p, AVG, tol);
+}
+
+void adaptivePrism::recurError(adaptivePrism *p, double AVG, double tol)
+{
+  if(!p->e[0])
+    p->visible = true;
+  else {
+    double vr;
+    if(!p->e[0]->e[0]) {
+      double v1 = p->e[0]->V();
+      double v2 = p->e[1]->V();
+      double v3 = p->e[2]->V();
+      double v4 = p->e[3]->V();
+      double v5 = p->e[4]->V();
+      double v6 = p->e[5]->V();
+      double v7 = p->e[6]->V();
+      double v8 = p->e[7]->V();
+      vr = (v1 + v2 + v3 + v4/2 +v5 +v6 +v7 +v8/2) / 7;
+      double v = p->V();
+      if(fabs(v - vr) > AVG * tol){
+        p->visible = false;
+        recurError(p->e[0], AVG, tol);
+        recurError(p->e[1], AVG, tol);
+        recurError(p->e[2], AVG, tol);
+        recurError(p->e[3], AVG, tol);
+        recurError(p->e[4], AVG, tol);
+        recurError(p->e[5], AVG, tol);
+        recurError(p->e[6], AVG, tol);
+        recurError(p->e[7], AVG, tol);
+      }
+      else
+        p->visible = true;
+    }
+    else {
+      bool err=false;
+      double ve[8];
+      for(int i=0; i<8; i++){
+        double v1 = p->e[i]->e[0]->V();
+        double v2 = p->e[i]->e[1]->V();
+        double v3 = p->e[i]->e[2]->V();
+        double v4 = p->e[i]->e[3]->V();
+        double v5 = p->e[i]->e[4]->V();
+        double v6 = p->e[i]->e[5]->V();
+        double v7 = p->e[i]->e[6]->V();
+        double v8 = p->e[i]->e[7]->V();
+        double vr = (v1 + v2 + v3 + v4/2 +v5 +v6 +v7 +v8/2) / 7;
+        ve[i] = p->e[i]->V();
+        err |= (fabs((ve[i] - vr)) > AVG*tol);
+      }
+      double vr = (ve[0] + ve[1] + ve[2] + ve[3] / 2 + 
+                   ve[4] + ve[5] + ve[6] + ve[7] / 2) / 7;
+      err |= (fabs((p->V() - vr))>AVG*tol);
+      if(err) {
+        p->visible = false;
+        for(int i=0;i<8;i++)
+          recurError(p->e[i], AVG, tol);
+      }
+      else
+        p->visible = true;
+    }
+  }
+}
+
 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());
-  Double_Vector valz(N), resz(adaptivePoint::all.size());
+  Double_Vector val(N),  res(T::allPoints.size());
+  Double_Vector valx(N), resx(T::allPoints.size());
+  Double_Vector valy(N), resy(T::allPoints.size());
+  Double_Vector valz(N), resz(T::allPoints.size());
   Double_Matrix xyz(_posX->size2(), 3);
-  Double_Matrix XYZ(adaptivePoint::all.size(), 3);
+  Double_Matrix XYZ(T::allPoints.size(), 3);
 
   for(int k = 0; k < _posX->size2(); ++k){
     xyz(k, 0) = (*_posX)(ielem, k);
@@ -873,8 +867,8 @@ int adaptiveElements<T>::_zoomElement(int ielem, int level, GMSH_Post_Plugin *pl
   _geometry->mult(xyz, XYZ);
 
   int k = 0;
-  for(std::set<adaptivePoint>::iterator it = adaptivePoint::all.begin();
-      it != adaptivePoint::all.end(); ++it){
+  for(std::set<adaptivePoint>::iterator it = T::allPoints.begin();
+      it != T::allPoints.end(); ++it){
     adaptivePoint *p = (adaptivePoint*)&(*it);
     p->val = res(k);
     if(_valX){
@@ -933,30 +927,34 @@ 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)
 {
-  T::create(level, _coeffs, _eexps);
+  T::create(level);
 
-  if(!adaptivePoint::all.size()) return;
+  if(!T::allPoints.size()) return;
 
   const int N = _coeffs->size1();
   if(_interpolate) delete _interpolate;
-  _interpolate = new Double_Matrix(adaptivePoint::all.size(), N);
+  _interpolate = new Double_Matrix(T::allPoints.size(), N);
 
   if(_geometry) delete _geometry;
-  _geometry = new Double_Matrix(adaptivePoint::all.size(), _posX->size2());
+  _geometry = new Double_Matrix(T::allPoints.size(), _posX->size2());
 
   double sf[100];
   int kk = 0;
-  for(std::set<adaptivePoint>::iterator it = adaptivePoint::all.begin(); 
-      it != adaptivePoint::all.end(); ++it) {
+  for(std::set<adaptivePoint>::iterator it = T::allPoints.begin(); 
+      it != T::allPoints.end(); ++it) {
     adaptivePoint *p = (adaptivePoint*)&(*it);
+
+    computeShapeFunctions(_coeffs, _eexps, p->x, p->y, p->z, sf);
     for(int k = 0; k < N; ++k)
-      (*_interpolate)(kk, k) = p->shapeFunctions[k];
+      (*_interpolate)(kk, k) = sf[k];
+
     if(_coeffsGeom)
       computeShapeFunctions(_coeffsGeom, _eexpsGeom, p->x, p->y, p->z, sf);
     else
       T::GSF(p->x, p->y, p->z, sf);
     for(int k = 0; k < _posX->size2(); k++)
       (*_geometry)(kk, k) = sf[k];
+
     kk++;
   }
 
@@ -997,7 +995,7 @@ void adaptiveElements<T>::initData(PViewData *data, int step)
   }
   if(!numEle) return;
 
-  int numNodes = getNumNodes();
+  int numNodes = _coeffsGeom ? _coeffsGeom->size1() : T::numNodes;
   int numVal = _coeffs->size1() * numComp;
   if(!numNodes || !numVal) return;
 
@@ -1083,9 +1081,8 @@ void adaptiveElements<T>::changeResolution(int level, double tol, GMSH_Post_Plug
 
   List_Reset(_listEle);
   *_numEle = 0;
-  
-  std::vector<int> done(_posX->size1(), 0);
 
+  std::vector<int> done(_posX->size1(), 0);
   // We first do the adaptive stuff at level 2 and will only process
   // elements that have reached the maximal recursion level
   int level_act = (level > 2) ? 2 : level;
@@ -1097,6 +1094,7 @@ void adaptiveElements<T>::changeResolution(int level, double tol, GMSH_Post_Plug
     if(level_act >= level) break;
     level_act++;
   }
+  //_changeResolution(level, plug, &done[0]);
 }
 
 adaptiveData::adaptiveData(PViewData *data)
@@ -1188,3 +1186,4 @@ void adaptiveData::changeResolution(int step, int level, double tol, GMSH_Post_P
   _level = level;
   _tol = tol;
 }
+
diff --git a/Post/adaptiveData.h b/Post/adaptiveData.h
index 90fc4cc4a08d0f766036d137eb6c4ca56f6479fb..200339a221b5c2a71afc282b78881c03443d4e77 100644
--- a/Post/adaptiveData.h
+++ b/Post/adaptiveData.h
@@ -19,11 +19,9 @@ class adaptivePoint {
  public:
   double x, y, z, X, Y, Z;
   double val, valx, valy, valz;
-  double shapeFunctions[128];
-  static std::set<adaptivePoint> all;
  public:
-  static adaptivePoint *create(double x, double y, double z, 
-			       Double_Matrix *coeffs, Double_Matrix *eexps);
+  static adaptivePoint *add(double x, double y, double z, 
+                            std::set<adaptivePoint> &allPoints);
   bool operator < (const adaptivePoint &other) const
   {
     if(other.x < x) return true;
@@ -41,6 +39,7 @@ class adaptiveLine {
   adaptivePoint *p[2];
   adaptiveLine *e[2];
   static std::list<adaptiveLine*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptiveLine(adaptivePoint *p1, adaptivePoint *p2)
@@ -59,9 +58,8 @@ class adaptiveLine {
     sf[0] = (1 - u) / 2.;
     sf[1] = (1 + u) / 2.;
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptiveLine *e, int maxlevel, int level, 
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptiveLine *e, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptiveLine *e, double AVG, double tol);
 };
@@ -72,6 +70,7 @@ class adaptiveTriangle {
   adaptivePoint *p[3];
   adaptiveTriangle *e[4];
   static std::list<adaptiveTriangle*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptiveTriangle(adaptivePoint *p1, adaptivePoint *p2, adaptivePoint *p3)
@@ -92,9 +91,8 @@ class adaptiveTriangle {
     sf[1] = u;
     sf[2] = v;
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptiveTriangle *t, int maxlevel, int level,
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptiveTriangle *t, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptiveTriangle *t, double AVG, double tol);
 };
@@ -105,6 +103,7 @@ class adaptiveQuadrangle {
   adaptivePoint *p[4];
   adaptiveQuadrangle *e[4];
   static std::list<adaptiveQuadrangle*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptiveQuadrangle(adaptivePoint *p1, adaptivePoint *p2, 
@@ -128,9 +127,8 @@ class adaptiveQuadrangle {
     sf[2] = 0.25 * (1. + u) * (1. + v);
     sf[3] = 0.25 * (1. - u) * (1. + v);
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptiveQuadrangle *q, int maxlevel, int level,
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptiveQuadrangle *q, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptiveQuadrangle *q, double AVG, double tol);
 };
@@ -141,6 +139,7 @@ class adaptivePrism {
   adaptivePoint *p[6];
   adaptivePrism *e[12];
   static std::list<adaptivePrism*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptivePrism(adaptivePoint *p1, adaptivePoint *p2, 
@@ -171,9 +170,8 @@ class adaptivePrism {
     sf[4] = u*(1+w)/2;
     sf[5] = v*(1+w)/2;
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptivePrism *p, int maxlevel, int level, 
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptivePrism *p, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptivePrism *p, double AVG, double tol);
 };
@@ -184,6 +182,7 @@ class adaptiveTetrahedron {
   adaptivePoint *p[4];
   adaptiveTetrahedron *e[8];
   static std::list<adaptiveTetrahedron*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptiveTetrahedron(adaptivePoint *p1, adaptivePoint *p2, 
@@ -208,9 +207,8 @@ class adaptiveTetrahedron {
     sf[2] = v;
     sf[3] = w;
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptiveTetrahedron *t, int maxlevel, int level, 
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptiveTetrahedron *t, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptiveTetrahedron *t, double AVG, double tol);
 };
@@ -221,6 +219,7 @@ class adaptiveHexahedron {
   adaptivePoint *p[8];
   adaptiveHexahedron *e[8];
   static std::list<adaptiveHexahedron*> all;
+  static std::set<adaptivePoint> allPoints;
   static int numNodes, numEdges;
  public:
   adaptiveHexahedron(adaptivePoint *p1, adaptivePoint *p2, adaptivePoint *p3, 
@@ -255,9 +254,8 @@ class adaptiveHexahedron {
     sf[6] = 0.125 * (1 + u) * (1 + v) * (1 + w);
     sf[7] = 0.125 * (1 - u) * (1 + v) * (1 + w);
   }
-  static void create(int maxlevel, Double_Matrix *coeffs, Double_Matrix *eexps);
-  static void recurCreate(adaptiveHexahedron *h, int maxlevel, int level,
-			  Double_Matrix *coeffs, Double_Matrix *eexps);
+  static void create(int maxlevel);
+  static void recurCreate(adaptiveHexahedron *h, int maxlevel, int level);
   static void error(double AVG, double tol);
   static void recurError(adaptiveHexahedron *h, double AVG, double tol);
 };
diff --git a/demos/primitives.pos b/demos/primitives.pos
index 0c7547de6a413941b7d7a8fc70b9f01b3e0081b3..2a291b6c20c27def267aacae70f0847e0780e15d 100644
--- a/demos/primitives.pos
+++ b/demos/primitives.pos
@@ -67,7 +67,6 @@ View "Primitives"{
 };
 
 Alias View[0];
-View[1].Name = "Primitives (dupl. 1)" ;
 View[1].OffsetX = 3 ;
 View[1].IntervalsType = 2 ;
 View[1].ShowElement = 1 ;
@@ -75,7 +74,6 @@ View[1].GlyphLocation = 2 ;
 View[1].ColorTable = {Red,Green,Magenta,Cyan,Brown,Pink} ;
 
 Alias View[0];
-View[2].Name = "Primitives (dupl. 2)" ;
 View[2].OffsetX = 6 ;
 View[2].IntervalsType = 4 ;
 View[2].ShowElement = 1 ;