D:/demona/2.0/manual/doc/2-31set/matlab/myF.m

Go to the documentation of this file.

function y=myF(x,n)

global lambda
global re
global unk


y=x; %so corners are zero.

h=1/(n-1.0);


j=1;
for i=2:n-1
   index=((j-1)*n+i-1)*unk+1;
   sc=x(index);
   vc=x(index+1);
   sn=x(index+n*unk);
   y(index)=sc-0.0;
   y(index+1)=vc-2*(sc-sn)/h^2;
end

j=n;
for i=2:n-1
   index=((j-1)*n+i-1)*unk+1;
   sc=x(index);
   vc=x(index+1);
   ss=x(index-n*unk);
   y(index)=sc-0.0;
   y(index+1)=vc-2*(sc-ss-h)/h^2;
end

i=1;
for j=2:n-1
   index=((j-1)*n+i-1)*unk+1;
   sc=x(index);
   vc=x(index+1);
   se=x(index+unk);
   y(index)=sc-0.0;
   y(index+1)=vc-2*(sc-se)/h^2;
end

i=n;
for j=2:n-1
   index=((j-1)*n+i-1)*unk+1;
   sc=x(index);
   vc=x(index+1);
   sw=x(index-unk);
   y(index)=sc-0.0;
   y(index+1)=vc-2*(sc-sw)/h^2;
end


for j=2:n-1
    
    for i=2:n-1
        
        index=((j-1)*n+i-1)*unk+1;
        
        sc=x(index);
        sw=x(index-unk);
        se=x(index+unk);
        sn=x(index+n*unk);
        ss=x(index-n*unk);
        vc=x(index+1);
        vw=x(index+1-unk);
        ve=x(index+1+unk);
        vn=x(index+1+n*unk);
        vs=x(index+1-n*unk);
                       
        dSTFdX=(se-sw)/(2.0*h);

                dSTFdY=(sn-ss)/(2.0*h);
        
%        dSTFdX=(se-sc)/(1.0*h);

%               dSTFdY=(sc-ss)/(1.0*h);

                myU=dSTFdY;
                myV=-dSTFdX;
                        


                        d2VORdX2=(vw+ve-2.0*vc)/(h^2.d0);

                        d2VORdY2=(vn+vs-2.0*vc)/(h^2.d0);
                                
                        dVORdX=(ve-vw)/(2.0*h);

                        dVORdY=(vn-vs)/(2.0*h);
            
%                       dVORdX=(vc-vw)/(1.0*h);

%                       dVORdY=(vc-vs)/(1.0*h);

                        UdVORdX=myU*dVORdX;
                        VdVORdY=myV*dVORdY;
                        
       
        y(index)=(ss+sn+sw+se-4*sc)/h^2+vc;
        y(index+1)=re*(UdVORdX+VdVORdY)-(d2VORdX2+d2VORdY2);
        
                
    end
    
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

% 
% y=x;
% 
% h=1/(n-1.0);
% 
% for j=2:n-1
%     
%     for i=2:n-1
%         
%         index=(j-1)*n+i;
%         
%         Tc=x(index);
%         Tw=x(index-1);
%         Te=x(index+1);
%         Tn=x(index+n);
%         Ts=x(index-n);
%                        
%         y(index)=(Ts+Tn+Tw+Te-4*Tc)/h^2+lambda*exp(Tc);
%  
%     end
%    
% end
% 

---