//********************************************************************************
//*
//*  C++ finite element method for heat equation
//*  James sandham
//*  15 April 2015
//*
//********************************************************************************

//********************************************************************************
//
// HeatFE is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.
//
//********************************************************************************

#include<iostream>
#include"Shape.h"


//--------------------------------------------------------------------------------
// Shape function Constructor 
//--------------------------------------------------------------------------------
Shape::Shape(int typ)
{
  Type = typ;

  switch (typ)
  {
  case 1:     //linear 1D lines
    Npe = 2; nDim =1;
    break;
  case 2:     //linear 2D triangles
    Npe = 3; nDim = 2;
    break;
  case 3:     //linear 2D squares
    Npe = 4; nDim = 2;
    break;
  case 4:     //linear 3D tetrahedra
    Npe = 4; nDim = 3;
    break;
  case 8:     //quadratic 1D lines
    Npe = 3; nDim = 1;
    break;
  case 9:     //quadratic 2D triangles
    Npe = 6; nDim = 2;
    break;
  case 10:    //quadratic 2D squares
    Npe = 8; nDim = 2;
    break;
  case 11:    //quadratic 3D tetrahedra
    Npe = 10; nDim = 3;
    break;
  case 15:    //point
    Npe = 1; nDim = 0;
    break;
  }
}



//--------------------------------------------------------------------------------
// Shape functions 
//--------------------------------------------------------------------------------
void Point0D::shape(double xi, double et, double ze)
{
  //shape functions for 0D points at node (xi=1,et=1,ze=1)
  N[0] = 1.0;
}


void LineLin1D::shape(double xi, double et, double ze)
{
  //shape functions for linear 1D lines at node (xi,et=1,ze=1)
  N[0] = 0.5*(1-xi);
  N[1] = 0.5*(1+xi);
}


void TriLin2D::shape(double xi, double et, double ze)
{
  //shape functions for linear 2D triangles at node (xi,et,ze=1)
  N[0] = xi;
  N[1] = et;
  N[2] = 1-xi-et;
}


void SqrLin2D::shape(double xi, double et, double ze)
{
  //shape functions for linear 2D squares at node (xi,et,ze=1)
  N[0] = 0.25*(1-xi)*(1-et);
  N[1] = 0.25*(1+xi)*(1-et);
  N[2] = 0.25*(1+xi)*(1+et);
  N[3] = 0.25*(1-xi)*(1+et);
}


void TetLin3D::shape(double xi, double et, double ze)
{
  //shape functions for linear 3D tetrahedra at node (xi,et,ze)
  N[0] = xi;
  N[1] = et;
  N[2] = ze;
  N[3] = 1-xi-et-ze;
}


void LineQuad1D::shape(double xi, double et, double ze)
{
  //shape functions for quadratic 1D lines at node (xi,et=1,ze=1)
  N[0] = -0.5*xi*(1-xi);
  N[1] = 0.5*xi*(1+xi);
  N[2] = 1-xi*xi;
}


void TriQuad2D::shape(double xi, double et, double ze)
{
  //shape functions for quadratic 2D triangles at node (xi,et,ze=1)
  N[0] = xi*(2*xi-1);
  N[1] = et*(2*et-1);
  N[2] = (1-xi-et)*(2*(1-xi-et)-1);
  N[3] = 4*xi*et;
  N[4] = 4*et*(1-xi-et);
  N[5] = 4*xi*(1-xi-et);
}


void SqrQuad2D::shape(double xi, double et, double ze)
{
  //shape functions for quadratic 2D squares at node (xi,et,ze=1)
  N[0] = 0.25*(1-xi)*(et-1)*(xi+et+1);
  N[1] = 0.25*(1+xi)*(et-1)*(et-xi+1);
  N[2] = 0.25*(1+xi)*(1+et)*(xi+et-1);
  N[3] = 0.25*(xi-1)*(et+1)*(xi-et+1);
  N[4] = 0.5*(1-et)*(1-xi*xi);
  N[5] = 0.5*(1+xi)*(1-et*et);
  N[6] = 0.5*(1+et)*(1-xi*xi);
  N[7] = 0.5*(1-xi)*(1-et*et);
}


void TetQuad3D::shape(double xi, double et, double ze)
{
  //shape functions for quadratic 3D tetrahedra at node (xi,et,ze)
  N[0] = (1-xi-et-ze)*(2*(1-xi-et-ze)-1);
  N[1] = xi*(2*xi-1);
  N[2] = et*(2*et-1);
  N[3] = ze*(2*ze-1);
  N[4] = 4*xi*(1-xi-et-ze);
  N[5] = 4*xi*et;
  N[6] = 4*et*(1-xi-et-ze);
  N[7] = 4*ze*(1-xi-et-ze);
  N[8] = 4*xi*ze;
  N[9] = 4*et*ze;
}




//--------------------------------------------------------------------------------
// Derivatives of shape functions (local coordinates)
//--------------------------------------------------------------------------------
void Point0D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for 0D points at node (xi=1,et=1,ze=1)
  dN[0][0] = 0.0;
}


void LineLin1D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for linear 1D lines at node (xi,et=1,ze=1)
  dN[0][0] = -0.5; dN[0][1] = 0.5;
}


void TriLin2D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for linear 2D triangles at node (xi,et,ze=1)
  dN[0][0] = 1; dN[0][1] = 0; dN[0][2] = -1;
  dN[1][0] = 0; dN[1][1] = 1; dN[1][2] = -1;
}


void SqrLin2D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for linear 2D squares at node (xi,et,ze=1)
  dN[0][0] = -0.25*(1-et); dN[0][1] = 0.25*(1-et); dN[0][2] = 0.25*(1+et); dN[0][3] = -0.25*(1+et);
  dN[1][0] = -0.25*(1-xi); dN[1][1] = -0.25*(1+xi); dN[1][2] = 0.25*(1+xi); dN[1][3] = 0.25*(1-xi);
}


void TetLin3D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for linear 3D tetrahedra at node (xi,et,ze)
  dN[0][0] = 1; dN[0][1] = 0; dN[0][2] = 0; dN[0][3] = -1;
  dN[1][0] = 0; dN[1][1] = 1; dN[1][2] = 0; dN[1][3] = -1;
  dN[2][0] = 0; dN[2][1] = 0; dN[2][2] = 1; dN[2][3] = -1;
}


void LineQuad1D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for linear 1D lines at node (xi,et=1,ze=1)
  dN[0][0] = -0.5+xi; dN[0][1] = 0.5+xi; dN[0][2] = -2.0*xi;
}


void TriQuad2D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for quadratic 2D triangles at node (xi,et,ze=1)
  dN[0][0] = 4*xi-1; dN[0][1] = 0; dN[0][2] = -4*(1-xi-et)+1; dN[0][3] = 4*et; dN[0][4] = -4*et; dN[0][5] = 4*(1-et)-8*xi;
  dN[1][0] = 0; dN[1][1] = 4*et-1; dN[1][2] = -4*(1-xi-et)+1; dN[1][3] = 4*xi; dN[1][4] = 4*(1-xi)-8*et; dN[1][5] = -4*xi;
}


void SqrQuad2D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for quadratic 2D squares at node (xi,et,ze=1)
}


void TetQuad3D::dshape(double xi, double et, double ze)
{
  //derivatives of shape functions for quadratic 3D tetrahedra at node (xi,et,ze)
}



