Warning, /utils/matlab/Graphix/GraphixUtility/calc_Bolus_CS.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 892fe9e2 on 2005-10-19 20:09:14 UTC
892fe9e2f3 Dani*0001 function [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
0002 % [ub,vb]=calc_Bolus_CS(rAc,dXc,dYc,dXg,dYg,delR,kwx,kwy,hFacW,hFacS);
0003 %-- GM-Bolus transp. :
0004 % from Kwx & Kwy (Skew-Flx => /2),
0005 % compute Volume Stream function psiX,psiY above uVel.vVel
0006 % (at interface between 2 levels), units=m^3/s :
0007 % psiX=(rAc*kwx)_i / dXc ; psiY=(rAc*kwy)_j / dYc ;
0008 % and then the bolus velocity (m/s):
0009 % ub = d_k(psiX)/dYg/drF ; vb = d_k(psiY)/dXg/drF ;
0010 %---------------------------------------------------------------------
0011 % set_axis
0012 [ncx,nc]=size(rAc); nr=length(delR);
0013 nPg=6*nc*nc; ncp=nc+1;
0014
0015 rAu=dXc.*dYg; rAu=reshape(rAu,nPg,1);
0016 rAv=dYc.*dXg; rAv=reshape(rAv,nPg,1);
0017
0018 %-- K*rAc + add 1 overlap :
0019 tmp=reshape(rAc,nPg,1)*ones(1,nr);
0020 kwx=reshape(kwx,nPg,nr);
0021 kwy=reshape(kwy,nPg,nr);
0022 tmp=0.5*tmp;
0023 kwx=tmp.*kwx;
0024 kwy=tmp.*kwy;
0025 kwx=reshape(kwx,ncx,nc,nr);
0026 kwy=reshape(kwy,ncx,nc,nr);
0027 v6X=split_C_cub(kwx,1);
0028 v6Y=split_C_cub(kwy,1);
0029 k6x=v6X(:,[2:ncp],:,:);
0030 k6y=v6Y([2:ncp],:,:,:);
0031
0032 %--- recip_hFac & mask :
0033 hFacW=reshape(hFacW,nPg,nr);
0034 hFacS=reshape(hFacS,nPg,nr);
0035 rhFcW=1./hFacW; rhFcW(find(hFacW==0))=0;
0036 rhFcS=1./hFacS; rhFcS(find(hFacS==0))=0;
0037 hFacW=ceil(hFacW); hFacW=min(1,hFacW);
0038 hFacS=ceil(hFacS); hFacS=min(1,hFacS);
0039
0040 %-----------------
0041
0042 v6X=zeros(nc,nc,nr,6);
0043 v6X([1:nc],:,:,:)=k6x([2:ncp],:,:,:)+k6x([1:nc],:,:,:);
0044 v6X=v6X/2;
0045 psiX=zeros(ncx,nc,nr+1);
0046 for n=1:6
0047 is=1+nc*(n-1);ie=nc*n;
0048 psiX([is:ie],[1:nc],[1:nr])=v6X([1:nc],[1:nc],[1:nr],n);
0049 end
0050 psiX=reshape(psiX,nPg,nr+1);
0051 psiX(:,[1:nr])=hFacW.*psiX(:,[1:nr]);
0052 ub=zeros(nPg,nr);
0053 ub=psiX(:,[2:nr+1])-psiX(:,[1:nr]);
0054 delR = reshape(delR,[1,length(delR)]);
0055 tmp=rAu*delR;
0056 ub=ub./tmp; ub=rhFcW.*ub;
0057 ub=reshape(ub,ncx,nc,nr);
0058
0059 %-----------------
0060
0061 v6Y=zeros(nc,nc,nr,6);
0062 v6Y(:,[1:nc],:,:)=k6y(:,[2:ncp],:,:)+k6y(:,[1:nc],:,:);
0063 v6Y=v6Y/2;
0064 psiY=zeros(ncx,nc,nr+1);
0065 for n=1:6
0066 is=1+nc*(n-1);ie=nc*n;
0067 psiY([is:ie],[1:nc],[1:nr])=v6Y([1:nc],[1:nc],[1:nr],n);
0068 end
0069 psiY=reshape(psiY,nPg,nr+1);
0070 psiY(:,[1:nr])=hFacS.*psiY(:,[1:nr]);
0071 vb=zeros(nPg,nr);
0072 vb=psiY(:,[2:nr+1])-psiY(:,[1:nr]);
0073 tmp=rAv*delR;
0074 vb=vb./tmp; vb=rhFcS.*vb;
0075 vb=reshape(vb,ncx,nc,nr);
0076
0077 %-----------------
0078 return
