diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 24cb8009a15d43167e643f2614383872722bcc80..71af67d7f994f2119052e10f95e90fbf7308c294 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -846,6 +846,31 @@ void MElement::xyz2uvw(double xyz[3], double uvw[3]) const
   // routines are implemented for simplices, where the basis functions
   // are linear)
   uvw[0] = uvw[1] = uvw[2] = 0.;
+
+  // For high order elements, start from the nearer point
+  if (getPolynomialOrder() > 2) {
+    int numNearer = 0;
+    const MVertex *v = getShapeFunctionNode(0);
+    double distNearer = (v->x()-xyz[0])*(v->x()-xyz[0]) +
+                        (v->y()-xyz[1])*(v->y()-xyz[1]) +
+                        (v->z()-xyz[2])*(v->z()-xyz[2]);
+    for (int i = 1; i < getNumShapeFunctions(); i++) {
+      const MVertex *v = getShapeFunctionNode(i);
+      double dist = (v->x()-xyz[0])*(v->x()-xyz[0]) +
+                    (v->y()-xyz[1])*(v->y()-xyz[1]) +
+                    (v->z()-xyz[2])*(v->z()-xyz[2]);
+      if (dist < distNearer) {
+        numNearer = i;
+        distNearer = dist;
+      }
+    }
+    const nodalBasis *nb = getFunctionSpace();
+    fullMatrix<double> refpnts = nb->getReferenceNodes();
+    uvw[0] = refpnts(numNearer, 0);
+    uvw[1] = refpnts(numNearer, 1);
+    uvw[2] = refpnts(numNearer, 2);
+  }
+
   int iter = 1, maxiter = 20;
   double error = 1., tol = 1.e-6;
 
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index ec61eff8cbd70f5c4b76257298543f4b0dbbe616..472aed9c16556e501faf1ff9e26a555c0b088fe1 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -23,6 +23,9 @@ class nodalBasis {
   void getReferenceNodes(fullMatrix<double> &nodes) const {
     nodes = points;
   }
+  const fullMatrix<double>& getReferenceNodes() const {
+    return points;
+  }
   void getReferenceNodesForBezier(fullMatrix<double> &nodes) const;
 
   // Basis functions & gradients evaluation