|
0
|
1 % Sparse function test routing |
|
|
2 % Copyright (C) 1998 Andy Adler. |
|
|
3 % This code has no warranty whatsoever. |
|
|
4 % You may do what you like with this code as long as you leave this copyright |
|
|
5 % in place. If you modify the code then include a notice saying so. |
|
|
6 % |
|
|
7 % This code comes from my thesis work. Its a finite element model |
|
|
8 % of the voltage distribution in a cylindrical body with current applied |
|
|
9 % at the boundary. The complete model was published in |
|
|
10 % "Electrical Impedance Tomography: Regularized Imaging and Contrast |
|
|
11 % Detection", A.Adler, R.Guardo, IEEE Trans Med Img, 15(2),170-179 |
|
|
12 % |
|
|
13 % $Id$ |
|
|
14 |
|
|
15 % octave commands, but OK in matlab too |
|
|
16 do_fortran_indexing= 1; |
|
|
17 empty_list_elements_ok = 1; |
|
|
18 |
|
|
19 TESTTIMES= 5; |
|
|
20 dotheplot=0; |
|
|
21 |
|
|
22 %increase N by multiples of 4 to do a larger test |
|
|
23 N= 8; |
|
|
24 %increase niveaux to do a larger test |
|
|
25 niveaux= [-.2:.1:.2]' ; |
|
|
26 |
|
|
27 pos_i= [0 1]; |
|
|
28 elec=16; |
|
|
29 ELEM=[]; |
|
|
30 NODE= [0;0]; |
|
|
31 int=1; |
|
|
32 for k=1:N |
|
|
33 phi= (0:4*k-1)*pi/2/k; |
|
|
34 NODE= [NODE k/N*[sin(phi);cos(phi)]]; |
|
|
35 |
|
|
36 ext= 2*(k*k-k+1); |
|
|
37 idxe=[0:k-1; 1:k]; |
|
|
38 idxi=[0:k-1]; |
|
|
39 elem= [ ext+idxe ext+2*k+[-idxe idxe] ext+rem(4*k-idxe,4*k) ... |
|
|
40 ext+idxe ext+2*k+[-idxe idxe] ext+rem(4*k-idxe,4*k); |
|
|
41 int+idxi int+2*(k-1)+[-idxi idxi] ... |
|
|
42 int+rem(4*(k-1)-idxi, 4*(k-1)+(k==1) ) ... |
|
|
43 ext+4*k+1+idxi ext+6*k+[1-idxi 3+idxi] ext+8*k+3-idxi ]; |
|
|
44 for j=1:k |
|
|
45 r1= rem(j+k-1,3)+1; |
|
|
46 r2= rem(j+k,3)+1; |
|
|
47 r3= 6-r1-r2; |
|
|
48 elem([r1 r2 r3],j+k*(0:7) )= elem(:,j+k*(0:7)); |
|
|
49 end |
|
|
50 |
|
|
51 ELEM=[ ELEM elem(:,1:(8-4*(k==N))*k) ]; |
|
|
52 int=ext; |
|
|
53 end %for k=1:N |
|
|
54 |
|
|
55 MES= ext+N*4/elec*([0:elec-1]); |
|
|
56 |
|
|
57 j=pos_i(1); k=floor(j); |
|
|
58 MES=[MES( 1+rem( (0:elec-1)+k,elec) )+floor(MES(1:2)*[-1;1]*(j-k)); |
|
|
59 MES( 1+rem( (1:elec)+k,elec) )+floor(MES(1:2)*[-1;1]*(j-k)) ]; |
|
|
60 |
|
|
61 cour= [ MES(1,:)' MES(1,1+rem(elec+(0:elec-1)+pos_i*[-1;1],elec) )'; |
|
|
62 ones(elec,1)*[-1 1] ]; |
|
|
63 |
|
|
64 QQCirC= [zeros(ext-1,1);cos(2*pi*(0:4/N:elec-.01)'/elec)]; |
|
|
65 QQCirC= 2*QQCirC/sum(abs(QQCirC)); |
|
|
66 if exist('niveaux') |
|
|
67 QQCirC= QQCirC*([-1 1]*niveaux([1 length(niveaux)]))/length(niveaux); |
|
|
68 end |
|
|
69 volt= [1; zeros(elec,1)]; |
|
|
70 |
|
|
71 ELS=rem(rem(0:elec^2-1,elec)-floor((0:elec^2-1)/elec)+elec,elec)'; |
|
|
72 ELS=~any(rem( elec+[-1 0 [-1 0]+pos_i*[-1;1] ] ,elec)' ... |
|
|
73 *ones(1,elec^2)==ones(4,1)*ELS')'; |
|
|
74 |
|
|
75 |
|
|
76 d= size(ELEM,1); %dimentions+1 |
|
|
77 n= size(NODE,2); %NODEs |
|
|
78 e= size(ELEM,2); %ELEMents |
|
|
79 |
|
|
80 if dotheplot |
|
|
81 % octave commands, we use eval so it doesn't fail in Matlab |
|
|
82 eval('gset nokey;','1'); |
|
|
83 fleche= [1.02 0;1.06 .05;1.06 .02;1.1 .02; ... |
|
|
84 1.1 -.02; 1.06 -.02;1.06 -.05;1.02 0]; |
|
|
85 jjj= .95; |
|
|
86 if ~isempty(MES) |
|
|
87 xy=NODE(:,MES(1,:)); |
|
|
88 end |
|
|
89 xxx=zeros(3,e); xxx(:)=NODE(1,ELEM(:)); |
|
|
90 xxx= jjj*xxx+ (1-jjj)*ones(3,1)*mean(xxx); |
|
|
91 yyy=zeros(3,e); yyy(:)=NODE(2,ELEM(:)); |
|
|
92 yyy= jjj*yyy+ (1-jjj)*ones(3,1)*mean(yyy); |
|
|
93 plot([xxx;xxx(1,:)],[yyy;yyy(1,:)],'b', ... |
|
|
94 fleche*xy,fleche*[0 1;-1 0]*xy,'r'); |
|
|
95 |
|
|
96 axis([ [-1.1 1.1]*max(NODE(1,:)) [-1.1 1.1]*max(NODE(2,:)) ]) |
|
|
97 end % plot mesh |
|
|
98 if exist('niveaux') |
|
|
99 node=NODE; |
|
|
100 elem= [ELEM([1 1 2 3 ],:) ... |
|
|
101 ELEM([2 1 2 3],:) ELEM([3 2 1 3],:)]; |
|
|
102 NODE= [node; niveaux(1)*ones(1,n) ]; |
|
|
103 ELEM= []; |
|
|
104 |
|
|
105 for k=2:length(niveaux); |
|
|
106 NODE=[NODE ,[node; niveaux(k)*ones(1,n)] ]; |
|
|
107 ELEM= [ELEM,(elem + ... |
|
|
108 [[(k-1)*n*ones(1,e);(k-2)*n*ones(3,e)] ... |
|
|
109 [(k-1)*n*ones(2,e);(k-2)*n*ones(2,e)] ... |
|
|
110 [(k-1)*n*ones(3,e);(k-2)*n*ones(1,e)]] ) ]; |
|
|
111 end %for k |
|
|
112 |
|
|
113 MES= MES + floor(length(niveaux)/2)*n; |
|
|
114 QQCirC= QQCirC*ones(1, length(niveaux)); QQCirC= QQCirC(:); |
|
|
115 cour(1:elec,:)= cour(1:elec,:)+ floor(length(niveaux)/2)*n; |
|
|
116 |
|
|
117 end %if exist('niveaux') |
|
|
118 |
|
|
119 d= size(ELEM,1); %dimentions+1 |
|
|
120 n= size(NODE,2); %NODEs |
|
|
121 e= size(ELEM,2); %ELEMents |
|
|
122 p= size(volt,1)-1; |
|
|
123 |
|
|
124 CC= sparse((1:d*e),ELEM(:),ones(d*e,1), d*e, n); |
|
|
125 |
|
|
126 sa= zeros(d*e,d); |
|
|
127 for j=1:e |
|
|
128 a= inv([ ones(d,1) NODE( :, ELEM(:,j) )' ]); |
|
|
129 sa(d*(j-1)+1:d*j,:)= a(2:d,:)'*a(2:d,:)/(d-1)/(d-2)/abs(det(a)); |
|
|
130 end %for j=1:ELEMs |
|
|
131 ridx= ones(d,1)*(1:e)*d; |
|
|
132 ridx= ridx(:)*ones(1,d) + ones(d*e,1)*[-d+1:0]; |
|
|
133 cidx= 1:d*e; |
|
|
134 SS= sparse( cidx'*ones(1,d), ridx, sa, d*e, d*e); |
|
|
135 |
|
|
136 QQ=sparse(cour(1:p,:),(1:p)'*ones(1,size(cour,2)), ... |
|
|
137 cour(p+1:2*p,:),n,p ); |
|
|
138 |
|
|
139 |
|
|
140 if 0 %OCTAVE |
|
|
141 # this code uses the first syntax for |
|
|
142 # octave sparse |
|
|
143 CCt= spfun(CC,'trans'); |
|
|
144 CCtSS= spfun(CCt,'mul',SS); |
|
|
145 ZZ= spfun(CCtSS,'mul',CC); |
|
|
146 ZZs= spfun(ZZ,'extract',2,n,2,n); |
|
|
147 QF= full(spfun(QQ,'extract',2,n,1,p)); |
|
|
148 fprintf('Solving Finite Element sparse eqn, n=%d nnz=%d density=%f\n', ... |
|
|
149 n, nnz(ZZ), nnz(ZZ)/n^2 ); |
|
|
150 t0= clock; |
|
|
151 for i=1:TESTTIMES |
|
|
152 VV= spfun(ZZs,'solve',QF,2); |
|
|
153 end |
|
|
154 tasktime= etime(clock,t0); |
|
|
155 else |
|
|
156 ZZ= CC'*SS*CC; |
|
|
157 ZZs= ZZ(2:n,2:n); |
|
|
158 QF= full(QQ(2:n,:)); |
|
|
159 % QF= (QQ(2:n,:)); |
|
|
160 fprintf('Solving Finite Element sparse eqn, n=%d nnz=%d density=%f\n', ... |
|
|
161 n, nnz(ZZ), nnz(ZZ)/n^2 ); |
|
|
162 t0= clock; |
|
|
163 for i=1:TESTTIMES |
|
|
164 VV= ZZs\QF; |
|
|
165 end |
|
|
166 tasktime= etime(clock,t0); |
|
|
167 end |
|
|
168 |
|
|
169 |
|
|
170 checkval= full( VV(MES(1,:),1)-VV(MES(2,:),1) ); |
|
|
171 goldvalue= [ -0.918243; 0.555779; 0.148452; 0.078826; |
|
|
172 0.052132; 0.040116; 0.034141; 0.031872; |
|
|
173 0.031279; 0.033594; 0.039429; 0.052189; |
|
|
174 0.078377; 0.149426; 0.596195; -1.003565]; |
|
|
175 |
|
|
176 rdif= sqrt(mean((checkval-goldvalue).^2)) / sqrt(mean(checkval.^2)); |
|
|
177 if rdif > 1e-5 |
|
|
178 fprintf('Sparse FEM solution fails'); |
|
|
179 end |
|
|
180 fprintf('Time per iteration= %f s\n', tasktime/ TESTTIMES); |
|
|
181 fprintf('Your machine is %f faster than a 486dx100!\n', ... |
|
|
182 4.45*TESTTIMES/tasktime ); |
|
|
183 |
|
|
184 #ZZt=ZZ'; tt=time; for i = 1:100; x= ZZ*ZZt; end ; (time-tt)/100 |
|
|
185 % ans = 0.039929 % 0.52245 |
|
|
186 |
|
|
187 % |
|
|
188 % Results: |
|
|
189 % 11 Nov 00: Your machine is 16.501822 faster than a 486dx100! |
|
|
190 % |
|
|
191 % $Log$ |
|
|
192 % Revision 1.1 2001/10/10 19:54:49 pkienzle |
|
|
193 % Initial revision |
|
|
194 % |
|
|
195 % Revision 1.4 2001/04/04 02:13:46 aadler |
|
|
196 % complete complex_sparse, templates, fix memory leaks |
|
|
197 % |
|
|
198 % Revision 1.3 2001/03/30 04:36:30 aadler |
|
|
199 % added multiply, solve, and sparse creation |
|
|
200 % |
|
|
201 % Revision 1.2 2000/12/18 03:31:16 aadler |
|
|
202 % Split code to multiple files |
|
|
203 % added sparse inverse |
|
|
204 % |
|
|
205 % Revision 1.1 2000/11/11 02:47:11 aadler |
|
|
206 % DLD functions for sparse support in octave |
|
|
207 % |
|
|
208 % |
