#ifndef _2D_LINEAR_CC
#define _2D_LINEAR_CC

//#include "2d_FE.h"

// Linear functions on triangles in 2D. The reference triangle
//   has vertices at (0,0), (1,0) and (0,1).

//mapping from a reference element to a given FE
Point map_to_FE(Point p, int ele)
{
   // The mapping from the reference element to a given element is the linear
   //   function x = Bx_ + a, here x is the point on the element, x_ is the
   //   point on the reference element, B is 2x2 matrix and a is a vector.

   Point result;

   // the transformation:
   double a1 = nodes[elements[ele].vertices[0]].p.x;
   double a2 = nodes[elements[ele].vertices[0]].p.y;
   double B11 = nodes[elements[ele].vertices[1]].p.x - a1;
   double B12 = nodes[elements[ele].vertices[2]].p.x - a1;   
   double B21 = nodes[elements[ele].vertices[1]].p.y - a2;
   double B22 = nodes[elements[ele].vertices[2]].p.y - a2;

   result.x = B11*p.x + B12*p.y + a1;
   result.y = B21*p.x + B22*p.y + a2;
   
   return result;
}

//mapping from a a given FE to a reference element
Point map_from_FE(Point p, int ele)
{
   // Here the mapping is the inverse of the above
   //   i.e. x_ = B^{-1}(x-a)

   Point result;

   // the transformation:
   double a1 = nodes[elements[ele].vertices[0]].p.x;
   double a2 = nodes[elements[ele].vertices[0]].p.y;
   double B11 = nodes[elements[ele].vertices[1]].p.x - a1;
   double B12 = nodes[elements[ele].vertices[2]].p.x - a1;   
   double B21 = nodes[elements[ele].vertices[1]].p.y - a2;
   double B22 = nodes[elements[ele].vertices[2]].p.y - a2;

   // the inverse matrix is 1/det (B22 -B12; -B21 B11)
   double det = (B11*B22 - B12*B21);

   result.x = (B22*(p.x-a1) - B12*(p.y-a2))/det;
   result.y = (-B21*(p.x-a1) + B11*(p.y-a2))/det;
   
   return result;
}

// the number of local degrees of freedom
int dofs_per_element()
{
   return 3;
}

// the mapping of local to global degrees of freedom
void fill_element_dofs(int* array, int ele)
{
   array[0] = elements[ele].vertices[0];
   array[1] = elements[ele].vertices[1];
   array[2] = elements[ele].vertices[2];   
}

//the i-th shape function on the reference element
double reference_shape_funct(int index, Point p)
{
   if (index==0)
      return 1-p.x-p.y;
   if (index==1)
      return p.x;
   if (index==2)
      return p.y;
   
}

//the x-derivative of the i-th shape function on the reference element
double reference_shape_funct_x(int index, Point p)
{
   if (index==0)
      return -1;
   if (index==1)
      return 1;
   if (index==2)
      return 0;
}

//the y-derivative of the i-th shape function on the reference element
double reference_shape_funct_y(int index, Point p)
{
   if (index==0)
      return -1;
   if (index==1)
      return 0;
   if (index==2)
      return 1;
}


// the jacobian of the transformation
//   (i.e. the change of variables in the integral over an element
//      to the reference element)
double jacobian(Point q_point, int ele)
{
   // the transformation:
   double a1 = nodes[elements[ele].vertices[0]].p.x;
   double a2 = nodes[elements[ele].vertices[0]].p.y;
   double B11 = nodes[elements[ele].vertices[1]].p.x - a1;
   double B12 = nodes[elements[ele].vertices[2]].p.x - a1;   
   double B21 = nodes[elements[ele].vertices[1]].p.y - a2;
   double B22 = nodes[elements[ele].vertices[2]].p.y - a2;

   // the inverse matrix is 1/det (B22 -B12; -B21 B11)
   double det = (B11*B22 - B12*B21);

   return det;
}


/*********
          NOTE
**********/          
// The following two functions are not needed in general,
//   since we are only working on the reference element.
// I've included them here for clarity.
double shape_funct(Point p, int ele, int index)
{
   return reference_shape_funct(index, map_from_FE(p,ele));
}

double shape_funct_x(Point p, int ele, int index)
{
   double dx;
   double dy;

   // the transformation:
   double a1 = nodes[elements[ele].vertices[0]].p.x;
   double a2 = nodes[elements[ele].vertices[0]].p.y;
   double B11 = nodes[elements[ele].vertices[1]].p.x - a1;
   double B12 = nodes[elements[ele].vertices[2]].p.x - a1;   
   double B21 = nodes[elements[ele].vertices[1]].p.y - a2;
   double B22 = nodes[elements[ele].vertices[2]].p.y - a2;

   // the inverse matrix is 1/det (B22 -B12; -B21 B11)
   double det = (B11*B22 - B12*B21);

   dx = B22/det;
   dy = -B21/det;

   return reference_shape_funct_x(index, map_from_FE(p,ele))*dx
      + reference_shape_funct_y(index, map_from_FE(p,ele))*dy;
}

double shape_funct_y(Point p, int ele, int index)
{
   double dx;
   double dy;

   // the transformation:
   double a1 = nodes[elements[ele].vertices[0]].p.x;
   double a2 = nodes[elements[ele].vertices[0]].p.y;
   double B11 = nodes[elements[ele].vertices[1]].p.x - a1;
   double B12 = nodes[elements[ele].vertices[2]].p.x - a1;   
   double B21 = nodes[elements[ele].vertices[1]].p.y - a2;
   double B22 = nodes[elements[ele].vertices[2]].p.y - a2;

   // the inverse matrix is 1/det (B22 -B12; -B21 B11)
   double det = (B11*B22 - B12*B21);

   dx = -B12/det;
   dy = B11/det;
   
   return reference_shape_funct_x(index, map_from_FE(p,ele))*dx
      + reference_shape_funct_y(index, map_from_FE(p,ele))*dy;
}

#endif