quicwind/cell_sizes.m

   1 function s=cell_sizes(X)
   2 % s=cell_sizes(X)
   3 % compute averaged sizes and areas projected on cartesian planes of cells in a deformed
   4 % and vertical stretched mesh
   5 %
   6 % input:
   7 %   X   cell array, nodal coordinatees x y z
   8 % output
   9 %   s structure contains arrays
  10 % same shape as u, v, w:
  11 %    dz_at_u   averaged increment in z on u-sides etc
  12 %    dy_at_u
  13 %     area_u   area of u-side
  14 %    dz_at_v   etc
  15 %    dx_at_v
  16 %     area_v
  17 %    dx_at_w
  18 %    dy_at_w
  19 %     area_w
  20 %   weight_u
  21 %   weight_v
  22 %   weight_w
  23 %  same shape as x,y,z:
  24 %    depth_x  increment in x from x side to next
  25 %    depth_y  etc
  26 %    depth_z
  27
  28 check_mesh(X);
  29
  30 x = X{1}; y = X{2}; z = X{3};
  31 [nx1,ny1,nz1] = size(x);
  32 nx = nx1-1; ny = ny1-1; nz = nz1-1;
  33
  34
  35 %                                             ^ z,w,k
  36 %     (i,j+1,k+1)---------(i+1,j+1,k+1        |
  37 %     /  |               / |                  |    /
  38 %    /   |              /  |                  |   / y,v,j
  39 %   /    |             /   |                  |  /
  40 %(i,j,k+1)----------(i+1,j,k+1)               | /
  41 %  |    /|            |    |                  |--------> x,u,i
  42 %  |  dy |            |    |
  43 %  | /   |            |    |
  44 %  |/    |            |    |
  45 %  |   (i,j+1,k)---------(i+1,j+1,k)
  46 %  |   /              |  /
  47 %  |  /               | /
  48 %  | /                |/
  49 %(i,j,k)----------(i+1,j,k)
  50
  51 % cell dimensions - distances between midpoints
  52 s.dz_at_u = zeros(nx+1,ny,nz);  % dz for u etc.
  53 s.dy_at_u = zeros(nx+1,ny,nz);  % dz for u etc.
  54 s.area_u  = zeros(nx+1,ny,nz);  % dz for u etc.
  55 for k=1:nz
  56     for j=1:ny
  57         for i=1:nx+1
  58             s.dz_at_u(i,j,k) = 0.5*(z(i,j,k+1)-z(i,j,k)+z(i,j+1,k+1)-z(i,j+1,k));
  59             s.dy_at_u(i,j,k) = 0.5*(y(i,j+1,k)-y(i,j,k)+y(i,j+1,k+1)-y(i,j,k+1));
  60             s.area_u (i,j,k) = s.dy_at_u(i,j,k)*s.dz_at_u(i,j,k);
  61         end
  62     end
  63 end
  64
  65
  66 s.dz_at_v = zeros(nx,ny+1,nz);
  67 s.dx_at_v = zeros(nx,ny+1,nz);
  68 s.area_v  = zeros(nx,ny+1,nz);
  69 for k=1:nz
  70     for j=1:ny+1
  71         for i=1:nx
  72             s.dz_at_v(i,j,k) = 0.5*(z(i,j,k+1)-z(i,j,k)+z(i+1,j,k+1)-z(i+1,j,k));
  73             s.dx_at_v(i,j,k) = 0.5*(x(i+1,j,k)-x(i,j,k)+x(i+1,j,k+1)-x(i,j,k+1));
  74             s.area_v(i,j,k)=s.dx_at_v(i,j,k)*s.dz_at_v(i,j,k);
  75         end
  76     end
  77 end
  78
  79 s.dx_at_w = zeros(nx,ny,nz+1);
  80 s.dy_at_w = zeros(nx,ny,nz+1);
  81 s.area_w = zeros(nx,ny,nz+1);
  82 for k=1:nz+1
  83     for j=1:ny
  84         for i=1:nx
  85             s.dx_at_w(i,j,k) = 0.5*(x(i+1,j,k)-x(i,j,k)+x(i+1,j,k)-x(i,j+1,k));
  86             s.dy_at_w(i,j,k) = 0.5*(y(i,j+1,k)-y(i,j,k)+y(i+1,j+1,k)-y(i+1,j,k));
  87             s.area_w(i,j,k) = s.dx_at_w(i,j,k) *  s.dy_at_w(i,j,k);
  88         end
  89     end
  90 end
  91
  92 % averaged depth of cells in the 3 directions
  93 s.depth_x=zeros(nx,ny,nz);
  94 s.depth_y=zeros(nx,ny,nz);
  95 s.depth_z=zeros(nx,ny,nz);
  96 for l1=0:1
  97     for l2=0:1
  98         for k=1:nz
  99             for j=1:ny
 100                 for i=1:nx
 101                     s.depth_x(i,j,k)=s.depth_x(i,j,k)+0.25*(x(i+1,j+l1,k+l2)-x(i,j+l1,k+l2));
 102                     s.depth_y(i,j,k)=s.depth_y(i,j,k)+0.25*(y(i+l1,j+1,k+l2)-y(i+l1,j,k+l2));
 103                     s.depth_z(i,j,k)=s.depth_z(i,j,k)+0.25*(z(i+l1,j+l2,k+1)-z(i+l1,j+l2,k));
 104                 end
 105             end
 106         end
 107     end
 108 end
 109
 110 s.weight_u=zeros(nx+1,ny,nz);
 111 s.weight_v=zeros(nx,ny+1,nz);
 112 s.weight_w=zeros(nx,ny,nz+1);
 113 for k=1:nz
 114     for j=1:ny
 115         for i=1:nx
 116             s.weight_u(i,j,k)=0.5*s.depth_x(i,j,k)*(s.area_u(i+1,j,k)+s.area_u(i,j,k));
 117         end
 118         i=nx+1;
 119         s.weight_u(i,j,k)=s.depth_x(i-1,j,k)*s.area_u(i,j,k);
 120     end
 121 end
 122 for k=1:nz
 123     for i=1:nx
 124         for j=1:ny
 125             s.weight_v(i,j,k)=0.5*s.depth_y(i,j,k)*(s.area_v(i,j+1,k)+s.area_v(i,j,k));
 126         end
 127         j=ny+1;
 128         s.weight_v(i,j,k)=s.depth_y(i,j-1,k)*s.area_v(i,j,k);
 129     end
 130 end
 131 for i=1:nx
 132    for j=1:ny
 133        for k=1:nz
 134             s.weight_w(i,j,k)=0.5*s.depth_z(i,j,k)*(s.area_w(i,j,k+1)+s.area_w(i,j,k));
 135        end
 136        k=nz+1;
 137        s.weight_w(i,j,k)=s.depth_z(i,j,k-1)*s.area_w(i,j,k);
 138    end
 139 end
 140
 141
 142 end