function pe=CreatePVector(load,elements,grids,materials,inputs,analysis_type);
%CREATEPVECTOR - create the load vector for an icefront
%
%   works for the icefront elements, segments or quad elements
%   this routine can be used for MacAyeal, Pattyn ans Stokes's models
%
%   Usage:
%      pe=CreatePVector(load,elements,grids,materials,inputs,analysis_type)
% 
%   See also PENALTYCREATEPVECTOR, PENALTYCREATEKMATRIX

if strcmpi(analysis_type,'diagnostic_horiz'),

	pe=CreatePVectorHoriz(load,elements,grids,materials,inputs);

elseif strcmpi(analysis_type,'diagnostic_stokes'),

	pe=CreatePVectorStokes(load,elements,grids,materials,inputs);

else
	error('CreatePVector error message: only horizontal diagnostic should use icefront objects');

end
end %end function

function pe=CreatePVectorStokes(load,elements,grids,materials,inputs);

%some variables
numgrids=4;
NDOF4=4;

numgridstria=4; %number of grids for the triangles
NDOF4=4;


%Create elementary vector.
pe=elemvector(numgrids*NDOF4);

%recover element object
pentaelem=elements(load.eid).element;
matpar=materials(end).constants;

%recover material parameters
gravity=matpar.g;
rho_ice=matpar.rho_ice;
rho_water=matpar.rho_water;

%recover extra inputs
[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
[bed_param bed_is_present]=recover_input(inputs,'bed');

%initialize extra inputs
thickness_list=zeros(6,1);
bed_list=zeros(6,1);

%Get all element grid data:
xyz_list=getgriddata(pentaelem,grids);

%Identify which grids are comprised in the segment: 
for i=1:6,
	if pentaelem.g(i)==load.g(1),
		grid1=i;
	end
	if pentaelem.g(i)==load.g(2),
		grid2=i;
	end
	if pentaelem.g(i)==load.g(3),
		grid3=i;
	end
	if pentaelem.g(i)==load.g(4),
		grid4=i;
	end
end
quad_grids=[grid1 grid2 grid3 grid4];

%Build row indices for elementary vector.
for i=1:numgrids,
	doflist=grids(load.g(i)).grid.doflist;
	for j=1:NDOF4,
		pe.row_indices((i-1)*NDOF4+j)=doflist(j);
	end
	dof=doflist(1);
end

for i=1:6,
	doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
	if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
	if(bed_is_present) bed_list(i)=bed_param(dof);end;
end


%Recover thicknesses and bed at all grids
if thickness_is_present,
	thickness=thickness_list(quad_grids);
	bed=bed_list(quad_grids);
else
	thickness=pentaelem.h(quad_grids);
	bed=pentaelem.b(quad_grids);
end

%build xyz list for [grid1, grid2,grid3,grid3]
xyz_list_quad=xyz_list(quad_grids,:);

%Create a new grid in the midle of the quad and add it at the end of the list
x5=mean(xyz_list_quad(:,1));
y5=mean(xyz_list_quad(:,2));
z5=mean(xyz_list_quad(:,3));

xyz_list_quad=[xyz_list_quad; x5 y5 z5];

%Compute four normals (the quad is divided into four triangles)
V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
normal1=1/norm(V1)*V1';

V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
normal2=1/norm(V2)*V2';

V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
normal3=1/norm(V3)*V3';

V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
normal4=1/norm(V4)*V4';

%Create a matrix with the four normals
normal=[normal1 normal2 normal3 normal4];

%Compute load contribution for this segment:
if pentaelem.shelf==1,
	%icefront ends in the water
	pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
else
	%icefront ends in the air
	pe.terms=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
end

end %End of function

function pe=QuadPressureLoadStokes(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);

%The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
%
%   grid 4 +-----------------+ grid 3
%          |\2            1 /|
%          |1\    tria3    /2|
%          |  \           /  |
%          |   \         /   |
%          |    \       /    |
%          |     \     /     |
%          |      \ 3 /      |
%          |tria4  \ / 3     |
%          |      3 \grid5   |
%          |       / \       | 
%          |      / 3 \ tria2|
%          |     /     \     |
%          |    /       \    |
%          |   /         \   |
%          |  /   tria1   \  |
%          |2/1           2\1|
%   grid1  +-----------------+ grid 2
%
%

%Build the four normal vectors 
nx=normal(1,:);
ny=normal(2,:);
nz=normal(3,:);

% Get gaussian points and weights. order 2 since we have a product of 2 nodal functions (We use GaussTria since we will build two triangles)
[num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);

%Recover the surface of the four nodes
s=thickness+bed;

%Add surface sor the fifth point (average of the surfaces)
s(5)=mean(s);

%Recover grid coordinates
x=xyz_list_quad(:,1);
y=xyz_list_quad(:,2);
z=xyz_list_quad(:,3);

%Build triangles in a 2D plan before using reference elements

%Create four triangle elements for each quad
tria=triaelem;

l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));

%First triangle : nodes 1, 2 and 5
x1tria1=0;
y1tria1=0;
x2tria1=sqrt(l12);
y2tria1=0;
x3tria1=cos_theta_triangle1*sqrt(l15);
y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);

xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];

%Second triangle : nodes 2, 3 and 5
x1tria2=0;
y1tria2=0;
x2tria2=sqrt(l23);
y2tria2=0;
x3tria2=cos_theta_triangle2*sqrt(l25);
y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);

xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];

%Third triangle : nodes 3, 4 and 5
x1tria3=0;
y1tria3=0;
x2tria3=sqrt(l34);
y2tria3=0;
x3tria3=cos_theta_triangle3*sqrt(l35);
y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);

xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];

%Fourth triangle : nodes 4, 1 and 5
x1tria4=0;
y1tria4=0;
x2tria4=sqrt(l14);
y2tria4=0;
x3tria4=cos_theta_triangle4*sqrt(l45);
y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);

xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];

%Create a list xyz_tria with coordinates of every triangle

xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];

%Initialize load vector
pe=zeros(16,1);

%Start  looping on the triangle's gaussian points: 
for ig=1:num_gauss,
	
	 %Get the value of nodal functions for each gauss point (nodal functions of the triangle to calculate z_g)
	 gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
	 l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
	 
	 %Get the coordinates of gauss point for each triangle in penta/quad
	 r_tria=gauss_coord(2)-gauss_coord(1);
	 s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
 
	 %Coordinates of gauss points in the reference triangle
	 r_quad(1,1)=r_tria;
	 s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
	 r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
	 s_quad(2,1)=r_tria;
	 r_quad(3,1)=-r_tria;
	 s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
	 r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
	 s_quad(4,1)=-r_tria;

	 %Get the nodal function of the quad for the gauss points of each triangle
	 l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];

	 %Compute jacobian of each triangle
	 for i=1:4
	complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero
	complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
	J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
	 end

	 %Calculation of the z coordinate for the gaussian point ig for each triangle
	 z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
	 z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
	 z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
	 z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);

	 %Loop on the triangles
	 for i=1:4

	%Loop on the grids of the quad
	air_pressure_tria=0;

	%Now deal with water pressure
	if fill==1, %icefront ends in water
		water_level_above_g_tria=min(0,z_g(i));              % 0 if the gaussian point is above water level
		water_pressure_tria=rho_water*gravity*water_level_above_g_tria;
	elseif fill==0,
		water_pressure_tria=0;
	else
		error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
	end

	%Add pressure from triangle i
	%Loop on the grids of the quad
	for j=1:4
		pressure_tria(j) = water_pressure_tria + air_pressure_tria;
	end

	pe=pe+J(i)*gauss_weights(ig)*...
	[ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
	  pressure_tria(1)*l1l4_tria(i,1)*ny(i)
	  pressure_tria(1)*l1l4_tria(i,1)*nz(i)
	  0
	  pressure_tria(2)*l1l4_tria(i,2)*nx(i)
	  pressure_tria(2)*l1l4_tria(i,2)*ny(i)
	  pressure_tria(2)*l1l4_tria(i,2)*nz(i)
	  0
	  pressure_tria(3)*l1l4_tria(i,3)*nx(i)
	  pressure_tria(3)*l1l4_tria(i,3)*ny(i)
	  pressure_tria(3)*l1l4_tria(i,3)*nz(i)
	  0
	  pressure_tria(4)*l1l4_tria(i,4)*nx(i)
	  pressure_tria(4)*l1l4_tria(i,4)*ny(i)
	  pressure_tria(4)*l1l4_tria(i,4)*nz(i)
	  0 ];

   end

end

end %End of function

function pe=CreatePVectorHoriz(load,elements,grids,materials,inputs);

	global element_debug  element_debugid


	if strcmpi(load.type,'segment'),
		%check that the element is onbed (collapsed formulation) otherwise:pe=0
		if strcmpi(elements(load.eid).element.type,'pentaelem'),
			if  ~elements(load.eid).element.onbed
				pe=elemvector(0);
				return;
			end
		end

		%some variables
		numgrids=2;
		NDOF2=2;

		%Create elementary vector.
		pe=elemvector(numgrids*NDOF2);

		%recover element object
		triaelem=elements(load.eid).element;
		matpar=materials(end).constants;

		if (element_debug & triaelem.id==element_debugid),
			disp(sprintf('%s%i','Inside icefront load  connected to element ',triaelem.id));
		end

		%recover material parameters
		gravity=matpar.g;
		rho_ice=matpar.rho_ice;
		rho_water=matpar.rho_water;
		
		%recover extra inputs
		[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
		[bed_param bed_is_present]=recover_input(inputs,'bed');
			
		%Get all element grid data:
		xyz_list=getgriddata(triaelem,grids);

		%initialize extra inputs
		thickness_list=zeros(3,1);
		bed_list=zeros(3,1);

		%Identify which grids are comprised in the segment: 
		for i=1:3,
			if triaelem.g(i)==load.g(1),
				grid1=i;
			end
			if triaelem.g(i)==load.g(2),
				grid2=i;
			end
		end

		%Build row indices for elementary vector.
		for i=1:numgrids,
			doflist=grids(load.g(i)).grid.doflist;
			for j=1:NDOF2,
				pe.row_indices((i-1)*NDOF2+j)=doflist(j);
			end
		end
		
		for i=1:3,
			doflist=grids(triaelem.g(i)).grid.doflist; dof=doflist(1);
			if(thickness_is_present),
				thickness_list(i)=thickness_param(dof);
			else
				thickness_list(i)=triaelem.h(i);
			end
			if(bed_is_present),
				bed_list(i)=bed_param(dof);
			else
				bed_list(i)=triaelem.b(i);
			end;
		end

		%Recover thicknesses and bed at grid1 and grid2: 
		if thickness_is_present,
			thickness=[thickness_list(grid1) thickness_list(grid2)];
			bed=[bed_list(grid1) bed_list(grid2)];
		else
			thickness=[triaelem.h(grid1) triaelem.h(grid2)];
			bed=[triaelem.b(grid1) triaelem.b(grid2)];
		end

		%Recover grid coordinates
		x1=xyz_list(grid1,1);
		y1=xyz_list(grid1,2);
		x2=xyz_list(grid2,1);
		y2=xyz_list(grid2,2);

		%Compute length and normal of segment
		normal=zeros(2,1);
		normal(1)=cos(atan2(x1-x2,y2-y1));
		normal(2)=sin(atan2(x1-x2,y2-y1));
		length=sqrt((x2-x1)^2+(y2-y1)^2);

		%Compute load contribution for this segment:
		if triaelem.shelf==1,
			%icefront ends in the water
			fill=1;
		else
			%icefront ends in the air
			fill=0;
		end
		pe.terms=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);
		
		if (element_debug & triaelem.id==element_debugid),

			disp(sprintf('\nicefront load'));
			disp(sprintf('grids %i %i',grid1,grid2));
			disp(sprintf('rho_water %g',rho_water));
			disp(sprintf('rho_ice %g',rho_ice));
			disp(sprintf('gravity %g',gravity));
			disp(sprintf('thickness (%g,%g)',thickness(1),thickness(2)));
			disp(sprintf('bed (%g,%g)',bed(1),bed(2)));
			disp(sprintf('normal (%g,%g)',normal(1),normal(2)));
			disp(sprintf('length %g',length));
			disp(sprintf('fill %i',fill));
			
			disp(sprintf('      pe_g->terms\n'));
			pe.terms
		end

	elseif strcmpi(load.type,'quad'),

		%some variables
		numgrids=4;
		NDOF2=2;

		numgridstria=4; %number of grids for the triangles
		NDOF2=2;


		%Create elementary vector.
		pe=elemvector(numgrids*NDOF2);

		%recover element object
		pentaelem=elements(load.eid).element;
		matpar=materials(end).constants;

		%recover material parameters
		gravity=matpar.g;
		rho_ice=matpar.rho_ice;
		rho_water=matpar.rho_water;
		
		%recover extra inputs
		[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
		[bed_param bed_is_present]=recover_input(inputs,'bed');

		%initialize extra inputs
		thickness_list=zeros(6,1);
		bed_list=zeros(6,1);

		%Get all element grid data:
		xyz_list=getgriddata(pentaelem,grids);

		%Identify which grids are comprised in the segment: 
		for i=1:6,
			if pentaelem.g(i)==load.g(1),
				grid1=i;
			end
			if pentaelem.g(i)==load.g(2),
				grid2=i;
			end
			if pentaelem.g(i)==load.g(3),
				grid3=i;
			end
			if pentaelem.g(i)==load.g(4),
				grid4=i;
			end
		end
		quad_grids=[grid1 grid2 grid3 grid4];

		%Build row indices for elementary vector.
		for i=1:numgrids,
			doflist=grids(load.g(i)).grid.doflist;
			for j=1:NDOF2,
				pe.row_indices((i-1)*NDOF2+j)=doflist(j);
			end
			dof=doflist(1);
		end

		for i=1:6,
			doflist=grids(pentaelem.g(i)).grid.doflist; dof=doflist(1);
			if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
			if(bed_is_present) bed_list(i)=bed_param(dof);end;
		end


		%Recover thicknesses and bed at all grids
		if thickness_is_present,
			thickness=thickness_list(quad_grids);
			bed=bed_list(quad_grids);
		else
			thickness=pentaelem.h(quad_grids);
			bed=pentaelem.b(quad_grids);
		end

		%build xyz list for [grid1, grid2,grid3,grid3]
		xyz_list_quad=xyz_list(quad_grids,:);

		%Create a new grid in the midle of the quad and add it at the end of the list
		x5=mean(xyz_list_quad(:,1));
		y5=mean(xyz_list_quad(:,2));
		z5=mean(xyz_list_quad(:,3));

		xyz_list_quad=[xyz_list_quad; x5 y5 z5];

		%Compute four normals (the quad is divided into four triangles)
		V1=cross(xyz_list_quad(1,:)-xyz_list_quad(5,:),xyz_list_quad(2,:)-xyz_list_quad(5,:));
		normal1=1/norm(V1)*V1';
		
		V2=cross(xyz_list_quad(2,:)-xyz_list_quad(5,:),xyz_list_quad(3,:)-xyz_list_quad(5,:));
		normal2=1/norm(V2)*V2';

		V3=cross(xyz_list_quad(3,:)-xyz_list_quad(5,:),xyz_list_quad(4,:)-xyz_list_quad(5,:));
		normal3=1/norm(V3)*V3';
		
		V4=cross(xyz_list_quad(4,:)-xyz_list_quad(5,:),xyz_list_quad(1,:)-xyz_list_quad(5,:));
		normal4=1/norm(V4)*V4';

		%Create a matrix with the four normals
		normal=[normal1 normal2 normal3 normal4];

		%Compute load contribution for this segment:
		if pentaelem.shelf==1,
			%icefront ends in the water
			pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,1);
		else
			%icefront ends in the air
			pe.terms=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,0);
		end

		if (element_debug & pentaelem.id==element_debugid),
			disp(sprintf('      pe_g->terms\n'));
			pe.terms
		end

	else
		error('unsupported element type')
	end
	
	end %End of function

	function pe=QuadPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,xyz_list_quad,fill);

	%The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows
	%
	%   grid 4 +-----------------+ grid 3
	%          |\2            1 /|
	%          |1\    tria3    /2|
	%          |  \           /  |
	%          |   \         /   |
	%          |    \       /    |
	%          |     \     /     |
	%          |      \ 3 /      |
	%          |tria4  \ / 3     |
	%          |      3 \grid5   |
	%          |       / \       | 
	%          |      / 3 \ tria2|
	%          |     /     \     |
	%          |    /       \    |
	%          |   /         \   |
	%          |  /   tria1   \  |
	%          |2/1           2\1|
	%   grid1  +-----------------+ grid 2
	%
	%

	%Build the four normal vectors 
	nx=normal(1,:);
	ny=normal(2,:);

	% Get gaussian points and weights. order 2 since we have a product of 2 nodal functions (We use GaussTria since we will build two triangles)
	[num_gauss,first_gauss_area,second_gauss_area,third_gauss_area,gauss_weights]=GaussTria(2);

	%Recover the surface of the four nodes
	s=thickness+bed;

	%Add surface sor the fifth point (average of the surfaces)
	s(5)=mean(s);

	%Recover grid coordinates
	x=xyz_list_quad(:,1);
	y=xyz_list_quad(:,2);
	z=xyz_list_quad(:,3);

	%Build triangles in a 2D plan before using reference elements

	%Create four triangle elements for each quad
	tria=triaelem;

	l12=(x(2)-x(1))^2+(y(2)-y(1))^2+(z(2)-z(1))^2;
	l14=(x(4)-x(1))^2+(y(4)-y(1))^2+(z(4)-z(1))^2;
	l15=(x(5)-x(1))^2+(y(5)-y(1))^2+(z(5)-z(1))^2;
	l23=(x(3)-x(2))^2+(y(3)-y(2))^2+(z(3)-z(2))^2;
	l25=(x(5)-x(2))^2+(y(5)-y(2))^2+(z(5)-z(2))^2;
	l34=(x(4)-x(3))^2+(y(4)-y(3))^2+(z(4)-z(3))^2;
	l35=(x(5)-x(3))^2+(y(5)-y(3))^2+(z(5)-z(3))^2;
	l45=(x(5)-x(4))^2+(y(5)-y(4))^2+(z(5)-z(4))^2;
	cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15));
	cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25));
	cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35));
	cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45));

	%First triangle : nodes 1, 2 and 5
	x1tria1=0;
	y1tria1=0;
	x2tria1=sqrt(l12);
	y2tria1=0;
	x3tria1=cos_theta_triangle1*sqrt(l15);
	y3tria1=sqrt(1-cos_theta_triangle1^2)*sqrt(l15);

	xyz_tria1=[x1tria1 y1tria1; x2tria1 y2tria1; x3tria1 y3tria1];

	%Second triangle : nodes 2, 3 and 5
	x1tria2=0;
	y1tria2=0;
	x2tria2=sqrt(l23);
	y2tria2=0;
	x3tria2=cos_theta_triangle2*sqrt(l25);
	y3tria2=sqrt(1-cos_theta_triangle2^2)*sqrt(l25);

	xyz_tria2=[x1tria2 y1tria2; x2tria2 y2tria2; x3tria2 y3tria2];

	%Third triangle : nodes 3, 4 and 5
	x1tria3=0;
	y1tria3=0;
	x2tria3=sqrt(l34);
	y2tria3=0;
	x3tria3=cos_theta_triangle3*sqrt(l35);
	y3tria3=sqrt(1-cos_theta_triangle3^2)*sqrt(l35);

	xyz_tria3=[x1tria3 y1tria3; x2tria3 y2tria3; x3tria3 y3tria3];

	%Fourth triangle : nodes 4, 1 and 5
	x1tria4=0;
	y1tria4=0;
	x2tria4=sqrt(l14);
	y2tria4=0;
	x3tria4=cos_theta_triangle4*sqrt(l45);
	y3tria4=sqrt(1-cos_theta_triangle4^2)*sqrt(l45);

	xyz_tria4=[x1tria4 y1tria4; x2tria4 y2tria4; x3tria4 y3tria4];

	%Create a list xyz_tria with coordinates of every triangle

	xyz_tria=[xyz_tria1; xyz_tria2; xyz_tria3; xyz_tria4];

	%Initialize load vector
	pe=zeros(8,1);

	%Start  looping on the triangle's gaussian points: 
	for ig=1:num_gauss,
		
		 %Get the value of nodal functions for each gauss point (nodal functions of the triangle to calculate z_g)
		 gauss_coord=[first_gauss_area(ig) second_gauss_area(ig) third_gauss_area(ig)];
		 l1l2l3_tria=GetNodalFunctions(tria,gauss_coord);
		 
		 %Get the coordinates of gauss point for each triangle in penta/quad
		 r_tria=gauss_coord(2)-gauss_coord(1);
		 s_tria=-3/sqrt(3)*(gauss_coord(1)+gauss_coord(2)-2/3);
	 
		 %Coordinates of gauss points in the reference triangle
		 r_quad(1,1)=r_tria;
		 s_quad(1,1)=1/sqrt(3)*s_tria-2/3;
		 r_quad(2,1)=-1/sqrt(3)*s_tria+2/3;
		 s_quad(2,1)=r_tria;
		 r_quad(3,1)=-r_tria;
		 s_quad(3,1)=-1/sqrt(3)*s_tria+2/3;
		 r_quad(4,1)=1/sqrt(3)*s_tria-2/3;
		 s_quad(4,1)=-r_tria;

		 %Get the nodal function of the quad for the gauss points of each triangle
		 l1l4_tria=1/4*[(r_quad-1).*(s_quad-1) -(r_quad+1).*(s_quad-1) (r_quad+1).*(s_quad+1) -(r_quad-1).*(s_quad+1)];

		 %Compute jacobian of each triangle
		 for i=1:4
		complete_list=zeros(3,3); %a third coordinate is needed for the jacobian determinant calculation, here it is zero
		complete_list(1:3,1:2)=xyz_tria(3*i-2:3*i,:);
		J(i)=GetJacobianDeterminant(tria,complete_list,l1l2l3_tria);
		 end

		 %Calculation of the z coordinate for the gaussian point ig for each triangle
		 z_g(1)=z(1)*l1l2l3_tria(1)+z(2)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
		 z_g(2)=z(2)*l1l2l3_tria(1)+z(3)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
		 z_g(3)=z(3)*l1l2l3_tria(1)+z(4)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);
		 z_g(4)=z(4)*l1l2l3_tria(1)+z(1)*l1l2l3_tria(2)+z(5)*l1l2l3_tria(3);

		 %Loop on the triangles
		 for i=1:4

		%Loop on the grids of the quad
		%Calculate the ice pressure
		for j=1:4
			ice_pressure_tria(j)=gravity*rho_ice*(s(j)-z_g(i));
		end
		air_pressure_tria=0;

		%Now deal with water pressure
		if fill==1, %icefront ends in water
		water_level_above_g_tria=min(0,z_g(i));              % 0 if the gaussian point is above water level
		water_pressure_tria=rho_water*gravity*water_level_above_g_tria;

		elseif fill==0,
		water_pressure_tria=0;

		else
		error('QuadPressureLoad error message: unknow fill type for icefront boundary condition');
		end

		%Add pressure from triangle i
		%Loop on the grids of the quad
		for j=1:4
			pressure_tria(j) = ice_pressure_tria(j) + water_pressure_tria + air_pressure_tria;
		end

		pe=pe+J(i)*gauss_weights(ig)*...
		[ pressure_tria(1)*l1l4_tria(i,1)*nx(i)
		  pressure_tria(1)*l1l4_tria(i,1)*ny(i)
		  pressure_tria(2)*l1l4_tria(i,2)*nx(i)
		  pressure_tria(2)*l1l4_tria(i,2)*ny(i)
		  pressure_tria(3)*l1l4_tria(i,3)*nx(i)
		  pressure_tria(3)*l1l4_tria(i,3)*ny(i)
		  pressure_tria(4)*l1l4_tria(i,4)*nx(i)
		  pressure_tria(4)*l1l4_tria(i,4)*ny(i)];

	   end

	end

	end



	function pe=SegmentPressureLoad(rho_water,rho_ice,gravity,thickness,bed,normal,length,fill);

	%Build the normal vector
	nx=normal(1,1);
	ny=normal(2,1);

	% Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
	num_gauss=3;
	[segment_gauss_coord, gauss_weights]=GaussSegment(num_gauss);

	%Recover the thickness and z_b of the two nodes
	h1=thickness(1);
	b1=bed(1);
	h2=thickness(2);
	b2=bed(2);

	%Compute jacobian of segment
	Jdet=1/2*length;

	%Initialize load vector
	pe=zeros(4,1);

	%Start  looping on segments's gaussian points: 
	for ig=1:num_gauss,

		thickness=h1*(1+segment_gauss_coord(ig))/2+h2*(1-segment_gauss_coord(ig))/2;
		bed=b1*(1+segment_gauss_coord(ig))/2+b2*(1-segment_gauss_coord(ig))/2;
		
		if fill==1, %icefront ends in water
			ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
			air_pressure=0;

			%Now deal with water pressure
			surface_under_water=min(0,thickness+bed); % 0 if the top of the glacier is above water level
			base_under_water=min(0,bed);              % 0 if the bottom of the glacier is above water level
			water_pressure=1.0/2.0*gravity*rho_water*(surface_under_water^2 - base_under_water^2);

		elseif fill==0,
			ice_pressure=1.0/2.0*gravity*rho_ice*thickness^2;
			air_pressure=0;
			water_pressure=0;

		else
			error('SegmentPressureLoad error message: unknow fill type for icefront boundary condition');
		end

		pressure = ice_pressure + water_pressure + air_pressure;
		

		pe=pe+pressure*Jdet*gauss_weights(ig)*...
		[ (1+segment_gauss_coord(ig))/2*nx
		  (1+segment_gauss_coord(ig))/2*ny
		  (1-segment_gauss_coord(ig))/2*nx
		  (1-segment_gauss_coord(ig))/2*ny];

	end

end %End of function
