diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 5881e05096c70058b2b9e693aa3b42ce555f2486..4c6dbae50818ee3018ad5d9090894ef982c24abf 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -158,6 +158,23 @@ SPoint3 MElement::barycenter()
   return p;
 }
 
+SPoint3 MElement::barycenterUVW()
+{
+  SPoint3 p(0., 0., 0.);
+  int n = getNumVertices();
+  for(int i = 0; i < n; i++) {
+    double x, y, z;
+    getNode(i, x, y, z);
+    p[0] += x;
+    p[1] += y;
+    p[2] += z;
+  }
+  p[0] /= (double)n;
+  p[1] /= (double)n;
+  p[2] /= (double)n;
+  return p;
+}
+
 double MElement::getVolume()
 {
   int npts;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 146d116803403e7c4fcfa56e08669b374d0c9ae2..4d3cc54a502a1bafc19e30bce668637eecfc259b 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -194,6 +194,7 @@ class MElement
 
   // compute the barycenter
   virtual SPoint3 barycenter();
+  virtual SPoint3 barycenterUVW();
 
   // revert the orientation of the element
   virtual void revert(){}
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index a60ae836a5bfd53d0688d08419f47a93f4bdc8eb..9c552acbab7138ad1f1aec74c09e062fa9ada2a7 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -137,6 +137,10 @@ class MHexahedron : public MElement {
     default: u =  0.; v =  0.; w =  0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 53e77dd72a663b72c21b461c91ae83ea7de525cf..a26fbaf26846034eba562dfd8e27b706abb22498 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -88,6 +88,10 @@ class MLine : public MElement {
     default: u =  0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0, 0, 0);
+  }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
 };
 
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index e294bd5813d665fae7db88b7f648b96668098992..1b094360514b10001631c690f77729a866c8a3fd 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -46,6 +46,10 @@ class MPoint : public MElement {
   {
     u = v = w = 0.;
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
   virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
   {
     s[0] = 1.;
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 33c81597fc52d0aa06083ffbd0e6aea10c920e05..d6e14731ca620190d024c153eb826a55d2b9bd81 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -142,6 +142,10 @@ class MPrism : public MElement {
     default: u = 0.; v = 0.; w =  0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(1/3., 1/3., 0.);
+  }
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 90676b6f21986a2ca5961f45865e55e2108373d7..065d1186dc74602c09b4c6d561e01984d6ebb580 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -176,6 +176,10 @@ class MPyramid : public MElement {
     default: u =  0.; v =  0.; w = 0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., .2);
+  }
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index a0c934faf3db9695d2a8263f25420fa36d43a1f4..f3affef3cc185ac5030ce1036534fdc273a59cf6 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -128,6 +128,10 @@ void projectInMeanPlane(double *xn, double *yn);
     default: u =  0.; v =  0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
   virtual void revert()
   {
     MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
@@ -240,6 +244,10 @@ class MQuadrangle8 : public MQuadrangle {
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
 };
 
 /*
@@ -316,6 +324,10 @@ class MQuadrangle9 : public MQuadrangle {
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
 };
 
 /*
@@ -412,6 +424,10 @@ class MQuadrangleN : public MQuadrangle {
   {
     num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(0., 0., 0.);
+  }
 };
 
 template <class T>
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index ffae493d343ead7f779d2185c8b45b3098134936..3f017d8e8fb382c5699e0f1950921ee98db60478 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -140,6 +140,10 @@ class MTetrahedron : public MElement {
     default: u = 0.; v = 0.; w = 0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(.25, .25, .25);
+  }
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 2f78ab4fc346ba126e8c57043d61e3af80cc1bd2..dbcb47e433c797deb5b72ce71b6d33cdf1026d45 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -136,6 +136,10 @@ class MTriangle : public MElement {
     default: u = 0.; v = 0.; break;
     }
   }
+  virtual SPoint3 barycenterUVW()
+  {
+    return SPoint3(1/3., 1/3., 0.);
+  }
   virtual bool isInside(double u, double v, double w)
   {
     double tol = _isInsideTolerance;