#include <cstdio>
#include <string.h>
#include <cmath>

#include "../../include/include.h"
#include "../../shared/Exceptions/exceptions.h"

#include "GeometricalEdge.h"
#include "Geometry.h"

using namespace std;

namespace bamg {

	/*Constructor/Destructor*/

	/*Methods*/
	/*FUNCTION GeometricalEdge::Cracked{{{1*/
	int    GeometricalEdge::Cracked() const  {
		return type &1;  
	}/*}}}*/
	/*FUNCTION GeometricalEdge::F{{{1*/
	R2 GeometricalEdge::F(double theta) const{
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/F)*/
		// parametrization of the curve edge

	   R2 A=v[0]->r,B=v[1]->r;
		double ca,cb,cta,ctb;

		//Check that theta is in [0 1]
		ISSMASSERT(theta>-1e-12 && theta<1+1e-12);

		if (TgA()) 
		 if (TgB()) // interpolation d'hermite
			{ cb =  theta*theta*(3-2*theta);
			 ca =  1-cb;     
			 cta = (1-theta)*(1-theta)*theta;
			 ctb = (theta-1)*theta*theta ;
			}
		 else { // 1-t*t, t-t*t, t*t
			 double t = theta;
			 cb = t*t;
			 ca = 1-cb;
			 cta= t-cb;
			 ctb=0;    
		 }    
		else
		 if (TgB()){
			 double t = 1-theta;
			 ca = t*t;
			 cb = 1-ca;
			 ctb= -t+ca;
			 cta=0;    
		 }
		 else {
			 ca =(1-theta),cb = theta,cta=ctb=0; // lagrange P1
		 }
		return A*ca + B*cb + tg[0]* cta + tg[1] * ctb;
	  }
	/*}}}1*/
	/*FUNCTION GeometricalEdge::Mark{{{1*/
	int    GeometricalEdge::Mark()    const  {
		return type &16; 
	}/*}}}*/
	/*FUNCTION GeometricalEdge::R1tg{{{1*/
	double GeometricalEdge::R1tg(double theta,R2 & t) const{
		/*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/
		// 1/R of radius of cuvature

		R2 A=v[0]->r,B=v[1]->r;
		double dca,dcb,dcta,dctb;
		double ddca,ddcb,ddcta,ddctb;
		double tt = theta*theta;

		//check theta
		ISSMASSERT(theta>=0 && theta<=1);

		if (TgA()){ 
			if (TgB()){
				// Tangent A and B provided:
				// interpolation d'hermite
				dcb = 6*theta*(1-theta);
				ddcb = 6*(1-2*theta);
				dca = -dcb;
				ddca = -ddcb;
				dcta =  (3*theta - 4)*theta + 1;
				ddcta=6*theta-4;
				dctb = 3*tt - 2*theta;
				ddctb = 6*theta-2;
			}
			else {
				//Tangent A provided but tangent B not provided
				// 1-t*t, t-t*t, t*t
				double t = theta;
				dcb = 2*t;
				ddcb = 2;
				dca = -dcb;
				ddca = -2;
				dcta = 1-dcb;
				ddcta = -ddcb;
				dctb=0;    
				ddctb=0;    
			}
		}
		else{
			if (TgB()){
				//Tangent B provided but tangent A not provided
				double t = 1-theta;
				dca = -2*t;
				ddca = 2;
				dcb = -dca;
				ddcb = -2;
				dctb = 1+dca;
				ddctb= ddca;
				dcta =0;
				ddcta =0;
			}
			else {
				//Neither thangent A nor tangent B provided
				// lagrange P1
				t=B-A;
				return 0;
			} 
		}
		R2 d  =  A*dca  + B*dcb  + tg[0]* dcta  + tg[1] * dctb;
		R2 dd =  A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
		double d2=(d,d);
		double sd2 = sqrt(d2);
		t=d;
		if(d2>1.0e-20){
			t/=sd2;
			return Abs(Det(d,dd))/(d2*sd2);
		}
		else return 0;
	}
	/*}}}1*/
	/*FUNCTION GeometricalEdge::Required{{{1*/
	int    GeometricalEdge::Required()       {
		return type &64; 
	}/*}}}*/
	/*FUNCTION GeometricalEdge::Set {{{1*/
	void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew){ 
		*this = rec;
		v[0] = GhNew.vertices + Gh.GetId(v[0]);    
		v[1] = GhNew.vertices + Gh.GetId(v[1]); 
		if (Adj[0]) Adj[0] =  GhNew.edges + Gh.GetId(Adj[0]);     
		if (Adj[1]) Adj[1] =  GhNew.edges + Gh.GetId(Adj[1]);     
	}
	/*}}}*/
	/*FUNCTION GeometricalEdge::SetCracked{{{1*/
	void   GeometricalEdge::SetCracked()     { 
		type |= 1;/*=>1st digit to 1*/
	}/*}}}*/
	/*FUNCTION GeometricalEdge::SetTgA{{{1*/
	void   GeometricalEdge::SetTgA()         { 
		type |=4; /*=>2d digit to 1*/
	}/*}}}*/
	/*FUNCTION GeometricalEdge::SetTgB{{{1*/
	void   GeometricalEdge::SetTgB()         { 
		type |=8; /*=> 3d digit to 1*/
	}/*}}}*/
	/*FUNCTION GeometricalEdge::SetMark{{{1*/
	void   GeometricalEdge::SetMark()        { 
		type |=16;/*=> 4th digiy to 1*/
	}/*}}}*/
	/*FUNCTION GeometricalEdge::SetUnMark{{{1*/
	void   GeometricalEdge::SetUnMark()      { 
		type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/;
	}/*}}}*/
	/*FUNCTION GeometricalEdge::SetRequired{{{1*/
	void   GeometricalEdge::SetRequired()    { 
		type |= 64;/*=>6th digit to 1*/ 
	}/*}}}*/
	  /*FUNCTION GeometricalEdge::Tg{{{1*/
	int    GeometricalEdge::Tg(int i) const  {
		return i==0 ? TgA() : TgB();
	}/*}}}*/
	/*FUNCTION GeometricalEdge::TgA{{{1*/
	int    GeometricalEdge::TgA()     const  {
		return type &4;  
	}/*}}}*/
	/*FUNCTION GeometricalEdge::TgB{{{1*/
	int    GeometricalEdge::TgB()     const  {
		return type &8;  
	}/*}}}*/
}