From: <cd...@us...> - 2008-03-14 13:16:49
|
Revision: 4748 http://octave.svn.sourceforge.net/octave/?rev=4748&view=rev Author: cdf Date: 2008-03-14 06:16:53 -0700 (Fri, 14 Mar 2008) Log Message: ----------- equalize mesh Modified Paths: -------------- trunk/octave-forge/extra/msh/INDEX Added Paths: ----------- trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m Modified: trunk/octave-forge/extra/msh/INDEX =================================================================== --- trunk/octave-forge/extra/msh/INDEX 2008-03-14 11:02:08 UTC (rev 4747) +++ trunk/octave-forge/extra/msh/INDEX 2008-03-14 13:16:53 UTC (rev 4748) @@ -8,3 +8,6 @@ MSH2Mgeomprop MSH2Mtopprop MSH2Mnodesonsides +Mesh adaptation + MSH2Mequalizemesh + MSH2Mdisplacementsmoothing Added: trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m =================================================================== --- trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m (rev 0) +++ trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m 2008-03-14 13:16:53 UTC (rev 4748) @@ -0,0 +1,195 @@ +function [Ax,Ay,Bx,By] = MSH2Mequalizemesh(msh) + + + ## -*- texinfo -*- + ## @deftypefn {Function File} {[@var{Ax}, @var{Ay}, @var{bx}, @ + ## @var{by}]} = MSH2Mequalizemesh(@var{msh}) + ## + ## To equalize the size of triangle edges, set a spring of equilibrium + ## length @var{factor}*@var{area} along each edge of the mesh and solve for + ## equilibrium. + ## + ## This function builds matrices containing the resulting + ## (linearized) equation for x and y coordinates of each mesh node. + ## Boundary conditions enforcing the displacement (Dirichlet type + ## problem) or the force (Neumann type) at the boundary must be added + ## separately to make the system solvable. + ## + ## May be useful when distorting a mesh, e.g.: + ## + ## @example + ## + ## msh = MSH2Mstructmesh(linspace(0,1,10),@ + ## linspace(0,1,10),@ + ## 1,1:4,"random"); + ## dnodes = MSH2Mnodesonsides(msh,1:4); + ## varnodes = setdiff([1:columns(msh.p)],dnodes); + ## xd = msh.p(1,dnodes)'; + ## yd = msh.p(2,dnodes)'; + ## dx = dy = zeros(columns(msh.p),1); + ## dytot = dxtot = -.4*sin(xd.*yd*pi/2); + ## Nsteps = 10; + ## for ii=1:Nsteps + ## dx(dnodes) = dxtot; + ## dy(dnodes) = dytot; + ## [Ax,Ay] = MSH2Mdisplacementsmoothing(msh,1); + ## dx(varnodes) = Ax(varnodes,varnodes) \ ... + ## (-Ax(varnodes,dnodes)*dx(dnodes)); + ## dy(varnodes) = Ay(varnodes,varnodes) \ ... + ## (-Ay(varnodes,dnodes)*dy(dnodes)); + ## msh.p(1,:) += dx'/Nsteps; + ## msh.p(2,:) += dy'/Nsteps; + ## dx(dnodes) = 0; + ## dy(dnodes) = 0; + ## [AX, AY, BX, BY] = MSH2Mequalizemesh(msh, .1); + ## dx(varnodes) = Ax(varnodes,varnodes) \ ... + ## (BX(varnodes)-Ax(varnodes,dnodes)*dx(dnodes)); + ## dy(varnodes) = Ay(varnodes,varnodes) \ ... + ## (BY(varnodes)-Ay(varnodes,dnodes)*dy(dnodes)); + ## msh.p(1,:) += dx'; + ## msh.p(2,:) += dy'; + ## triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); + ## pause(.01) + ## endfor + ## @end example + ## @seealso{MSH2Mdisplacementsmoothing} + ## + ## @end deftypefn + + ## This file is part of + ## + ## MSH - Meshing Software Package for Octave + ## ------------------------------------------------------------------- + ## Copyright (C) 2007 Carlo de Falco and Culpo Massimiliano + ## + ## MSH is free software; you can redistribute it and/or modify + ## it under the terms of the GNU General Public License as published by + ## the Free Software Foundation; either version 2 of the License, or + ## (at your option) any later version. + ## + ## MSH is distributed in the hope that it will be useful, + ## but WITHOUT ANY WARRANTY; without even the implied warranty of + ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + ## GNU General Public License for more details. + ## + ## You should have received a copy of the GNU General Public License + ## along with MSH; if not, write to the Free Software + ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 + ## USA + ## + ## + ## AUTHORS: + ## Carlo de Falco + ## Dublin City University + ## School of Mathemetical Sciences + ## Ireland + ## + ## Culpo Massimiliano + ## Bergische Universitaet Wuppertal + ## Fachbereich C - Mathematik und Naturwissenschaften + ## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 + ## D-42119 Wuppertal, Germany + + nel= columns(msh.t); + + x = msh.p(1,:); + y = msh.p(2,:); + +# dx = (x(msh.t([1 2 3],:))-x(msh.t([2 3 1],:))); +# dy = (y(msh.t([1 2 3],:))-y(msh.t([2 3 1],:))); +# dx2 = dx.^2; +# dy2 = dy.^2; + +# l2 = dx2 + dy2; +# l = sqrt(l2); +# area = MSH2Mgeomprop(msh,'area'); +# l0 = factor * sqrt(area) ([1 1 1],:); + + + Ax = spalloc(length(x),length(x),1); + Ay = spalloc(length(x),length(x),1); + + ax = zeros(3,3,nel); + ay = zeros(3,3,nel); + bx = spalloc(3,nel,0); + by = spalloc(3,nel,0); + + for inode=1:3 + giinode(inode,:)=msh.t(inode,:); + for jnode=1:3 + ginode(inode,jnode,:)=msh.t(inode,:); + gjnode(inode,jnode,:)=msh.t(jnode,:); + end + end + + for ii=1:3 + +# bx (ii,:) = ((l0(ii,:)-l(ii,:))./l(ii,:)) .* dx(ii,:); +# by (ii,:) = ((l0(ii,:)-l(ii,:))./l(ii,:)) .* dy(ii,:); + + for jj=ii+1:3 + + ax(ii,jj,:) = ax(jj,ii,:) = -ones(1,1,nel); + ay(ii,jj,:) = ay(jj,ii,:) = -ones(1,1,nel); + + ax(ii,ii,:) -= ax(ii,jj,:); + ax(jj,jj,:) -= ax(ii,jj,:); + ay(ii,ii,:) -= ay(ii,jj,:); + ay(jj,jj,:) -= ay(ii,jj,:); + + endfor + endfor + + Ax = sparse(ginode(:),gjnode(:),ax(:)); + Ay = sparse(ginode(:),gjnode(:),ay(:)); + Bx = sparse(giinode(:),1,bx(:)); + By = sparse(giinode(:),1,by(:)); + + +%!demo +%! ### equalize a structured mesh without moving boundary nodes +%! msh = MSH2Mstructmesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"random"); +%! dnodes = MSH2Mnodesonsides(msh,1:4); +%! varnodes = setdiff([1:columns(msh.p)],dnodes); +%! x = msh.p(1,:)'; +%! y = msh.p(2,:)'; +%! [AX, AY, BX, BY] = MSH2Mequalizemesh(msh); +%! x(varnodes) = AX(varnodes,varnodes) \ (BX(varnodes)-AX(varnodes,dnodes)*x(dnodes)); +%! y(varnodes) = AY(varnodes,varnodes) \ (BY(varnodes)-AY(varnodes,dnodes)*y(dnodes)); +%! msh.p(1,:) = x'; +%! msh.p(2,:) = y'; +%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); +%! pause(.01) + +%!demo +%! ### distort a mesh on a square equalizing at each step +%! msh = MSH2Mstructmesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"random"); +%! dnodes = MSH2Mnodesonsides(msh,1:4); +%! varnodes = setdiff([1:columns(msh.p)],dnodes); +%! x = msh.p(1,:)'; +%! y = msh.p(2,:)'; +%! dx = dy = zeros(columns(msh.p),1); +%! dytot = dxtot = -.7*sin(x(dnodes).*y(dnodes)*pi/2); +%! Nsteps = 30; +%! for ii=1:Nsteps +%! dx(dnodes) = dxtot; +%! dy(dnodes) = dytot; +%! [Ax,Ay] = MSH2Mdisplacementsmoothing(msh,1); +%! dx(varnodes) = Ax(varnodes,varnodes) \ ... +%! (-Ax(varnodes,dnodes)*dx(dnodes)); +%! dy(varnodes) = Ay(varnodes,varnodes) \ ... +%! (-Ay(varnodes,dnodes)*dy(dnodes)); +%! msh.p(1,:) += dx'/Nsteps; +%! msh.p(2,:) += dy'/Nsteps; +%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)','r'); +%! pause(.5) +%! x = msh.p(1,:)'; +%! y = msh.p(2,:)'; +%! [AX, AY, BX, BY] = MSH2Mequalizemesh(msh); +%! x(varnodes) = AX(varnodes,varnodes) \ (BX(varnodes)-AX(varnodes,dnodes)*x(dnodes)); +%! y(varnodes) = AY(varnodes,varnodes) \ (BY(varnodes)-AY(varnodes,dnodes)*y(dnodes)); +%! msh.p(1,:) = x'; +%! msh.p(2,:) = y'; +%! hold on;triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)');hold off +%! pause(.5) +%! endfor \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |