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

Go to the documentation of this file.

function y=fv(s,v,n)

y=zeros(n*n,1);

y=v;

h=1/(n-1.0);
RE=100;


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

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

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

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


for j=2:n-1
    
    for i=2:n-1
        
        index=(j-1)*n+i;
        
        sc=s(index);
        sw=s(index-1);
        se=s(index+1);
        sn=s(index+n);
        ss=s(index-n);
        vc=v(index);
        vw=v(index-1);
        ve=v(index+1);
        vn=v(index+n);
        vs=v(index-n);
                       
        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)=RE*(UdVORdX+VdVORdY)-(d2VORdX2+d2VORdY2);
        
                
    end
    
    
end

---