diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 5091bc18cc74e96bf111e436de4e25720a975e2d..3ffe7f1e106e9efdb02c53b42d3a0f6c0740cf57 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -65,3 +65,8 @@ double MLine::getLength()
 {
   return _v[0]->distance(_v[1]);
 }
+
+double MLine::getVolume()
+{
+  return getLength();
+}
diff --git a/Geo/MLine.h b/Geo/MLine.h
index a26fbaf26846034eba562dfd8e27b706abb22498..8e60f50a6af679f62974b488c7a1ace7892755b5 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -39,6 +39,7 @@ class MLine : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual double getInnerRadius(); // half-length of segment line
   virtual double getLength(); // length of segment line
+  virtual double getVolume();
   virtual void getVertexInfo(const MVertex * vertex, int &ithVertex) const
   {
     ithVertex = _v[0] == vertex ? 0 : 1;
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index c448fabab4ba94d2d3cb60acf486d6c7e737f0ae..6f41fe2dc36669ea45aae70c8ccfc090f4bd4626 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -94,7 +94,20 @@ int MQuadrangleN::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges
 int MQuadrangle8::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
 int MQuadrangle9::getNumEdgesRep(){ return 4 * CTX::instance()->mesh.numSubEdges; }
 
-
+double MQuadrangle::getVolume()
+{
+  if(getNumVertices() > 4)
+    return MElement::getVolume();
+  double a = _v[0]->distance(_v[1]);
+  double b = _v[1]->distance(_v[2]);
+  double c = _v[2]->distance(_v[3]);
+  double d = _v[3]->distance(_v[0]);
+  double m = _v[0]->distance(_v[2]);
+  double n = _v[1]->distance(_v[3]);
+  double mn = 2. * m * n;
+  double abcd = a*a - b*b + c*c - d*d;
+  return sqrt( mn*mn - abcd*abcd ) / 4.;
+}
 
 static void _myGetEdgeRep(MQuadrangle *q, int num, double *x, double *y, double *z,
                           SVector3 *n, int numSubEdges)
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index f3affef3cc185ac5030ce1036534fdc273a59cf6..4484c1dad8edc2a2dac8354e0e23b3c436fc5a53 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -38,7 +38,7 @@ class MQuadrangle : public MElement {
     v[2] = _v[2];
     v[3] = _v[3];
   }
-void projectInMeanPlane(double *xn, double *yn);
+  void projectInMeanPlane(double *xn, double *yn);
  public :
   MQuadrangle(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, int num=0, int part=0)
     : MElement(num, part)
@@ -132,6 +132,7 @@ void projectInMeanPlane(double *xn, double *yn);
   {
     return SPoint3(0., 0., 0.);
   }
+  virtual double getVolume();
   virtual void revert()
   {
     MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index a15d74f8ef7611f7f4b96ca7b9193aabe09ccc66..98f977d4f1b8207420a7d009bc6e35f88790d5d9 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -25,6 +25,17 @@ SPoint3 MTriangle::circumcenter()
   return SPoint3(res[0], res[1], res[2]);
 }
 
+double MTriangle::getVolume()
+{
+  if(getNumVertices() > 3)
+    return MElement::getVolume();
+  SPoint3 p0(_v[0]->x(), _v[0]->y(), _v[0]->z());
+  SPoint3 p1(_v[1]->x(), _v[1]->y(), _v[1]->z());
+  SPoint3 p2(_v[2]->x(), _v[2]->y(), _v[2]->z());
+  SVector3 v1(p0, p1), v2(p0, p2);
+  return norm(crossprod(v1, v2)) / 2.;
+}
+
 double MTriangle::distoShapeMeasure()
 {
 #if defined(HAVE_MESH)
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index dbcb47e433c797deb5b72ce71b6d33cdf1026d45..54a298904e0ce0f9e51fd11625759fae9b4d85fe 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -149,6 +149,7 @@ class MTriangle : public MElement {
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
   virtual SPoint3 circumcenter();
+  virtual double getVolume();
   static int edges_tri(const int edge, const int vert)
   {
     static const int e[3][2] = {