0001
0002
0003
0004
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
0080 for l=1:2
0081 for j=1:2
0082
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
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
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
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