function Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type);
%CREATEKMATRIX - create the stiffness matrix for acceleratedtriaelem
%
%   global finite element stiffness matrix of the whole model
%   works only for MacAyeal's model
%
%   Usage:
%      Ke=CreateKMatrix(element,grids,materials,inputs,analysis_type)
%
%   See also CREATEPVECTOR


%keep permanent pieces of stiffness matrix in memory
global permanent_pieces_of_Ke alpha beta row_location col_location aire

%Ok, this is the accelerated MacAyeal's element. Very particular methodology here! 

%Some variables
numdof=3*2; %3 grids, 2 degrees of freedom per grid.

%recover extra inputs from users, at current convergence iteration.
[flow_law_param flow_law_is_present]=recover_input(inputs,'B');
[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');
[velocity_param velocity_is_present]=recover_input(inputs,'velocity');
[oldvelocity_param oldvelocity_is_present]=recover_input(inputs,'oldvelocity');

if velocity_is_present,
	vx=velocity_param(1:6:numdof*element.nods);
	vy=velocity_param(2:6:numdof*element.nods);
else
	vx={};
	vy={};
	permanent_pieces_of_Ke={};
end
if oldvelocity_is_present,
	oldvx=oldvelocity_param(1:6:numdof*element.nods);
	oldvy=oldvelocity_param(2:6:numdof*element.nods);
else
	oldvx={};
	oldvy={};
end
if flow_law_is_present,
	B_list=flow_law_param(1:6:6*element.nods);
	B_bar=1/3*(B_list(element.index(:,1))+B_list(element.index(:,2))+B_list(element.index(:,3)));
else
	B_bar=element.B_bar;
end

if basal_drag_is_present,
	drag=basal_drag_param(1:6:6*element.nods);
else
	drag=element.drag;
end

if thickness_is_present,
	thickness=thickness_param(1:6:6*element.nods);
else
	thickness=element.thickness;
end

thickness_el=thickness(element.index)*[1;1;1]/3;

if surface_is_present,
	surface=surface_param(1:6:6*element.nods);
else
	surface=element.surface;
end

if bed_is_present,
	bed=bed_param(1:6:6*element.nods);
else
	bed=element.bed;
end

if isempty(permanent_pieces_of_Ke),
	%Start building areas
	aire=zeros(element.nel,1);

	for n=1:element.nel
		aire(n)=1/2 * det([1 1 1;element.x(element.index(n,:))';element.y(element.index(n,:))']);
	end

	aire=abs(aire); % if index is sorted from its original value, then aire could be negative

	alpha=zeros(element.nel,3);
	beta=zeros(element.nel,3);
	gamma=zeros(element.nel,3);

	for n=1:element.nel
		X=inv([element.x(element.index(n,:)) element.y(element.index(n,:)) ones(3,1)]);
		alpha(n,:)=X(1,:);
		beta(n,:)=X(2,:);
		gamma(n,:)=X(3,:);
	end

	clear X


	%  Do once and for all the initial computation of matrix-locations:
	row_location=zeros(4*element.nel*9,1);
	col_location=zeros(4*element.nel*9,1);

	count=-element.nel+1;

	for i=1:3
		for j=1:3
			for k=1:2,
				for l=1:2,
					count=count+element.nel;
					row_location(count:count+element.nel-1)=(element.index(:,i)-1)*numdof+k;
					col_location(count:count+element.nel-1)=(element.index(:,j)-1)*numdof+l;
				end
			end
		end
	end


	% Create permanent part of stiffness matrix.
	permanent_pieces_of_Ke=zeros(2*element.nel*9,1);

	count=-element.nel+1;

	for i=1:3
		for j=1:3
			for k=1:2,
				for l=1:2,
					count=count+element.nel;

					permanent_pieces_of_Ke(count:count+element.nel-1)= thickness_el .* aire.*(...
					((k==1) & (l==1)).*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j)) +...
					((k==2) & (l==1)).*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i)) + ...
					((k==1) & (l==2)).*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)) + ...
					((k==2) & (l==2)).*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j)) ...
					);
				end
			end
		end
	end
	
end %if isempty(permanent_pieces_of_Ke)
	
%Compute viscosity
viscosity_overshoot=materials(end).constants.viscosity_overshoot;
newviscosity=viscosity(element.index,element.nel,alpha,beta,vx,vy,B_bar,element.glen_coeff);
oldviscosity=viscosity(element.index,element.nel,alpha,beta,oldvx,oldvy,B_bar,element.glen_coeff);
nu_bar=newviscosity+viscosity_overshoot*(newviscosity-oldviscosity);

%Multiply viscosity by permanent_pieces_of_Ke
count=-element.nel+1;

%Create sparse global matrix terms: 
Ke_terms=zeros(4*element.nel*9,1);

for i=1:3
	for j=1:3
		for k=1:2,
			for l=1:2,
				count=count+element.nel;
				Ke_terms(count:(count+element.nel-1))=nu_bar.*permanent_pieces_of_Ke(count:(count+element.nel-1));
			end
		end
	end
end

Drag_operator=sparse2(6*element.nods,6*element.nods);

%Create basal stiffness
if (element.friction_type==2) & velocity_is_present

	%average of p and q over the grids (size nel->nods)
	pcoeff_grid=zeros(element.nods,1);
	qcoeff_grid=zeros(element.nods,1);
	for i=1:element.nods
		%1: find the elements that contain the grid i
		neighbors_gridi=[];
		for j=1:3
			neighbors_gridi=[neighbors_gridi find(element.index(:,j)==i)'];
		end
		numberofneighbors_gridi=length(neighbors_gridi);
		%2 retrieve the value of p and q over each of these elements. The average is
		%plugged into dbx_grid
		qcoeff_grid(i)=sum(element.q(neighbors_gridi))/numberofneighbors_gridi;
		pcoeff_grid(i)=sum(element.p(neighbors_gridi))/numberofneighbors_gridi;
	end
    
	%Retrieve some parameters needed to compute alpha2 (taub=alpha2*ub)
	frictionparameters=struct();
	frictionparameters.element_type='2d';
	frictionparameters.rho_ice=materials(end).constants.rho_ice;
	frictionparameters.rho_water=materials(end).constants.rho_water;
	frictionparameters.g=materials(end).constants.g;	
	frictionparameters.p=pcoeff_grid;
	frictionparameters.q=qcoeff_grid;
	frictionparameters.velocities=[vx vy];
	frictionparameters.h=thickness;
	frictionparameters.b=bed;
	frictionparameters.k=drag;

	alpha2=Getalpha2(frictionparameters);

	drag_operator_value=zeros(2*element.nel*27,1);

	row_location_drag=zeros(2*element.nel*27,1);
	col_location_drag=zeros(2*element.nel*27,1);

	count=-element.nel+1;
	    
	for i=1:3
		for j=1:3
			for m=1:3
				for k=1:2,

					count=count+element.nel;
			 
					row_location_drag(count:count+element.nel-1)=(element.index(:,i)-1)*numdof+k;
					col_location_drag(count:count+element.nel-1)=(element.index(:,j)-1)*numdof+k;
					drag_operator_value(count:count+element.nel-1)=aire.*(...
					((i==j) & (j==m)).*alpha2(element.index(:,m))*1/10+...
					(((i==j) | (j==m) | (i==m)) & ~((i==j) & (j==m))).*alpha2(element.index(:,m))*1/30+...
					((i~=j) & (j~=m) & (i~=m)).*alpha2(element.index(:,m))*1/60);

				end
			end
		end
	end  

	Drag_operator=sparse2(row_location_drag,col_location_drag,drag_operator_value,6*element.nods,6*element.nods);
end

%Create final global stiffness matrix
Ke=globalmatrix;
Ke.nrows=numdof*element.nods;
Ke.matrix=Drag_operator+sparse2(row_location,col_location,Ke_terms,6*element.nods,6*element.nods);
