/*!\file TriaRef.c
 * \brief: implementation of the TriaRef object
 */

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"
/*}}}*/

/*Element macros*/
#define NUMNODES 3

/*Object constructors and destructor*/
/*FUNCTION TriaRef::TriaRef(){{{*/
TriaRef::TriaRef(){
	this->element_type_list=NULL;
}
/*}}}*/
/*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{*/

TriaRef::TriaRef(const int nummodels){

	/*Only allocate pointer*/
	element_type_list=xNew<int>(nummodels);

}
/*}}}*/
/*FUNCTION TriaRef::~TriaRef(){{{*/
TriaRef::~TriaRef(){
	xDelete<int>(element_type_list);
}
/*}}}*/

/*Management*/
/*FUNCTION TriaRef::SetElementType{{{*/
void TriaRef::SetElementType(int type,int type_counter){

	_assert_(type==P1Enum || type==P1DGEnum);

	/*initialize element type*/
	this->element_type_list[type_counter]=type;
}
/*}}}*/

/*Reference Element numerics*/
/*FUNCTION TriaRef::GetBMacAyeal {{{*/
void TriaRef::GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dh/dx           0    ]
	 *          [   0           dh/dy  ]
	 *          [ 1/2*dh/dy  1/2*dh/dx ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
	 */

	int i;
	IssmDouble dbasis[NDOF2][NUMNODES];

	/*Get dh1dh2dh3 in actual coordinate system: */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build B: */
	for (i=0;i<NUMNODES;i++){
		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
		*(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dbasis[1][i]; 
		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetBMacAyealStokes {{{*/
void TriaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){

	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[   dh/dx         0     ]
	 *          [       0       dh/dy   ]
	 *          [  1/2*dh/dy  1/2*dh/dx ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*NUMNODES)
	 */

	/*Same thing in the actual coordinate system: */
	IssmDouble dbasis[NDOF2][NUMNODES];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build B': */
	for (int i=0;i<NUMNODES;i++){
		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0; 
		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
		*(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dbasis[1][i]; 
		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dbasis[0][i]; 
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentBFlux{{{*/
void TriaRef::GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2){
	/*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
	 *
	 * and phi1=phi3 phi2=phi4
	 *
	 * We assume B has been allocated already, of size: 1x4
	 */

	IssmDouble l1l3[NUMNODES];

	GetNodalFunctions(&l1l3[0],gauss);

	B[0] = +l1l3[index1];
	B[1] = +l1l3[index2];
	B[2] = -l1l3[index1];
	B[3] = -l1l3[index2];
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentBprimeFlux{{{*/
void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2){
	/*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
	 *
	 * and phi1=phi3 phi2=phi4
	 *
	 * We assume Bprime has been allocated already, of size: 1x4
	 */

	IssmDouble l1l3[NUMNODES];

	GetNodalFunctions(&l1l3[0],gauss);

	Bprime[0] = l1l3[index1];
	Bprime[1] = l1l3[index2];
	Bprime[2] = l1l3[index1];
	Bprime[3] = l1l3[index2];
}
/*}}}*/
/*FUNCTION TriaRef::GetBPrognostic{{{*/
void TriaRef::GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ h ]
	 *          [ h ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*NUMNODES)
	 */

	IssmDouble basis[NUMNODES];

	/*Get dh1dh2dh3 in actual coordinate system: */
	GetNodalFunctions(&basis[0],gauss);

	/*Build B_prog: */
	for (int i=0;i<NUMNODES;i++){
		*(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=basis[i];
		*(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=basis[i];
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimeMacAyeal {{{*/
void TriaRef::GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){

	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_prime=[ 2*dh/dx    dh/dy ]
	 *                [   dh/dx  2*dh/dy ]
	 *                [   dh/dy    dh/dx ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
	 */

	/*Same thing in the actual coordinate system: */
	IssmDouble dbasis[NDOF2][NUMNODES];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build B': */
	for (int i=0;i<NUMNODES;i++){
		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dbasis[0][i]; 
		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dbasis[1][i]; 
		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dbasis[0][i]; 
		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dbasis[1][i]; 
		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/
void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){

	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
	 * For node i, Bprimei can be expressed in the actual coordinate system
	 * by: 
	 *       Bprimei=[  dh/dx    0   ]
	 *               [    0    dh/dy ]
	 *               [  dh/dy  dh/dx ]
	 *               [  dh/dx  dh/dy ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*NUMNODES)
	 */

	/*Same thing in the actual coordinate system: */
	IssmDouble dbasis[NDOF2][NUMNODES];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build Bprime: */
	for (int i=0;i<NUMNODES;i++){
		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=0; 
		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dbasis[0][i]; 
		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dbasis[1][i]; 
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimePrognostic{{{*/
void TriaRef::GetBprimePrognostic(IssmDouble* Bprime_prog, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_prime=[ dh/dx ]
	 *                [ dh/dy ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B' has been allocated already, of size: 3x(NDOF2*NUMNODES)
	 */

	/*Same thing in the actual coordinate system: */
	IssmDouble dbasis[NDOF2][NUMNODES];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build B': */
	for (int i=0;i<NUMNODES;i++){
		*(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dbasis[0][i]; 
		*(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dbasis[1][i]; 
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetL{{{*/
void TriaRef::GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof){
	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
	 * For node i, Li can be expressed in the actual coordinate system
	 * by: 
	 *       numdof=1: 
	 *                 Li=h;
	 *       numdof=2:
	 *                 Li=[ h   0 ]
	 *                    [ 0   h ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume L has been allocated already, of size: NUMNODES (numdof=1), or numdofx(numdof*NUMNODES) (numdof=2)
	 */

	int i;
	IssmDouble basis[3];

	/*Get basis in actual coordinate system: */
	GetNodalFunctions(basis,gauss);

	/*Build L: */
	if(numdof==1){
		for (i=0;i<NUMNODES;i++){
			L[i]=basis[i]; 
		}
	}
	else{
		for (i=0;i<NUMNODES;i++){
			*(L+numdof*NUMNODES*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
			*(L+numdof*NUMNODES*0+numdof*i+1)=0;
			*(L+numdof*NUMNODES*1+numdof*i)=0;
			*(L+numdof*NUMNODES*1+numdof*i+1)=basis[i];
		}
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetJacobian{{{*/
void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
	IssmDouble x1,y1,x2,y2,x3,y3;

	x1=*(xyz_list+NUMNODES*0+0);
	y1=*(xyz_list+NUMNODES*0+1);
	x2=*(xyz_list+NUMNODES*1+0);
	y2=*(xyz_list+NUMNODES*1+1);
	x3=*(xyz_list+NUMNODES*2+0);
	y3=*(xyz_list+NUMNODES*2+1);


	*(J+NDOF2*0+0)=0.5*(x2-x1);
	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
	*(J+NDOF2*0+1)=0.5*(y2-y1);
	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{*/
void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated*/
	IssmDouble x1,y1,x2,y2;

	x1=*(xyz_list+3*0+0);
	y1=*(xyz_list+3*0+1);
	x2=*(xyz_list+3*1+0);
	y2=*(xyz_list+3*1+1);

	*Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.));
	if(*Jdet<0) _error2_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION TriaRef::GetJacobianDeterminant2d{{{*/
void TriaRef::GetJacobianDeterminant2d(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
	IssmDouble J[2][2];

	/*Get Jacobian*/
	GetJacobian(&J[0][0],xyz_list,gauss);

	/*Get Determinant*/
	Matrix2x2Determinant(Jdet,&J[0][0]);
	if(*Jdet<0) _error2_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION TriaRef::GetJacobianDeterminant3d {{{*/
void TriaRef::GetJacobianDeterminant3d(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	IssmDouble x1,x2,x3,y1,y2,y3,z1,z2,z3;

	x1=*(xyz_list+3*0+0);
	y1=*(xyz_list+3*0+1);
	z1=*(xyz_list+3*0+2);
	x2=*(xyz_list+3*1+0);
	y2=*(xyz_list+3*1+1);
	z2=*(xyz_list+3*1+2);
	x3=*(xyz_list+3*2+0);
	y3=*(xyz_list+3*2+1);
	z3=*(xyz_list+3*2+2);

	*Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
	if(*Jdet<0) _error2_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION TriaRef::GetJacobianInvert{{{*/
void TriaRef::GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss){

	/*Jacobian*/
	IssmDouble J[2][2];

	/*Call Jacobian routine to get the jacobian:*/
	GetJacobian(&J[0][0], xyz_list, gauss);

	/*Invert Jacobian matrix: */
	Matrix2x2Invert(Jinv,&J[0][0]);

}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctions{{{*/
void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	basis[0]=gauss->coord1;
	basis[1]=gauss->coord2;
	basis[2]=gauss->coord3;

}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	IssmDouble BasisFunctions[3];

	GetNodalFunctions(&BasisFunctions[0],gauss);

	_assert_(index1>=0 && index1<3);
	_assert_(index2>=0 && index2<3);
	basis[0]=BasisFunctions[index1];
	basis[1]=BasisFunctions[index2];
}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{*/
void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * actual coordinate system): */
	int       i;
	IssmDouble    dbasis_ref[NDOF2][NUMNODES];
	IssmDouble    Jinv[NDOF2][NDOF2];

	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss); 

	/*Get Jacobian invert: */
	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);

	/*Build dbasis: 
	 *
	 * [dhi/dx]= Jinv*[dhi/dr]
	 * [dhi/dy]       [dhi/ds]
	 */
	for (i=0;i<NUMNODES;i++){
		dbasis[NUMNODES*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
		dbasis[NUMNODES*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
	}

}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dl1dl3,GaussTria* gauss){
	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */

	/*First nodal function: */
	*(dl1dl3+NUMNODES*0+0)=-0.5; 
	*(dl1dl3+NUMNODES*1+0)=-1.0/(2.0*SQRT3);

	/*Second nodal function: */
	*(dl1dl3+NUMNODES*0+1)=0.5;
	*(dl1dl3+NUMNODES*1+1)=-1.0/(2.0*SQRT3);

	/*Third nodal function: */
	*(dl1dl3+NUMNODES*0+2)=0;
	*(dl1dl3+NUMNODES*1+2)=1.0/SQRT3;

}
/*}}}*/
/*FUNCTION TriaRef::GetInputDerivativeValue{{{*/
void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){

	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
	 * point specified by gauss_basis:
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *
	 * p is a vector of size 2x1 already allocated.
	 */

	/*Nodal Derivatives*/
	IssmDouble dbasis[2][3]; //nodal derivative functions in actual coordinate system.

	/*Get dh1dh2dh3 in actual coordinate system: */
	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss);

	/*Assign values*/
	*(p+0)=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2];
	*(p+1)=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2];

}
/*}}}*/
/*FUNCTION TriaRef::GetInputValue{{{*/
void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){

	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 
	 * point specifie by gauss: */

	/*nodal functions annd output: */
	IssmDouble basis[3];

	/*Get nodal functions*/
	GetNodalFunctions(basis, gauss);

	/*Get parameter*/
	*p=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2];
}
/*}}}*/
