Warning, /utils/matlab/cs_grid/calcBolusVelCube.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 81aa3ab1 on 2006-08-18 00:46:40 UTC
7ab08d6865 Dani*0001 function [ub,vb]=calcBolusVelCube(d,g);
0002
0003 % [ub,vb] = calcEulerPsiCube(d,g);
0004 %
0005 % Input arguements:
0006 % The incoming field data (d) and grid data (g) must be in a structured
0007 % array format (which is the format that comes from rdmnc):
0008 % d [Field data] Kwx,Kwy
0009 % g [Grid data ] drF,rA,dxC,dyC,dxG,dyG,HFacW,HFacS
0010 %
0011 % Original author: Jean-Michel Campin
0012 % Modifications: Daniel Enderton
0013 %
0014 % JMCs old comments:
0015 %-- GM-Bolus transp. :
0016 % from Kwx & Kwy (Skew-Flx => /2),
0017 % compute Volume Stream function psiX,psiY above uVel.vVel
0018 % (at interface between 2 levels), units=m^3/s :
0019 % psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ;
0020 % and then the bolus velocity (m/s):
0021 % ub = d_k(psiX)/dYg/drF ; vb = d_k(psiY)/dXg/drF ;
0022 %---------------------------------------------------------------------
0023
0024
0025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0026 % Prepare / reform incoming data %
0027 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0028
0029 nc = size(g.XC,2);
0030 nr = length(g.drF);
81aa3ab116 Davi*0031 nt = size(d.GM_Kwx,4);
7ab08d6865 Dani*0032
0033 dr = g.drF;
0034 hw = reshape(g.HFacW(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0035 hs = reshape(g.HFacS(1:6*nc,1:nc,1:nr),[6*nc*nc,nr]);
0036 ra = reshape(g.rA(1:6*nc,1:nc) ,[6*nc*nc,1]);
0037 dxc = reshape(g.dxC(1:6*nc,1:nc),[6*nc*nc,1]);
0038 dyc = reshape(g.dyC(1:6*nc,1:nc),[6*nc*nc,1]);
0039 dxg = reshape(g.dxG(1:6*nc,1:nc),[6*nc*nc,1]);
0040 dyg = reshape(g.dyG(1:6*nc,1:nc),[6*nc*nc,1]);
81aa3ab116 Davi*0041 kwx_all = reshape(d.GM_Kwx,[6*nc*nc,nr,nt]);
0042 kwy_all = reshape(d.GM_Kwy,[6*nc*nc,nr,nt]);
7ab08d6865 Dani*0043
0044 rAu=dxc.*dyg;
0045 rAv=dyc.*dxg;
81aa3ab116 Davi*0046 %--- recip_hFac & mask :
0047 hw_recip=1./hw; hw_recip(find(hw==0))=0;
0048 hs_recip=1./hs; hs_recip(find(hs==0))=0;
0049 hw=ceil(hw); hw=min(1,hw);
0050 hs=ceil(hs); hs=min(1,hs);
7ab08d6865 Dani*0051
0052 for it = 1:nt
0053 kwx = kwx_all(:,:,it);
0054 kwy = kwy_all(:,:,it);
0055
0056 %-- K*ra + add 1 overlap :
0057 kwx = 0.5*(ra*ones(1,nr)).*kwx;
0058 kwy = 0.5*(ra*ones(1,nr)).*kwy;
0059 kwx = reshape(kwx,[6*nc,nc,nr]);
0060 kwy = reshape(kwy,[6*nc,nc,nr]);
0061 v6X = split_C_cub(kwx,1);
0062 v6Y = split_C_cub(kwy,1);
0063 k6x = v6X(:,[2:nc+1],:,:);
0064 k6y = v6Y([2:nc+1],:,:,:);
0065
0066 %-----------------
0067 v6X = zeros(nc,nc,nr,6);
0068 v6Y = zeros(nc,nc,nr,6);
0069
0070 v6X([1:nc],:,:,:) = k6x([2:nc+1],:,:,:) + k6x([1:nc],:,:,:);
0071 v6Y(:,[1:nc],:,:) = k6y(:,[2:nc+1],:,:) + k6y(:,[1:nc],:,:);
0072
0073 v6X = v6X/2;
0074 v6Y = v6Y/2;
0075
0076 psiX = zeros(6*nc,nc,nr+1);
0077 psiY = zeros(6*nc,nc,nr+1);
0078
0079 for n = 1:6
0080 is = 1+nc*(n-1);ie=nc*n;
0081 psiX([is:ie],[1:nc],[1:nr]) = v6X([1:nc],[1:nc],[1:nr],n);
0082 psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n);
0083 end
0084
0085 psiX = reshape(psiX,6*nc*nc,nr+1);
0086 psiY = reshape(psiY,6*nc*nc,nr+1);
0087
0088 psiX(:,[1:nr]) = hw.*psiX(:,[1:nr]);
0089 psiY(:,[1:nr]) = hs.*psiY(:,[1:nr]);
0090 ub = psiX(:,[2:nr+1]) - psiX(:,[1:nr]);
0091 vb = psiY(:,[2:nr+1]) - psiY(:,[1:nr]);
0092
0093 dr = reshape(dr,[1,length(dr)]);
0094 ub = reshape(hw_recip.*ub./(rAu*dr),[6*nc,nc,nr]);
0095 vb = reshape(hs_recip.*vb./(rAv*dr),[6*nc,nc,nr]);
0096
0097 ub_all(:,:,:,it) = ub;
0098 vb_all(:,:,:,it) = vb;
0099 end
0100
0101 ub = ub_all;
81aa3ab116 Davi*0102 vb = vb_all;
