Mercurial > octave-nkf
view doc/interpreter/images/sparseimages.m @ 5926:ddfe04062467
[project @ 2006-08-14 19:16:07 by jwe]
author | jwe |
---|---|
date | Mon, 14 Aug 2006 19:16:07 +0000 |
parents | b2a802aa0cda |
children |
line wrap: on
line source
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