%C_VELFINDER - compute a velocity field
%
%   This routine is used by macayealcontrol, it calculates the velocity field of with
%   MacAyeal's model and computes the value of some matrices used by macayealcontrol

summer=[1;1;1];

rowD=zeros(nel,9);
colD=zeros(nel,9);
rowDshort=zeros(nel,6);
colDshort=zeros(nel,6);

valueuu=zeros(nel,6);
valueuv=zeros(nel,9);
valuevv=zeros(nel,6);


nu_bar_uu=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
nu_bar_vv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
nu_bar_uv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];


count=0;
for m=1:3
	for k=1:3
	  count=count+1;
	  rowD(:,count)=index(:,m);
     colD(:,count)=index(:,k);
   end
end

count=0;
for m=1:3
	for k=m:3
	  count=count+1;
	  rowDshort(:,count)=index(:,m);
     colDshort(:,count)=index(:,k);
   end
end

valuex=zeros(nel,27);
valuey=zeros(nel,27);
row=zeros(nel,27);

count=0;

for m=1:3
   for k=1:3
	   for l=1:3
		   count=count+1;
		   row(:,count)=index(:,m);
		   
		   if ( (m==k) + (m==l) + (l==k) )==3
			   fac=1/10;
		   elseif ( (m==k) + (m==l) + (l==k) )==1
			   fac=1/30;
		   else
			   fac=1/60;
		   end
		   
         valuex(:,count) = -fac*area(:) ...
	  						.*z_thick(index(:,l)).* ...
	                       z_surf_xbar(:);
         valuey(:,count) = -fac*area(:) ...
	  					   .*z_thick(index(:,l)).* ...
	                  z_surf_ybar(:);
      end
   end
end


Rhs=rho_ice*g*[sparse(row,ones(size(row)),valuex,nods,1)
sparse(row,ones(size(row)),valuey,nods,1)];
Rhs=full(Rhs);
Rhs_a_nofront=Rhs;



for k=1:2 %for each of a segment's  two nodes
	for l=1:2
		for j=1:2
%x-direction			
				Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ...
(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
- rho_water*g/2*z_bed(index_icefront(:,l)).*z_bed(index_icefront(:,j))) ...
.*normal_icefront(:,1).*length_icefront(:) ...
/(4*(l==k & j==k) + 12*(l~=k | j~=k) );

%y-direction
				Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ...
(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
- rho_water*g/2*z_bed(index_icefront(:,l)).*z_bed(index_icefront(:,j))) ...
.*normal_icefront(:,2).*length_icefront(:) ...
/(4*(l==k & j==k) + 12*(l~=k | j~=k) );

		end
	end
end

Rhs_a_front=Rhs-Rhs_a_nofront;






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%create P matrix%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_specified= 2*sum(nodes_on_dirichlet);
row=zeros(2*nods-num_specified,1);
col=zeros(2*nods-num_specified,1);
value=zeros(2*nods-num_specified,1);
count=0;
for n=1:nods
   if(~(nodes_on_dirichlet(n)))
   		 count=count+1;
		row(count)=count;
		col(count)=n;
		value(count)=1;
	end
end
for n=1:nods
   if (~(nodes_on_dirichlet(n)))
         count=count+1;
		row(count)=count;
		col(count)=nods+n;
		value(count)=1;
	end
end

P=sparse(row,col,value,2*nods-num_specified,2*nods);


specified_velocity=zeros(2*nods,1);
pos=find(nodes_on_dirichlet);

specified_velocity(pos)=vx_obs(pos);
specified_velocity(nods+pos)=vy_obs(pos);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%create P0 matrix%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

num_specified= 2*sum(nodes_on_boundary);
row=zeros(2*nods-num_specified,1);
col=zeros(2*nods-num_specified,1);
value=zeros(2*nods-num_specified,1);
count=0;
for n=1:nods
   if(~(nodes_on_boundary(n)))
   		 count=count+1;
		row(count)=count;
		col(count)=n;
		value(count)=1;
	end
end
for n=1:nods
   if (~(nodes_on_boundary(n)))
         count=count+1;
		row(count)=count;
		col(count)=nods+n;
		value(count)=1;
	end
end

P0=sparse(row,col,value,2*nods-num_specified,2*nods);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


count=0;
   for m=1:3
		for k=m:3
         count=count+1;
         valueuu(:,count)=area(:).*z_thick_bar.* ...
            ((2*alpha(:,m).* ...
            alpha(:,k)+1/2*beta(:,m).* ...
            beta(:,k)));  


         valuevv(:,count)=area(:).*z_thick_bar.* ...
            ((2*beta(:,m).* ...
            beta(:,k)+1/2*alpha(:,m).* ...
            alpha(:,k)));
    end
end

count=0;
   for m=1:3
		for k=1:3
         count=count+1;
         valueuv(:,count)=area(:).*z_thick_bar.* ...
            (1/2*beta(:,m).* ...
            alpha(:,k)+alpha(:,m).* ...
            beta(:,k));
      end
   end

Duu=sparse(rowDshort,colDshort,nu_bar_uu.*valueuu,nods,nods);
Duu=Duu+triu(Duu,1)';

Dvv=sparse(rowDshort,colDshort,nu_bar_vv.*valuevv,nods,nods);
Dvv=Dvv+triu(Dvv,1)';

F=[Duu sparse(rowD,colD,nu_bar_uv.*valueuv,nods,nods)
   sparse(colD,rowD,nu_bar_uv.*valueuv,nods,nods) Dvv]; 


Rhs_parsed=P*(Rhs - F*specified_velocity);
Fprime=P*F*P';


solution=Fprime\Rhs_parsed;
solution=P'*solution + specified_velocity;
u=solution(1:nods);
v=solution(nods+1:2*nods);

