Triangle.h

Go to the documentation of this file.

#ifndef TRIANGLE
#define TRIANGLE

#include "AuxDefs.h"
#include "ElementGeometry.h"
#include "Segment.h"

#include <cmath>
#include <cassert>



template<size_t SPD> 
class Triangle: public AbstractGeom<SPD> {
 public:
  Triangle (const std::vector<double>& globalCoordVec, const std::vector<GlobalNodalIndex>& connectivity)
    :AbstractGeom<SPD> (globalCoordVec, connectivity) {
      assert (connectivity.size () == 3);
  }


  inline virtual ~Triangle(){}

  Triangle(const Triangle<SPD> & that) : AbstractGeom<SPD>(that) {
  }

  virtual Triangle<SPD>* clone() const {
    return new Triangle<SPD>(*this);
  }


  inline const size_t getNumVertices() const { return 3; }

  inline const std::string getPolytopeName() const { return "TRIANGLE"; }

  inline const size_t getParametricDimension() const { return 2; }

  inline const size_t getEmbeddingDimension() const { return SPD; }

  void map(const double * X, double *Y) const; 

  void dMap(const double * X, double *DY, double &Jac) const; 

  inline size_t getNumFaces() const { return 3; }

  virtual const double getInRadius(void) const;

  virtual const double getOutRadius(void) const;

  virtual Segment<SPD> * getFaceGeometry(size_t e) const; 

  virtual void computeNormal (size_t e, std::vector<double>& vNormal) const;

private:
  static const size_t    SegmentNodes[];
  static const double ParamCoord[];
};


// Class implementation

template <size_t SPD>
const size_t    Triangle<SPD>::SegmentNodes[] = {0,1,1,2,2,0};

template <size_t SPD>
const double Triangle<SPD>::ParamCoord[] = {1,0,0,1,0,0};


template<size_t SPD> 
void Triangle<SPD>::map(const double * X, double *Y) const
{
  for(size_t i=0; i<SPD; i++)
    Y[i] = X[0]*AbstractGeom<SPD>::getCoordinate(0,i) + X[1]*AbstractGeom<SPD>::getCoordinate(1,i) + (1-X[0]-X[1])*AbstractGeom<SPD>::getCoordinate(2,i);

  return;
}




template<size_t SPD> 
void Triangle<SPD>::dMap(const double * X, double *DY, double &Jac) const
{
  for(size_t i=0; i<SPD; i++)
    {
      DY[                  i] = AbstractGeom<SPD>::getCoordinate(0,i) - AbstractGeom<SPD>::getCoordinate(2,i);
      DY[SPD+i] = AbstractGeom<SPD>::getCoordinate(1,i) - AbstractGeom<SPD>::getCoordinate(2,i);
    }
  
  double g11=0;
  double g22=0;
  double g12=0;

  for(size_t i=0; i<SPD; i++)
    {
      g11 += DY[i]*DY[i];
      g22 += DY[SPD+i]*DY[SPD+i];
      g12 += DY[SPD+i]*DY[i];
    }
  
  Jac=sqrt(g11*g22-g12*g12);
  
  return;
}


template<size_t SPD>   
Segment<SPD> * Triangle<SPD>::getFaceGeometry(size_t e) const
{
  std::vector<GlobalNodalIndex> conn(2);
  switch(e)
  {
    case 0:
      conn[0] = AbstractGeom<SPD>::getConnectivity ()[0];
      conn[1] = AbstractGeom<SPD>::getConnectivity ()[1];
      break;

    case 1:
      conn[0] = AbstractGeom<SPD>::getConnectivity ()[1];
      conn[1] = AbstractGeom<SPD>::getConnectivity ()[2];
      break;

    case 2:
      conn[0] = AbstractGeom<SPD>::getConnectivity ()[2];
      conn[1] = AbstractGeom<SPD>::getConnectivity ()[0];
      break;

    default:
      return 0;
  }

  return new Segment<SPD> (AbstractGeom<SPD>::getGlobalCoordVec (), conn);
}

template<size_t SPD> 
const double Triangle<SPD>:: getInRadius(void) const {
  double a,b,c,s;
  a = b = c = s = 0.0;
  for(size_t i=0; i<SPD; i++) {
    a += (AbstractGeom<SPD>::getCoordinate(1,i) - AbstractGeom<SPD>::getCoordinate(0,i))*
      (AbstractGeom<SPD>::getCoordinate(1,i) - AbstractGeom<SPD>::getCoordinate(0,i)) ;
    b += (AbstractGeom<SPD>::getCoordinate(2,i) - AbstractGeom<SPD>::getCoordinate(1,i))*
      (AbstractGeom<SPD>::getCoordinate(2,i) - AbstractGeom<SPD>::getCoordinate(1,i)) ;
    c += (AbstractGeom<SPD>::getCoordinate(0,i) - AbstractGeom<SPD>::getCoordinate(2,i))*
      (AbstractGeom<SPD>::getCoordinate(0,i) - AbstractGeom<SPD>::getCoordinate(2,i)) ;
  }
  a = sqrt(a);
  b = sqrt(b);
  c = sqrt(c);
  s = (a + b + c)/2.0;
  return(2.0*sqrt(s*(s-a)*(s-b)*(s-c))/(a+b+c));
};


template<size_t SPD> 
const double Triangle<SPD>:: getOutRadius(void) const {
  double a,b,c;
  a = b = c = 0.0;
  for(size_t i=0; i<SPD; i++) {
    a += (AbstractGeom<SPD>::getCoordinate(1,i) - AbstractGeom<SPD>::getCoordinate(0,i))*
      (AbstractGeom<SPD>::getCoordinate(1,i) - AbstractGeom<SPD>::getCoordinate(0,i)) ;
    b += (AbstractGeom<SPD>::getCoordinate(2,i) - AbstractGeom<SPD>::getCoordinate(1,i))*
      (AbstractGeom<SPD>::getCoordinate(2,i) - AbstractGeom<SPD>::getCoordinate(1,i)) ;
    c += (AbstractGeom<SPD>::getCoordinate(0,i) - AbstractGeom<SPD>::getCoordinate(2,i))*
      (AbstractGeom<SPD>::getCoordinate(0,i) - AbstractGeom<SPD>::getCoordinate(2,i)) ;
  }
  a = sqrt(a);
  b = sqrt(b);
  c = sqrt(c);
  return(a*b*c/sqrt((a+b+c)*(b+c-a)*(c+a-b)*(a+b-c)));
};

template <size_t SPD>
void Triangle<SPD>::computeNormal(size_t e, std::vector<double> &VNormal) const {
  double NodalCoord[4];

  size_t    n[2];
  double v[2];

  n[0] = SegmentNodes[e*2];
  n[1] = SegmentNodes[e*2+1];

  Triangle<SPD>::map(&Triangle<SPD>::ParamCoord[2*n[0]], NodalCoord );
  Triangle<SPD>::map(&Triangle<SPD>::ParamCoord[2*n[1]], NodalCoord+2);

  v[0] = NodalCoord[2]-NodalCoord[0];
  v[1] = NodalCoord[3]-NodalCoord[1];

  double norm = sqrt(v[0]*v[0]+v[1]*v[1]);

  if(norm<=0) {
    std::cerr << 
      "The normal cannot be computed. Two vertices of a polytope seem to coincide\n";
  }

  VNormal.push_back( v[1]/norm);
  VNormal.push_back(-v[0]/norm);
}


#endif

---