Mercurial > hg > octave-lyh
changeset 6003:7f8e4f5e2eb7
[project @ 2006-09-26 21:09:31 by jwe]
author | jwe |
---|---|
date | Tue, 26 Sep 2006 21:09:31 +0000 |
parents | c5d0490aad55 |
children | 01556febbaaf |
files | doc/interpreter/sparseimages.m |
diffstat | 1 files changed, 221 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/doc/interpreter/sparseimages.m @@ -0,0 +1,221 @@ +function sparseimages(nm,typ) + if (strcmp(typ,"txt")) + txtimages(nm,15,typ); + else + if (strcmp (nm, "gplot")) + gplotimages("gplot",typ); + elseif (strcmp (nm, "grid")) + femimages("grid",typ); + else + otherimages(nm,200,typ); + endif + endif + ## Kluge to give gnuplot enough time to process last figure before we + ## exit. Otherwise, Octave will delete the temporary data files when + ## it exits and gnuplot will fail... + sleep (1); +endfunction + +## Use this function before plotting commands and after every call to +## print since print() resets output to stdout (unfortunately, gnpulot +## can't pop output as it can the terminal type). +function bury_output () + automatic_replot(0); + __gnuplot_set__ term dumb + [status, dummy] = fileattrib("/dev/null"); + if (status) + __gnuplot_raw__ ("set output \"/dev/null\"\n"); + endif +endfunction + +function gplotimages(nm,typ) + bury_output (); + A = sparse([2,6,1,3,2,4,3,5,4,6,1,5], + [1,1,2,2,3,3,4,4,5,5,6,6],1,6,6); + xy = [0,4,8,6,4,2;5,0,5,7,5,7]'; + gplot(A,xy) + print(strcat(nm,".",typ),strcat("-d",typ)) + bury_output (); +endfunction + +function txtimages(nm,n,typ) + a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... + sparse(ceil([1:n]/2),1:n,1,n,n); + if (strcmp (nm, "gplot") || strcmp (nm, "grid")) + fid = fopen (sprintf ("%s.txt", nm), "wt"); + fputs (fid, "+---------------------------------+\n"); + fputs (fid, "| Image unavailable in text mode. |\n"); + fputs (fid, "+---------------------------------+\n"); + fclose (fid); + elseif (strcmp (nm, "spmatrix")) + printsparse(a,strcat("spmatrix.",typ)); + else + if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && + !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD"))) + if (strcmp (nm, "spchol")) + r1 = chol(a); + printsparse(r1,strcat("spchol.",typ)); + elseif (strcmp (nm, "spcholperm")) + [r2,p2,q2]=chol(a); + printsparse(r2,strcat("spcholperm.",typ)); + endif + ## printf("Text NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); + endif + endif +endfunction + +function otherimages(nm,n,typ) + bury_output (); + a = 10*speye(n) + sparse(1:n,ceil([1:n]/2),1,n,n) + ... + sparse(ceil([1:n]/2),1:n,1,n,n); + if (strcmp (nm, "spmatrix")) + spy(a); + axis("ij") + print(strcat("spmatrix.",typ),strcat("-d",typ)) + bury_output (); + else + if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && + !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD"))) + if (strcmp (nm, "spchol")) + r1 = chol(a); + spy(r1); + axis("ij") + print(strcat("spchol.",typ),strcat("-d",typ)) + bury_output (); + elseif (strcmp (nm, "spcholperm")) + [r2,p2,q2]=chol(a); + spy(r2); + axis("ij") + print(strcat("spcholperm.",typ),strcat("-d",typ)) + bury_output (); + endif + ## printf("Image NNZ: Matrix %d, Chol %d, PermChol %d\n",nnz(a),nnz(r1),nnz(r2)); + endif + endif +endfunction + +function printsparse(a,nm) + fid = fopen (nm,"wt"); + for i = 1:size(a,1) + if (rem(i,5) == 0) + fprintf (fid," %2d - ", i); + else + fprintf (fid," | "); + endif + for j = 1:size(a,2) + if (a(i,j) == 0) + fprintf(fid," ") + else + fprintf(fid," *") + endif + endfor + fprintf(fid,"\n") + endfor + fprintf(fid," |-"); + for j=1:size(a,2) + if (rem(j,5)==0) + fprintf(fid,"-|"); + else + fprintf(fid,"--"); + endif + endfor + fprintf(fid,"\n") + fprintf(fid," "); + for j=1:size(a,2) + if (rem(j,5)==0) + fprintf(fid,"%2d",j); + else + fprintf(fid," "); + endif + endfor + fclose(fid); +endfunction + +function femimages (nm,typ) + bury_output (); + if (!isempty(findstr(octave_config_info ("DEFS"),"HAVE_COLAMD")) && + !isempty(findstr(octave_config_info ("DEFS"),"HAVE_CHOLMOD")) && + !isempty(findstr(octave_config_info ("DEFS"),"HAVE_UMFPACK"))) + ## build a rectangle + node_y = [1;1.2;1.5;1.8;2]*ones(1,11); + 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]; + nodes = [node_x(:), node_y(:)]; + + [h,w] = size(node_x); + elems = []; + for idx = 1:w-1 + widx = (idx-1)*h; + elems = [elems; widx+[(1:h-1);(2:h);h+(1:h-1)]']; + elems = [elems; widx+[(2:h);h+(2:h);h+(1:h-1)]']; + endfor + + E = size(elems,1); #No. of elements + N = size(nodes,1); #No. of elements + D = size(elems,2); #dimentions+1 + + ## Plot FEM Geometry + elemx = elems(:,[1,2,3,1])'; + xelems = reshape( nodes(elemx, 1), 4, E); + yelems = reshape( nodes(elemx, 2), 4, E); + + ## Set element conductivity + conductivity = [1*ones(1,16),2*ones(1,48),1*ones(1,16)]; + + ## Dirichlet boundary conditions + D_nodes = [1:5, 51:55]; + D_value = [10*ones(1,5), 20*ones(1,5)]; + + ## Neumann boundary conditions + ## Note that N_value must be normalized by the boundary + ## length and element conductivity + N_nodes = []; + N_value = []; + + ## Calculate connectivity matrix + C = sparse((1:D*E), reshape(elems',D*E,1),1, D*E, N); + + ## Calculate stiffness matrix + Siidx = floor([0:D*E-1]'/D)*D*ones(1,D) + ones(D*E,1)*(1:D) ; + Sjidx = [1:D*E]'*ones(1,D); + Sdata = zeros(D*E,D); + dfact = prod(2:(D-1)); + for j = 1:E + a = inv([ ones(D,1), nodes( elems(j,:), : ) ]); + const = conductivity(j)*2/dfact/abs(det(a)); + Sdata(D*(j-1)+(1:D),:)= const * a(2:D,:)'*a(2:D,:); + endfor + + ## Element-wise system matrix + SE = sparse(Siidx,Sjidx,Sdata); + ## Global system matrix + S = C'* SE *C; + + ## Set Dirichlet boundary + V = zeros(N,1); + V(D_nodes) = D_value; + idx = 1:N; + idx(D_nodes) = []; + + ## Set Neumann boundary + Q = zeros(N,1); + Q(N_nodes) = N_value; # FIXME + + V(idx) = S(idx,idx)\( Q(idx) - S(idx,D_nodes)*V(D_nodes) ); + + velems = reshape( V(elemx), 4, E); + + sz = size(xelems,2); + ## FIXME How can I do this without a gnuplot specific commands? plot3 anyone? + unwind_protect + __gnuplot_set__ parametric; + __gnuplot_raw__ ("set nohidden3d;\n"); + tmp = [([xelems; NaN*ones(1,sz)])(:), ([yelems; NaN*ones(1,sz)])(:), ([velems; NaN*ones(1,sz)])(:)]; + __gnuplot_splot__(tmp); + __gnuplot_raw__ ("set view 80,10;\n") + print(strcat(nm,".",typ),strcat("-d",typ)) + bury_output (); + unwind_protect_cleanup + __gnuplot_set__ noparametric; + end_unwind_protect + endif +endfunction