quicwind/saddlepoint/cell_A_B.m

   1 function [A,B] = cell_A_B(dx,dy,z,i,j,k,moduli)
   2 % B=matrix_B(X)
   3 % in
   4 %   dx      spacing of the grid in the x direction
   5 %   dy      spacing of the grid in the y direction
   6 %   z       matrix of the height of the nodes
   7 %   i,j,k   cell indices
   8 %   moduli  moduli penalizations for each dimension (size 3)
   9 % out
  10 %   B     multiplication by B transforms fluxes through the faces of the
  11 %   cell to 2 wind vectors in Cartesian coordinates
  12
  13 % A is a diagonal matrix with volume and moduli
  14 % volum
  15 v = dx*dy*(z(i,j,k+1)-z(i,j,k)+z(i+1,j,k+1)-z(i+1,j,k)+z(i+1,j+1,k+1)-z(i+1,j+1,k)+z(i,j+1,k+1)-z(i,j+1,k))/4;
  16 A = v*diag([moduli(1),moduli(1),moduli(2),moduli(2),moduli(3),moduli(3)]);
  17
  18 % v = B*u
  19 % u = flux through faces left x1, right x1, front x2, back x2, bottom x3, top x3 - always low to high coordinates
  20 % v = v1(1),v2(1),v1(2),v2(2),v1(3),v2(3) cartesian components of the wind vectors
  21 % X{1} = x1 components of node coordinates, size 2x2x2, low-high x1,x2,x3
  22 % X{2} = x2 components of node coordinates, size 2x2x2, low-high x1,x2,x3
  23 % X{3} = x3 components of node coordinates, size 2x2x2, low-high x1,x2,x3
  24
  25 % tangents
  26 txk = (z(i+1,j,k)-z(i,j,k)+z(i+1,j+1,k)-z(i,j+1,k))/(2*dx); % tangent of thetax at k
  27 tyk = (z(i,j+1,k)-z(i,j,k)+z(i+1,j+1,k)-z(i+1,j,k))/(2*dy); % tangent of thetay at k
  28 txk1 = (z(i+1,j,k+1)-z(i,j,k+1)+z(i+1,j+1,k+1)-z(i,j+1,k+1))/(2*dx); % tangent of thetax at k+1
  29 tyk1 = (z(i,j+1,k+1)-z(i,j,k+1)+z(i+1,j+1,k+1)-z(i+1,j,k+1))/(2*dy); % tangent of thetay at k+1
  30 % secants
  31 sck = sqrt((1+txk*txk)*(1+tyk*tyk)); % sec(thetax)*sec(thetay) at k
  32 sck1 = sqrt((1+txk1*txk1)*(1+tyk1*tyk1)); % sec(thetax)*sec(thetay) at k+1
  33 % areas at the faces
  34 a1 = dy*.5*(z(i,j,k+1)-z(i,j,k)+z(i,j+1,k+1)-z(i,j+1,k)); % left face area
  35 a2 = dy*.5*(z(i+1,j,k+1)-z(i+1,j,k)+z(i+1,j+1,k+1)-z(i+1,j+1,k)); % right face area
  36 a3 = dx*.5*(z(i,j,k+1)-z(i,j,k)+z(i+1,j,k+1)-z(i+1,j,k)); % front face area
  37 a4 = dx*.5*(z(i,j+1,k+1)-z(i,j+1,k)+z(i+1,j+1,k+1)-z(i+1,j+1,k)); % front face area
  38 a5 = (dx*dy)*sck; % bottom face area
  39 a6 = (dx*dy)*sck1; % top face area
  40 a = [a1,a2,a3,a4,a5,a6];
  41 % columns are: left, right, front, back, bottom, and top faces
  42 % rows are: v1(1), v2(1), v1(2), v2(2), v1(3), and v2(3) wind coordinates
  43 B =  ...
  44         [-1,0,0,0,0,0;          % 1: v1(1) from flux through the left face
  45          0,1,0,0,0,0;           % 2: v2(1) from flux through the right face
  46          0,0,-1,0,0,0;          % 3: v1(2) from flux through the front face
  47          0,0,0,1,0,0;           % 4: v2(2) from flux through the back face
  48          -txk,0,-tyk,0,-1,0;    % 5: v1(3) from flux through the bottom face
  49          0,txk1,0,tyk1,0,1]...  % 6: v2(3) from flux through the top face
  50          * diag(1./a);
  51 if false,
  52     Binv = diag(a)* ...
  53             [-1,0,0,0,0,0;          % 1: flux through the left face from v1
  54              0,1,0,0,0,0;           % 2: flux through the right face from v2
  55              0,0,-1,0,0,0;          % 3: flux through the front face from v1
  56              0,0,0,1,0,0;           % 4: flux through the back face from v2
  57              txk,0,tyk,0,-1,0;      % 5: flux through the bottom face from v1
  58              0,-txk1,0,-tyk1,0,1];  % 6: flux through the top face from v2
  59     err_B_inv=norm(B*Binv-eye(6),'fro')
  60 end
  61
  62 end