MTriangle.h

// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#ifndef _MTRIANGLE_H_
#define _MTRIANGLE_H_

#include "MElement.h"

/*
 * MTriangle
 *
 *   v
 *   ^
 *   |
 *   2
 *   |`\
 *   |  `\
 *   |    `\
 *   |      `\
 *   |        `\
 *   0----------1 --> u
 *
 */
class MTriangle : public MElement {
 protected:
  MVertex *_v[3];
  void _getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v[0] = _v[edges_tri(num, 0)];
    v[1] = _v[edges_tri(num, 1)];
  }
  void _getFaceVertices(std::vector<MVertex*> &v) const
  {
    v[0] = _v[0];
    v[1] = _v[1];
    v[2] = _v[2];
  }
 public :
  MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
    : MElement(num, part)
  {
    _v[0] = v0; _v[1] = v1; _v[2] = v2;
  }
  MTriangle(std::vector<MVertex*> &v, int num=0, int part=0)
    : MElement(num, part)
  {
    for(int i = 0; i < 3; i++) _v[i] = v[i];
  }
  ~MTriangle(){}
  virtual int getDim() const { return 2; }
  virtual double etaShapeMeasure();
  virtual double gammaShapeMeasure();
  virtual double distoShapeMeasure();
  virtual double getInnerRadius();
  virtual double angleShapeMeasure();
  virtual int getNumVertices() const { return 3; }
  virtual MVertex *getVertex(int num){ return _v[num]; }
  virtual MVertex *getVertexMED(int num)
  {
    static const int map[3] = {0, 2, 1};
    return getVertex(map[num]);
  }
  virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
  {
    if(_v[0] != v1 && _v[0] != v2) return _v[0];
    if(_v[1] != v1 && _v[1] != v2) return _v[1];
    if(_v[2] != v1 && _v[2] != v2) return _v[2];
    return 0;
  }
  virtual int getNumEdges(){ return 3; }
  virtual MEdge getEdge(int num)
  {
    return MEdge(_v[edges_tri(num, 0)], _v[edges_tri(num, 1)]);
  }
  virtual void getEdgeInfo (const MEdge & edge, int &ithEdge, int &sign) const
  {
    for (ithEdge = 0; ithEdge < 3; ithEdge++){
      const MVertex *v0 = _v[edges_tri(ithEdge, 0)];
      const MVertex *v1 = _v[edges_tri(ithEdge, 1)];
      if (v0 == edge.getVertex(0) && v1 == edge.getVertex(1)){
        sign = 1; return;
      }
      if (v1 == edge.getVertex(0) && v0 == edge.getVertex(1)){
        sign = -1; return;
      }
    }
    Msg::Error("Could not get edge information for triangle %d", getNum());
  }
  virtual int getNumEdgesRep(){ return 3; }
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
  {
    MEdge e(getEdge(num));
    _getEdgeRep(e.getVertex(0), e.getVertex(1), x, y, z, n, 0);
  }
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(2);
    _getEdgeVertices(num, v);
  }
  virtual int getNumFaces(){ return 1; }
  virtual MFace getFace(int num)
  {
    return MFace(_v[0], _v[1], _v[2]);
  }
  virtual int getNumFacesRep(){ return 1; }
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
  {
    _getFaceRep(_v[0], _v[1], _v[2], x, y, z, n);
  }
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    _getFaceVertices(v);
  }
  virtual int getType() const { return TYPE_TRI; }
  virtual int getTypeForMSH() const { return MSH_TRI_3; }
  virtual int getTypeForUNV() const { return 91; } // thin shell linear triangle
  virtual int getTypeForVTK() const { return 5; }
  virtual const char *getStringForPOS() const { return "ST"; }
  virtual const char *getStringForBDF() const { return "CTRIA3"; }
  virtual const char *getStringForDIFF() const { return "ElmT3n2D"; }
  virtual const char *getStringForINP() const { return "C2D3"; }
  virtual void revert()
  {
    MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
  }
  virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
  virtual bool isInside(double u, double v, double w)
  {
    double tol = _isInsideTolerance;
    if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v))
      return false;
    return true;
  }
  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
  virtual SPoint3 circumcenter();
 private:
  int edges_tri(const int edge, const int vert) const
  {
    static const int e[3][2] = {
      {0, 1},
      {1, 2},
      {2, 0}
    };
    return e[edge][vert];
  }
  public:
  static void registerBindings(binding *b);
};

/*
 * MTriangle6
 *
 *   2
 *   |`\
 *   |  `\
 *   5    `4
 *   |      `\
 *   |        `\
 *   0-----3----1
 *
 */
class MTriangle6 : public MTriangle {
 protected:
  MVertex *_vs[3];
 public :
  MTriangle6(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
             MVertex *v5, int num=0, int part=0)
    : MTriangle(v0, v1, v2, num, part)
  {
    _vs[0] = v3; _vs[1] = v4; _vs[2] = v5;
    for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
  }
  MTriangle6(std::vector<MVertex*> &v, int num=0, int part=0)
    : MTriangle(v, num, part)
  {
    for(int i = 0; i < 3; i++) _vs[i] = v[3 + i];
    for(int i = 0; i < 3; i++) _vs[i]->setPolynomialOrder(2);
  }
  ~MTriangle6(){}
  virtual int getPolynomialOrder() const { return 2; }
  virtual int getNumVertices() const { return 6; }
  virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
  virtual MVertex *getVertexUNV(int num)
  {
    static const int map[6] = {0, 3, 1, 4, 2, 5};
    return getVertex(map[num]);
  }
  virtual MVertex *getVertexMED(int num)
  {
    static const int map[6] = {0, 2, 1, 5, 4, 3};
    return getVertex(map[num]);
  }
  virtual int getNumEdgeVertices() const { return 3; }
  virtual int getNumEdgesRep();
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3);
    MTriangle::_getEdgeVertices(num, v);
    v[2] = _vs[num];
  }
  virtual int getNumFacesRep();
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(6);
    MTriangle::_getFaceVertices(v);
    v[3] = _vs[0];
    v[4] = _vs[1];
    v[5] = _vs[2];
  }
  virtual int getTypeForMSH() const { return MSH_TRI_6; }
  virtual int getTypeForUNV() const { return 92; } // thin shell parabolic triangle
  //virtual int getTypeForVTK() const { return 22; }
  virtual const char *getStringForPOS() const { return "ST2"; }
  virtual const char *getStringForBDF() const { return "CTRIA6"; }
  virtual const char *getStringForDIFF() const { return "ElmT6n2D"; }
  virtual const char *getStringForINP() const { return "C2D6"; }
  virtual void revert()
  {
    MVertex *tmp;
    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
    tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
  }
};

/*
 * MTriangleN  FIXME: check the plot
 *
 *   2
 *   |`\                E = order - 1;
 *   |  `\              N = total number of vertices
 * 3+2E   2+2E
 *   |      `\          Interior vertex numbers
 *  ...       ...         for edge 0 <= i <= 2: 3+i*E to 2+(i+1)*E
 *   |          `\        in volume           : 3+3*E to N-1
 * 2+3E           3+E
 *   |  3+3E to N-1 `\
 *   |                `\
 *   0---3--...---2+E---1
 *
 */
class MTriangleN : public MTriangle {
 protected:
  std::vector<MVertex *> _vs;
  const char _order;
 public:
  MTriangleN(MVertex *v0, MVertex *v1, MVertex *v2,
             std::vector<MVertex*> &v, char order, int num=0, int part=0)
    : MTriangle(v0, v1, v2, num, part), _vs(v), _order(order)
  {
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
  MTriangleN(std::vector<MVertex*> &v, char order, int num=0, int part=0)
    : MTriangle(v[0], v[1], v[2], num, part), _order(order)
  {
    for(unsigned int i = 3; i < v.size(); i++) _vs.push_back(v[i]);
    for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
  }
  ~MTriangleN(){}
  virtual int getPolynomialOrder() const { return _order; }
  virtual int getNumVertices() const { return 3 + _vs.size(); }
  virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
  virtual int getNumFaceVertices() const
  {
    if(_order == 3 && _vs.size() == 6) return 0;
    if(_order == 3 && _vs.size() == 7) return 1;
    if(_order == 4 && _vs.size() == 9) return 0;
    if(_order == 4 && _vs.size() == 12) return 3;
    if(_order == 5 && _vs.size() == 12) return 0;
    if(_order == 5  && _vs.size() == 18) return 6;
    if(_order == 6  && _vs.size() == 25) return 10;
    if(_order == 7  && _vs.size() == 33) return 12;
    if(_order == 8  && _vs.size() == 42) return 15;
    if(_order == 9  && _vs.size() == 52) return 21;
    if(_order == 10 && _vs.size() == 63) return 28;
    if(_order == 6  && _vs.size() == 15) return 0;
    if(_order == 7  && _vs.size() == 18) return 0;
    if(_order == 8  && _vs.size() == 21) return 0;
    if(_order == 9  && _vs.size() == 24) return 0;
    if(_order == 10  && _vs.size() == 27) return 0;
    return 0;
  }
  virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
  virtual int getNumEdgesRep();
  virtual int getNumFacesRep();
  virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(_order + 1);
    MTriangle::_getEdgeVertices(num, v);
    int j = 2;
    const int ie = (num + 1) * (_order - 1);
    for(int i = num * (_order-1); i != ie; ++i) v[j++] = _vs[i];
  }
  virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
  virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
  {
    v.resize(3 + _vs.size());
    MTriangle::_getFaceVertices(v);
    for(unsigned int i = 0; i != _vs.size(); ++i) v[i + 3] = _vs[i];
  }
  virtual int getTypeForMSH() const
  {
    if(_order == 2 && _vs.size() == 3) return MSH_TRI_6;
    if(_order == 3 && _vs.size() == 6) return MSH_TRI_9;
    if(_order == 3 && _vs.size() == 7) return MSH_TRI_10;
    if(_order == 4 && _vs.size() == 9) return MSH_TRI_12;
    if(_order == 4 && _vs.size() == 12) return MSH_TRI_15;
    if(_order == 5 && _vs.size() == 12) return MSH_TRI_15I;
    if(_order == 5 && _vs.size() == 18) return MSH_TRI_21;
    if(_order == 6 && _vs.size() == 25) return MSH_TRI_28;
    if(_order == 7 && _vs.size() == 33) return MSH_TRI_36;
    if(_order == 8 && _vs.size() == 42) return MSH_TRI_45;
    if(_order == 9 && _vs.size() == 52) return MSH_TRI_55;
    if(_order ==10 && _vs.size() == 63) return MSH_TRI_66;
    return 0;
  }
  virtual void revert()
  {
    MVertex *tmp;
    tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
    std::vector<MVertex*> inv;
    inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
    _vs = inv;
  }
};

template <class T>
void sort3(T *t[3])
{
  T *temp;
  if(t[0] > t[1]){
    temp = t[1];
    t[1] = t[0];
    t[0] = temp;
  }
  if(t[1] > t[2]){
    temp = t[2];
    t[2] = t[1];
    t[1] = temp;
  }
  if(t[0] > t[1]){
    temp = t[1];
    t[1] = t[0];
    t[0] = temp;
  }
}

struct compareMTriangleLexicographic
{
  bool operator () (MTriangle *t1, MTriangle *t2) const
  {
    MVertex *_v1[3] = {t1->getVertex(0), t1->getVertex(1), t1->getVertex(2)};
    MVertex *_v2[3] = {t2->getVertex(0), t2->getVertex(1), t2->getVertex(2)};
    sort3(_v1);
    sort3(_v2);
    if(_v1[0] < _v2[0]) return true;
    if(_v1[0] > _v2[0]) return false;
    if(_v1[1] < _v2[1]) return true;
    if(_v1[1] > _v2[1]) return false;
    if(_v1[2] < _v2[2]) return true;
    return false;
  }
};

#endif