neumannbc

 NEUMANNBC  Adds Neumann data to r.h.s.
   Updates array f (r.h.s) with Neumann data when ifro(i)>30

 [f]=neumannbc(f,h,jacx,jacy,wx,wy,xy,ifro);

 Input: f = column array of r.h.s
        h = @(xi,)[...] unction handle to the expression of Neumann
              boundary data. It is a vector of 4 functions, each for any
              side.
         jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
         jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
         xy = column array with global mesh, length: noe=nov(npdx*npdy,ne)
         ifro = column array of length noe=nov(npdx*npdy,ne):
            if (x_i,y_i) is internal to Omega then ifro(i)=0,
            if (x_i,y_i) is on \partial\Omega and Dirichlet boundary
            condition is imposed there, then ifro(i)=1,
            if (x_i,y_i) is on \partial\Omega and Neumann boundary
            condition is imposed there, then ifro(i)=31,
        nov = local -global map, previously generated by cosnov_2d

 Output: f = updated column array of r.h.s

0001 function [f]=neumannbc(f,h,jacx,jacy,wx,wy,xy,ifro,nov);
0002 % NEUMANNBC  Adds Neumann data to r.h.s.
0003 %   Updates array f (r.h.s) with Neumann data when ifro(i)>30
0004 %
0005 % [f]=neumannbc(f,h,jacx,jacy,wx,wy,xy,ifro);
0006 %
0007 % Input: f = column array of r.h.s
0008 %        h = @(xi,)[...] unction handle to the expression of Neumann
0009 %              boundary data. It is a vector of 4 functions, each for any
0010 %              side.
0011 %         jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
0012 %         jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
0013 %         xy = column array with global mesh, length: noe=nov(npdx*npdy,ne)
0014 %         ifro = column array of length noe=nov(npdx*npdy,ne):
0015 %            if (x_i,y_i) is internal to Omega then ifro(i)=0,
0016 %            if (x_i,y_i) is on \partial\Omega and Dirichlet boundary
0017 %            condition is imposed there, then ifro(i)=1,
0018 %            if (x_i,y_i) is on \partial\Omega and Neumann boundary
0019 %            condition is imposed there, then ifro(i)=31,
0020 %        nov = local -global map, previously generated by cosnov_2d
0021 %
0022 % Output: f = updated column array of r.h.s
0023 %
0024 
0025 %   Written by Paola Gervasio
0026 %   $Date: 2007/04/01$
0027 
0028 
0029 [ldnov,ne]=size(nov);
0030 npdx=length(wx);
0031 npdy=length(wy);
0032 nm1=ldnov-npdx;
0033 for ie=1:ne
0034 
0035 % side 1 of Omega_ie
0036 for i=1:npdx
0037 ip=nov(i,ie);
0038 if ifro(ip)>30
0039 hh=h(xy(ip,1),xy(ip,2));
0040 f(ip)=f(ip)+hh(1)*wx(i)*jacx(ie);
0041 end
0042 end
0043 % side 2 of Omega_ie
0044 for j=1:npdy
0045 ip=nov(j*npdx,ie);
0046 if ifro(ip)>30
0047 hh=h(xy(ip,1),xy(ip,2));
0048 f(ip)=f(ip)+hh(2)*wy(j)*jacy(ie);
0049 end
0050 end
0051 % side 3 of Omega_ie
0052 for i=1:npdx
0053 ip=nov(nm1+i,ie);
0054 if ifro(ip)>30
0055 hh=h(xy(ip,1),xy(ip,2));
0056 f(ip)=f(ip)+hh(3)*wx(i)*jacx(ie);
0057 end
0058 end
0059 % side 4 of Omega_ie
0060 for j=1:npdy
0061 ip=nov((j-1)*npdx+1,ie);
0062 if ifro(ip)>30
0063 hh=h(xy(ip,1),xy(ip,2));
0064 f(ip)=f(ip)+hh(4)*wy(j)*jacy(ie);
0065 end
0066 end
0067 
0068 end
0069 return
0070