Warning, /utils/matlab/cs_grid/split_UV_cub.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 6b4f212f on 2024-06-05 16:55:17 UTC
896b20e91e Jean*0001 function [u6t,v6t] = split_UV_cub(u3d,v3d,ksign,kad)
0002 % [u6t,v6t] = split_UV_cub(u3d,v3d,[ksign,kad])
0003 %---------------------------------------------
3b83d20fc0 Jean*0004 % Split 2d/3d vector field u,v (on C-grid if kad > 0, on A-grid if kad < 0)
0005 % into 3d x 6 faces and:
0006 % kad=1: (C-grid) add 1 column to U and one row to V <== at the end !!!
0007 % => output is u6t(nc+1,nc,[nr],6) & v6t(nc,nc+1,[nr],6)
6b4f212f67 Jean*0008 % kad=2: add 1 column + 2 rows to U and 2 columns + 1 row to V
3b83d20fc0 Jean*0009 % => output is u6t(nc+1,nc+2,[nr],6) & v6t(nc+2,nc+1,[nr],6)
6b4f212f67 Jean*0010 % kad=3: add 2 columns + 2 rows to U and V (i.e, 1 on both side)
0011 % => output is u6t(nc+2,nc+2,[nr],6) & v6t(nc+2,nc+2,[nr],6)
0012 % kad=-1: assuming input u3d & v3d are on A-grid,
3b83d20fc0 Jean*0013 % add 1 column to U and one row to V <== at the begining !!!
0014 % => output is u6t(nc+1,nc,[nr],6) & v6t(nc,nc+1,[nr],6)
6b4f212f67 Jean*0015 % ksign=0 : --> no sign change ; ksign=1 : --> change the sign where needed
0016 % note: ksign used only if kad=2 (or =3)
896b20e91e Jean*0017 %----------------------------------------------
4ec37fd829 Jean*0018 % Written by jmc@mit.edu, 2005.
0019
896b20e91e Jean*0020 if nargin < 3 , ksign = 0; end
0021 if nargin < 4 , kad = 1; end
6b4f212f67 Jean*0022 if rem(kad,1) ~= 0 | kad < -1 | kad > 3, fprintf('kad= %f => Bad value\n',kad); return; end
3b83d20fc0 Jean*0023 plmn = 1 - 2*ksign ;
896b20e91e Jean*0024
4ec37fd829 Jean*0025 %------------------------------
0026 if ~isequal(size(u3d),size(v3d))
0027 error(['Error in Input-array dimensions: ',...
0028 'u = ',mat2str(size(u3d)),', ','v = ',mat2str(size(v3d))]);
0029 end
0030 dims=size(u3d); nDim=length(dims); n1h=dims(1); n2h=dims(2);
896b20e91e Jean*0031 %fprintf(' nDim= %i , dims:',nDim);fprintf(' %i',dims);fprintf('\n');
4ec37fd829 Jean*0032 if nDim == 2, nr=1; else nr=prod(dims(3:end)); end
0033 if n1h == 6*n2h, nc=n2h;
0034 u3d=permute(reshape(u3d,[nc 6 nc nr]),[1 3 4 2]);
0035 v3d=permute(reshape(v3d,[nc 6 nc nr]),[1 3 4 2]);
0036 elseif n1h*6 == n2h, nc=n1h;
0037 u3d=permute(reshape(u3d,[nc nc 6 nr]),[1 2 4 3]);
0038 v3d=permute(reshape(v3d,[nc nc 6 nr]),[1 2 4 3]);
0039 else
0040 error([' input var size: ',int2str(n1h),' x ',int2str(n2h),' does not fit regular cube !']);
0041 end
0042 ncp=nc+1 ; n2p=nc+2; %nye=nc+kad;
0043 %--
896b20e91e Jean*0044 %=================================================================
0045
0046 if kad == 0,
0047 %- split on to 6 faces with no overlap:
6b4f212f67 Jean*0048 dims(1)=nc; dims(2)=nc;
896b20e91e Jean*0049 u6t=u3d;
0050 v6t=v3d;
0051
0052 elseif kad == 1,
0053
6b4f212f67 Jean*0054 dims(1)=ncp; dims(2)=nc;
896b20e91e Jean*0055 %- split on to 6 faces with overlap in i+1 for u and j+1 for v :
0056 u6t=zeros(ncp,nc,nr,6);
0057 v6t=zeros(nc,ncp,nr,6);
0058 u6t([1:nc],:,:,:)=u3d(:,:,:,:);
0059 v6t(:,[1:nc],:,:)=v3d(:,:,:,:);
0060
0061 u6t(ncp,[1:nc],:,1)=u3d(1,[1:nc],:,2);
0062 u6t(ncp,[1:nc],:,3)=u3d(1,[1:nc],:,4);
0063 u6t(ncp,[1:nc],:,5)=u3d(1,[1:nc],:,6);
6b4f212f67 Jean*0064 u6t(ncp,[1:nc],:,2)=v3d([nc:-1:1],1,:,4);
0065 u6t(ncp,[1:nc],:,4)=v3d([nc:-1:1],1,:,6);
896b20e91e Jean*0066 u6t(ncp,[1:nc],:,6)=v3d([nc:-1:1],1,:,2);
0067
0068 v6t([1:nc],ncp,:,1)=u3d(1,[nc:-1:1],:,3);
0069 v6t([1:nc],ncp,:,3)=u3d(1,[nc:-1:1],:,5);
0070 v6t([1:nc],ncp,:,5)=u3d(1,[nc:-1:1],:,1);
6b4f212f67 Jean*0071 v6t([1:nc],ncp,:,2)=v3d([1:nc],1,:,3);
0072 v6t([1:nc],ncp,:,4)=v3d([1:nc],1,:,5);
3b83d20fc0 Jean*0073 v6t([1:nc],ncp,:,6)=v3d([1:nc],1,:,1);
896b20e91e Jean*0074
6b4f212f67 Jean*0075 elseif kad == 2 | kad == 3,
896b20e91e Jean*0076
6b4f212f67 Jean*0077 %- fill for kad=2 case:
0078 dims(1)=ncp; dims(2)=n2p;
896b20e91e Jean*0079 %- split on to 6 faces:
0080 u6t=zeros(ncp,n2p,nr,6);
0081 v6t=zeros(n2p,ncp,nr,6);
0082 u6t([1:nc],[2:ncp],:,:)=u3d(:,:,:,:);
0083 v6t([2:ncp],[1:nc],:,:)=v3d(:,:,:,:);
0084
0085 %- add overlap in i=nc+1 for u and j=nc+1 for v :
0086 u6t(ncp,[2:ncp],:,1)=u3d(1,[1:nc],:,2);
0087 u6t(ncp,[2:ncp],:,3)=u3d(1,[1:nc],:,4);
0088 u6t(ncp,[2:ncp],:,5)=u3d(1,[1:nc],:,6);
6b4f212f67 Jean*0089 u6t(ncp,[2:ncp],:,2)=v3d([nc:-1:1],1,:,4);
0090 u6t(ncp,[2:ncp],:,4)=v3d([nc:-1:1],1,:,6);
896b20e91e Jean*0091 u6t(ncp,[2:ncp],:,6)=v3d([nc:-1:1],1,:,2);
0092
0093 v6t([2:ncp],ncp,:,1)=u3d(1,[nc:-1:1],:,3);
0094 v6t([2:ncp],ncp,:,3)=u3d(1,[nc:-1:1],:,5);
0095 v6t([2:ncp],ncp,:,5)=u3d(1,[nc:-1:1],:,1);
6b4f212f67 Jean*0096 v6t([2:ncp],ncp,:,2)=v3d([1:nc],1,:,3);
0097 v6t([2:ncp],ncp,:,4)=v3d([1:nc],1,:,5);
3b83d20fc0 Jean*0098 v6t([2:ncp],ncp,:,6)=v3d([1:nc],1,:,1);
896b20e91e Jean*0099
0100 %- add overlap in j=0,nc+1 for u and i=0,nc+1 for v :
0101 u6t([1:ncp], 1 ,:,1)=u6t([1:ncp],ncp,:,6);
0102 u6t([1:ncp], 1 ,:,3)=u6t([1:ncp],ncp,:,2);
0103 u6t([1:ncp], 1 ,:,5)=u6t([1:ncp],ncp,:,4);
6b4f212f67 Jean*0104 u6t([1:ncp],n2p,:,1)=v6t( 2 ,[ncp:-1:1],:,3)*plmn;
0105 u6t([1:ncp],n2p,:,3)=v6t( 2 ,[ncp:-1:1],:,5)*plmn;
896b20e91e Jean*0106 u6t([1:ncp],n2p,:,5)=v6t( 2 ,[ncp:-1:1],:,1)*plmn;
0107
0108 u6t([1:ncp], 1 ,:,2)=v6t(ncp,[ncp:-1:1],:,6)*plmn;
0109 u6t([1:ncp], 1 ,:,4)=v6t(ncp,[ncp:-1:1],:,2)*plmn;
0110 u6t([1:ncp], 1 ,:,6)=v6t(ncp,[ncp:-1:1],:,4)*plmn;
6b4f212f67 Jean*0111 u6t([1:ncp],n2p,:,2)=u6t([1:ncp], 2 ,:,3);
0112 u6t([1:ncp],n2p,:,4)=u6t([1:ncp], 2 ,:,5);
896b20e91e Jean*0113 u6t([1:ncp],n2p,:,6)=u6t([1:ncp], 2 ,:,1);
0114
0115 v6t( 1 ,[1:ncp],:,1)=u6t([ncp:-1:1],ncp,:,5)*plmn;
0116 v6t( 1 ,[1:ncp],:,3)=u6t([ncp:-1:1],ncp,:,1)*plmn;
0117 v6t( 1 ,[1:ncp],:,5)=u6t([ncp:-1:1],ncp,:,3)*plmn;
6b4f212f67 Jean*0118 v6t(n2p,[1:ncp],:,1)=v6t( 2 ,[1:ncp],:,2);
0119 v6t(n2p,[1:ncp],:,3)=v6t( 2 ,[1:ncp],:,4);
896b20e91e Jean*0120 v6t(n2p,[1:ncp],:,5)=v6t( 2 ,[1:ncp],:,6);
0121
0122 v6t( 1 ,[1:ncp],:,2)=v6t(ncp,[1:ncp],:,1);
0123 v6t( 1 ,[1:ncp],:,4)=v6t(ncp,[1:ncp],:,3);
0124 v6t( 1 ,[1:ncp],:,6)=v6t(ncp,[1:ncp],:,5);
6b4f212f67 Jean*0125 v6t(n2p,[1:ncp],:,2)=u6t([ncp:-1:1], 2 ,:,4)*plmn;
0126 v6t(n2p,[1:ncp],:,4)=u6t([ncp:-1:1], 2 ,:,6)*plmn;
896b20e91e Jean*0127 v6t(n2p,[1:ncp],:,6)=u6t([ncp:-1:1], 2 ,:,2)*plmn;
0128
6b4f212f67 Jean*0129 if kad == 3,
0130 %- finish the kad=3 case:
0131 dims(1)=n2p;
0132 var=u6t; u6t=zeros(n2p,n2p,nr,6); u6t(2:n2p,:,:,:)=var;
0133 var=v6t; v6t=zeros(n2p,n2p,nr,6); v6t(:,2:n2p,:,:)=var;
0134 clear var
0135
0136 u6t(1,[2:ncp],:,1)=v6t([ncp:-1:2],ncp,:,5);
0137 u6t(1,[2:ncp],:,3)=v6t([ncp:-1:2],ncp,:,1);
0138 u6t(1,[2:ncp],:,5)=v6t([ncp:-1:2],ncp,:,3);
0139 u6t(1,[2:ncp],:,2)=u6t(ncp,[2:ncp],:,1);
0140 u6t(1,[2:ncp],:,4)=u6t(ncp,[2:ncp],:,3);
0141 u6t(1,[2:ncp],:,6)=u6t(ncp,[2:ncp],:,5);
0142
0143 v6t([2:ncp],1,:,1)=v6t([2:ncp],ncp,:,6);
0144 v6t([2:ncp],1,:,3)=v6t([2:ncp],ncp,:,2);
0145 v6t([2:ncp],1,:,5)=v6t([2:ncp],ncp,:,4);
0146 v6t([2:ncp],1,:,2)=u6t(ncp,[ncp:-1:2],:,6);
0147 v6t([2:ncp],1,:,4)=u6t(ncp,[ncp:-1:2],:,2);
0148 v6t([2:ncp],1,:,6)=u6t(ncp,[ncp:-1:2],:,4);
0149 end
0150
3b83d20fc0 Jean*0151 elseif kad == -1,
0152
6b4f212f67 Jean*0153 dims(1)=ncp; dims(2)=nc;
3b83d20fc0 Jean*0154 %- split on to 6 faces with overlap in i-1 for u and j-1 for v :
0155 u6t=zeros(ncp,nc,nr,6);
0156 v6t=zeros(nc,ncp,nr,6);
0157 u6t([2:ncp],:,:,:)=u3d(:,:,:,:);
0158 v6t(:,[2:ncp],:,:)=v3d(:,:,:,:);
0159
0160 u6t(1,[1:nc],:,1)=v6t([nc:-1:1],ncp,:,5);
0161 u6t(1,[1:nc],:,3)=v6t([nc:-1:1],ncp,:,1);
0162 u6t(1,[1:nc],:,5)=v6t([nc:-1:1],ncp,:,3);
0163 u6t(1,[1:nc],:,2)=u6t(ncp,[1:nc],:,1);
0164 u6t(1,[1:nc],:,4)=u6t(ncp,[1:nc],:,3);
0165 u6t(1,[1:nc],:,6)=u6t(ncp,[1:nc],:,5);
0166
0167 v6t([1:nc],1,:,1)=v6t([1:nc],ncp,:,6);
0168 v6t([1:nc],1,:,3)=v6t([1:nc],ncp,:,2);
0169 v6t([1:nc],1,:,5)=v6t([1:nc],ncp,:,4);
0170 v6t([1:nc],1,:,2)=u6t(ncp,[nc:-1:1],:,6);
0171 v6t([1:nc],1,:,4)=u6t(ncp,[nc:-1:1],:,2);
0172 v6t([1:nc],1,:,6)=u6t(ncp,[nc:-1:1],:,4);
0173
896b20e91e Jean*0174 end
0175
0176 %- restaure the right shape:
0177 if nDim == 2,
0178 u6t=squeeze(u6t);
0179 v6t=squeeze(v6t);
0180 else
0181 u6t=reshape(u6t,[dims 6]);
0182 v6t=reshape(v6t,[dims(2) dims(1) dims(3:end) 6]);
0183 end
0184
0185 return
