libgmsh-dev 2.8.5+dfsg-1.1+b1 / usr / include / gmsh / MTriangle.h

// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <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(const 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 getInnerRadius();
  virtual double getOuterRadius();
  virtual double angleShapeMeasure();
  virtual int getNumVertices() const { return 3; }
  virtual MVertex *getVertex(int num){ return _v[num]; }
  virtual const MVertex *getVertex(int num) const { return _v[num]; }
  virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
  virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
  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) const
  {
    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(bool curved){ return 3; }
  virtual void getEdgeRep(bool curved, 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(bool curved){ return 1; }
  virtual void getFaceRep(bool curved, 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 "CPS3"/*"STRI3"*//*"C2D3"*/; }
  virtual void reverse()
  {
    MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
  }
  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    w = 0.;
    switch(num) {
    case 0 : u = 0.; v = 0.; break;
    case 1 : u = 1.; v = 0.; break;
    case 2 : u = 0.; v = 1.; break;
    default: u = 0.; v = 0.; break;
    }
  }
  virtual SPoint3 barycenterUVW() const
  {
    return SPoint3(1/3., 1/3., 0.);
  }
  virtual bool isInside(double u, double v, double w) const
  {
    double tol = _isInsideTolerance;
    if(u < (-tol) || v < (-tol) || u > ((1. + tol) - v) || fabs(w) > tol)
      return false;
    return true;
  }
  virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
  virtual SPoint3 circumcenter();
  virtual double getVolume();
  static int edges_tri(const int edge, const int vert)
  {
    static const int e[3][2] = {
      {0, 1},
      {1, 2},
      {2, 0}
    };
    return e[edge][vert];
  }
};

/*
 * 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(const 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 const MVertex *getVertex(int num) const { 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 void xyz2uvw(double xyz[3], double uvw[3]) const{ MElement::xyz2uvw(xyz, uvw); }
  virtual int getNumEdgeVertices() const { return 3; }
  virtual int getNumEdgesRep(bool curved);
  virtual void getEdgeRep(bool curved, 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(bool curved);
  virtual void getFaceRep(bool curved, 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 "CPS6"/*"STRI65"*//*"C2D6"*/; }
  virtual void reverse()
  {
    MVertex *tmp;
    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) const
  {
    num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

/*
 * 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,
             const 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(const 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 const MVertex *getVertex(int num) const { return num < 3 ? _v[num] : _vs[num - 3]; }
  virtual int getNumFaceVertices() const
  {
    if (ElementType::SerendipityFromTag(getTypeForMSH()) > 0)
      return 0;
    else
      return  (_order - 1) * (_order - 2) / 2;
  }
  virtual void xyz2uvw(double xyz[3], double uvw[3]) const { MElement::xyz2uvw(xyz, uvw); }
  virtual int getNumEdgeVertices() const { return 3 * (_order - 1); }
  virtual int getNumEdgesRep(bool curved);
  virtual int getNumFacesRep(bool curved);
  virtual void getEdgeRep(bool curved, 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(bool curved, 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 == 1 && _vs.size() == 0) return MSH_TRI_3;
    if(_order == 2 && _vs.size() == 3) return MSH_TRI_6;
    if(_order == 3 && _vs.size() == 7) return MSH_TRI_10;
    if(_order == 4 && _vs.size() == 12) return MSH_TRI_15;
    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;

    if(_order == 3 && _vs.size() == 6) return MSH_TRI_9;
    if(_order == 4 && _vs.size() == 9) return MSH_TRI_12;
    if(_order == 5 && _vs.size() == 12) return MSH_TRI_15I;
    if(_order == 6 && _vs.size() == 15) return MSH_TRI_18;
    if(_order == 7 && _vs.size() == 18) return MSH_TRI_21I;
    if(_order == 8 && _vs.size() == 21) return MSH_TRI_24;
    if(_order == 9 && _vs.size() == 24) return MSH_TRI_27;
    if(_order ==10 && _vs.size() == 27) return MSH_TRI_30;
    Msg::Error("no tag matches a p%d triangle with %d vertices", _order, 3+_vs.size());
    return 0;
  }
  virtual int getTypeForVTK() const { return (_order==2) ? 22 : MTriangle::getTypeForVTK(); }
  virtual void reverse()
  {
    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;
  }
  virtual void getNode(int num, double &u, double &v, double &w) const
  {
    num < 3 ? MTriangle::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
  }
};

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