spetialize xyz2uvw for p1Triangle

Geo/MTriangle.cpp

@@ -91,43 +91,24 @@ double MTriangle::gammaShapeMeasure()
 #endif
 }
 
-// void MTriangle::xyz2uvw(double xyz[3], double uvw[3]){
-
-//   double X[3] = {_v[0]->x(), _v[1]->x(), _v[2]->x()};
-//   double Y[3] = {_v[0]->y(), _v[1]->y(), _v[2]->y()};
-//   double Z[3] = {_v[0]->z(), _v[1]->z(), _v[2]->z()};
-
-  // uvw[0] = uvw[1] = uvw[2] = 0.;
-  
-  // double jac[3][3] = {{X[1]-X[0], X[2]-X[Ø], 0},{}, {}};
-
-
-  // if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
-
-  // double xn = 0., yn = 0., zn = 0.;
-  // double sf[1256];
-  // getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
-  // for (int i = 0; i < getNumShapeFunctions(); i++) {
-  //   MVertex *v = getShapeFunctionNode(i);
-  //   xn += v->x() * sf[i];
-  //   yn += v->y() * sf[i];
-  //   zn += v->z() * sf[i];
-  // }
-  // double inv[3][3];
-  // inv3x3(jac, inv);
-  
-  // double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
-  //   inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
-  // double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) +
-  //   inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn);
-  // double wn = uvw[2] + inv[0][2] * (xyz[0] - xn) +
-  //   inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn);
-  
-  // uvw[0] = un;
-  // uvw[1] = vn;
-  // uvw[2] = wn;
-  
-  //}
+ void MTriangle::xyz2uvw(double xyz[3], double uvw[3])
+ {
+   const double O[3] = {_v[0]->x(), _v[0]->y(), _v[0]->z()};
+   const double d[3] = {xyz[0] - O[0], xyz[1] - O[1], xyz[2] - O[2]};
+   const double d1[3] = {_v[1]->x() - O[0], _v[1]->y() - O[1], _v[1]->z() - O[2]};
+   const double d2[3] = {_v[2]->x() - O[0], _v[2]->y() - O[1], _v[2]->z() - O[2]};
+   const double Jxy = d1[0] * d2[1] - d1[1] * d2[0];
+   const double Jxz = d1[0] * d2[2] - d1[2] * d2[0];
+   if (fabs(Jxy) > fabs(Jxz)) {
+     uvw[0] = (d[0] * d2[1] - d[1] * d2[0]) / Jxy;
+     uvw[1] = (d[1] * d1[0] - d[0] * d1[1]) / Jxy;
+   } 
+   else {
+     uvw[0] = (d[0] * d2[2] - d[2] * d2[0]) / Jxz;
+     uvw[1] = (d[2] * d1[0] - d[0] * d1[2]) / Jxz;
+   }
+   uvw[2] = 0.;
+}
 
 const polynomialBasis* MTriangle::getFunctionSpace(int o) const
 {

Geo/MTriangle.h

@@ -63,7 +63,7 @@ class MTriangle : public MElement {
     static const int map[3] = {0, 2, 1};
     return getVertex(map[num]);
   }
-  //virtual void xyz2uvw(double xyz[3], double uvw[3]);
+  virtual void xyz2uvw(double xyz[3], double uvw[3]);
   virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
   {
     if(_v[0] != v1 && _v[0] != v2) return _v[0];