1## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo 2## 3## This file is part of: 4## MSH - Meshing Software Package for Octave 5## 6## MSH is free software; you can redistribute it and/or modify 7## it under the terms of the GNU General Public License as published by 8## the Free Software Foundation; either version 2 of the License, or 9## (at your option) any later version. 10## 11## MSH is distributed in the hope that it will be useful, 12## but WITHOUT ANY WARRANTY; without even the implied warranty of 13## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14## GNU General Public License for more details. 15## 16## You should have received a copy of the GNU General Public License 17## along with MSH; If not, see <http://www.gnu.org/licenses/>. 18## 19## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> 20## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> 21 22## -*- texinfo -*- 23## @deftypefn {Function File} {[@var{Ax},@var{Ay}]} = @ 24## msh2m_displacement_smoothing(@var{msh},@var{k}) 25## 26## Displace the boundary of a 2D mesh setting a spring with force/length 27## constant @var{k} along each edge and enforcing equilibrium. 28## 29## This function builds matrices containing the resulting (linearized) 30## equation for x and y coordinates of each mesh node. Boundary 31## conditions enforcing the displacement (Dirichlet type problem) or the 32## force (Neumann type) at the boundary must be added to make the system 33## solvable, e.g.: 34## 35## @example 36## msh = msh2m_structured_mesh(linspace(0,1,10),@ 37## linspace(0,1,10),@ 38## 1,1:4,"left"); 39## 40## dnodes = msh2m_nodes_on_sides(msh,1:4); 41## varnodes = setdiff([1:columns(msh.p)],dnodes); 42## xd = msh.p(1,dnodes)'; 43## yd = msh.p(2,dnodes)'; 44## dx = dy = zeros(columns(msh.p),1); 45## dxtot = dytot = -.5*sin(xd.*yd*pi/2); 46## Nsteps = 10; 47## 48## for ii = 1:Nsteps 49## dx(dnodes) = dxtot; 50## dy(dnodes) = dytot; 51## [Ax,Ay] = msh2m_displacement_smoothing(msh,1); 52## dx(varnodes) = Ax(varnodes,varnodes) \ ... 53## (-Ax(varnodes,dnodes)*dx(dnodes)); 54## dy(varnodes) = Ay(varnodes,varnodes) \ ... 55## (-Ay(varnodes,dnodes)*dy(dnodes)); 56## msh.p += [ dx'/Nsteps; dy'/Nsteps ] ; 57## triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); 58## pause(.01) 59## endfor 60## @end example 61## 62## @seealso{msh2m_jiggle_mesh} 63## 64## @end deftypefn 65 66function [Ax,Ay] = msh2m_displacement_smoothing(msh, k) 67 68 ## Check input 69 if nargin != 2 # Number of input parameters 70 error("msh2m_displacement_smoothing: wrong number of input parameters."); 71 elseif !(isstruct(msh) && isfield(msh,"p") && 72 isfield(msh,"t") && isfield(msh,"e")) 73 error("msh2m_displacement_smoothing: first input is not a valid mesh structure."); 74 elseif !isscalar(k) 75 error("msh2m_displacement_smoothing: k must be a valid scalar"); 76 endif 77 78 ## Construct matrices 79 x = msh.p(1,:); 80 y = msh.p(2,:); 81 82 dx2 = (x(msh.t([1 2 3],:))-x(msh.t([2 3 1],:))).^2; 83 dy2 = (y(msh.t([1 2 3],:))-y(msh.t([2 3 1],:))).^2; 84 85 l2 = dx2 + dy2; 86 87 Ax = spalloc(length(x),length(x),1); 88 Ay = spalloc(length(x),length(x),1); 89 90 ax = zeros(3,3,columns(msh.t)); 91 ay = zeros(3,3,columns(msh.t)); 92 93 for inode=1:3 94 for jnode=1:3 95 ginode(inode,jnode,:)=msh.t(inode,:); 96 gjnode(inode,jnode,:)=msh.t(jnode,:); 97 endfor 98 endfor 99 100 for ii=1:3 101 for jj=ii+1:3 102 103 ax(ii,jj,:) = ax(jj,ii,:) = reshape(-k * dx2(ii,:)./l2(ii,:),1,1,[]); 104 ay(ii,jj,:) = ay(jj,ii,:) = reshape(-k * dy2(ii,:)./l2(ii,:),1,1,[]); 105 106 ax(ii,ii,:) -= ax(ii,jj,:); 107 ax(jj,jj,:) -= ax(ii,jj,:); 108 ay(ii,ii,:) -= ay(ii,jj,:); 109 ay(jj,jj,:) -= ay(ii,jj,:); 110 111 endfor 112 endfor 113 114 Ax = sparse(ginode(:),gjnode(:),ax(:)); 115 Ay = sparse(ginode(:),gjnode(:),ay(:)); 116 117endfunction 118 119%!demo 120%! msh = msh2m_structured_mesh(linspace(0,1,10), 121%! linspace(0,1,10), 122%! 1,1:4,"left"); 123%! dnodes = msh2m_nodes_on_sides(msh,1:4); 124%! varnodes = setdiff([1:columns(msh.p)],dnodes); 125%! 126%! xd = msh.p(1,dnodes)'; 127%! yd = msh.p(2,dnodes)'; 128%! 129%! dy = zeros(columns(msh.p),1); 130%! dx = dy; 131%! 132%! dxtot = -.5*sin(xd.*yd*pi/2); 133%! dytot = -.5*sin(xd.*yd*pi/2); 134%! 135%! Nsteps = 5; 136%! for ii=1:Nsteps 137%! 138%! dx(dnodes) = dxtot; 139%! dy(dnodes) = dytot; 140%! 141%! [Ax,Ay] = msh2m_displacement_smoothing(msh,1); 142%! 143%! dx(varnodes) = Ax(varnodes,varnodes) \ ... 144%! (-Ax(varnodes,dnodes)*dx(dnodes)); 145%! dy(varnodes) = Ay(varnodes,varnodes) \ ... 146%! (-Ay(varnodes,dnodes)*dy(dnodes)); 147%! 148%! msh.p(1,:) += dx'/Nsteps; 149%! msh.p(2,:) += dy'/Nsteps; 150%! 151%! if mod(ii,2)==0 152%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); 153%! pause(.01) 154%! endif 155%! endfor 156