From 716f2d351fa7ab0e8bcaf7372a677d1fa4d1dc43 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Tue, 29 Mar 2011 05:27:10 +0000
Subject: [PATCH] A crude getNode function for MElements without a lookup
 table.

---
 Geo/MElement.cpp   | 34 ++++++++++++++++++++++++++++------
 Geo/MElement.h     |  8 ++------
 Geo/MHexahedron.h  | 14 +++++++++++++-
 Geo/MLine.h        | 10 +++++++++-
 Geo/MPoint.h       |  2 +-
 Geo/MPrism.h       | 10 +++++++++-
 Geo/MPyramid.h     | 10 +++++++++-
 Geo/MQuadrangle.h  | 14 +++++++++++++-
 Geo/MTetrahedron.h | 10 +++++++++-
 Geo/MTriangle.h    | 10 +++++++++-
 10 files changed, 102 insertions(+), 20 deletions(-)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 7b46f6ebdc..a87c799e6f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -107,6 +107,19 @@ double MElement::rhoShapeMeasure()
     return 0.;
 }
 
+void MElement::getNode(int num, double &u, double &v, double &w)
+{
+  // only for MElements that don't have a lookup table for this
+  // (currently only 1st order elements have)
+  double uvw[3];
+  MVertex* ver = getVertex(num);
+  double xyz[3] = {ver->x(), ver->y(), ver->z()};
+  xyz2uvw(xyz, uvw);
+  u = uvw[0];
+  v = uvw[1];
+  w = uvw[2];
+}
+
 void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
 {
   const polynomialBasis* fs = getFunctionSpace(o);
@@ -565,12 +578,21 @@ double MElement::integrateFlux(double val[], int face, int pOrder, int order)
   getFaceVertices(face, v);
   MElementFactory f;
   int type = 0;
-  if(getType() == TYPE_TRI || getType() == TYPE_TET) type = getTriangleType(getPolynomialOrder());
-  else if(getType() == TYPE_QUA || getType() == TYPE_HEX) type = getQuadType(getPolynomialOrder());
-  else if(getType() == TYPE_PYR && face < 4) type = getTriangleType(getPolynomialOrder());
-  else if(getType() == TYPE_PYR && face >= 4) type = getQuadType(getPolynomialOrder());
-  else if(getType() == TYPE_PRI && face < 2) type = getTriangleType(getPolynomialOrder());
-  else if(getType() == TYPE_PRI && face >= 2) type = getQuadType(getPolynomialOrder());
+  switch(getType()) {
+  case TYPE_TRI : type = getTriangleType(getPolynomialOrder()); break;
+  case TYPE_TET : type = getTriangleType(getPolynomialOrder()); break;
+  case TYPE_QUA : type = getQuadType(getPolynomialOrder());     break;
+  case TYPE_HEX : type = getQuadType(getPolynomialOrder());     break;
+  case TYPE_PYR : 
+    if(face < 4) type = getTriangleType(getPolynomialOrder());
+    else type = getQuadType(getPolynomialOrder());
+    break;
+  case TYPE_PRI :
+    if(face < 2) type = getTriangleType(getPolynomialOrder());
+    else type = getQuadType(getPolynomialOrder());
+    break;
+  default: type = 0; break;
+  }
   MElement* fe = f.create(type, v);
 
   double intv[3];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index aa649f97d1..b4836a4dc1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -220,12 +220,8 @@ class MElement
   // get the function space for the jacobian of the element
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const { return 0; }
 
-  // return parametric coordinates (u, v, w) of a vertex
-  virtual void getNode(int num, double &u, double &v, double &w)
-  {
-    u = v = w = 0.;
-    Msg::Error("Node get not available for this element");
-  }
+  // return parametric coordinates (u,v,w) of a vertex
+  virtual void getNode(int num, double &u, double &v, double &w);
 
   // return the interpolating nodal shape functions evaluated at point
   // (u,v,w) in parametric coordinates (if order == -1, use the
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index f015c66e3f..4336efe7cc 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -161,7 +161,7 @@ class MHexahedron : public MElement {
     s[7][1] =  0.125 * (1. - u) * (1. + w);
     s[7][2] =  0.125 * (1. - u) * (1. + v);
   }
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     switch(num) {
     case 0 : u = -1.; v = -1.; w = -1.; break;
@@ -352,6 +352,10 @@ class MHexahedron20 : public MHexahedron {
     _vs[8] = old[10]; _vs[10] = old[8];
     _vs[9] = old[11]; _vs[11] = old[9];
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -493,6 +497,10 @@ class MHexahedron27 : public MHexahedron {
     _vs[13] = old[15]; _vs[15] = old[13];
     _vs[14] = old[16]; _vs[16] = old[14];
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -632,6 +640,10 @@ class MHexahedronN : public MHexahedron {
   {
     Msg::Error("not done yet reverse hexN");
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 8 ? MHexahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 #endif
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 207aaceae6..bfbcde1eac 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -78,7 +78,7 @@ class MLine : public MElement {
       return false;
     return true;
   }
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     v = w = 0.;
     switch(num) {
@@ -142,6 +142,10 @@ class MLine3 : public MLine {
   //virtual int getTypeForVTK() const { return 21; }
   virtual const char *getStringForPOS() const { return "SL2"; }
   virtual const char *getStringForINP() const { return "C1D3"; }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -198,6 +202,10 @@ class MLineN : public MLine {
     if(_vs.size() == 9) return MSH_LIN_11;
     return 0;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 struct compareMLinePtr {
diff --git a/Geo/MPoint.h b/Geo/MPoint.h
index 1d0e201f82..01402a5553 100644
--- a/Geo/MPoint.h
+++ b/Geo/MPoint.h
@@ -42,7 +42,7 @@ class MPoint : public MElement {
   virtual int getTypeForMSH() const { return MSH_PNT; }
   virtual int getTypeForVTK() const { return 1; }
   virtual const char *getStringForPOS() const { return "SP"; }
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     u = v = w = 0.;
   }
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index cd3b22c2f2..8702faa1c1 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -159,7 +159,7 @@ class MPrism : public MElement {
     s[5][1] =  0.5 * (1. + w)    ;
     s[5][2] =  0.5 * v           ;
   }*/
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     switch(num) {
     case 0 : u = 0.; v = 0.; w = -1.; break;
@@ -329,6 +329,10 @@ class MPrism15 : public MPrism {
     tmp = _vs[2]; _vs[2] = _vs[4]; _vs[4] = tmp;
     tmp = _vs[7]; _vs[7] = _vs[8]; _vs[8] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -451,6 +455,10 @@ class MPrism18 : public MPrism {
     // quad face vertices
     tmp = _vs[10]; _vs[10] = _vs[11]; _vs[11] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 6 ? MPrism::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 #endif
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 423352f07e..6ee6e36ea5 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -165,7 +165,7 @@ class MPyramid : public MElement {
     s[4][1] = 0.;
     s[4][2] = 1.;
   }
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     switch(num) {
     case 0 : u = -1.; v = -1.; w = 0.; break;
@@ -321,6 +321,10 @@ class MPyramid13 : public MPyramid {
     tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp;
     tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -437,6 +441,10 @@ class MPyramid14 : public MPyramid {
     tmp = _vs[1]; _vs[1] = _vs[5]; _vs[5] = tmp;
     tmp = _vs[2]; _vs[2] = _vs[6]; _vs[6] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 5 ? MPyramid::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 #endif
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 28c2fa5cb1..11e1bc18dc 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -117,7 +117,7 @@ void projectInMeanPlane(double *xn, double *yn);
   virtual const char *getStringForINP() const { return "C2D4"; }
   virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     w = 0.;
     switch(num) {
@@ -237,6 +237,10 @@ class MQuadrangle8 : public MQuadrangle {
     tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -309,6 +313,10 @@ class MQuadrangle9 : public MQuadrangle {
     tmp = _vs[0]; _vs[0] = _vs[3]; _vs[3] = tmp;
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -399,6 +407,10 @@ class MQuadrangleN : public MQuadrangle {
     inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
     _vs = inv;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 4 ? MQuadrangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 template <class T>
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 6d4b761bed..f666db4a18 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -130,7 +130,7 @@ class MTetrahedron : public MElement {
   void xyz2uvw(double xyz[3], double uvw[3]);
   virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     switch(num) {
     case 0 : u = 0.; v = 0.; w = 0.; break;
@@ -269,6 +269,10 @@ class MTetrahedron10 : public MTetrahedron {
     tmp = _vs[1]; _vs[1] = _vs[2]; _vs[2] = tmp;
     tmp = _vs[5]; _vs[5] = _vs[3]; _vs[3] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -394,6 +398,10 @@ class MTetrahedronN : public MTetrahedron {
   virtual int getNumEdgesRep();
   virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
   virtual int getNumFacesRep();
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 4 ? MTetrahedron::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 #endif
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 0e64bb3e27..6b675fe4ac 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -125,7 +125,7 @@ class MTriangle : public MElement {
   }
   virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
   virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
-  void getNode(int num, double &u, double &v, double &w)
+  virtual void getNode(int num, double &u, double &v, double &w)
   {
     w = 0.;
     switch(num) {
@@ -226,6 +226,10 @@ class MTriangle6 : public MTriangle {
     tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
     tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 /*
@@ -328,6 +332,10 @@ class MTriangleN : public MTriangle {
     inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
     _vs = inv;
   }
+  virtual void getNode(int num, double &u, double &v, double &w)
+  {
+    num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
+  }
 };
 
 template <class T>
-- 
GitLab