Description of HADC

0001 function [ADC_DE,ADC_DE_allcmpts,elapsed_time] ...
0002     = HADC(experiment,mymesh,DIFF_cmpts,IC_cmpts)
0003 
0004 % diffusion equation (zero IC) to get the time-dependent diffusion coefficient
0005 %
0006 % Input:
0007 %     1. experiment is a structure with 8 elements:
0008 %         ngdir_total
0009 %         gdir
0010 %         sdeltavec
0011 %         bdeltavec
0012 %         seqvec
0013 %         npervec
0014 %         rtol
0015 %         atol
0016 %     2. mymesh is a structure with 10 elements:
0017 %         Nnode
0018 %         Nele
0019 %         Nface
0020 %         Pts_cmpt_reorder
0021 %         Ele_cmpt_reorder
0022 %         Pts_ind
0023 %         Pts_boundary_reorder
0024 %         Fac_boundary_reorder
0025 %         Nboundary
0026 %         Ncmpt
0027 %     3. DIFF_cmpts
0028 %     4. IC_cmpts
0029 %
0030 % Output:
0031 %     1. ADC_DE
0032 %     2. ADC_DE_allcmpts
0033 %     3. elapsed_time
0034 
0035 gdir = experiment.gdir;
0036 sdeltavec = experiment.sdeltavec;
0037 bdeltavec = experiment.bdeltavec;
0038 seqvec = experiment.seqvec;
0039 npervec = experiment.npervec;
0040 rtol = experiment.rtol;
0041 atol = experiment.atol;
0042 
0043 global FEM_M FEM_K FEM_A FEM_Q FEM_G
0044 global QVAL UG
0045 global BDELTA SDELTA SEQ OGSEPER
0046 
0047 disp(['H-ADC Model: solving DE']);
0048 
0049 UG = gdir';
0050 UG = UG/norm(UG);
0051         
0052 for icmpt = 1:mymesh.Ncmpt
0053     [VOL(icmpt)] ...
0054         = get_volume_mesh(mymesh.Pts_cmpt_reorder{icmpt},mymesh.Ele_cmpt_reorder{icmpt});
0055 end
0056 
0057 VOL_allcmpts = 0;
0058 
0059 for icmpt = 1:mymesh.Ncmpt
0060     VOL_allcmpts  = VOL_allcmpts + VOL(icmpt);
0061 end
0062 
0063 for icmpt = 1:mymesh.Ncmpt
0064     VOL_frac(icmpt) = VOL(icmpt)/VOL_allcmpts;
0065 end
0066 
0067 nexperi = length(sdeltavec);
0068 elapsed_time=zeros(mymesh.Ncmpt, nexperi);
0069 
0070 
0071 for icmpt = 1:mymesh.Ncmpt
0072     
0073 
0074         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0075         % To replace the pde matrices
0076         coordinates = mymesh.Pts_cmpt_reorder{icmpt}; 
0077         elements = mymesh.Ele_cmpt_reorder{icmpt};    
0078         [FEM_K,volumes]=stiffness_matrixP1_3D(elements',coordinates',DIFF_cmpts(icmpt));
0079         FEM_M=mass_matrixP1_3D(elements',volumes);   
0080         FEM_G=sparse(size(FEM_K,1),1);
0081         FEM_A=sparse(size(FEM_K,1),size(FEM_K,2));
0082         FEM_Q=sparse(size(FEM_K,1),size(FEM_K,2)); 
0083         coordinates_t=coordinates';
0084         one = sparse(size(FEM_K,1),1);
0085         one(:,:) = 1;
0086 
0087         for iboundary = 1:mymesh.Nboundary
0088             neumann = mymesh.Fac_boundary_reorder{icmpt}{iboundary}';
0089             if sum(size(neumann))>0
0090                 [FacA,FacC,FacN] = get_surfacenormal_mesh(mymesh.Pts_cmpt_reorder{icmpt},mymesh.Ele_cmpt_reorder{icmpt},neumann');
0091                 mycoeff = (FacN(1,:)*UG(1) + FacN(2,:)*UG(2)+FacN(3,:)*UG(3));
0092                 GG = flux_matrixP1_3D(neumann,coordinates', DIFF_cmpts(icmpt)*mycoeff');     
0093                 FEM_G = FEM_G + GG*one;
0094             end;
0095         end;
0096         
0097         % MyM = model_FEM_matrices{icmpt}.M;
0098         % MyK = model_FEM_matrices{icmpt}.K;
0099         % MyA = model_FEM_matrices{icmpt}.A;
0100         % MyQ = model_FEM_matrices{icmpt}.Q;
0101         % MyG = DIFF_cmpts(icmpt)*model_FEM_matrices{icmpt}.G;
0102         % errorG=sum((MyG - FEM_G).^2);
0103         % errorK=sum(sum((MyK - FEM_K).^2));
0104         % errorM=sum(sum((MyM - FEM_M).^2));
0105         % if errorG+errorK+errorM>1e-7
0106             % disp('Matrices are different!');
0107             % disp(['Error M:', num2str(errorM),', Error K:', num2str(errorK), ', Error G:', num2str(errorG)]);
0108             % % [full(MyG),  full(FEM_G)]
0109             % stop;
0110         % end;
0111         % disp(['Error M:', num2str(errorM),', Error K:', num2str(errorK), ', Error G:', num2str(errorG)]);
0112         % To replace the pde matrices
0113         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0114    
0115     
0116     ODEsolve_atol = atol;
0117     ODEsolve_rtol = rtol;
0118     
0119     options = odeset('Mass',FEM_M,'AbsTol',ODEsolve_atol,'RelTol',...
0120         ODEsolve_rtol,'Vectorized','on','Stats','off',...
0121         'Jacobian',@odejac_bt_includeb);
0122     disp(['    Compartment ',num2str(icmpt)]);
0123     nexperi = length(sdeltavec);  % Is it correct to be here???
0124     for iexperi = 1:nexperi
0125         disp(['      Experiment ',num2str(iexperi)]);
0126         iex_start_time = clock;
0127         SDELTA = sdeltavec(iexperi);
0128         BDELTA = bdeltavec(iexperi);
0129         SEQ = seqvec(iexperi);% for choosing case PGSE, OGSEcos or OGSEsin
0130         omega = 2*pi*npervec(iexperi)/SDELTA;
0131         OGSEPER = 1./omega*2*pi;%% set up number for OGSE
0132         QVAL = 0;
0133         TLIST = [0,SDELTA+BDELTA];
0134         ICC = zeros(size(FEM_M,1),1);
0135         sol = ode23t(@odefun_bt_includeb,TLIST,ICC,options);
0136        
0137         %deff_PDE_formulation_src{iexperi}{icmpt} = FEM_G.'*sol.y/VOL(icmpt)/VOL(icmpt)/DIFF_cmpts(icmpt);
0138         deff_PDE_formulation_src{iexperi}{icmpt} = FEM_G.'*sol.y/VOL(icmpt)/VOL(icmpt);
0139 
0140         deff_PDE_formulation_src_time{iexperi}{icmpt} = sol.x;
0141 
0142         hvec = deff_PDE_formulation_src{iexperi}{icmpt};
0143         tvec11 = deff_PDE_formulation_src_time{iexperi}{icmpt};
0144         Ftvec11 = seqintprofile(tvec11);
0145         a = trapz(tvec11,Ftvec11.*hvec*VOL(icmpt))/trapz(tvec11,Ftvec11.^2);
0146         
0147         ADC_DE(icmpt,iexperi) = DIFF_cmpts(icmpt)-a;
0148         elapsed_time(icmpt, iexperi) = etime(clock, iex_start_time);
0149     end
0150 end
0151 ADC_DE_allcmpts = nan*ones(nexperi,1);
0152 for iexperi = 1:nexperi
0153     ADC_DE_allcmpts(iexperi,1) = sum((IC_cmpts.*VOL_frac)'.*ADC_DE(:,iexperi))./sum((IC_cmpts.*VOL_frac)');
0154 end
HADC

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
