c_velfinder

PURPOSE ^

C_VELFINDER - compute a velocity field

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %C_VELFINDER - compute a velocity field
0002 %
0003 %   This routine is used by macayealcontrol, it calculates the velocity field of with
0004 %   MacAyeal's model and computes the value of some matrices used by macayealcontrol
0005 
0006 summer=[1;1;1];
0007 
0008 rowD=zeros(nel,9);
0009 colD=zeros(nel,9);
0010 rowDshort=zeros(nel,6);
0011 colDshort=zeros(nel,6);
0012 
0013 valueuu=zeros(nel,6);
0014 valueuv=zeros(nel,9);
0015 valuevv=zeros(nel,6);
0016 
0017 
0018 nu_bar_uu=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0019 nu_bar_vv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0020 nu_bar_uv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0021 
0022 
0023 count=0;
0024 for m=1:3
0025     for k=1:3
0026       count=count+1;
0027       rowD(:,count)=index(:,m);
0028      colD(:,count)=index(:,k);
0029    end
0030 end
0031 
0032 count=0;
0033 for m=1:3
0034     for k=m:3
0035       count=count+1;
0036       rowDshort(:,count)=index(:,m);
0037      colDshort(:,count)=index(:,k);
0038    end
0039 end
0040 
0041 valuex=zeros(nel,27);
0042 valuey=zeros(nel,27);
0043 row=zeros(nel,27);
0044 
0045 count=0;
0046 
0047 for m=1:3
0048    for k=1:3
0049        for l=1:3
0050            count=count+1;
0051            row(:,count)=index(:,m);
0052            
0053            if ( (m==k) + (m==l) + (l==k) )==3
0054                fac=1/10;
0055            elseif ( (m==k) + (m==l) + (l==k) )==1
0056                fac=1/30;
0057            else
0058                fac=1/60;
0059            end
0060            
0061          valuex(:,count) = -fac*area(:) ...
0062                               .*z_thick(index(:,l)).* ...
0063                            z_surf_xbar(:);
0064          valuey(:,count) = -fac*area(:) ...
0065                              .*z_thick(index(:,l)).* ...
0066                       z_surf_ybar(:);
0067       end
0068    end
0069 end
0070 
0071 
0072 Rhs=rho_ice*g*[sparse(row,ones(size(row)),valuex,nods,1)
0073 sparse(row,ones(size(row)),valuey,nods,1)];
0074 Rhs=full(Rhs);
0075 Rhs_a_nofront=Rhs;
0076 
0077 
0078 
0079 for k=1:2 %for each of a segment's  two nodes
0080     for l=1:2
0081         for j=1:2
0082 %x-direction
0083                 Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ...
0084 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
0085 - rho_water*g/2*z_bed(index_icefront(:,l)).*z_bed(index_icefront(:,j))) ...
0086 .*normal_icefront(:,1).*length_icefront(:) ...
0087 /(4*(l==k & j==k) + 12*(l~=k | j~=k) );
0088 
0089 %y-direction
0090                 Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ...
0091 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
0092 - rho_water*g/2*z_bed(index_icefront(:,l)).*z_bed(index_icefront(:,j))) ...
0093 .*normal_icefront(:,2).*length_icefront(:) ...
0094 /(4*(l==k & j==k) + 12*(l~=k | j~=k) );
0095 
0096         end
0097     end
0098 end
0099 
0100 Rhs_a_front=Rhs-Rhs_a_nofront;
0101 
0102 
0103 
0104 
0105 
0106 
0107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%create P matrix%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0108 
0109 num_specified= 2*sum(nodes_on_dirichlet);
0110 row=zeros(2*nods-num_specified,1);
0111 col=zeros(2*nods-num_specified,1);
0112 value=zeros(2*nods-num_specified,1);
0113 count=0;
0114 for n=1:nods
0115    if(~(nodes_on_dirichlet(n)))
0116             count=count+1;
0117         row(count)=count;
0118         col(count)=n;
0119         value(count)=1;
0120     end
0121 end
0122 for n=1:nods
0123    if (~(nodes_on_dirichlet(n)))
0124          count=count+1;
0125         row(count)=count;
0126         col(count)=nods+n;
0127         value(count)=1;
0128     end
0129 end
0130 
0131 P=sparse(row,col,value,2*nods-num_specified,2*nods);
0132 
0133 
0134 specified_velocity=zeros(2*nods,1);
0135 pos=find(nodes_on_dirichlet);
0136 
0137 specified_velocity(pos)=vx_obs(pos);
0138 specified_velocity(nods+pos)=vy_obs(pos);
0139 
0140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0141 
0142 
0143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%create P0 matrix%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0144 
0145 num_specified= 2*sum(nodes_on_boundary);
0146 row=zeros(2*nods-num_specified,1);
0147 col=zeros(2*nods-num_specified,1);
0148 value=zeros(2*nods-num_specified,1);
0149 count=0;
0150 for n=1:nods
0151    if(~(nodes_on_boundary(n)))
0152             count=count+1;
0153         row(count)=count;
0154         col(count)=n;
0155         value(count)=1;
0156     end
0157 end
0158 for n=1:nods
0159    if (~(nodes_on_boundary(n)))
0160          count=count+1;
0161         row(count)=count;
0162         col(count)=nods+n;
0163         value(count)=1;
0164     end
0165 end
0166 
0167 P0=sparse(row,col,value,2*nods-num_specified,2*nods);
0168 
0169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0170 
0171 
0172 count=0;
0173    for m=1:3
0174         for k=m:3
0175          count=count+1;
0176          valueuu(:,count)=area(:).*z_thick_bar.* ...
0177             ((2*alpha(:,m).* ...
0178             alpha(:,k)+1/2*beta(:,m).* ...
0179             beta(:,k)));  
0180 
0181 
0182          valuevv(:,count)=area(:).*z_thick_bar.* ...
0183             ((2*beta(:,m).* ...
0184             beta(:,k)+1/2*alpha(:,m).* ...
0185             alpha(:,k)));
0186     end
0187 end
0188 
0189 count=0;
0190    for m=1:3
0191         for k=1:3
0192          count=count+1;
0193          valueuv(:,count)=area(:).*z_thick_bar.* ...
0194             (1/2*beta(:,m).* ...
0195             alpha(:,k)+alpha(:,m).* ...
0196             beta(:,k));
0197       end
0198    end
0199 
0200 Duu=sparse(rowDshort,colDshort,nu_bar_uu.*valueuu,nods,nods);
0201 Duu=Duu+triu(Duu,1)';
0202 
0203 Dvv=sparse(rowDshort,colDshort,nu_bar_vv.*valuevv,nods,nods);
0204 Dvv=Dvv+triu(Dvv,1)';
0205 
0206 F=[Duu sparse(rowD,colD,nu_bar_uv.*valueuv,nods,nods)
0207    sparse(colD,rowD,nu_bar_uv.*valueuv,nods,nods) Dvv]; 
0208 
0209 
0210 Rhs_parsed=P*(Rhs - F*specified_velocity);
0211 Fprime=P*F*P';
0212 
0213 
0214 solution=Fprime\Rhs_parsed;
0215 solution=P'*solution + specified_velocity;
0216 u=solution(1:nods);
0217 v=solution(nods+1:2*nods);
0218

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003