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

#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 "../../DataSet/DataSet.h"
#include "../../include/include.h"

/*Object constructors and destructor*/
/*FUNCTION PentaVertexInput::PentaVertexInput(){{{1*/
PentaVertexInput::PentaVertexInput(){
	return;
}
/*}}}*/
/*FUNCTION PentaVertexInput::PentaVertexInput(int in_enum_type,double* values){{{1*/
PentaVertexInput::PentaVertexInput(int in_enum_type,double* in_values){

	enum_type=in_enum_type;
	values[0]=in_values[0];
	values[1]=in_values[1];
	values[2]=in_values[2];
	values[3]=in_values[3];
	values[4]=in_values[4];
	values[5]=in_values[5];
}
/*}}}*/
/*FUNCTION PentaVertexInput::~PentaVertexInput(){{{1*/
PentaVertexInput::~PentaVertexInput(){
	return;
}
/*}}}*/

/*Object management*/
/*FUNCTION PentaVertexInput::copy{{{1*/
Object* PentaVertexInput::copy() {
	
	return new PentaVertexInput(this->enum_type,this->values);

}
/*}}}*/
/*FUNCTION PentaVertexInput::DeepEcho{{{1*/
void PentaVertexInput::DeepEcho(void){

	printf("PentaVertexInput:\n");
	printf("   enum: %i (%s)\n",this->enum_type,EnumAsString(this->enum_type));
	printf("   values: [%g %g %g %g %g %g]\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]);
}
/*}}}*/
/*FUNCTION PentaVertexInput::Demarshall{{{1*/
void  PentaVertexInput::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   i;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/
	memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	memcpy(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values);

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
/*FUNCTION PentaVertexInput::Echo {{{1*/
void PentaVertexInput::Echo(void){
	this->DeepEcho();
}
/*}}}*/
/*FUNCTION PentaVertexInput::Enum{{{1*/
int PentaVertexInput::Enum(void){

	return PentaVertexInputEnum;

}
/*}}}*/
/*FUNCTION PentaVertexInput::EnumType{{{1*/
int PentaVertexInput::EnumType(void){

	return this->enum_type;

}
/*}}}*/
/*FUNCTION PentaVertexInput::Id{{{1*/
int    PentaVertexInput::Id(void){ return -1; }
/*}}}*/
/*FUNCTION PentaVertexInput::Marshall{{{1*/
void  PentaVertexInput::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_value=0;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum value of PentaVertexInput: */
	enum_value=PentaVertexInputEnum;
	
	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_value,sizeof(enum_value));marshalled_dataset+=sizeof(enum_value);
	
	/*marshall PentaVertexInput data: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values);

	*pmarshalled_dataset=marshalled_dataset;
}
/*}}}*/
/*FUNCTION PentaVertexInput::MarshallSize{{{1*/
int   PentaVertexInput::MarshallSize(){
	
	return sizeof(values)+
		+sizeof(enum_type)+
		+sizeof(int); //sizeof(int) for enum value
}
/*}}}*/
/*FUNCTION PentaVertexInput::MyRank{{{1*/
int    PentaVertexInput::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
/*FUNCTION PentaVertexInput::SpawnSingInput{{{1*/
Input* PentaVertexInput::SpawnSingInput(int index){

	/*output*/
	SingVertexInput* outinput=NULL;

	/*Create new Sing input (copy of current input)*/
	ISSMASSERT(index<6 && index>=0);
	outinput=new SingVertexInput(this->enum_type,this->values[index]);

	/*Assign output*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaVertexInput::SpawnBeamInput{{{1*/
Input* PentaVertexInput::SpawnBeamInput(int* indices){

	/*output*/
	BeamVertexInput* outinput=NULL;
	double newvalues[2];

	/*Loop over the new indices*/
	for(int i=0;i<2;i++){

		/*Check index value*/
		ISSMASSERT(indices[i]>=0 && indices[i]<6);

		/*Assign value to new input*/
		newvalues[i]=this->values[indices[i]];
	}

	/*Create new Beam input*/
	outinput=new BeamVertexInput(this->enum_type,&newvalues[0]);

	/*Assign output*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaVertexInput::SpawnTriaInput{{{1*/
Input* PentaVertexInput::SpawnTriaInput(int* indices){

	/*output*/
	TriaVertexInput* outinput=NULL;
	double newvalues[3];

	/*Loop over the new indices*/
	for(int i=0;i<3;i++){

		/*Check index value*/
		ISSMASSERT(indices[i]>=0 && indices[i]<6);

		/*Assign value to new input*/
		newvalues[i]=this->values[indices[i]];
	}

	/*Create new Tria input*/
	outinput=new TriaVertexInput(this->enum_type,&newvalues[0]);

	/*Assign output*/
	return outinput;

}
/*}}}*/
/*FUNCTION PentaVertexInput::SpawnResult{{{1*/
Result* PentaVertexInput::SpawnResult(int step, double time){

	return new PentaVertexResult(this->enum_type,this->values,step,time);

}
/*}}}*/

/*Object functions*/
/*FUNCTION PentaVertexInput::GetParameterValue(bool* pvalue) {{{1*/
void PentaVertexInput::GetParameterValue(bool* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(int* pvalue){{{1*/
void PentaVertexInput::GetParameterValue(int* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,Node* node){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss){
	/*P1 interpolation on Gauss point*/

	/*intermediary*/
	double l1l6[6];

	/*nodal functions: */
	GetNodalFunctionsP1(&l1l6[0],gauss);

	/*Assign output pointers:*/
	*pvalue=l1l6[0]*values[0]+l1l6[1]*values[1]+l1l6[2]*values[2]+l1l6[3]*values[3]+l1l6[4]*values[4]+l1l6[5]*values[5];

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
void PentaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){ISSMERROR(" not supported yet!");}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){{{1*/
void PentaVertexInput::GetParameterValues(double* values,double* gauss_pointers, int numgauss){
	/*It is assumed that output values has been correctly allocated*/

	int i,j;
	double gauss[4];

	for (i=0;i<numgauss;i++){

		/*Get current Gauss point coordinates*/
		for (j=0;j<4;j++) gauss[j]=gauss_pointers[i*4+j];

		/*Assign parameter value*/
		GetParameterValue(&values[i],&gauss[0]);
	}
}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss){{{1*/
void PentaVertexInput::GetParameterDerivativeValue(double* p, double* xyz_list, double* gauss){
	/*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_coord:
	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
	 *
	 *   p is a vector of size 3x1 already allocated.
	 */

	const int NDOF3=3;
	const int numgrids=6;
	double dh1dh6[NDOF3][numgrids];

	/*Get nodal funnctions derivatives in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);

	p[0]=this->values[0]*dh1dh6[0][0]+this->values[1]*dh1dh6[0][1]+this->values[2]*dh1dh6[0][2]+this->values[3]*dh1dh6[0][3]+this->values[4]*dh1dh6[0][4]+this->values[5]*dh1dh6[0][5];
	p[1]=this->values[0]*dh1dh6[1][0]+this->values[1]*dh1dh6[1][1]+this->values[2]*dh1dh6[1][2]+this->values[3]*dh1dh6[1][3]+this->values[4]*dh1dh6[1][4]+this->values[5]*dh1dh6[1][5];
	p[2]=this->values[0]*dh1dh6[2][0]+this->values[1]*dh1dh6[2][1]+this->values[2]*dh1dh6[2][2]+this->values[3]*dh1dh6[2][3]+this->values[4]*dh1dh6[2][4]+this->values[5]*dh1dh6[2][5];

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
void PentaVertexInput::GetVxStrainRate3d(double* epsilonvx,double* xyz_list, double* gauss){
	int i,j;

	const int numgrids=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numgrids];
	double velocity[numgrids][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

	/*Here, we are computing the strain rate of (vx,0,0)*/
	for(i=0;i<numgrids;i++){
		velocity[i][0]=this->values[i];
		velocity[i][1]=0.0;
		velocity[i][2]=0.0;
	}
	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
void PentaVertexInput::GetVyStrainRate3d(double* epsilonvy,double* xyz_list, double* gauss){
	int i,j;

	const int numgrids=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numgrids];
	double velocity[numgrids][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

	/*Here, we are computing the strain rate of (0,vy,0)*/
	for(i=0;i<numgrids;i++){
		velocity[i][0]=0.0;
		velocity[i][1]=this->values[i];
		velocity[i][2]=0.0;
	}
	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss) {{{1*/
void PentaVertexInput::GetVzStrainRate3d(double* epsilonvz,double* xyz_list, double* gauss){
	int i,j;

	const int numgrids=6;
	const int DOFVELOCITY=3;
	double B[8][27];
	double B_reduced[6][DOFVELOCITY*numgrids];
	double velocity[numgrids][DOFVELOCITY];

	/*Get B matrix: */
	GetBStokes(&B[0][0], xyz_list, gauss);
	/*Create a reduced matrix of B to get rid of pressure */
	for (i=0;i<6;i++){
		for (j=0;j<3;j++){
			B_reduced[i][j]=B[i][j];
		}
		for (j=4;j<7;j++){
			B_reduced[i][j-1]=B[i][j];
		}
		for (j=8;j<11;j++){
			B_reduced[i][j-2]=B[i][j];
		}
		for (j=12;j<15;j++){
			B_reduced[i][j-3]=B[i][j];
		}
		for (j=16;j<19;j++){
			B_reduced[i][j-4]=B[i][j];
		}
		for (j=20;j<23;j++){
			B_reduced[i][j-5]=B[i][j];
		}
	}

	/*Here, we are computing the strain rate of (0,0,vz)*/
	for(i=0;i<numgrids;i++){
		velocity[i][0]=0.0;
		velocity[i][1]=0.0;
		velocity[i][2]=this->values[i];
	}

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply(&B_reduced[0][0],6,DOFVELOCITY*numgrids,0,&velocity[0][0],DOFVELOCITY*numgrids,1,0,epsilonvz,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss) {{{1*/
void PentaVertexInput::GetVxStrainRate3dPattyn(double* epsilonvx,double* xyz_list, double* gauss){

	int i;
	const int numgrids=6;
	const int NDOF2=2;
	double B[5][NDOF2*numgrids];
	double velocity[numgrids][NDOF2];

	/*Get B matrix: */
	GetBPattyn(&B[0][0], xyz_list, gauss);

	/*Here, we are computing the strain rate of (vx,0)*/
	for(i=0;i<numgrids;i++){
		velocity[i][0]=this->values[i];
		velocity[i][1]=0.0;
	}

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
				&velocity[0][0],NDOF2*numgrids,1,0,
				epsilonvx,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss) {{{1*/
void PentaVertexInput::GetVyStrainRate3dPattyn(double* epsilonvy,double* xyz_list, double* gauss){

	int i;
	const int numgrids=6;
	const int NDOF2=2;
	double B[5][NDOF2*numgrids];
	double velocity[numgrids][NDOF2];

	/*Get B matrix: */
	GetBPattyn(&B[0][0], xyz_list, gauss);

	/*Here, we are computing the strain rate of (0,vy)*/
	for(i=0;i<numgrids;i++){
		velocity[i][0]=0.0;
		velocity[i][1]=this->values[i];
	}

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
				&velocity[0][0],NDOF2*numgrids,1,0,
				epsilonvy,0);

}
/*}}}*/
/*FUNCTION PentaVertexInput::ChangeEnum(int newenumtype){{{1*/
void PentaVertexInput::ChangeEnum(int newenumtype){
	this->enum_type=newenumtype;
}
/*}}}*/
/*FUNCTION PentaVertexInput::GetParameterAverage(double* pvalue){{{1*/
void PentaVertexInput::GetParameterAverage(double* pvalue){
	*pvalue=1./6.*(values[0]+values[1]+values[2]+values[3]+values[4]+values[5]);
}
/*}}}*/

/*Intermediary*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsP1 {{{1*/
void PentaVertexInput::GetNodalFunctionsP1(double* l1l6, double* gauss_coord){

	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	l1l6[0]=gauss_coord[0]*(1-gauss_coord[3])/2.0;

	l1l6[1]=gauss_coord[1]*(1-gauss_coord[3])/2.0;

	l1l6[2]=gauss_coord[2]*(1-gauss_coord[3])/2.0;

	l1l6[3]=gauss_coord[0]*(1+gauss_coord[3])/2.0;

	l1l6[4]=gauss_coord[1]*(1+gauss_coord[3])/2.0;

	l1l6[5]=gauss_coord[2]*(1+gauss_coord[3])/2.0;

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsMINI{{{1*/
void PentaVertexInput::GetNodalFunctionsMINI(double* l1l7, double* gauss_coord){

	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	/*First nodal function: */
	l1l7[0]=gauss_coord[0]*(1.0-gauss_coord[3])/2.0;

	/*Second nodal function: */
	l1l7[1]=gauss_coord[1]*(1.0-gauss_coord[3])/2.0;

	/*Third nodal function: */
	l1l7[2]=gauss_coord[2]*(1.0-gauss_coord[3])/2.0;

	/*Fourth nodal function: */
	l1l7[3]=gauss_coord[0]*(1.0+gauss_coord[3])/2.0;

	/*Fifth nodal function: */
	l1l7[4]=gauss_coord[1]*(1.0+gauss_coord[3])/2.0;

	/*Sixth nodal function: */
	l1l7[5]=gauss_coord[2]*(1.0+gauss_coord[3])/2.0;

	/*Seventh nodal function: */
	l1l7[6]=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(1.0+gauss_coord[3])*(1.0-gauss_coord[3]);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsP1Derivatives {{{1*/
void PentaVertexInput::GetNodalFunctionsP1Derivatives(double* dh1dh6,double* xyz_list, double* gauss_coord){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the actual coordinate system: */
	int i;
	const int NDOF3=3;
	const int numgrids=6;

	double dh1dh6_ref[NDOF3][numgrids];
	double Jinv[NDOF3][NDOF3];

	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsP1DerivativesReference(&dh1dh6_ref[0][0], gauss_coord); 

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

	/*Build dh1dh3: 
	 *
	 * [dhi/dx]= Jinv*[dhi/dr]
	 * [dhi/dy]       [dhi/ds]
	 * [dhi/dz]       [dhi/dn]
	 */

	for (i=0;i<numgrids;i++){
		*(dh1dh6+numgrids*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
		*(dh1dh6+numgrids*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
		*(dh1dh6+numgrids*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
	}

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivatives{{{1*/
void PentaVertexInput::GetNodalFunctionsMINIDerivatives(double* dh1dh7,double* xyz_list, double* gauss_coord){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * actual coordinate system: */

	int i;

	const  int numgrids=7;
	double dh1dh7_ref[3][numgrids];
	double Jinv[3][3];


	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsMINIDerivativesReference(&dh1dh7_ref[0][0], gauss_coord); 

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

	/*Build dh1dh6: 
	 *
	 * [dhi/dx]= Jinv'*[dhi/dr]
	 * [dhi/dy]        [dhi/ds]
	 * [dhi/dz]        [dhi/dzeta]
	 */

	for (i=0;i<numgrids;i++){
		*(dh1dh7+numgrids*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
		*(dh1dh7+numgrids*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
		*(dh1dh7+numgrids*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
	}

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsP1DerivativesReference {{{1*/
void PentaVertexInput::GetNodalFunctionsP1DerivativesReference(double* dl1dl6,double* gauss_coord){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */

	const int numgrids=6;
	double A1,A2,A3,z;

	A1=gauss_coord[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
	A2=gauss_coord[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
	A3=gauss_coord[2]; //third area coordinate value  In term of xi and eta: A3=y/SQRT3;
	z=gauss_coord[3]; //fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2


	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+0)=-0.5*(1.0-z)/2.0;
	*(dl1dl6+numgrids*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+numgrids*2+0)=-0.5*A1;

	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+1)=0.5*(1.0-z)/2.0;
	*(dl1dl6+numgrids*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+numgrids*2+1)=-0.5*A2;

	/*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+2)=0.0;
	*(dl1dl6+numgrids*1+2)=1.0/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+numgrids*2+2)=-0.5*A3;

	/*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+3)=-0.5*(1.0+z)/2.0;
	*(dl1dl6+numgrids*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+numgrids*2+3)=0.5*A1;

	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+4)=0.5*(1.0+z)/2.0;
	*(dl1dl6+numgrids*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+numgrids*2+4)=0.5*A2;

	/*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+numgrids*0+5)=0.0;
	*(dl1dl6+numgrids*1+5)=1.0/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+numgrids*2+5)=0.5*A3;
}
/*}}}*/
/*FUNCTION PentaVertexInput::GetNodalFunctionsMINIDerivativesReference{{{1*/
void PentaVertexInput::GetNodalFunctionsMINIDerivativesReference(double* dl1dl7,double* gauss_coord){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */

	int    numgrids=7; //six plus bubble grids

	double r=gauss_coord[1]-gauss_coord[0];
	double s=-3.0/SQRT3*(gauss_coord[0]+gauss_coord[1]-2.0/3.0);
	double zeta=gauss_coord[3];

	/*First nodal function: */
	*(dl1dl7+numgrids*0+0)=-0.5*(1.0-zeta)/2.0;
	*(dl1dl7+numgrids*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
	*(dl1dl7+numgrids*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Second nodal function: */
	*(dl1dl7+numgrids*0+1)=0.5*(1.0-zeta)/2.0;
	*(dl1dl7+numgrids*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
	*(dl1dl7+numgrids*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Third nodal function: */
	*(dl1dl7+numgrids*0+2)=0;
	*(dl1dl7+numgrids*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
	*(dl1dl7+numgrids*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);

	/*Fourth nodal function: */
	*(dl1dl7+numgrids*0+3)=-0.5*(1.0+zeta)/2.0;
	*(dl1dl7+numgrids*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
	*(dl1dl7+numgrids*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Fith nodal function: */
	*(dl1dl7+numgrids*0+4)=0.5*(1.0+zeta)/2.0;
	*(dl1dl7+numgrids*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
	*(dl1dl7+numgrids*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Sixth nodal function: */
	*(dl1dl7+numgrids*0+5)=0;
	*(dl1dl7+numgrids*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
	*(dl1dl7+numgrids*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);

	/*Seventh nodal function: */
	*(dl1dl7+numgrids*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
	*(dl1dl7+numgrids*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
	*(dl1dl7+numgrids*2+6)=27*gauss_coord[0]*gauss_coord[1]*gauss_coord[2]*(-2.0*zeta);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetJacobian {{{1*/
void PentaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss_coord){

	const int NDOF3=3;
	int i,j;

	/*The Jacobian is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	double A1,A2,A3; //area coordinates
	double xi,eta,zi; //parametric coordinates

	double x1,x2,x3,x4,x5,x6;
	double y1,y2,y3,y4,y5,y6;
	double z1,z2,z3,z4,z5,z6;

	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
	A1=gauss_coord[0];
	A2=gauss_coord[1];
	A3=gauss_coord[2];

	xi=A2-A1;
	eta=SQRT3*A3;
	zi=gauss_coord[3];

	x1=*(xyz_list+3*0+0);
	x2=*(xyz_list+3*1+0);
	x3=*(xyz_list+3*2+0);
	x4=*(xyz_list+3*3+0);
	x5=*(xyz_list+3*4+0);
	x6=*(xyz_list+3*5+0);

	y1=*(xyz_list+3*0+1);
	y2=*(xyz_list+3*1+1);
	y3=*(xyz_list+3*2+1);
	y4=*(xyz_list+3*3+1);
	y5=*(xyz_list+3*4+1);
	y6=*(xyz_list+3*5+1);

	z1=*(xyz_list+3*0+2);
	z2=*(xyz_list+3*1+2);
	z3=*(xyz_list+3*2+2);
	z4=*(xyz_list+3*3+2);
	z5=*(xyz_list+3*4+2);
	z6=*(xyz_list+3*5+2);

	*(J+NDOF3*0+0)=0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
	*(J+NDOF3*1+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
	*(J+NDOF3*2+0)=SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);

	*(J+NDOF3*0+1)=0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
	*(J+NDOF3*1+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
	*(J+NDOF3*2+1)=SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);

	*(J+NDOF3*0+2)=0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
	*(J+NDOF3*1+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
	*(J+NDOF3*2+2)=SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetJacobianInvert {{{1*/
void PentaVertexInput::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_coord){

	double Jdet;
	const int NDOF3=3;

	/*Call Jacobian routine to get the jacobian:*/
	GetJacobian(Jinv, xyz_list, gauss_coord);

	/*Invert Jacobian matrix: */
	MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
}
/*}}}*/
/*FUNCTION PentaVertexInput::GetBPattyn {{{1*/
void PentaVertexInput::GetBPattyn(double* B, double* xyz_list, double* gauss_coord){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
	 * For grid 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  ]
	 *          [ 1/2*dh/dz      0      ]
	 *          [  0         1/2*dh/dz  ]
	 * where h is the interpolation function for grid i.
	 *
	 * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
	 */

	int i;
	const int numgrids=6;
	const int NDOF3=3;
	const int NDOF2=2;

	double dh1dh6[NDOF3][numgrids];

	/*Get dh1dh6 in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);

	/*Build B: */
	for (i=0;i<numgrids;i++){
		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh6[0][i]; 
		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0.0;

		*(B+NDOF2*numgrids*1+NDOF2*i)=0.0;
		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh6[1][i];

		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 

		*(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 
		*(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;

		*(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
		*(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 
	}

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetBStokes {{{1*/
void PentaVertexInput::GetBStokes(double* B, double* xyz_list, double* gauss_coord){

	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*DOFPERGRID. 
	 * For grid i, Bi can be expressed in the actual coordinate system
	 * by: 		Bi=[ dh/dx          0             0       0  ]
	 *					[   0           dh/dy           0       0  ]
	 *					[   0             0           dh/dy     0  ]
	 *					[ 1/2*dh/dy    1/2*dh/dx        0       0  ]
	 *					[ 1/2*dh/dz       0         1/2*dh/dx   0  ]
	 *					[   0          1/2*dh/dz    1/2*dh/dy   0  ]
	 *					[   0             0             0       h  ]
	 *					[ dh/dx         dh/dy         dh/dz     0  ]
	 *	where h is the interpolation function for grid i.
	 *	Same thing for Bb except the last column that does not exist.
	 */

	int i;
	const int calculationdof=3;
	const int numgrids=6;
	int DOFPERGRID=4;

	double dh1dh7[calculationdof][numgrids+1];
	double l1l6[numgrids];


	/*Get dh1dh7 in actual coordinate system: */
	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss_coord);
	GetNodalFunctionsP1(l1l6, gauss_coord);

	/*Build B: */
	for (i=0;i<numgrids+1;i++){
		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i)=dh1dh7[0][i]; //B[0][DOFPERGRID*i]=dh1dh6[0][i];
		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+1)=0;
		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+2)=0;
		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i)=0;
		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+1)=dh1dh7[1][i];
		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+2)=0;
		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i)=0;
		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+1)=0;
		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+2)=dh1dh7[2][i];
		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i)=(float).5*dh1dh7[1][i]; 
		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+1)=(float).5*dh1dh7[0][i]; 
		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+2)=0;
		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i)=(float).5*dh1dh7[2][i];
		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+1)=0;
		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+2)=(float).5*dh1dh7[0][i];
		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i)=0;
		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+1)=(float).5*dh1dh7[2][i];
		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+2)=(float).5*dh1dh7[1][i];
		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i)=0;
		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+1)=0;
		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+2)=0;
		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i)=dh1dh7[0][i];
		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+1)=dh1dh7[1][i];
		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+2)=dh1dh7[2][i];
	}

	for (i=0;i<numgrids;i++){ //last column not for the bubble function
		*(B+(DOFPERGRID*numgrids+3)*0+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*1+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*2+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*3+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*4+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*5+DOFPERGRID*i+3)=0;
		*(B+(DOFPERGRID*numgrids+3)*6+DOFPERGRID*i+3)=l1l6[i];
		*(B+(DOFPERGRID*numgrids+3)*7+DOFPERGRID*i+3)=0;
	}

}
/*}}}*/
/*FUNCTION PentaVertexInput::SquareMin(double* psquaremin, bool process_units){{{1*/
void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){

	int i;
	const int numnodes=6;
	double valuescopy[numnodes];
	double squaremin;

	/*First,  copy values, to process units if requested: */
	for(i=0;i<numnodes;i++)valuescopy[i]=this->values[i];

	/*Process units if requested: */
	if(process_units)NodalValuesUnitConversion(&valuescopy[0],numnodes,enum_type,parameters);

	/*Now, figure out minimum of valuescopy: */
	squaremin=pow(valuescopy[0],2);
	for(i=1;i<numnodes;i++){
		if(pow(valuescopy[i],2)<squaremin)squaremin=pow(valuescopy[i],2);
	}
	/*Assign output pointers:*/
	*psquaremin=squaremin;
}
/*}}}*/
/*FUNCTION PentaVertexInput::Scale(double scale_factor){{{1*/
void PentaVertexInput::Scale(double scale_factor){
	
	int i;
	const int numgrids=6;

	for(i=0;i<numgrids;i++)values[i]=values[i]*scale_factor;
}
/*}}}*/
/*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
void PentaVertexInput::AXPY(Input* xinput,double scalar){

	int i;
	const int numgrids=6;
	PentaVertexInput*  xpentavertexinput=NULL;

	/*xinput is of the same type, so cast it: */
	xpentavertexinput=(PentaVertexInput)xinput;

	/*Carry out the AXPY operation:*/
	for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];

}
/*}}}*/
/*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
void PentaVertexInput::Constrain(double cm_min, double cm_max){

	int i;
	const int numgrids=6;
		
	if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
	if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;

}
/*}}}*/
/*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){

	const int numvertices=6;
	VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);


}
/*}}}*/