function Ke=CreateKMatrix(pentaelem,grids,materials,inputs,analysis_type)
%CREATEKMATRIX - create the stiffmess matrix for pentaelem
%
%   this stiffness matrix works for MacAyeal, Pattyn's  and Stokes' model ,
%   thermal model (transient ans steady), prognostic, melting and 
%   for the computation of the bed slope and surface slope
%
%   Usage:
%      Ke=CreateKMatrix(pentaelem,grids,materials,inputs,analysis_type)
%
%   See also CREATEPVECTOR

if strcmpi(analysis_type,'diagnostic_horiz'),

	Ke=CreateKMatrixHoriz(pentaelem,grids,materials,inputs,analysis_type);

elseif strcmpi(analysis_type,'diagnostic_vert'),

	Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'diagnostic_basevert'),

	Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'prognostic'),

	Ke=CreateKMatrixPrognostic(pentaelem,grids,materials,inputs);

elseif (strcmpi(analysis_type,'thermaltransient') | strcmpi(analysis_type,'thermalsteady')),

	Ke=CreateKMatrixThermal(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'melting')

	Ke=CreateKMatrixMelting(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'diagnostic_stokes'),

	Ke=CreateKMatrixStokes(pentaelem,grids,materials,inputs);

elseif strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y') | strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),

	Ke=CreateKMatrixSlopeCompute(pentaelem,grids,materials,inputs,analysis_type);


end
end %end function

function Ke=CreateKMatrixSlopeCompute(pentaelem,grids,materials,inputs,analysis_type)

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.

	if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
		if pentaelem.onbed~=1,
			Ke=elemmatrix(0);
			return;
		end
	elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
		if pentaelem.onsurface~=1,
			Ke=elemmatrix(0);
			return;
		end
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

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

	%Just keep the first 3 grids
	if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
		xyz_list=xyz_list(1:3,:);
	elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
		xyz_list=xyz_list(4:6,:);
	end	

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		if strcmpi(analysis_type,'bed_slope_compute_x') | strcmpi(analysis_type,'bed_slope_compute_y'),
			doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		elseif strcmpi(analysis_type,'surface_slope_compute_x') | strcmpi(analysis_type,'surface_slope_compute_y'),
			doflist=grids(pentaelem.g(i+3)).grid.doflist; %list of dofs in the g-set: last three grids
		end
		dof=doflist(1);
		Ke.row_indices(i)=dof;
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L matrix if viscous basal drag present:
		L=GetL(triaelem,gauss_coord,DOFPERGRID);
   
		DL_scalar=gauss_weight*Jdettria;

		%Do the triple product tL*D*L. 
		Ke_gg_gaussian=L'*DL_scalar*L;
		
		%Add Ke_gg_drag_gaussian to Ke
		Ke.terms=Ke.terms+Ke_gg_gaussian;
	end

end %end function

function Ke=CreateKMatrixStokes(pentaelem,grids,materials,inputs,analysis_type)

%some variables
numgrids=6;
DOFPERGRID=4;
numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

numgrids2d=3;
numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.

%some parameters
MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT

%Create elementary stiffness matrix 
Ke=elemmatrix(numdof);

%recover material
matice=materials(pentaelem.matid).material;
matpar=materials(end).constants;

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

%recover conditioning number for the matrix
conditioning=pentaelem.reconditioning_number;

%recover extra inputs from users, at current convergence iteration.
[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
[bed_param bed_is_present]=recover_input(inputs,'bed');
[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
[velocity_param velocity_is_present]=recover_input(inputs,'velocity');

%initialize extra arrays if some inputs are present
if velocity_is_present,
	vxvyvz_list=zeros(numgrids,3);
end
B_list=zeros(numgrids,1);
bed_list=zeros(numgrids,1);
thickness_list=zeros(numgrids,1);
K_list=zeros(numgrids,1);
temperature_list=zeros(numgrids,1);

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

%Build linear indices for elementary stiffness matrix.
for i=1:numgrids,
	doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
	for j=1:DOFPERGRID,
		dof=doflist(j);
		Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
		if j<4, %the fourth dof corresponds to pressure and not velocity
			if velocity_is_present, vxvyvz_list(i,j)=velocity_param(dof); end;
		end
	end
	
	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;
	if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
	if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
	if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;
end

% Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
%get tria gaussian points as well as segment gaussian points. For tria gaussian 
%points, order of integration is 2, because we need to integrate the product tB*D*B' 
%which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
%points, same deal, which yields 3 gaussian points.
area_order=5;
num_vert_gauss=5;
[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

Ke_temp=zeros(27,27);

%Start  looping on the number of gaussian points:
for igarea=1:num_area_gauss,
	for igvert=1:num_vert_gauss,
		%Pick up the gaussian point and its weight:
		gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
		gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];

		%Get strain rate, if velocity has been supplied: 
		if velocity_is_present,
			epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord);
		else
			epsilon=zeros(6,1);
		end
		%disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
		
		%Update material parameter that deals with ice rigidity. The update depends either on the 
		%temperature (if we are running thermal, or transient), or on B (if we are running inverse 
		%control methods on B). The latter prevails. 
		
		%Update material if temperature is provided.
		if temperature_is_present,
			temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord);
			matice.B=paterson(temperature);
		end

		%Update material if flow law is specified. This will erase the previous change 
		%on B when temperature is provided. 
		if flow_law_is_present,
			B_param=GetParameterValueStokes(pentaelem,B_list,gauss_coord);
			matice.B=B_param; clear B_param.
		end

		%Get viscosity at last two iterations: 
		viscosity=GetViscosity3d(matice,epsilon);
		%disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));

		%Get B and Bprime matrices:
		B=GetBStokes(pentaelem,xyz_list,gauss_coord);
		Bprime=GetBprimeStokes(pentaelem,xyz_list,gauss_coord);

		%Get Jacobian determinant:
		Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
		%disp(sprintf('Jacobian determinant: %f\n', Jdet));

		% Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant 
		% onto this scalar matrix, so that we win some computational time: */
		D=gauss_weight*Jdet*diag([viscosity*ones(6,1);-conditioning*ones(2,1)]);

		% Do the triple product tB*D*Bprime: 
		Ke_gg_gaussian=B'*D*Bprime;

		%Add Ke_gg_gaussian to Ke_temp: 
		Ke_temp=Ke_temp+Ke_gg_gaussian;

	end %for igarea=1:num_area_gauss,
end %for igvert=1:num_vert_gauss,

%Deal with 2d friction at the bedrock interface
if (pentaelem.onbed==1 & (pentaelem.shelf==0))

	%build alpha^2 of each of the 3 grids of the triangle
	if (pentaelem.friction_type==2),  %friction_type=2 implies viscous, friction_type=1 implies plastic
	
		%Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
		frictionparameters=struct();
		frictionparameters.element_type='3d';
		frictionparameters.rho_ice=rho_ice;
		frictionparameters.rho_water=rho_water;
		frictionparameters.g=gravity;	
		frictionparameters.p=pentaelem.p;
		frictionparameters.q=pentaelem.q;

		if velocity_is_present,
			frictionparameters.velocities=vxvyvz_list(1:3,:);
		else
			frictionparameters.velocities=zeros(3,3);
		end
		if thickness_is_present,
			frictionparameters.h=thickness_list(1:3);
		else
			frictionparameters.h=pentaelem.h(1:3);
		end
		if bed_is_present,
			frictionparameters.b=bed_list(1:3);
		else
			frictionparameters.b=pentaelem.b(1:3);
		end
		if basal_drag_is_present,
			frictionparameters.k=K_list(1:3);
		else
			frictionparameters.k=pentaelem.k(1:3);
		end

		alpha2=Getalpha2(frictionparameters);
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord_tria=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		gauss_coord_penta=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig) -1];

		%Get the Jacobian determinant
		xyz_list_tria=xyz_list(1:3,:);
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord_tria);

		%Get L and Lprime matrix if viscous basal drag present:
		L=GetLStokes(triaelem,gauss_coord_tria);
		Lprime=GetLprimeStokes(triaelem,xyz_list,gauss_coord_tria,gauss_coord_penta);

		%Build normal vector to the surface pointing outward the glacier system
		n=zeros(3,1);
		normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(3,:)-xyz_list_tria(2,:)));
		n(1)=1/norm(normal)*normal(1);
		n(2)=1/norm(normal)*normal(2);
		n(3)=1/norm(normal)*normal(3);

		%Get strain rate, if velocity has been supplied: 
		if velocity_is_present,
			epsilon=GetStrainRateStokes(pentaelem,vxvyvz_list,xyz_list,gauss_coord_penta);
		else
			epsilon=zeros(6,1);
		end


		%Update material if temperature is provided.
		if temperature_is_present,
			temperature=GetParameterValueStokes(pentaelem,temperature_list,gauss_coord_penta);
			matice.B=paterson(temperature);
		end

		%Get viscosity at last two iterations: 
		viscosity=GetViscosity3d(matice,epsilon);
		%disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));


		%Get alpha2 on current gaussian point
		alpha2_g=GetParameterValue(triaelem,alpha2,gauss_coord_tria);

		%Create DL diagonal matrix
		DL=diag([alpha2_g*[1;1;-n(1)*n(3);-n(2)*n(3);-n(1)*n(3);-n(2)*n(3)]; -viscosity*[n(1);n(2);n(3);n(1)/2.0;n(2)/2.0];conditioning*[n(1);n(2);n(3)]])*gauss_weight*Jdettria;

		%Create stiffness matrix for drag
		Ke_gg_drag_gaussian=L'*DL*Lprime;

		%Add Ke_gg_drag_gaussian Ke 
		Ke_temp(1:12,1:12)=Ke_temp(1:12,1:12)+Ke_gg_drag_gaussian;
	end
end%ends friction

%Reduce the matrix
Ke_reduced=ReduceMatrixStokes(pentaelem,Ke_temp);

%Add the components of the matrix for the 6 grids 
Ke.terms=Ke.terms+Ke_reduced;

end%end function


function Ke=CreateKMatrixVert(pentaelem,grids,materials,inputs)

	%some variables
	numgrids=6;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	numgrids2d=3;
	numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

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

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		dof=doflist(3); %z degree of freedom
		Ke.row_indices(i)=dof;
	end

	% Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
	%get tria gaussian points as well as segment gaussian points. For tria gaussian 
	%points, order of integration is 2, because we need to integrate the product tB*D*B' 
	%which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
	%points, same deal, which yields 3 gaussian points.
	area_order=2;
	num_vert_gauss=2;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
		for igvert=1:num_vert_gauss,
			%Pick up the gaussian point and its weight:
			gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
			gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
			
			%Get B and Bprime matrices:
			B=GetB_vert(pentaelem,xyz_list,gauss_coord);
			Bprime=GetBprime_vert(pentaelem,gauss_coord);

			%Get Jacobian determinant:
			Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
			%disp(sprintf('Jacobian determinant: %f\n', Jdet));

			D=gauss_weight*Jdet;

			% Do the triple product tB*D*Bprime: 
			Ke_gg_gaussian=B'*D*Bprime;

			%Add Ke_gg_gaussian to Ke: 
			Ke.terms=Ke.terms+Ke_gg_gaussian;
			
		end %for ig=1:num_area_gauss,
	  end %for ig=1:num_vert_gauss,

	%Deal with upper surface integration: 

	if (pentaelem.onsurface==1),
		
		%Get the coordinates of the three grids
		xyz_list_tria=xyz_list(4:6,:);

		% Get gaussian points and weights.
		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

		%Build normal vector to the surface
		normal=cross((xyz_list_tria(1,:)-xyz_list_tria(3,:)),(xyz_list_tria(2,:)-xyz_list_tria(3,:)));
		nz=1/norm(normal)*normal(3);

		for ig=1:num_gauss2D,

			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			
			%Get the Jacobian determinant
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);

			%Get L matrix if viscous basal drag present:
			L=GetL(triaelem,gauss_coord,DOFPERGRID);
	   
			DL_scalar=gauss_weight*Jdettria*nz;

			%Do the triple product tL*D*L. Do not forget the '-' sign!
			Ke_gg_vert_gaussian=-L'*DL_scalar*L;
			

			%Add Ke_gg_drag_gaussian to Ke
			Ke.terms(4:6,4:6)=Ke.terms(4:6,4:6)+Ke_gg_vert_gaussian;
		end
	end%ends friction

end %end function

function Ke=CreateKMatrixBaseVert(pentaelem,grids,materials,inputs)

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	if pentaelem.onbed~=1,
		Ke=elemmatrix(0);
		return;
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

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

	%Just keep the first 3 grids
	xyz_list=xyz_list(1:3,:);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		dof=doflist(3); %z degree of freedom
		Ke.row_indices(i)=dof;
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L matrix if viscous basal drag present:
		L=GetL(triaelem,gauss_coord,DOFPERGRID);
   
		DL_scalar=gauss_weight*Jdettria;

		%Do the triple product tL*D*L. 
		Ke_gg_vert_gaussian=L'*DL_scalar*L;
		
		%Add Ke_gg_drag_gaussian to Ke
		Ke.terms=Ke.terms+Ke_gg_vert_gaussian;
	end

end %end function




function Ke=CreateKMatrixThermal(pentaelem,grids,materials,inputs)


	%some variables
	numgrids=6;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%recover material
	matice=materials(pentaelem.matid).material;
	matpar=materials(end).constants;

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

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

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

	%recover extra inputs from users, at current convergence iteration.
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	[dt dt_is_present]=recover_input(inputs,'dt');
	if (~velocity_is_present )
		error('CreateKMatrixThermal error message: missing velocity parameter!');
	end
	if (~dt_is_present & ~pentaelem.thermal_steadystate)
		error('CreateKMatrixThermal error message: transient thermal solution requires present of time step as parameter!');
	end

	%initialize vxvyvz_list 
	vxvyvz_list=zeros(numgrids,3);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,

		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		for j=1:3,
			dof=doflist(j);
			vxvyvz_list(i,j)=velocity_param(dof);
		end	
		dof=doflist(1);
		Ke.row_indices(i)=dof;
		
	end

	% Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
	%get tria gaussian points as well as segment gaussian points. For tria gaussian 
	%points, order of integration is 2, because we need to integrate the product tB*D*B' 
	%which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
	%points, same deal, which yields 3 gaussian points.
	area_order=2;
	num_vert_gauss=2;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
		for igvert=1:num_vert_gauss,
			%Pick up the gaussian point and its weight:
			gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
			gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];

			%Get Jacobian determinant:
			Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
			%disp(sprintf('Jacobian determinant: %f\n', Jdet));
			    
			%conduction

			%Get B and Bprime matrices:
			Bconduct=GetB_conduct(pentaelem,xyz_list,gauss_coord);

			% Build the D matrix
			D=gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity))*eye(3);

			% Do the triple product tB*D*B: 
			if pentaelem.thermal_steadystate,
				Kconduct=Bconduct'*D*Bconduct;
			else
				Kconduct=dt*Bconduct'*D*Bconduct;
			end

			%advection 

			%Get B and Bprime matrices:
			Badvec=GetB_advec(pentaelem,gauss_coord);
			Bprimeadvec=GetBprime_advec(pentaelem,xyz_list,gauss_coord);

			% Build the D matrix
			vx=GetParameterValue(pentaelem,vxvyvz_list(:,1),gauss_coord);
			vy=GetParameterValue(pentaelem,vxvyvz_list(:,2),gauss_coord);
			vz=GetParameterValue(pentaelem,vxvyvz_list(:,3),gauss_coord);

			D=gauss_weight*Jdet*[vx,0,0; 0,vy,0; 0,0,vz];

			% Do the triple product tB*D*Bprime: 
			if pentaelem.thermal_steadystate,
				Kadvec=Badvec'*D*Bprimeadvec;
			else
				Kadvec=dt*Badvec'*D*Bprimeadvec;
			end

			%transient 
			if pentaelem.thermal_steadystate,
				Ktransient=0;
			else
				L=GetNodalFunctions(pentaelem,gauss_coord);
				Ktransient=L'*L*gauss_weight*Jdet;
			end

			%Add all contributions to Ke: 
			Ke.terms=Ke.terms+Kadvec+Kconduct+Ktransient;
			%Ke.terms=Ke.terms+Kconduct+Ktransient;

		end %for ig=1:vert_gauss,
	end %for ig=1:num_area_gauss,
			
	%Ice/ocean heat exchange flux on ice shelf base 
	if (pentaelem.onbed & pentaelem.shelf),

		%recover material parameters
		mixed_layer_capacity=matpar.mixed_layer_capacity;
		thermal_exchange_velocity=matpar.thermal_exchange_velocity;

		% Get gaussian points and weights.
		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

		for ig=1:num_gauss2D,

			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			
			%Get the Jacobian determinant
			xyz_list_tria=xyz_list(1:3,:);
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);
			
			%Get nodal functions value
			l1l2l3=GetNodalFunctions(triaelem,gauss_coord);

			%Compute the elementary load vector
			if pentaelem.thermal_steadystate,
				Ke_ocean_gaussian=l1l2l3'*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice)*l1l2l3;
			else
				Ke_ocean_gaussian=l1l2l3'*dt*gauss_weight*Jdettria*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice)*l1l2l3;
			end

			%Add Ke_gg_drag_gaussian to Ke
			Ke.terms(1:3,1:3)=Ke.terms(1:3,1:3)+Ke_ocean_gaussian;
		end
	end%ends ice/ocean exchange flux 


end %end function







function Ke=CreateKMatrixHoriz(pentaelem,grids,materials,inputs,analysis_type);

global element_debug  element_debugid

%Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the 
%bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 
%the stiffness matrix. 
if pentaelem.collapse & ~pentaelem.onbed

	Ke=elemmatrix(0);
	return;

elseif pentaelem.collapse & pentaelem.onbed
        

	Ke=PentaCollapseIntoTriaMatrix(pentaelem, grids, materials, inputs, analysis_type);

	%That's it, element will get destroyed upon return, and we have the correct 
	%collapsed stiffness matrix. 
	return;

else 

	%implement standard penta element.

	%some variables
	numgrids=6;
	DOFPERGRID=2;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	numgrids2d=3;
	numdof2d=numgrids2d*DOFPERGRID; %number of dof for friction element triaelem.


	%some parameters
	MAXSLOPE=.06;  %any element with slope>MAXSLOPE is considered a "rock" element, with infinite stiffness.
	MOUNTAINKEXPONENT=10; % "infinite" stiffness is going to be  10^MOUNTAINKEXPONENT


	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

	%recover material
	matice=materials(pentaelem.matid).material;
	matpar=materials(end).constants;

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

	%recover extra inputs from users, at current convergence iteration.
	[basal_drag_param basal_drag_is_present]=recover_input(inputs,'drag');
	[thickness_param thickness_is_present]=recover_input(inputs,'thickness');
	[surface_param surface_is_present]=recover_input(inputs,'surface');
	[bed_param bed_is_present]=recover_input(inputs,'bed');
	[temperature_param temperature_is_present]=recover_input(inputs,'temperature');
	[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
	[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
	[oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');

	%initialize extra arrays if some inputs are present
	if velocity_is_present,
		vxvy_list=zeros(numgrids,2);
	end
	if oldvelocity_is_present,
		oldvxvy_list=zeros(numgrids,2);
	end

	B_list=zeros(numgrids,1);
	K_list=zeros(numgrids,1);
	thickness_list=zeros(numgrids,1);
	surface_list=zeros(numgrids,1);
	bed_list=zeros(numgrids,1);
	temperature_list=zeros(numgrids,1);

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

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		for j=1:DOFPERGRID,
			dof=doflist(j);
			Ke.row_indices((i-1)*DOFPERGRID+j)=dof;
			if velocity_is_present, vxvy_list(i,j)=velocity_param(dof); end;
			if oldvelocity_is_present, oldvxvy_list(i,j)=oldvelocity_param(dof); end;
		end
		
		dof=doflist(1);
		if(flow_law_is_present), B_list(i) = flow_law_param(dof);end;
		if(basal_drag_is_present), K_list(i)=basal_drag_param(dof);end;
		if(thickness_is_present) thickness_list(i)=thickness_param(dof);end;
		if(surface_is_present) surface_list(i)=surface_param(dof);end;
		if(bed_is_present) bed_list(i)=bed_param(dof);end;
		if(temperature_is_present) temperature_list(i)=temperature_param(dof);end;

	end

		% Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 
	%get tria gaussian points as well as segment gaussian points. For tria gaussian 
	%points, order of integration is 2, because we need to integrate the product tB*D*B' 
	%which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian 
	%points, same deal, which yields 3 gaussian points.
	area_order=5;
	num_vert_gauss=5;
	[num_area_gauss,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,area_gauss_weights, vert_gauss_coord,vert_gauss_weights]=GaussPenta(area_order,num_vert_gauss);

	%Start  looping on the number of gaussian points:
	for igarea=1:num_area_gauss,
		for igvert=1:num_vert_gauss,
		%Pick up the gaussian point and its weight:
		gauss_weight=area_gauss_weights(igarea)*vert_gauss_weights(igvert);
		gauss_coord=[first_area_gauss_coord(igarea) second_area_gauss_coord(igarea) third_area_gauss_coord(igarea) vert_gauss_coord(igvert)];
	
		%Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
		if thickness_is_present,
			thickness=GetParameterValue(pentaelem,thickness_list,gauss_coord);
		else
			thickness=GetParameterValue(pentaelem,pentaelem.h,gauss_coord);
		end

		%Get strain rate, if velocity has been supplied: 
		if velocity_is_present,
			epsilon=GetStrainRate(pentaelem,vxvy_list,xyz_list,gauss_coord);
		else
			epsilon=zeros(5,1);
		end
		if oldvelocity_is_present,
			oldepsilon=GetStrainRate(pentaelem,oldvxvy_list,xyz_list,gauss_coord);
		else
			oldepsilon=zeros(5,1);
		end

		%disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
		
		%Update material parameter that deals with ice rigidity. The update depends either on the 
		%temperature (if we are running thermal, or transient), or on B (if we are running inverse 
		%control methods on B). The latter prevails. 

		%Update material if temperature is provided.
		if temperature_is_present,
			temperature=GetParameterValue(pentaelem,temperature_list,gauss_coord);
			matice.B=paterson(temperature);
		end

		%Update material if flow law is specified. This will erase the previous change 
		%on B when temperature is provided. 
		if flow_law_is_present,
			B_param=GetParameterValue(pentaelem,B_list,gauss_coord);
			matice.B=B_param; clear B_param.
		end

		%Get viscosity at last two iterations: 
		newviscosity=GetViscosity3d(matice,epsilon);
		oldviscosity=GetViscosity3d(matice,oldepsilon);
		viscosity=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);
		%disp(sprintf('Element id %i Viscosity: %g \n',pentaelem.id,viscosity));

		%Get B and Bprime matrices:
		B=GetB_horiz(pentaelem,xyz_list,gauss_coord);
		Bprime=GetBprime_horiz(pentaelem,xyz_list,gauss_coord);
		
		%Get Jacobian determinant:
		Jdet=GetJacobianDeterminant(pentaelem,xyz_list,gauss_coord);
		%disp(sprintf('Jacobian determinant: %f\n', Jdet));

		% Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
		% onto this scalar matrix, so that we win some computational time: */
		D=viscosity*gauss_weight*Jdet*diag(ones(5,1));
	  

		% Do the triple product tB*D*Bprime: 
		Ke_gg_gaussian=B'*D*Bprime;

		%Add Ke_gg_gaussian to Ke: 
		Ke.terms=Ke.terms+Ke_gg_gaussian;
		
		end %for ig=1:num_area_gauss,
	  end %for ig=1:num_vert_gauss,

		
	  

	%Deal with 2d friction at the bedrock interface

	if (pentaelem.onbed==1 & (pentaelem.shelf==0))

		%build alpha^2 of each of the 3 grids of the triangle
		if (pentaelem.friction_type==2),  %friction_type=2 implies viscous, friction_type=1 implies plastic
		
			%Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
			frictionparameters=struct();
			frictionparameters.element_type='3d';
			frictionparameters.rho_ice=rho_ice;
			frictionparameters.rho_water=rho_water;
			frictionparameters.g=gravity;	
			frictionparameters.p=pentaelem.p;
			frictionparameters.q=pentaelem.q;

			if velocity_is_present,
				frictionparameters.velocities=vxvy_list(1:3,:);
			else
				frictionparameters.velocities=zeros(3,2);
			end
			if thickness_is_present,
				frictionparameters.h=thickness_list(1:3);
			else
				frictionparameters.h=pentaelem.h(1:3);
			end
			if bed_is_present,
				frictionparameters.b=bed_list(1:3);
			else
				frictionparameters.b=pentaelem.b(1:3);
			end
			if basal_drag_is_present,
				frictionparameters.k=K_list(1:3);
			else
				frictionparameters.k=pentaelem.k(1:3);
			end

			alpha2=Getalpha2(frictionparameters);
		end

		% Get gaussian points and weights.
		[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

		for ig=1:num_gauss2D,

			%Pick up the gaussian point and its weight:
			gauss_weight=gauss_weights(ig);
			gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
			gauss_coord_penta=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig) -1];
			
			%Get the Jacobian determinant
			xyz_list_tria=xyz_list(1:3,:);
			Jdettria=GetJacobianDeterminant(triaelem,xyz_list_tria,gauss_coord);

			%Get L if viscous basal drag present:
			L=GetL(triaelem,gauss_coord,DOFPERGRID);
	   
			%Get alpha2 on current gaussian point
			alpha2_g=GetParameterValue(triaelem,alpha2,gauss_coord);

			%Build D
			DL_scalar=alpha2_g*gauss_weight*Jdettria;
			DL=diag(DL_scalar*ones(2,1));

			%Do the triple product tL*D*L
			Ke_gg_drag_gaussian=L'*DL*L;

			%Add Ke_gg_drag_gaussian to Ke
			Ke.terms(1:6,1:6)=Ke.terms(1:6,1:6)+Ke_gg_drag_gaussian;
		end
	end%ends friction

	
end %elseif pentaelem.collapse & pentaelem.onbed
end%ends function

function Ke=CreateKMatrixPrognostic(pentaelem,grids,materials,inputs)

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	
	if pentaelem.onbed~=1,
		Ke=elemmatrix(0);
		return;
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

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

	%Just keep the first 3 grids
	xyz_list=xyz_list(1:3,:);

	%recover extra inputs from users, at current convergence iteration.
	[dt dt_is_present]=recover_input(inputs,'dt');
	[velocity_average_param velocity_average_is_present]=recover_input(inputs,'velocity_average');

	%check on all parameters
	if ((~dt_is_present) | (~velocity_average_is_present)),
		error('CreateKMatrixPrognostic error message: missing input parameters!');
	end
	
	%initialize parameter lists (for our 3 grids)
	vxvy_list=zeros(numgrids,2);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		
		Ke.row_indices(i)=doflist(1);
	
		vxvy_list(i,:)=velocity_average_param(doflist(1:2));
	end

	%Create Artificial diffusivity if asked
	if pentaelem.artificial_diffusivity,
		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,[1/3 1/3 1/3]);

		%Build K matrix (artificial diffusivity matrix)
		v_gauss=1/3*(vxvy_list(1,:)+vxvy_list(2,:)+vxvy_list(3,:));
		K=Jdettria^(1/2)/2*[abs(v_gauss(1)) 0 ; 0 abs(v_gauss(2))];
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];

		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L matrix if viscous basal drag present:
		L=GetL(triaelem,gauss_coord,DOFPERGRID);

		DL_scalar=gauss_weight*Jdettria;

		%Do the triple product tL*D*L. 
		Ke_gg_gaussian=L'*DL_scalar*L;

		%Get B matrix 
		B=GetB_prog(triaelem,xyz_list,gauss_coord);
		
		%Get Bprime matrix 
		Bprime=GetBprime_prog(triaelem,xyz_list,gauss_coord);

		%Get u, v and their derivatives at gauss point
		u_g=GetParameterValue(triaelem,vxvy_list(:,1),gauss_coord);
		v_g=GetParameterValue(triaelem,vxvy_list(:,2),gauss_coord);
		du_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,1),xyz_list,gauss_coord);
		dv_g=GetParameterDerivativeValue(triaelem,vxvy_list(:,2),xyz_list,gauss_coord);
		dux_g=du_g(1);
		dvy_g=dv_g(2);
		
		DL_scalar=dt*gauss_weight*Jdettria;

		%Create DL and DLprime matrix
		DL=DL_scalar*[dux_g, 0; 0, dvy_g];
		DLprime=DL_scalar*[u_g, 0; 0, v_g];

		%Do the triple product tL*D*L. 
		Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;	

		%Add Ke_gg_drag_gaussian and Ke_gg_gaussian to Ke
		Ke.terms=Ke.terms+Ke_gg_gaussian+Ke_gg_thickness;

		%Add artificial diffusivity  if requested
		if pentaelem.artificial_diffusivity,

			%Do the triple product tL*D*L. 
			Ke_gg_gaussian=Bprime'*DL_scalar*K*Bprime;
			
			%Add Ke_gg_drag_gaussian to Ke
			Ke.terms=Ke.terms+Ke_gg_gaussian;
		end
	end

end %end function

function Ke=CreateKMatrixMelting(pentaelem,grids,materials,inputs)

	%We are dealing with a collapsed formulation on the lower triangle of the penta, 
	%only on the bedrock.
	if pentaelem.onbed~=1,
		Ke=elemmatrix(0);
		return;
	end

	%some variables
	numgrids=3;
	DOFPERGRID=1;
	numdof=numgrids*DOFPERGRID; %number of dof for element pentaelem.

	%Create elementary stiffness matrix 
	Ke=elemmatrix(numdof);

	%recover material
	matpar=materials(end).constants;

	%recover material parameters
	latentheat=matpar.latentheat;
	heatcapacity=matpar.heatcapacity;

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

	%Just keep the first 3 grids
	xyz_list=xyz_list(1:3,:);

	%Build linear indices for elementary stiffness matrix.
	for i=1:numgrids,
		doflist=grids(pentaelem.g(i)).grid.doflist; %list of dofs in the g-set
		dof=doflist(1); %melting is 1st degree of freedom
		Ke.row_indices(i)=dof;
	end

	% Get gaussian points and weights.
	[num_gauss2D,first_area_gauss_coord,second_area_gauss_coord,third_area_gauss_coord,gauss_weights]=GaussTria(2);

	for ig=1:num_gauss2D,

		%Pick up the gaussian point and its weight:
		gauss_weight=gauss_weights(ig);
		gauss_coord=[first_area_gauss_coord(ig) second_area_gauss_coord(ig) third_area_gauss_coord(ig)];
		
		%Get the Jacobian determinant
		Jdettria=GetJacobianDeterminant(triaelem,xyz_list,gauss_coord);

		%Get L matrix if viscous basal drag present:
		L=GetL(triaelem,gauss_coord,DOFPERGRID);
   
		DL_scalar=latentheat/heatcapacity*gauss_weight*Jdettria;

		%Do the triple product tL*D*L. 
		Ke_gg_vert_gaussian=L'*DL_scalar*L;
		
		%Add Ke_gg_drag_gaussian to Ke
		Ke.terms=Ke.terms+Ke_gg_vert_gaussian;
	end

end %end function
