0001 function [J,u,v]=objectivefunction_C(search_scalar,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,...
0002 nods,nel,vx_obs,vy_obs,u_objective,v_objective,index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,...
0003 valueuv,matrix_bar,criterion)
0004
0005
0006
0007
0008
0009
0010
0011
0012 B_new=B+search_scalar*direction;
0013
0014
0015 B_bar=(B_new(index(:,1))+B_new(index(:,2))+B_new(index(:,3)))/3;
0016
0017
0018
0019 nu_bar=viscosity(index,nel,alpha,beta,u_objective,v_objective,B_bar,glen_coeff);
0020
0021
0022 converged_yet=0;
0023 convergence_count=1;
0024 loop=0;
0025
0026
0027 while (~converged_yet)
0028 clear u v
0029
0030
0031 nu_bar_uu=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0032 nu_bar_vv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0033 nu_bar_uv=[nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar nu_bar];
0034
0035 if (loop>100)
0036 u=u_old; v=v_old;
0037 break;
0038 end;
0039 nu_bar_oldvalue=nu_bar;
0040 loop=loop+1;
0041
0042 Duu=sparse(rowDshort,colDshort,nu_bar_uu.*valueuu,nods,nods);
0043 Duu=Duu+triu(Duu,1)';
0044
0045 Dvv=sparse(rowDshort,colDshort,nu_bar_vv.*valuevv,nods,nods);
0046 Dvv=Dvv+triu(Dvv,1)';
0047
0048 F=[Duu sparse(rowD,colD,nu_bar_uv.*valueuv,nods,nods)
0049 sparse(colD,rowD,nu_bar_uv.*valueuv,nods,nods) Dvv];
0050
0051
0052 Rhs_parsed=P*(Rhs - F*specified_velocity);
0053 Fprime=P*F*P';
0054
0055
0056 solution=Fprime\Rhs_parsed;
0057 solution=P'*solution + specified_velocity;
0058 u=solution(1:nods);
0059 v=solution(nods+1:2*nods);
0060 nu_bar=viscosity(index,nel,alpha,beta, u,v,B_bar,glen_coeff);
0061
0062
0063 if convergence_count>1,
0064
0065 dug=[u; v]-[u_old; v_old];
0066 ndug=norm(dug,2);
0067 ug=[u_old;v_old];
0068 nug=norm(ug,2);
0069
0070 change=ndug/nug;
0071
0072
0073 if change>criterion,
0074 converged_yet=0;
0075 else
0076 converged_yet=1;
0077 end
0078 else
0079 converged_yet=0;
0080 end
0081
0082 u_old=u;
0083 v_old=v;
0084
0085 convergence_count=convergence_count+1;
0086
0087 end
0088
0089
0090
0091 save F_file F
0092
0093 Du=[(u-vx_obs).*weighting
0094 (v-vy_obs).*weighting];
0095
0096
0097 J=Du'*S*Du;