|
6126
|
1 function sparseimages (nm, typ) |
|
6040
|
2 if (! isempty (findstr (octave_config_info ("DEFS"), "HAVE_COLAMD")) |
|
|
3 && ! isempty (findstr (octave_config_info ("DEFS"), "HAVE_CHOLMOD")) |
|
|
4 && ! isempty (findstr (octave_config_info ("DEFS"), "HAVE_UMFPACK"))) |
|
6126
|
5 if (strcmp(typ,"txt")) |
|
|
6 txtimages (nm, 15, typ); |
|
6003
|
7 else |
|
6126
|
8 if (strcmp (nm, "gplot")) |
|
|
9 gplotimages ("gplot", typ); |
|
|
10 elseif (strcmp (nm, "grid")) |
|
|
11 femimages ("grid", typ); |
|
|
12 else |
|
|
13 otherimages (nm, 200, typ); |
|
|
14 endif |
|
6003
|
15 endif |
|
6040
|
16 else ## There is no sparse matrix implementation available because |
|
|
17 ## of missing libraries, plot sombreros instead |
|
|
18 sombreroimage (nm, typ); |
|
|
19 endif |
|
6003
|
20 endfunction |
|
|
21 |
|
|
22 ## Use this function before plotting commands and after every call to |
|
|
23 ## print since print() resets output to stdout (unfortunately, gnpulot |
|
|
24 ## can't pop output as it can the terminal type). |
|
|
25 function bury_output () |
|
6040
|
26 automatic_replot (0); |
|
6003
|
27 __gnuplot_set__ term dumb |
|
6040
|
28 [status, dummy] = fileattrib ("/dev/null"); |
|
6003
|
29 if (status) |
|
|
30 __gnuplot_raw__ ("set output \"/dev/null\"\n"); |
|
|
31 endif |
|
|
32 endfunction |
|
|
33 |
|
6040
|
34 function gplotimages (nm, typ) |
|
6003
|
35 bury_output (); |
|
6040
|
36 A = sparse ([2,6,1,3,2,4,3,5,4,6,1,5], |
|
|
37 [1,1,2,2,3,3,4,4,5,5,6,6], 1, 6, 6); |
|
6003
|
38 xy = [0,4,8,6,4,2;5,0,5,7,5,7]'; |
|
6040
|
39 gplot (A, xy) |
|
|
40 print (strcat (nm, ".", typ), strcat ("-d", typ)) |
|
6003
|
41 bury_output (); |
|
|
42 endfunction |
|
|
43 |
|
|
44 function txtimages(nm,n,typ) |
|
|
45 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... |
|
|
46 sparse(ceil([1:n]/2),1:n,1,n,n); |
|
|
47 if (strcmp (nm, "gplot") || strcmp (nm, "grid")) |
|
|
48 fid = fopen (sprintf ("%s.txt", nm), "wt"); |
|
|
49 fputs (fid, "+---------------------------------+\n"); |
|
|
50 fputs (fid, "| Image unavailable in text mode. |\n"); |
|
|
51 fputs (fid, "+---------------------------------+\n"); |
|
|
52 fclose (fid); |
|
|
53 elseif (strcmp (nm, "spmatrix")) |
|
|
54 printsparse(a,strcat("spmatrix.",typ)); |
|
|
55 else |
|
|
56 if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && |
|
|
57 !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD"))) |
|
|
58 if (strcmp (nm, "spchol")) |
|
|
59 r1 = chol(a); |
|
|
60 printsparse(r1,strcat("spchol.",typ)); |
|
|
61 elseif (strcmp (nm, "spcholperm")) |
|
|
62 [r2,p2,q2]=chol(a); |
|
|
63 printsparse(r2,strcat("spcholperm.",typ)); |
|
|
64 endif |
|
|
65 ## printf("Text NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); |
|
|
66 endif |
|
|
67 endif |
|
|
68 endfunction |
|
|
69 |
|
|
70 function otherimages(nm,n,typ) |
|
|
71 bury_output (); |
|
|
72 a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... |
|
|
73 sparse(ceil([1:n]/2),1:n,1,n,n); |
|
|
74 if (strcmp (nm, "spmatrix")) |
|
|
75 spy(a); |
|
|
76 axis("ij") |
|
|
77 print(strcat("spmatrix.",typ),strcat("-d",typ)) |
|
|
78 bury_output (); |
|
|
79 else |
|
|
80 if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && |
|
|
81 !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD"))) |
|
|
82 if (strcmp (nm, "spchol")) |
|
|
83 r1 = chol(a); |
|
|
84 spy(r1); |
|
|
85 axis("ij") |
|
|
86 print(strcat("spchol.",typ),strcat("-d",typ)) |
|
|
87 bury_output (); |
|
|
88 elseif (strcmp (nm, "spcholperm")) |
|
|
89 [r2,p2,q2]=chol(a); |
|
|
90 spy(r2); |
|
|
91 axis("ij") |
|
|
92 print(strcat("spcholperm.",typ),strcat("-d",typ)) |
|
|
93 bury_output (); |
|
|
94 endif |
|
|
95 ## printf("Image NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); |
|
|
96 endif |
|
|
97 endif |
|
|
98 endfunction |
|
|
99 |
|
|
100 function printsparse(a,nm) |
|
|
101 fid = fopen (nm,"wt"); |
|
|
102 for i = 1:size(a,1) |
|
|
103 if (rem(i,5) == 0) |
|
|
104 fprintf (fid," %2d - ", i); |
|
|
105 else |
|
|
106 fprintf (fid," | "); |
|
|
107 endif |
|
|
108 for j = 1:size(a,2) |
|
|
109 if (a(i,j) == 0) |
|
|
110 fprintf(fid," ") |
|
|
111 else |
|
|
112 fprintf(fid," *") |
|
|
113 endif |
|
|
114 endfor |
|
|
115 fprintf(fid,"\n") |
|
|
116 endfor |
|
|
117 fprintf(fid," |-"); |
|
|
118 for j=1:size(a,2) |
|
|
119 if (rem(j,5)==0) |
|
|
120 fprintf(fid,"-|"); |
|
|
121 else |
|
|
122 fprintf(fid,"--"); |
|
|
123 endif |
|
|
124 endfor |
|
|
125 fprintf(fid,"\n") |
|
|
126 fprintf(fid," "); |
|
|
127 for j=1:size(a,2) |
|
|
128 if (rem(j,5)==0) |
|
|
129 fprintf(fid,"%2d",j); |
|
|
130 else |
|
|
131 fprintf(fid," "); |
|
|
132 endif |
|
|
133 endfor |
|
|
134 fclose(fid); |
|
|
135 endfunction |
|
|
136 |
|
|
137 function femimages (nm,typ) |
|
|
138 bury_output (); |
|
|
139 if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && |
|
|
140 !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD")) && |
|
|
141 !isempty(findstr(octave_config_info ("DEFS"),"HAVE_UMFPACK"))) |
|
|
142 ## build a rectangle |
|
|
143 node_y = [1;1.2;1.5;1.8;2]*ones(1,11); |
|
|
144 node_x = ones(5,1)*[1,1.05,1.1,1.2,1.3,1.5,1.7,1.8,1.9,1.95,2]; |
|
|
145 nodes = [node_x(:), node_y(:)]; |
|
|
146 |
|
|
147 [h,w] = size(node_x); |
|
|
148 elems = []; |
|
|
149 for idx = 1:w-1 |
|
|
150 widx = (idx-1)*h; |
|
|
151 elems = [elems; widx+[(1:h-1);(2:h);h+(1:h-1)]']; |
|
|
152 elems = [elems; widx+[(2:h);h+(2:h);h+(1:h-1)]']; |
|
|
153 endfor |
|
|
154 |
|
|
155 E = size(elems,1); #No. of elements |
|
|
156 N = size(nodes,1); #No. of elements |
|
|
157 D = size(elems,2); #dimentions+1 |
|
|
158 |
|
|
159 ## Plot FEM Geometry |
|
|
160 elemx = elems(:,[1,2,3,1])'; |
|
|
161 xelems = reshape( nodes(elemx, 1), 4, E); |
|
|
162 yelems = reshape( nodes(elemx, 2), 4, E); |
|
|
163 |
|
|
164 ## Set element conductivity |
|
|
165 conductivity = [1*ones(1,16),2*ones(1,48),1*ones(1,16)]; |
|
|
166 |
|
|
167 ## Dirichlet boundary conditions |
|
|
168 D_nodes = [1:5, 51:55]; |
|
|
169 D_value = [10*ones(1,5), 20*ones(1,5)]; |
|
|
170 |
|
|
171 ## Neumann boundary conditions |
|
|
172 ## Note that N_value must be normalized by the boundary |
|
|
173 ## length and element conductivity |
|
|
174 N_nodes = []; |
|
|
175 N_value = []; |
|
|
176 |
|
|
177 ## Calculate connectivity matrix |
|
|
178 C = sparse((1:D*E), reshape(elems',D*E,1),1, D*E, N); |
|
|
179 |
|
|
180 ## Calculate stiffness matrix |
|
|
181 Siidx = floor([0:D*E-1]'/D)*D*ones(1,D) + ones(D*E,1)*(1:D) ; |
|
|
182 Sjidx = [1:D*E]'*ones(1,D); |
|
|
183 Sdata = zeros(D*E,D); |
|
|
184 dfact = prod(2:(D-1)); |
|
|
185 for j = 1:E |
|
|
186 a = inv([ ones(D,1), nodes( elems(j,:), : ) ]); |
|
|
187 const = conductivity(j)*2/dfact/abs(det(a)); |
|
|
188 Sdata(D*(j-1)+(1:D),:)= const * a(2:D,:)'*a(2:D,:); |
|
|
189 endfor |
|
|
190 |
|
|
191 ## Element-wise system matrix |
|
|
192 SE = sparse(Siidx,Sjidx,Sdata); |
|
|
193 ## Global system matrix |
|
|
194 S = C'* SE *C; |
|
|
195 |
|
|
196 ## Set Dirichlet boundary |
|
|
197 V = zeros(N,1); |
|
|
198 V(D_nodes) = D_value; |
|
|
199 idx = 1:N; |
|
|
200 idx(D_nodes) = []; |
|
|
201 |
|
|
202 ## Set Neumann boundary |
|
|
203 Q = zeros(N,1); |
|
|
204 Q(N_nodes) = N_value; # FIXME |
|
|
205 |
|
|
206 V(idx) = S(idx,idx)\( Q(idx) - S(idx,D_nodes)*V(D_nodes) ); |
|
|
207 |
|
|
208 velems = reshape( V(elemx), 4, E); |
|
|
209 |
|
|
210 sz = size(xelems,2); |
|
|
211 ## FIXME How can I do this without a gnuplot specific commands? plot3 anyone? |
|
|
212 unwind_protect |
|
|
213 __gnuplot_set__ parametric; |
|
|
214 __gnuplot_raw__ ("set nohidden3d;\n"); |
|
|
215 tmp = [([xelems; NaN*ones(1,sz)])(:), ([yelems; NaN*ones(1,sz)])(:), ([velems; NaN*ones(1,sz)])(:)]; |
|
|
216 __gnuplot_splot__(tmp); |
|
|
217 __gnuplot_raw__ ("set view 80,10;\n") |
|
|
218 print(strcat(nm,".",typ),strcat("-d",typ)) |
|
|
219 bury_output (); |
|
|
220 unwind_protect_cleanup |
|
|
221 __gnuplot_set__ noparametric; |
|
|
222 end_unwind_protect |
|
|
223 endif |
|
|
224 endfunction |
|
6040
|
225 |
|
|
226 ## There is no sparse matrix implementation available because of missing |
|
|
227 ## libraries, plot sombreros instead. Also plot a nice title that we are |
|
|
228 ## sorry about that. |
|
|
229 function sombreroimage (nm, typ) |
|
|
230 if (strcmp (typ, "txt")) |
|
|
231 fid = fopen (sprintf ("%s.txt", nm), "wt"); |
|
|
232 fputs (fid, "+---------------------------------------+\n"); |
|
|
233 fputs (fid, "| Image unavailable because of a |\n"); |
|
|
234 fputs (fid, "| missing sparse matrix implementation. |\n"); |
|
|
235 fputs (fid, "+---------------------------------------+\n"); |
|
|
236 fclose (fid); |
|
|
237 return; |
|
|
238 else ## if (!strcmp (typ, "txt")) |
|
|
239 |
|
|
240 bury_output (); |
|
|
241 |
|
|
242 x = y = linspace (-8, 8, 41)'; |
|
|
243 [xx, yy] = meshgrid (x, y); |
|
|
244 r = sqrt (xx .^ 2 + yy .^ 2) + eps; |
|
|
245 z = sin (r) ./ r; |
|
|
246 xlen = length (x); |
|
|
247 ylen = length (y); |
|
|
248 len = 3 * xlen; |
|
|
249 zz = zeros (ylen, len); |
|
|
250 k = 1; |
|
|
251 for i = 1:3:len |
|
|
252 zz(:,i) = x(k) * ones (ylen, 1); |
|
|
253 zz(:,i+1) = y; |
|
|
254 zz(:,i+2) = z(:,k); |
|
|
255 k++; |
|
|
256 endfor |
|
|
257 unwind_protect |
|
|
258 __gnuplot_raw__ ("set hidden3d;\n"); |
|
|
259 __gnuplot_raw__ ("set data style lines;\n"); |
|
|
260 __gnuplot_raw__ ("set surface;\n"); |
|
|
261 __gnuplot_raw__ ("set nocontour;\n"); |
|
|
262 __gnuplot_raw__ ("set nologscale;\n"); |
|
|
263 __gnuplot_set__ parametric; |
|
|
264 __gnuplot_raw__ ("set view 60, 30, 1, 1;\n"); |
|
|
265 __gnuplot_raw__ ("set nokey;\n"); |
|
|
266 __gnuplot_raw__ ("set nocolorbox;\n"); |
|
|
267 msg = strcat (""); |
|
|
268 __gnuplot_raw__ ("set title \"Sorry, graphics not available because octave was\\ncompiled without the sparse matrix implementation.\";\n"); |
|
|
269 __plt3__ (zz, "", ""); |
|
|
270 unwind_protect_cleanup |
|
|
271 __gnuplot_set__ noparametric; |
|
|
272 print (strcat (nm, ".", typ), strcat ("-d", typ)); |
|
|
273 bury_output (); |
|
|
274 end_unwind_protect |
|
|
275 endif |
|
|
276 endfunction |
