function Ke=CreateKMatrix(triaelem,grids,materials,inputs,analysis_type);
%CREATEKMATRIX - create the stiffmess matrix for triaelem
%
%   this stiffness matrix works for MacAyeal's model, prognostic
%   and for the computation of the bed slope and surface slope
%
%   Usage:
%      Ke=CreateKMatrix(triaelem,grids,materials,inputs,analysis_type)
%
%   See also CREATEPVECTOR

	if strcmpi(analysis_type,'diagnostic_horiz'),

		Ke=CreateKMatrixHoriz(triaelem,grids,materials,inputs);

	elseif strcmpi(analysis_type,'diagnostic_vert'),

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

	elseif strcmpi(analysis_type,'temperature'),

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

	elseif strcmpi(analysis_type,'prognostic'),

		Ke=CreateKMatrixPrognostic(triaelem,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(triaelem,grids,materials,inputs);

	end
end %end function


function Ke=CreateKMatrixVert(triaelem,grids,materials,inputs)
error('CreateKMatrix error message: Vertical analysis not implemented yet');
end %end function


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

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

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

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

	%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(triaelem.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 once for all if requested
	if triaelem.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 triaelem.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=CreateKMatrixThermal(triaelem,grids,materials,inputs)
error('CreateKMatrix error message: Thermal analysis not implemented yet');
end %end function


function Ke=CreateKMatrixSlopeCompute(triaelem,grids,materials,inputs)

	%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(triaelem,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(triaelem.g(i)).grid.doflist; %list of dofs in the g-set
		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=CreateKMatrixHoriz(triaelem,grids,materials,inputs);

global element_debug  element_debugid


%Check if MacAyeal's acceleration is not on? 
if triaelem.acceleration==1,
	Ke=elemmatrix(0);
	return;
end

%some variables
numgrids=3;
DOFPERGRID=2;
numdof=numgrids*DOFPERGRID; %number of dof for 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(triaelem.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_average_param temperature_average_is_present]=recover_input(inputs,'temperature_average');
[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
vxvy_list=zeros(numgrids,2);
oldvxvy_list=zeros(numgrids,2);

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_average_list=zeros(numgrids,1);


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


%Build linear indices for elementary stiffness matrix.
for i=1:numgrids,
	doflist=grids(triaelem.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);
	else
		K_list(i)=triaelem.k(i);
	end
	if(thickness_is_present),
		thickness_list(i)=thickness_param(dof);
	else
		thickness_list(i)=triaelem.h(i);
	end
	if(surface_is_present),
		surface_list(i)=surface_param(dof);
	else
		surface_list(i)=triaelem.s(i);
	end
	if(bed_is_present),
		bed_list(i)=bed_param(dof);
	else
		bed_list(i)=triaelem.b(i);
	end
	if(temperature_average_is_present) temperature_average_list(i)=temperature_average_param(dof);end;
end
if (element_debug & triaelem.id==element_debugid),
	disp(sprintf('El id %i TriaElemnet input list before gaussian loop: \n',element_debugid)); 
	disp(sprintf('   rho_ice: %g ',rho_ice));
	disp(sprintf('   gravity: %g',gravity));
	disp(sprintf('   rho_water: %g',rho_water));
	disp(sprintf('   Velocity: '));
	for i=1:3,
		disp(sprintf('      grid %i  [%g,%g]',i,vxvy_list(i,1),vxvy_list(i,2)));
	end
	disp(sprintf('   B [%g %g %g ]',B_list(1),B_list(2),B_list(3)));
	disp(sprintf('   K [%g %g %g ]',K_list(1),K_list(2),K_list(3)));
	disp(sprintf('   thickness [%g %g %g]',thickness_list(1),thickness_list(2),thickness_list(3)));
	disp(sprintf('   surface [%g %g %g ]',surface_list(1),surface_list(2),surface_list(3)));
	disp(sprintf('   bed [%g %g %g]',bed_list(1),bed_list(2),bed_list(3)));
	if(temperature_average_is_present)disp(sprintf('   temperature_average [%g %g %g]',temperature_average_list(1),temperature_average_list(2),temperature_average_list(3)));end;
end



alpha2_list=zeros(3,1);
%Build alpha2 used by drag staffness matrix
if (~triaelem.shelf & triaelem.friction_type==2)
    
	%Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
	frictionparameters=struct();
	frictionparameters.element_type='2d';
	frictionparameters.rho_ice=rho_ice;
	frictionparameters.rho_water=rho_water;
	frictionparameters.g=gravity;	
	frictionparameters.p=triaelem.p;
	frictionparameters.q=triaelem.q;

	 if velocity_is_present,
		frictionparameters.velocities=vxvy_list;
	else
		frictionparameters.velocities=zeros(3,2);
	end
	frictionparameters.h=thickness_list;
	frictionparameters.b=bed_list;
	frictionparameters.k=K_list;

	alpha2=Getalpha2(frictionparameters);
	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('   alpha2_list (%g %g %g )',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
	end
end

% Get gaussian points and weights. 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) 

[num_gauss,first_gauss_area_coord,second_gauss_area_coord,third_gauss_area_coord,gauss_weights]=GaussTria(2);
    
if (element_debug & triaelem.id==element_debugid),
	disp(sprintf('   gaussian points: '));
	for i=1:num_gauss,
		disp(sprintf('    %g %g %g : %g',first_gauss_area_coord(i),second_gauss_area_coord(i),third_gauss_area_coord(i),gauss_weights(i)));
	end
end
    
%Start  looping on the number of gaussian points:
for ig=1:num_gauss,
	%Pick up the gaussian point and its weight:
	gauss_weight=gauss_weights(ig);
	gauss_l1l2l3=[first_gauss_area_coord(ig) second_gauss_area_coord(ig) third_gauss_area_coord(ig)];

	%Compute thickness at gaussian point from t1,t2 and t3 fields in the element itself
	if thickness_is_present,
		thickness=GetParameterValue(triaelem,thickness_list,gauss_l1l2l3);
	else
		thickness=GetParameterValue(triaelem,triaelem.h,gauss_l1l2l3);
	end
	%disp(sprintf('Thickness: %f\n', thickness));

	%If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
	%velocity should be = 0. To achieve this result, we set k to a very low value.
	slopevector=zeros(2,1);
	if(~triaelem.shelf),
	
		if surface_is_present,
			slopevector=GetParameterDerivativeValue(triaelem,surface_list,xyz_list,gauss_l1l2l3);
		else
			slopevector=GetParameterDerivativeValue(triaelem,triaelem.s,xyz_list,gauss_l1l2l3);
		end
		slope_magnitude=sqrt(slopevector(1)^2+slopevector(2)^2);

		if (slope_magnitude>MAXSLOPE),
			K=10^(-MOUNTAINKEXPONENT);
		end
	end
	
	%Get strain rate, if velocity has been supplied: 
	if velocity_is_present,
 		epsilon=GetStrainRate(triaelem,vxvy_list,xyz_list,gauss_l1l2l3);
	else
 		epsilon=zeros(3,1);
	end
	if oldvelocity_is_present,
 		oldepsilon=GetStrainRate(triaelem,oldvxvy_list,xyz_list,gauss_l1l2l3);
	else
 		oldepsilon=zeros(3,1);
	end

	%disp(sprintf('Epsilon: %f %f %f\n',epsilon(1),epsilon(2),epsilon(3)));
	
	%Update material if temperature is provided.
	if temperature_average_is_present,
		temperature_average=GetParameterValue(triaelem,temperature_average_list,gauss_l1l2l3);
		matice.B=paterson(temperature_average);
	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(triaelem,B_list,gauss_l1l2l3);
		matice.B=B_param; clear B_param.
	end


	%Get viscosity at last two iterations: 
 	newviscosity=GetViscosity2d(matice,epsilon);
 	oldviscosity=GetViscosity2d(matice,oldepsilon);
	viscosity=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);

	%Get Jacobian determinant:
	Jdet=GetJacobianDeterminant(triaelem,xyz_list,gauss_l1l2l3);

	% 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*thickness*gauss_weight*Jdet*diag(ones(numgrids,1));

	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('   gaussian loop %i\n',ig));
		disp(sprintf('      thickness %g',thickness));
		disp(sprintf('      slope (%g,%g)',slopevector(1),slopevector(2)));
		disp(sprintf('      alpha2_list (%g,%g,%g)',alpha2_list(1),alpha2_list(2),alpha2_list(3)));
		disp(sprintf('      epsilon (%g,%g,%g)',epsilon(1),epsilon(2),epsilon(3)));
		disp(sprintf('      Matice: '));
		matice
		disp(sprintf('\n      viscosity: %g ',viscosity));
		disp(sprintf('      jacobian: %g ',Jdet));
		disp(sprintf('      gauss_weight: %g ',gauss_weight));
	end


	%Get B and Bprime matrices:
	B=GetB(triaelem,xyz_list,gauss_l1l2l3);
	Bprime=GetBprime(triaelem,xyz_list,gauss_l1l2l3);

	%Get L matrix if viscous basal drag present:
	L=zeros(2,DOFPERGRID*numgrids);
	if (triaelem.friction_type==2 && triaelem.shelf==0),
		L=GetL(triaelem,gauss_l1l2l3,DOFPERGRID);
	end

	
	% 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;

	%Now, take care of the basal friction if there is any
	alpha2_g=0;
	if (~triaelem.shelf) & (triaelem.friction_type==2),
 
		%compute alpha2 for the current gaussian point
		alpha2_g=GetParameterValue(triaelem,alpha2,gauss_l1l2l3);
        
		if velocity_is_present
			DL_scalar=alpha2_g*gauss_weight*Jdet;
		else
			DL_scalar=0;
		end
        
		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=Ke.terms+Ke_gg_drag_gaussian;
	end

	if (element_debug & triaelem.id==element_debugid),
		disp(sprintf('      alpha2 %g\n',alpha2_g));
		disp(sprintf('B:\n'));
		B
		disp(sprintf('Bprime:\n'));
		Bprime
		disp(sprintf('L:\n'));
		L
	end

end %for ig=1:num_gauss,
	
if (element_debug & triaelem.id==element_debugid),
	disp(sprintf('Ke_gg->terms:\n'));
	Ke.terms
	disp(sprintf('Ke_gg->row_indices:\n'));
	Ke.row_indices'
end

end %end function

