0001 function nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 if (isempty(u) & isempty(v)),
0012 nu_bar=10^14*ones(nel,1);
0013 else
0014 nu_bar=zeros(nel,1);
0015
0016 power=(glen_coeff-1)./(2*glen_coeff);
0017 [ux ,uy ,vx, vy]=shear(index,alpha,beta,u,v);
0018 second_inv=(ux.^2 + vy.^2 + ((uy+vx).^2)/4 + ux.*vy);
0019
0020 location=find(second_inv~=0);
0021 nu_bar(location)=B_bar(location)./(second_inv(location).^power(location));
0022
0023 location=find(second_inv==0 & power~=0);
0024 nu_bar(location)=10^18;
0025
0026 location=find(second_inv==0 & power==0);
0027 nu_bar(location)=B_bar(location);
0028 end