diff --git a/Common/gmshpy.i b/Common/gmshpy.i
index bcd1e06c1c551f66c4947ce701d35a251a96b27b..b924241e2af669343902a72f0cb7ff687bc720c3 100644
--- a/Common/gmshpy.i
+++ b/Common/gmshpy.i
@@ -59,12 +59,16 @@ namespace std {
    %template(IntVector) vector<int>;
    %template(DoubleVector) vector<double, std::allocator<double> >;
    %template(StringVector) vector<std::string, std::allocator<std::string> >;
+   %template(GEntityVector) vector<GEntity*, std::allocator<GEntity*> >;
    %template(GVertexVector) vector<GVertex*, std::allocator<GVertex*> >;
    %template(GEdgeVector) vector<GEdge*, std::allocator<GEdge*> >;
    %template(GFaceVector) vector<GFace*, std::allocator<GFace*> >;
    %template(GRegionVector) vector<GRegion*, std::allocator<GRegion*> >;
+   %template(MVertexVector) vector< MVertex *,std::allocator< MVertex * > >;
+   %template(MElementVector) vector< MElement *,std::allocator< MElement * > >;
    %template(VectorFunctionConst) vector<const function*, std::allocator<const function*> >;
    %template(GEdgeVectorVector) vector< std::vector< GEdge *,std::allocator< GEdge * > >,std::allocator< std::vector< GEdge *,std::allocator< GEdge * > > > >;
+   %template(GFaceList) list<GFace*, std::allocator<GFace*> >;
 }
 
 %include "fullMatrix.h"
diff --git a/Geo/MFace.h b/Geo/MFace.h
index 665484ec28828085befec06c6d4b8636dfa30e76..01cbfd46225475d1c4aab32c7296b473ce98281f 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -22,7 +22,7 @@ class MFace {
  public:
   MFace() {}
   MFace(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3=0);
-  MFace(std::vector<MVertex*> v);
+  MFace(const std::vector<MVertex*> v);
   inline int getNumVertices() const { return _v.size(); }
   inline MVertex *getVertex(const int i) const { return _v[i]; }
   inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index ed982b52df54032138dad746faa7acf007d15a3b..efb616abd9a203ab4efeac1a16d9b28152a09ac3 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -48,7 +48,7 @@ class MHexahedron : public MElement {
     _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
     _v[4] = v4; _v[5] = v5; _v[6] = v6; _v[7] = v7;
   }
-  MHexahedron(std::vector<MVertex*> &v, int num=0, int part=0)
+  MHexahedron(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 8; i++) _v[i] = v[i];
@@ -239,7 +239,7 @@ class MHexahedron20 : public MHexahedron {
     _vs[10] = v18; _vs[11] = v19;
     for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MHexahedron20(std::vector<MVertex*> &v, int num=0, int part=0)
+  MHexahedron20(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MHexahedron(v, num, part)
   {
     for(int i = 0; i < 12; i++) _vs[i] = v[8 + i];
@@ -385,7 +385,7 @@ class MHexahedron27 : public MHexahedron {
     _vs[15] = v23; _vs[16] = v24; _vs[17] = v25; _vs[18] = v26;
     for(int i = 0; i < 19; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MHexahedron27(std::vector<MVertex*> &v, int num=0, int part=0)
+  MHexahedron27(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MHexahedron(v, num, part)
   {
     for(int i = 0; i < 19; i++) _vs[i] = v[8 + i];
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 19319f8939e55bb9ffa53f41c43119288873aebf..d9dfa89a420db875103a87e370a21e0f9de9b67b 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -28,7 +28,7 @@ class MLine : public MElement {
   {
     _v[0] = v0; _v[1] = v1;
   }
-  MLine(std::vector<MVertex*> &v, int num=0, int part=0)
+  MLine(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 2; i++) _v[i] = v[i];
@@ -98,7 +98,7 @@ class MLine3 : public MLine {
     _vs[0] = v2;
     _vs[0]->setPolynomialOrder(2);
   }
-  MLine3(std::vector<MVertex*> &v, int num=0, int part=0)
+  MLine3(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MLine(v, num, part)
   {
     _vs[0] = v[2];
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index dae3bcb0b45be025205621d8702bb5e12d06891f..e81eaa95fb9eeb20f0c1bdee51542430dc1e3706 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -21,7 +21,7 @@ class MPoint : public MElement {
   {
     _v[0] = v0;
   }
-  MPoint(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPoint(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     _v[0] = v[0];
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index c16a81efe284ad3f2d4151e34c14e29c4cf8688e..3ff52cb36fad87ab433e3be1aea1f5eb6de612bf 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -55,7 +55,7 @@ class MPrism : public MElement {
     _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
     _v[4] = v4; _v[5] = v5;
   }
-  MPrism(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPrism(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 6; i++) _v[i] = v[i];
@@ -239,7 +239,7 @@ class MPrism15 : public MPrism {
     _vs[5] = v11; _vs[6] = v12; _vs[7] = v13; _vs[8] = v14;
     for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MPrism15(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPrism15(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MPrism(v, num, part)
   {
     for(int i = 0; i < 9; i++) _vs[i] = v[6 + i];
@@ -367,7 +367,7 @@ class MPrism18 : public MPrism {
     _vs[10] = v16; _vs[11] = v17;
     for(int i = 0; i < 12; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MPrism18(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPrism18(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MPrism(v, num, part)
   {
     for(int i = 0; i < 12; i++) _vs[i] = v[6 + i];
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index a17859360f6834b7800d595896d25708e0c1a03a..b955baad16b175448bbe8e1948f83539e7590fdb 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -58,7 +58,7 @@ class MPyramid : public MElement {
   {
     _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3; _v[4] = v4;
   }
-  MPyramid(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPyramid(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 5; i++) _v[i] = v[i];
@@ -238,7 +238,7 @@ class MPyramid13 : public MPyramid {
     _vs[5] = v10; _vs[6] = v11; _vs[7] = v12;
     for(int i = 0; i < 8; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MPyramid13(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPyramid13(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MPyramid(v, num, part)
   {
     for(int i = 0; i < 8; i++) _vs[i] = v[5 + i];
@@ -355,7 +355,7 @@ class MPyramid14 : public MPyramid {
     _vs[5] = v10; _vs[6] = v11; _vs[7] = v12; _vs[8] = v13;
     for(int i = 0; i < 9; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MPyramid14(std::vector<MVertex*> &v, int num=0, int part=0)
+  MPyramid14(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MPyramid(v, num, part)
   {
     for(int i = 0; i < 9; i++) _vs[i] = v[5 + i];
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 0ec0310863575b65d7fcc50308a4a8083d7f78a6..412224a8ff2043f9c0e3e21d7e7b26db358c9a41 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -45,7 +45,7 @@ class MQuadrangle : public MElement {
   {
     _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
   }
-  MQuadrangle(std::vector<MVertex*> &v, int num=0, int part=0)
+  MQuadrangle(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 4; i++) _v[i] = v[i];
@@ -180,7 +180,7 @@ class MQuadrangle8 : public MQuadrangle {
     _vs[0] = v4; _vs[1] = v5; _vs[2] = v6; _vs[3] = v7;
     for(int i = 0; i < 4; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MQuadrangle8(std::vector<MVertex*> &v, int num=0, int part=0)
+  MQuadrangle8(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MQuadrangle(v, num, part)
   {
     for(int i = 0; i < 4; i++) _vs[i] = v[4 + i];
@@ -263,7 +263,7 @@ class MQuadrangle9 : public MQuadrangle {
     _vs[0] = v4; _vs[1] = v5; _vs[2] = v6; _vs[3] = v7; _vs[4] = v8;
     for(int i = 0; i < 5; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MQuadrangle9(std::vector<MVertex*> &v, int num=0, int part=0)
+  MQuadrangle9(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MQuadrangle(v, num, part)
   {
     for(int i = 0; i < 5; i++) _vs[i] = v[4 + i];
@@ -334,12 +334,12 @@ class MQuadrangleN : public MQuadrangle {
   const char _order;
  public:
   MQuadrangleN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
-             std::vector<MVertex*> &v, char order, int num=0, int part=0)
+             const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MQuadrangle(v0, v1, v2, v3, num, part), _vs(v), _order(order)
   {
     for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
   }
-  MQuadrangleN(std::vector<MVertex*> &v, char order, int num=0, int part=0)
+  MQuadrangleN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MQuadrangle(v[0], v[1], v[2], v[3], num, part), _order(order)
   {
     for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index b7b12890105bec317de403e341eae27fc3fc26b0..10c6f8df36dc5b8f9c8d05d56a1bd603565b865d 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -51,7 +51,7 @@ class MTetrahedron : public MElement {
   {
     _v[0] = v0; _v[1] = v1; _v[2] = v2; _v[3] = v3;
   }
-  MTetrahedron(std::vector<MVertex*> &v, int num=0, int part=0)
+  MTetrahedron(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 4; i++) _v[i] = v[i];
@@ -198,7 +198,7 @@ class MTetrahedron10 : public MTetrahedron {
     _vs[0] = v4; _vs[1] = v5; _vs[2] = v6; _vs[3] = v7; _vs[4] = v8; _vs[5] = v9;
     for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MTetrahedron10(std::vector<MVertex*> &v, int num=0, int part=0)
+  MTetrahedron10(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MTetrahedron(v, num, part)
   {
     for(int i = 0; i < 6; i++) _vs[i] = v[4 + i];
@@ -314,12 +314,12 @@ class MTetrahedronN : public MTetrahedron {
   double _disto;
  public:
   MTetrahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
-                std::vector<MVertex*> &v, char order, int num=0, int part=0)
+                const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order),_disto(-1.e22)
   {
     for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
   }
-  MTetrahedronN(std::vector<MVertex*> &v, char order, int num=0, int part=0)
+  MTetrahedronN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order),_disto(-1.e22)
   {
     for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 3cb6b223310b2d1a2364bc18dc50335b258b4a7f..87150f25bc09891f0c6e82bd85f38c37cc9685d5 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -43,7 +43,7 @@ class MTriangle : public MElement {
   {
     _v[0] = v0; _v[1] = v1; _v[2] = v2;
   }
-  MTriangle(std::vector<MVertex*> &v, int num=0, int part=0)
+  MTriangle(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MElement(num, part)
   {
     for(int i = 0; i < 3; i++) _v[i] = v[i];
@@ -174,7 +174,7 @@ class MTriangle6 : public MTriangle {
     _vs[0] = v3; _vs[1] = v4; _vs[2] = v5;
     for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
   }
-  MTriangle6(std::vector<MVertex*> &v, int num=0, int part=0)
+  MTriangle6(const std::vector<MVertex*> &v, int num=0, int part=0)
     : MTriangle(v, num, part)
   {
     for(int i = 0; i < 3; i++) _vs[i] = v[3 + i];
@@ -250,12 +250,12 @@ class MTriangleN : public MTriangle {
   const char _order;
  public:
   MTriangleN(MVertex *v0, MVertex *v1, MVertex *v2,
-             std::vector<MVertex*> &v, char order, int num=0, int part=0)
+             const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MTriangle(v0, v1, v2, num, part), _vs(v), _order(order)
   {
     for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
   }
-  MTriangleN(std::vector<MVertex*> &v, char order, int num=0, int part=0)
+  MTriangleN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
     : MTriangle(v[0], v[1], v[2], num, part), _order(order)
   {
     for(unsigned int i = 3; i < v.size(); i++) _vs.push_back(v[i]);