Mercurial > forge
changeset 6632:b9ccab7c988b octave-forge
Updated OF package to new naming convention. Deleted old functions. Updated INDEX and DESCRIPTION.
line wrap: on
line diff
--- a/extra/msh/DESCRIPTION Mon Feb 01 03:57:21 2010 +0000 +++ b/extra/msh/DESCRIPTION Mon Feb 01 07:38:09 2010 +0000 @@ -1,10 +1,10 @@ -Name: MSH -Version: 0.1.1 -Date: 2009-02-25 -Author: Carlo de Falco, Culpo Massimiliano -Maintainer: Carlo de Falco, Culpo Massimiliano -Title: Meshing Software Package for Octave -Description: Package for creating and managing triangular and tetrahedral meshes for Finite Element or Finite Volume PDE solvers. Uses a mesh data structure compatible with pdetool. Relies on gmsh for unstructured mesh generation. +Name: msh +Version: 1.0.0 +Date: 2010-01-31 +Author: Carlo de Falco, Massimiliano Culpo +Maintainer: Carlo de Falco, Massimiliano Culpo +Title: MeSHing software package for octave +Description: Create and manage triangular and tetrahedral meshes for Finite Element or Finite Volume PDE solvers. Use a mesh data structure compatible with PDEtool. Rely on gmsh for unstructured mesh generation. Depends: octave (>= 3.0), splines SystemRequirements: gmsh (>= 1.6.5), awk Autoload: no
--- a/extra/msh/INDEX Mon Feb 01 03:57:21 2010 +0000 +++ b/extra/msh/INDEX Mon Feb 01 07:38:09 2010 +0000 @@ -1,21 +1,27 @@ -MSH >> MSH - Meshing Software Package for Octave -Mesh creation - MSH2Mstructmesh - MSH3Mstructmesh - MSH2Mgmsh - MSH3Mgmsh - MSH2Mjoinstructm - MSH3Mjoinstructm - MSH2Msubmesh - MSH3Msubmesh - MSH2Mmeshalongspline +MSH >> MSH - MeSHing software package for octave +Structured mesh creation + msh2m_structured_mesh + msh3m_structured_mesh + msh2m_mesh_along_spline +Unstructured mesh creation + msh2m_gmsh + msh3m_gmsh +Mesh manipulation + msh2m_join_structured_mesh + msh3m_join_structured_mesh Mesh properties - MSH2Mgeomprop - MSH3Mgeomprop - MSH2Mtopprop - MSH2Mnodesonsides - MSH3Mnodesonfaces + msh2m_geometrical_properties + msh3m_geometrical_properties + msh2m_topological_properties + msh2m_nodes_on_sides + msh3m_nodes_on_faces Mesh adaptation - MSH2Mequalizemesh - MSH2Mdisplacementsmoothing - MSH2Mjigglemesh \ No newline at end of file + msh2m_equalize_mesh + msh2m_displacement_smoothing + msh2m_jiggle_mesh +Mesh extraction + msh3e_surface_mesh + msh2m_submesh + msh3m_submesh +Mesh plotting + msh2p_mesh \ No newline at end of file
--- a/extra/msh/inst/MSH2Mdisplacementsmoothing.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,154 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{Ax},@var{Ay}]} = MSH2Mdisplacementsmoothing(@var{msh},@var{k}) -## -## To displace the boundary of a 2D mesh, set a spring with -## force/length constant @var{k} along each edge and enforce -## 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 -## to make the system solvable, e.g.: -## -## @example -## msh = MSH2Mstructmesh(linspace(0,1,10),@ -## linspace(0,1,10),@ -## 1,1:4, "left"); -## 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); -## dxtot = dytot = -.5*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 += [ dx'/Nsteps; dy'/Nsteps ] ; -## triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); -## pause(.01) -## endfor -## @end example -## -## @seealso{MSH2Mjigglemesh} -## -## @end deftypefn - -function [Ax,Ay] = MSH2Mdisplacementsmoothing(msh, k) - - x = msh.p(1,:); - y = msh.p(2,:); - - dx2 = (x(msh.t([1 2 3],:))-x(msh.t([2 3 1],:))).^2; - dy2 = (y(msh.t([1 2 3],:))-y(msh.t([2 3 1],:))).^2; - - l2 = dx2 + dy2; - - Ax = spalloc(length(x),length(x),1); - Ay = spalloc(length(x),length(x),1); - - ##for iel=1:columns(msh.t) - - ax = zeros(3,3,columns(msh.t)); - ay = zeros(3,3,columns(msh.t)); - - for inode=1:3 - for jnode=1:3 - ginode(inode,jnode,:)=msh.t(inode,:); - gjnode(inode,jnode,:)=msh.t(jnode,:); - endfor - endfor - - for ii=1:3 - for jj=ii+1:3 - - ax(ii,jj,:) = ax(jj,ii,:) = reshape(-k * dx2(ii,:)./l2(ii,:),1,1,[]); - ay(ii,jj,:) = ay(jj,ii,:) = reshape(-k * dy2(ii,:)./l2(ii,:),1,1,[]); - - 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(:)); - - ##endfor - -endfunction - -%!demo -%! msh = MSH2Mstructmesh(linspace(0,1,10), -%! linspace(0,1,10), -%! 1,1:4,"left"); -%! dnodes = MSH2Mnodesonsides(msh,1:4); -%! varnodes = setdiff([1:columns(msh.p)],dnodes); -%! -%! xd = msh.p(1,dnodes)'; -%! yd = msh.p(2,dnodes)'; -%! -%! dy = zeros(columns(msh.p),1); -%! dx = dy; -%! -%! dxtot = -.5*sin(xd.*yd*pi/2); -%! dytot = -.5*sin(xd.*yd*pi/2); -%! -%! Nsteps = 10; -%! for ii=1:Nsteps -%! ii -%! 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; -%! -%! if mod(ii,2)==0 -%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); -%! pause(.01) -%! end -%! end
--- a/extra/msh/inst/MSH2Mequalizemesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,132 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{Ax}, @var{Ay}, @var{bx}, @ -## @var{by}]} = MSH2Mequalizemesh(@var{msh}) -## -## To equalize the size of triangle edges, apply a baricentric@ -## regularization, i.e. move each node to the @ -## center of mass of the patch of triangles to which it belongs. -## -## May be useful when distorting a mesh, -## type @code{ demo MSH2Mequalizemesh } to see some examples. -## -## @seealso{MSH2Mdisplacementsmoothing} -## -## @end deftypefn - -function [msh] = MSH2Mequalizemesh(msh) - - nel= columns(msh.t); - - x = msh.p(1,:)'; - y = msh.p(2,:)'; - - dnodes = unique(msh.e(1:2,:)(:)); - varnodes = setdiff([1:columns(msh.p)],dnodes); - - Ax = spalloc(length(x),length(x),1); - Ay = spalloc(length(x),length(x),1); - - ax = zeros(3,3,nel); - ay = zeros(3,3,nel); - - 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,:); - endfor - endfor - - for ii=1:3 - - 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(:)); - - x(varnodes) = Ax(varnodes,varnodes) \ (-Ax(varnodes,dnodes)*x(dnodes)); - y(varnodes) = Ay(varnodes,varnodes) \ (-Ay(varnodes,dnodes)*y(dnodes)); - msh.p(1,:) = x'; - msh.p(2,:) = y'; - -endfunction - -%!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,:)'; -%! msh = MSH2Mequalizemesh(msh); -%! 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,:)'; -%! msh = MSH2Mequalizemesh(msh); -%! hold on;triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)');hold off -%! pause(.5) -%! endfor \ No newline at end of file
--- a/extra/msh/inst/MSH2Mgeomprop.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,429 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{varargout}]} = MSH2Mgeomprop(@var{mesh},[@var{string1},@var{string2},...]) -## -## Computes geometrical properties of the specified mesh -## -## Input: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{string1}, @var{string2},...: identifier of the property to compute. Possible choices are listed below. -## @itemize @bullet -## @item "bar" (center of mass): return a matrix with size 2 times number of elements containing the coordinates of the center of mass of every trg. -## @item "cir" (circumcenter): return a matrix with size 2 times number of elements containing the coordinates of the circumcenter of every trg. -## @item "emidp" (boundary edges midpoint): return a matrix with size 2 times number of b.edges containing the coordinates of the midpoint. -## @item "slength" (length of the sides): return a matrix with size 3 times number of elements containing the length of the sides. -## @item "cdist" (distance among circumcenters of neighbouring elements): return a matrix with size 3 times number of elements containing the -## distance among circumcenters of neighbouring elements. If the corresponding side lies on the edge, the distance between -## circumcenter and border edge is returned in the matrix. -## @item "wjacdet" : -## @item "shg": gradient of the P1 shape functions for BIM method -## @item "area" (trg area): return a row vector, with length equal to number of elements, containing the area of every trg in the mesh. -## @item "midedge" (midpoint coordinates of every edge): return a matrix with size 2(x and y coordinates) times 3(edge number) times n of elements -## containing the coordinates of the midpoint of every trg edge. -## @end itemize -## @end itemize -## -## The output will contain the geometrical properties requested in the input in the same order specified in the function call -## @seealso{MSH2Mtopprop} -## @end deftypefn - -function [varargout] = MSH2Mgeomprop(mesh,varargin) - - p = mesh.p; e = mesh.e; t = mesh.t; - ## Number of elements in the mesh - nelem = columns(t); - - ## Check if varargin is composed of strings as expected - if iscellstr(varargin) == 0 - warning("Unexpected input. See help for more information."); - print_usage; - endif - - ## Edges coefficients - [k,j,w] = coeflines(p,t,nelem); - - for nn = 1:length(varargin) - - request = varargin{nn}; - switch request - - case "bar"##Center of mass coordinates - if isfield(mesh,'bar') - varargout{nn} = mesh.bar; - else - [b] = coordinates(p,t,nelem,j,w,k,'bar'); - varargout{nn} = b; - clear b; - endif - - case "cir"##Circum center coordinates - if isfield(mesh,'cir') - varargout{nn} = mesh.cir; - else - [b] = coordinates(p,t,nelem,j,w,k,'cir'); - varargout{nn} = b; - clear b; - endif - - case "emidp"##Boundary edges midpoint coordinates - if isfield(mesh,'emidp') - varargout{nn} = mesh.emidp; - else - [b] = midpoint(p,e); - varargout{nn} = b; - clear b; - endif - - case "slength"##Length of every side - if isfield(mesh,'slength') - varargout{nn} = mesh.slength; - else - [b] = sidelength(p,t,nelem); - varargout{nn} = b; - clear b; - endif - - case "cdist"##Distaneleme among circumcenters of neighbouring elements - if isfield(mesh,'cdist') - varargout{nn} = mesh.cdist; - else - - if isfield(mesh,'cir') - cir = mesh.cir; - else - [cir] = coordinates(p,t,nelem,j,w,k,"cir"); - endif - - if isfield(mesh,'n') - n = mesh.n; - else - [n] = MSH2Mtopprop(mesh,'n'); - endif - - [b] = distance(cir,n,nelem); - [semib] = semidistance(cir,nelem,j,w,k); - border = isnan(n); - [index1] = find( border(1,:) ); - [index2] = find( border(2,:) ); - [index3] = find( border(3,:) ); - b(1,index1) = semib(1,index1); - b(2,index2) = semib(2,index2); - b(3,index3) = semib(3,index3); - varargout{nn} = b; - clear b semib index1 index2 index3 border; - endif - - case "wjacdet" - if isfield(mesh,'wjacdet') - varargout{nn} = mesh.wjacdet; - else - [b] = computearea(p,e,t,'wjac'); - varargout{nn} = b; - clear b - endif - - case "area"##Area of the elements - if isfield(mesh,'area') - varargout{nn} = mesh.area; - else - [b] = computearea(p,e,t,'area'); - varargout{nn} = b; - clear b - endif - - case "shg"##Gradient of the hat functions - if isfield(mesh,'shg') - varargout{nn} = mesh.shg; - else - [b] = shapegrad(p,t); - varargout{nn} = b; - clear b - endif - - case "midedge" - if isfield(mesh,'midedge') - varargout{nn} = mesh.midedge; - else - [b] = midedge(p,t,nelem); - varargout{nn} = b; - clear b; - endif - - otherwise - warning("Unexpected value in passed string. Empty vector passed as output.") - varargout{nn} = []; - endswitch - - endfor - -endfunction - -function [k,j,w] = coeflines(p,t,nelem) - ##The edges are described by the analytical expression: - ## - ## k*x + j*y + w = 0 - ## - ##The coefficients k,j,w are stored in matrixes - - ##i-th edge list, i =1,2,3 - s1 = sort(t(2:3,:),1); - s2 = sort(t([3,1],:),1); - s3 = sort(t(1:2,:),1); - ##Initialization of the matrix data-structure - k = ones(3,nelem); - j = ones(3,nelem); - w = ones(3,nelem); - ##Searching for lines parallel to x axis - [i1] = find( (p(2,s1(2,:)) - p(2,s1(1,:))) != 0 ); - noti1 = setdiff([1:nelem], i1); - [i2] = find( (p(2,s2(2,:)) - p(2,s2(1,:))) != 0 ); - noti2 = setdiff([1:nelem], i2); - [i3] = find( (p(2,s3(2,:)) - p(2,s3(1,:))) != 0 ); - noti3 = setdiff([1:nelem], i3); - ##Computation of the coefficients - ##Edge 1 - j(1,i1) = (p(1,s1(1,i1)) - p(1,s1(2,i1)))./(p(2,s1(2,i1)) - p(2,s1(1,i1))); - w(1,i1) = -(p(1,s1(1,i1)) + p(2,s1(1,i1)).*j(1,i1)); - k(1,noti1) = 0; - j(1,noti1) = 1; - w(1,noti1) = - p(2,s1(1,noti1)); - ##Edge 2 - j(2,i2) = (p(1,s2(1,i2)) - p(1,s2(2,i2)))./(p(2,s2(2,i2)) - p(2,s2(1,i2))); - w(2,i2) = -(p(1,s2(1,i2)) + p(2,s2(1,i2)).*j(2,i2)); - k(2,noti2) = 0; - j(2,noti2) = 1; - w(2,noti2) = - p(2,s2(1,noti2)); - ##Edge 3 - j(3,i3) = (p(1,s3(1,i3)) - p(1,s3(2,i3)))./(p(2,s3(2,i3)) - p(2,s3(1,i3))); - w(3,i3) = -(p(1,s3(1,i3)) + p(2,s3(1,i3)).*j(3,i3)); - k(3,noti3) = 0; - j(3,noti3) = 1; - w(3,noti3) = - p(2,s3(1,noti3)); -endfunction - -function [b] = coordinates(p,t,nelem,j,w,k,string) - ##Computes the coordinate of the geometrical entity specified by string - ##Initialization of the output vectors - b = zeros(2, nelem); - switch string - case "bar" - b(1,:) = ( p(1,t(1,:)) + p(1,t(2,:)) + p(1,t(3,:)) )/3; - b(2,:) = ( p(2,t(1,:)) + p(2,t(2,:)) + p(2,t(3,:)) )/3; - case "cir" - ##Computation of the midpoint of the first two edges - mid1 = zeros(2,nelem); - mid2 = zeros(2,nelem); - ##X coordinate - mid1(1,:) = ( p(1,t(2,:)) + p(1,t(3,:)) )/2; - mid2(1,:) = ( p(1,t(3,:)) + p(1,t(1,:)) )/2; - ##Y coordinate - mid1(2,:) = ( p(2,t(2,:)) + p(2,t(3,:)) )/2; - mid2(2,:) = ( p(2,t(3,:)) + p(2,t(1,:)) )/2; - ##Computation of the intersect between axis 1 and axis 2 - ##Searching for element with edge 1 parallel to x-axes - [parx] = find( j(1,:) == 0); - [notparx] = setdiff(1:nelem, parx); - coefy = zeros(1,nelem); - ##If it is not parallel - coefy(notparx) = ((j(2,notparx)./j(1,notparx)).*k(1,notparx) - k(2,notparx)).^(-1); - b(2,notparx) = coefy(1,notparx).*( j(2,notparx).*mid2(1,notparx) - k(2,notparx).*mid2(2,notparx) + k(1,notparx)./j(1,notparx).*j(2,notparx).*mid1(2,notparx) - j(2,notparx).*mid1(1,notparx) ); - b(1,notparx) =( k(1,notparx).*b(2,notparx) + j(1,notparx).*mid1(1,notparx) - k(1,notparx).*mid1(2,notparx) )./j(1,notparx); - ##If it is parallel - b(2,parx) = mid1(2,parx); - b(1,parx) = k(2,parx)./j(2,parx).*( b(2,parx) - mid2(2,parx) ) + mid2(1,parx); - endswitch -endfunction - -function [b] = midpoint(p,e) - ##Computes the coordinates of the midpoint on the boundary edges - b = zeros(2,columns(e)); - b(1,:) = ( p(1,e(1,:)) + p(1,e(2,:)) )./2; - b(2,:) = ( p(2,e(1,:)) + p(2,e(2,:)) )./2; -endfunction - -function [l] = sidelength(p,t,nelem) - ##Computes the length of every side - - l = zeros(3, nelem); - ##i-th edge list, i =1,2,3 - s1 = sort(t(2:3,:),1); - s2 = sort(t([3,1],:),1); - s3 = sort(t(1:2,:),1); - ##First side length - l(1,:) = sqrt( (p(1,s1(1,:)) - p(1,s1(2,:))).^2 + (p(2,s1(1,:)) - p(2,s1(2,:))).^2 ); - ##Second side length - l(2,:) = sqrt( (p(1,s2(1,:)) - p(1,s2(2,:))).^2 + (p(2,s2(1,:)) - p(2,s2(2,:))).^2 ); - ##Third side length - l(3,:) = sqrt( (p(1,s3(1,:)) - p(1,s3(2,:))).^2 + (p(2,s3(1,:)) - p(2,s3(2,:))).^2 ); -endfunction - -function [d] = semidistance(b,nelem,j,w,k) - ##Computes the distance to the sides of the nodes with coordinates b - ##The edges are described by the analytical expression: - ## - ## k*x + j*y + w = 0 - ## - ##The coefficients k,j,w are stored in matrixes - - ##Initialization of the output vector of distances - d = zeros(3, nelem); - ##Computation of the distances from the geometrical entity to the edges - d(1,:) = abs( k(1,:).*b(1,:) + j(1,:).*b(2,:) + w(1,:))./(sqrt( k(1,:).^2 + j(1,:).^2)); - d(2,:) = abs( k(2,:).*b(1,:) + j(2,:).*b(2,:) + w(2,:))./(sqrt( k(2,:).^2 + j(2,:).^2)); - d(3,:) = abs( k(3,:).*b(1,:) + j(3,:).*b(2,:) + w(3,:))./(sqrt( k(3,:).^2 + j(3,:).^2)); -endfunction - -function [d] = distance(b,n,nelem) - ##Computes the distance between two neighbouring entities - - ##Initialization of the output vector of distances - d = NaN(3, nelem); - ##Trg not on the geometrical border - border = isnan(n); - [index1] = find( border(1,:) == 0 ); - [index2] = find( border(2,:) == 0 ); - [index3] = find( border(3,:) == 0 ); - ##Computation of the distances between two neighboring geometrical entities - d(1,index1) = sqrt( ( b(1,index1) - b(1,n(1,index1)) ).^2 + ( b(2,index1) - b(2,n(1,index1)) ).^2 ); - d(2,index2) = sqrt( ( b(1,index2) - b(1,n(2,index2)) ).^2 + ( b(2,index2) - b(2,n(2,index2)) ).^2 ); - d(3,index3) = sqrt( ( b(1,index3) - b(1,n(3,index3)) ).^2 + ( b(2,index3) - b(2,n(3,index3)) ).^2 ); -endfunction - -function [b] = computearea(p,e,t,string) - ##Compute the area of every element in the mesh - weight = [1/3 1/3 1/3]; - areakk = 1/2; - Nelements = columns(t); - - jac([1,2],:) = [p(1,t(2,:))-p(1,t(1,:)); - p(1,t(3,:))-p(1,t(1,:))]; - jac([3,4],:) = [p(2,t(2,:))-p(2,t(1,:)); - p(2,t(3,:))-p(2,t(1,:))]; - jacdet = jac(1,:).*jac(4,:)-jac(2,:).*jac(3,:); - - degen=find(jacdet <= 0); - if ~isempty(degen) - ## XXX FIXME: there should be a -verbose option to allow to see this - ##fprintf(1,'invalid mesh element: %d fixing...\n',degen); - t(1:3,degen) = t([2,1,3],degen); - jac([1,2],degen) = [p(1,t(2,degen))-p(1,t(1,degen)); - p(1,t(3,degen))-p(1,t(1,degen))]; - jac([3,4],degen) = [p(2,t(2,degen))-p(2,t(1,degen)); - p(2,t(3,degen))-p(2,t(1,degen))]; - jacdet(degen) = jac(1,degen).*jac(4,degen)-jac(2,degen).*jac(3,degen); - endif - - for inode = 1:3 - wjacdet(inode,:) = areakk .* jacdet .* weight(inode); - endfor - - if string == 'wjac' - b = wjacdet(); - elseif string == 'area' - b = sum(wjacdet)'; - endif - -endfunction - -function [d] = midedge(p,t,nelem) - ##Computes the midpoint coordinates for every edge - s1 = t(2:3,:); s2 = t([3,1],:); s3 = t(1:2,:); - edge = cell(3,1); - edge(1) = s1; edge(2) = s2; edge(3) = s3; - d = zeros(2,3,nelem);#Lati * Coordinate * Elementi - for jj = 1:3 - tempx = (p(1,edge{jj}(1,:)) + p(1,edge{jj}(2,:)))/2; - tempy = (p(2,edge{jj}(1,:)) + p(2,edge{jj}(2,:)))/2; - temp = [tempx; tempy]; - d(:,jj,:) = temp; - endfor -endfunction - -function [shg] = shapegrad(p,t) - ##Computes the gradient of the hat functions - x0 = p(1,t(1,:)); - y0 = p(2,t(1,:)); - x1 = p(1,t(2,:)); - y1 = p(2,t(2,:)); - x2 = p(1,t(3,:)); - y2 = p(2,t(3,:)); - - denom = (-(x1.*y0) + x2.*y0 + x0.*y1 - x2.*y1 - x0.*y2 + x1.*y2); - shg(1,1,:) = (y1 - y2)./denom; - shg(2,1,:) = -(x1 - x2)./denom; - shg(1,2,:) = -(y0 - y2)./denom; - shg(2,2,:) = (x0 - x2)./denom; - shg(1,3,:) = (y0 - y1)./denom; - shg(2,3,:) = -(x0 - x1)./denom; -endfunction - -%!test -%! [mesh] = MSH2Mstructmesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); -%! [mesh.bar, mesh.cir, mesh.emidp, mesh.slength, mesh.cdist, mesh.area,mesh.midedge] = MSH2Mgeomprop(mesh,'bar','cir','emidp','slength','cdist','area','midedge'); -%! bar = [0.16667 0.16667 0.66667 0.66667 0.33333 0.33333 0.83333 0.83333 -%! 0.16667 0.66667 0.16667 0.66667 0.33333 0.83333 0.33333 0.83333]; -%! cir = [0.25000 0.25000 0.75000 0.75000 0.25000 0.25000 0.75000 0.75000 -%! 0.25000 0.75000 0.25000 0.75000 0.25000 0.75000 0.25000 0.75000]; -%! emidp =[0.25000 0.75000 1.00000 1.00000 0.25000 0.75000 0.00000 0.00000 -%! 0.00000 0.00000 0.25000 0.75000 1.00000 1.00000 0.25000 0.75000]; -%! slength =[0.70711 0.70711 0.70711 0.70711 0.50000 0.50000 0.50000 0.50000 -%! 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 -%! 0.50000 0.50000 0.50000 0.50000 0.70711 0.70711 0.70711 0.70711]; -%! cdist = [0.00000 0.00000 0.00000 0.00000 0.50000 0.50000 0.25000 0.25000 -%! 0.25000 0.25000 0.50000 0.50000 0.50000 0.25000 0.50000 0.25000 -%! 0.25000 0.50000 0.25000 0.50000 0.00000 0.00000 0.00000 0.00000]; -%! area = [ 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500]; -%! midedge = zeros(2,3,8); -%! midedge(:,:,1) = [0.25000 0.00000 0.25000 -%! 0.25000 0.25000 0.00000]; -%! midedge(:,:,2) = [0.25000 0.00000 0.25000 -%! 0.75000 0.75000 0.50000]; -%! midedge(:,:,3) = [0.75000 0.50000 0.75000 -%! 0.25000 0.25000 0.00000]; -%! midedge(:,:,4) = [0.75000 0.50000 0.75000 -%! 0.75000 0.75000 0.50000]; -%! midedge(:,:,5) = [0.50000 0.25000 0.25000 -%! 0.25000 0.50000 0.25000]; -%! midedge(:,:,6) = [0.50000 0.25000 0.25000 -%! 0.75000 1.00000 0.75000]; -%! midedge(:,:,7) = [1.00000 0.75000 0.75000 -%! 0.25000 0.50000 0.25000]; -%! midedge(:,:,8) = [1.00000 0.75000 0.75000 -%! 0.75000 1.00000 0.75000]; -%! toll = 1e-4; -%! assert(mesh.bar,bar,toll); -%! assert(mesh.cir,cir,toll); -%! assert(mesh.emidp,emidp,toll); -%! assert(mesh.slength,slength,toll); -%! assert(mesh.cdist,cdist,toll); -%! assert(mesh.area,area,toll); -%! assert(mesh.midedge,midedge,toll);
--- a/extra/msh/inst/MSH2Mgmsh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,151 +0,0 @@ -## Copyright (C) 2007,2009 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH2Mgmsh(@var{geometry},@var{option},@var{value},...) -## -## Construct an unstructured triangular 2D mesh making use of the free -## software gmsh. Return a PDE-tool like mesh structure. -## -## Input: -## @itemize @minus -## @item @var{geometry}: basename of the ".geo" file to be meshed. -## @item @var{option}: string of the option to pass gmsh. -## @item @var{value}: value of the option to pass gmsh. -## @end itemize -## -## For more information regarding the possible option to pass, refer to gmsh manual or gmsh site: -## http://www.geuz.org/gmsh/ -## -## @var{mesh} structure is composed of the following fields: -## -## @itemize @minus -## @item @var{p}: 2 X (# nodes) matrix. Contain mesh points coordinates. -## @item @var{e}: 7 X (# side edges) matrix. Contain mesh side -## edges information. -## @item @var{t}: 4 X (# triangles) matrix. Contain pointer to @var{p} -## field, as well as region number. -## @end itemize -## -## @seealso{MSH2Mstructmesh, MSH2Mjoinstructm, MSH2Msubmesh} -## @end deftypefn - -function [mesh] = MSH2Mgmsh(geometry,varargin) - - ## 1. Check input - ## 1a. Number of input - if !mod(nargin,2) - warning("WRONG NUMBER OF INPUT."); - print_usage; - endif - - ## 2. Build mesh - noptions = (nargin - 1) / 2; # Number of passed options - - ## 2a. Construct system command string - optstring = ""; - for ii = 1:noptions - option = varargin{2*(ii)-1}; - value = varargin{2*ii}; - if !ischar(value) - value = num2str(value); - endif - optstring = [optstring," -",option," ",value]; - endfor - - ## 2b. Invoke gmsh - printf("\n"); - printf("Generating mesh...\n"); - system(["gmsh -format msh -2 -o " geometry ".msh" optstring " " geometry ".geo"]); - - ## 2c. Build structure fields - printf("Processing gmsh data...\n"); - ## Points - com_p = "awk '/\\$Nodes/,/\\$EndNodes/ {print $2, $3 > ""msh_p.txt""}' "; - ## Side edges - com_e = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""1"") print $7, $8,$5 > ""msh_e.txt""}' "; - ## Triangles - com_t = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""2"") print $7, $8, $9, $5 > ""msh_t.txt""}' "; - - command = [com_p,geometry,".msh ; "]; - command = [command,com_e,geometry,".msh ; "]; - command = [command,com_t,geometry,".msh"]; - - system(command); - - ## 2d. Create PDE-tool like structure - printf("Creating PDE-tool like mesh...\n"); - p = load("msh_p.txt")'; # Mesh-points - tmp = load("msh_e.txt")'; # Mesh surface-edges - be = zeros(7,columns(tmp)); - be([1,2,5],:) = tmp; - t = load("msh_t.txt")'; # Mesh tetrahedra - - ## 3. Remove hanging nodes - printf("Check for hanging nodes...\n"); - nnodes = columns(p); - in_msh = intersect( 1:nnodes , t(1:3,:) ); - if length(in_msh) != nnodes - new_num(in_msh) = [1:length(in_msh)]; - t(1:3,:) = new_num(t(1:3,:)); - be(1:2,:) = new_num(be(1:2,:)); - p = p(:,in_msh); - endif - - ## 4. Set region numbers in edge structure - printf("Setting region number in edge structure...\n"); - msh = struct("p",p,"t",t,"e",be); - msh.be = MSH2Mtopprop(msh, "boundary"); - msh.e(6,:) = msh.t(4,msh.be(1,:)); - jj = find (sum(msh.be>0)==4); - msh.e(7,jj) = msh.t(4,msh.be(3,jj)); - - mesh = struct("p",p,"e",be,"t",t); - - ## 5. Delete temporary files - printf("Deleting temporary files...\n"); - system(["rm -f msh_p.txt msh_e.txt msh_t.txt msh_s.txt *.msh"]); - -endfunction - -%!test -%! fid = fopen("circle.geo","w"); -%! fprintf(fid,"Point(1) = {0, 0, 0, 1};\n"); -%! fprintf(fid,"Point(2) = {1, 0, 0, 1};\n"); -%! fprintf(fid,"Point(3) = {-1, 0, 0, 1};\n"); -%! fprintf(fid,"Circle(1) = {3, 1, 2};\n"); -%! fprintf(fid,"Circle(2) = {2, 1, 3};\n"); -%! fprintf(fid,"Line Loop(4) = {2, 1};\n"); -%! fprintf(fid,"Plane Surface(4) = {4};"); -%! fclose(fid); -%! mesh = MSH2Mgmsh("circle","v",0); -%! system("rm circle.geo"); -%! nnodest = length(unique(mesh.t)); -%! nnodesp = columns(mesh.p); -%! assert(nnodest,nnodesp); \ No newline at end of file
--- a/extra/msh/inst/MSH2Mjigglemesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,112 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[newmsh]} = MSH2Mjigglemesh(@var{msh}, @var{steps}) -## -## To equalize the size of triangle edges, set a spring of rest -## length @var{factor}*@var{area} along each edge of the mesh and -## solve for static equilibrium. -## -## The non-linear eqautions of the system obtained are soved via a -## non-linear Gass-Seidel method. @var{step} is the number of steps of -## the method to be applied. -## -## May be useful when distorting a mesh, type @code{demo -## MSH2Mjigglemesh} to see some examples. -## -## @seealso{MSH2Mdisplacementsmoothing, MSH2Mequalizemesh} -## -## @end deftypefn - -function [msh] = MSH2Mjigglemesh(msh, steps) - - nel= columns(msh.t); - nnodes = columns(msh.p); - - x = msh.p(1,:)'; - y = msh.p(2,:)'; - - dnodes = unique(msh.e(1:2,:)(:)); - vnodes = setdiff(1:nnodes,dnodes); - - ## find node neighbours XXX FIXME: should this go - ## into MSH2Mtopprop ? - sides = MSH2Mtopprop(msh,'sides'); - for inode = 1:nnodes - neig{inode} = (sides(:, sides(1,:) == inode | sides(2,:) == inode))(:); - neig{inode} (neig{inode} == inode) = []; - endfor - - for istep = 1:steps - for inode =vnodes - xx = x(neig{inode}) * ones(size(neig{inode}))'; - lx = abs ( xx - xx' )(:); - mx = ( xx + xx' )(:)/2; - x(inode) = sum(mx.*lx)/sum(lx); - - yy = y(neig{inode}) * ones(size(neig{inode}))'; - ly = abs ( yy - yy' )(:); - my = (yy + yy')(:)/2; - y(inode) = sum(my.*ly)/sum(ly); - endfor - endfor - - msh.p = [x';y']; - -endfunction - -%!demo -%! ### distort a mesh on a square equalizing at each step -%! msh = MSH2Mstructmesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"right"); -%! 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 = -.4*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,:)'; -%! msh = MSH2Mjigglemesh(msh,10); -%! hold on;triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)');hold off -%! pause(.5) -%! endfor \ No newline at end of file
--- a/extra/msh/inst/MSH2Mjoinstructm.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,158 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH2Mjoinstructm(@var{mesh1},@var{mesh2},@var{s1},@var{s2}) -## -## Join two structured meshes (created by MSH2Mstructmesh) into one -## mesh structure variable. -## -## Input: -## @itemize @minus -## @item @var{mesh1}, @var{mesh2}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{s1}, @var{s2}: number of the corresponding geometrical border edge for respectively mesh1 and mesh2. -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @end itemize -## -## WARNING: on the common edge the two meshes must share the same vertexes. -## -## @seealso{MSH2Mstructmesh,MSH2Mgmsh,MSH2Msubmesh} -## @end deftypefn - -function [mesh] = MSH2Mjoinstructm(mesh1,mesh2,s1,s2) - - ## make sure that the outside world is always - ## on the same side of the boundary of mesh1 - [mesh1.e(6:7,:),I] = sort(mesh1.e(6:7,:)); - for ic=1:size(mesh1.e,2) - mesh1.e(1:2,ic) = mesh1.e(I(:,ic),ic); - endfor - - intnodes1=[]; - intnodes2=[]; - - j1=[];j2=[]; - for is=1:length(s1) - side1 = s1(is); - side2 = s2(is); - [i,j] = find(mesh1.e(5,:)==side1); - j1=[j1 j]; - [i,j] = find(mesh2.e(5,:)==side2); - oldregion(side1) = max(max(mesh2.e(6:7,j))); - j2=[j2 j]; - endfor - - intnodes1=[mesh1.e(1,j1),mesh1.e(2,j1)]; - intnodes2=[mesh2.e(1,j2),mesh2.e(2,j2)]; - - intnodes1 = unique(intnodes1); - [tmp,I] = sort(mesh1.p(1,intnodes1)); - intnodes1 = intnodes1(I); - [tmp,I] = sort(mesh1.p(2,intnodes1)); - intnodes1 = intnodes1(I); - - intnodes2 = unique(intnodes2); - [tmp,I] = sort(mesh2.p(1,intnodes2)); - intnodes2 = intnodes2(I); - [tmp,I] = sort(mesh2.p(2,intnodes2)); - intnodes2 = intnodes2(I); - - ##delete redundant edges - mesh2.e(:,j2) = []; - - ## change edge numbers - indici=[]; - consecutivi=[]; - indici = unique(mesh2.e(5,:)); - consecutivi (indici) = [1:length(indici)]+max(mesh1.e(5,:)); - mesh2.e(5,:)=consecutivi(mesh2.e(5,:)); - - - ##change node indices in connectivity matrix - ##and edge list - indici=[];consecutivi=[]; - indici = 1:size(mesh2.p,2); - offint = setdiff(indici,intnodes2); - consecutivi (offint) = [1:length(offint)]+size(mesh1.p,2); - consecutivi (intnodes2) = intnodes1; - mesh2.e(1:2,:)=consecutivi(mesh2.e(1:2,:)); - mesh2.t(1:3,:)=consecutivi(mesh2.t(1:3,:)); - - - ##delete redundant points - mesh2.p(:,intnodes2) = []; - - ##set region numbers - regions = unique(mesh1.t(4,:)); - newregions(regions) = 1:length(regions); - mesh1.t(4,:) = newregions(mesh1.t(4,:)); - - ##set region numbers - regions = unique(mesh2.t(4,:)); - newregions(regions) = [1:length(regions)]+max(mesh1.t(4,:)); - mesh2.t(4,:) = newregions(mesh2.t(4,:)); - ##set adjacent region numbers in edge structure 2 - [i,j] = find(mesh2.e(6:7,:)); - i = i+5; - mesh2.e(i,j) = newregions(mesh2.e(i,j)); - ##set adjacent region numbers in edge structure 1 - mesh1.e(6,j1) = newregions(oldregion(mesh1.e(5,j1))); - - ##make the new p structure - mesh.p = [mesh1.p mesh2.p]; - mesh.e = [mesh1.e mesh2.e]; - mesh.t = [mesh1.t mesh2.t]; - -endfunction - -%!test -%! [mesh1] = MSH2Mstructmesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); -%! [mesh2] = MSH2Mstructmesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); -%! [mesh] = MSH2Mjoinstructm(mesh1,mesh2,2,4); -%! p = [0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 1.00000 1.00000 1.00000 1.50000 1.50000 1.50000 2.00000 2.00000 2.00000 -%! 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000]; -%! e = [1 4 7 8 3 6 1 2 7 10 13 14 9 12 -%! 4 7 8 9 6 9 2 3 10 13 14 15 12 15 -%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -%! 1 1 2 2 3 3 4 4 5 5 6 6 7 7 -%! 0 0 2 2 0 0 0 0 0 0 0 0 0 0 -%! 1 1 1 1 1 1 1 1 2 2 2 2 2 2]; -%! t = [1 2 4 5 2 3 5 6 7 8 10 11 8 9 11 12 -%! 4 5 7 8 4 5 7 8 10 11 13 14 10 11 13 14 -%! 2 3 5 6 5 6 8 9 8 9 11 12 11 12 14 15 -%! 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]; -%! toll = 1e-4; -%! assert(mesh.p,p,toll); -%! assert(mesh.p,p,toll); -%! assert(mesh.p,p,toll);
--- a/extra/msh/inst/MSH2Mmeshalongspline.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,92 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{msh}]} = -## MSH2Mmeshalongspline(@var{xc}, @var{yc}, @var{Nnx}, @var{Nny}, -## @var{sigma}) -## -## Generates a structured mesh in a thin layer of size @var{sigma} -## sitting on a natural Catmull-Rom type cubic spline with control -## points @var{xc}, @var{yc}. -## If @var{Nnx} and @var{Nny} are scalars, -## the mesh will have @var{Nnx} nodes in the direction along the -## spline and @var{Nny} in the normal direction. -## If @var{Nnx} and @var{Nny} are vectors they will indicate -## the curvilinear coordinates of the mesh nodes. -## Be aware that if @var{sigma} is not much smaller than the curvature -## of the line the resulting mesh may be invalid. -## -## @seealso{MSH2Mstructmesh} -## @end deftypefn - - -function msh2 = MSH2Mmeshalongspline(xc,yc,Nnx,Nny,sigma) - - if (nargin != 5) - print_usage(); - else - s = [0:length(xc)-1]; - xsPP = catmullrom ( s, xc ); - ysPP = catmullrom ( s, yc ); - - if (length(Nnx)>1) - ss = Nnx(:).'; - else - ss = linspace(0,s(end),Nnx); - endif - xs = ppval(xsPP,ss); - ys = ppval(ysPP,ss); - - dxsPP = fnder(xsPP,1); - dysPP = fnder(ysPP,1); - - - nx = -ppval(dysPP,ss)'; - ny = ppval(dxsPP,ss)'; - - nx = nx ./ sqrt(nx.^2+ny.^2); - ny = ny ./ sqrt(nx.^2+ny.^2); - - if (length(Nny)>1) - ssy = Nny(:).'; - else - ssy = linspace(0,1,Nny); - endif - - msh2 = MSH2Mstructmesh ([1:length(ss)], ssy, 1, 1:4); - - jj = (msh2.p(1,:)); - p(1,:) = xs(jj) + sigma*nx(jj)' .* msh2.p(2,:); - p(2,:) = ys(jj) + sigma*ny(jj)' .* msh2.p(2,:); - - msh2.p = p; - endif - -endfunction \ No newline at end of file
--- a/extra/msh/inst/MSH2Mnodesonsides.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,72 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{nodelist}]} =@ -## MSH2Mnodesonsides(@var{mesh}, @var{sidelist}) -## -## Returns a list of the nodes lying on the sides @var{sidelist} of -## the mesh @var{mesh}. -## -## Input: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{sidelist}: row vector containing the number of the sides (numbering referred to mesh.e(5,:)). -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{nodelist}: list of the nodes that lies on the specified sides. -## @end itemize -## -## @seealso{MSH2Mgeomprop,MSH2Mtopprop} -## @end deftypefn - -function [nodelist] = MSH2Mnodesonsides(mesh,sidelist) - - edgelist =[]; - - for ii = 1:length(sidelist) - edgelist=[edgelist,find(mesh.e(5,:)==sidelist(ii))]; - endfor - - ##Set list of nodes with Dirichelet BCs - nodelist = mesh.e(1:2,edgelist); - nodelist = [nodelist(1,:) nodelist(2,:)]; - nodelist = unique(nodelist); - -endfunction - -%!test -%! [mesh1] = MSH2Mstructmesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); -%! [mesh2] = MSH2Mstructmesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); -%! [mesh] = MSH2Mjoinstructm(mesh1,mesh2,2,4); -%! [nodelist] = MSH2Mnodesonsides(mesh,[1 2]); -%! reallist = [1 4 7 8 9]; -%! assert(nodelist,reallist);
--- a/extra/msh/inst/MSH2Mstructmesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,209 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH2Mstructmesh(@var{x},@var{y},@var{region},@var{sides},@var{string}) -## -## Constructs a structured triangular 2D mesh on a rectangular domain, -## and returns a PDEtool-like mesh structure. -## -## Input: -## @itemize @minus -## @item @var{x}: vector representing the 1D meshing of the side parallel to x axis. -## @item @var{y}: vector representing the 1D meshing of the side parallel to y axis. -## @item @var{region}: number assigned to the meshed region. -## @item @var{sides}: row vector containing the four numbers assigned to the geometrical edges. -## @item @var{string}: (optional) orientation of the diagonal edge of the structured mesh. -## Values: "right", "left", "random". Default is "right". -## @end itemize -## -## Output: mesh basic structure, composed of the following fields -## @itemize @minus -## @item @var{p}: matrix with size 2 times number of mesh point. -## @itemize @bullet -## @item 1st row: x-coordinates of the points. -## @item 2nd row: y-coordinates of the points. -## @end itemize -## @item @var{e}: matrix with size 7 times number of mesh border edges. -## @itemize @bullet -## @item 1st row: p-matrix column number of the first edge-vertex. -## @item 2nd row: p-matrix column number of the second edge-vertex. -## @item 3rd row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 4th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 5th row: number of the geometrical border upon which the referred mesh edge is lying on. -## @item 6th row: number of the region to the right of the referred mesh edge. -## @item 7th row: number of the region to the left of the referred mesh edge. -## @end itemize -## @item @var{t}: -## @itemize @bullet -## @item 1st row: p-matrix column number of the first trg-vertex. -## @item 2nd row: p-matrix column number of the second trg-vertex. -## @item 3rd row: p-matrix column number of the third trg-vertex. -## @item 4th row: number of the region upon which the referred trg is lying on. -## @end itemize -## @end itemize -## -## @seealso{MSH2Mgmsh,MSH2Mjoinstructm,MSH2Msubmesh} -## @end deftypefn - -function [mesh] = MSH2Mstructmesh(x,y,region,sides,varargin) - - default = 'right'; - - if length(varargin)==0 - string = default; - else - string = varargin{1}; - endif - - switch string - case 'right' - [mesh] = Ustructmesh_right(x, y, region, sides); - case 'left' - [mesh] = Ustructmesh_left(x, y, region, sides); - case 'random' - [mesh] = Ustructmesh_random(x, y, region, sides); - otherwise - error(['Passed string has not a valid value: ' string]); - endswitch - -endfunction - -##RIGHT DIAGONAL STRUCTURED MESH -function [mesh]=Ustructmesh_right(x,y,region,sides) - - x = sort(x); - y = sort(y); - - nx = length(x); - ny = length(y); - [XX,YY] = meshgrid(x,y); - p = [XX(:),YY(:)]'; - iiv (ny,nx)=0; - iiv(:)=1:nx*ny; - iiv(end,:)=[]; - iiv(:,end)=[]; - iiv=iiv(:)'; - t = [[iiv;iiv+ny;iiv+ny+1],[iiv;iiv+ny+1;iiv+1] ]; - t (4,:)=region; - - l1 = 1+ny*([1:nx]-1); - l4 = 1:ny; - l2 = ny*(nx-1)+1:nx*ny; - l3 = ny + l1 -1; - - e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) - l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region - ]; - mesh.p = p; mesh.e = e; mesh.t = t; - -endfunction - -##LEFT DIAGONAL STRUCTURED MESH -function [mesh]=Ustructmesh_left(x,y,region,sides) - - x = sort(x); - y = sort(y); - - nx = length(x); - ny = length(y); - [XX,YY] = meshgrid(x,y); - p = [XX(:),YY(:)]'; - iiv (ny,nx)=0; - iiv(:)=1:nx*ny; - iiv(end,:)=[]; - iiv(:,end)=[]; - iiv=iiv(:)'; - t = [[iiv;iiv+ny;iiv+1],[iiv+1;iiv+ny;iiv+ny+1] ]; - t (4,:)=region; - - l1 = 1+ny*([1:nx]-1); - l4 = 1:ny; - l2 = ny*(nx-1)+1:nx*ny; - l3 = ny + l1 -1; - - e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) - l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region - ]; - mesh.p = p; mesh.e = e; mesh.t = t; - -endfunction - -##RANDOM DIAGONAL STRUCTURED MESH -function [mesh]=Ustructmesh_random(x,y,region,sides) - - x = sort(x); - y = sort(y); - - nx = length(x); - ny = length(y); - [XX,YY] = meshgrid(x,y); - p = [XX(:),YY(:)]'; - iiv (ny,nx)=0; - iiv(:)=1:nx*ny; - iiv(end,:)=[]; - iiv(:,end)=[]; - iiv=iiv(:)'; - - niiv = length(iiv); - theperm = iiv(randperm(niiv)); - first = theperm(1:floor(niiv/2)); - second = theperm(floor(niiv/2)+1:end); - - t = [[first;first+ny;first+ny+1],[first;first+ny+1;first+1] ]; - t = [t,[second;second+ny;second+1],[second+ny;second+ny+1;second+1] ]; - - t (4,:)=region; - - l1 = 1+ny*([1:nx]-1); - l4 = 1:ny; - l2 = ny*(nx-1)+1:nx*ny; - l3 = ny + l1 -1; - - e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) - l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 - [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region - ]; - mesh.p = p; mesh.e = e; mesh.t = t; - -endfunction
--- a/extra/msh/inst/MSH2Msubmesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,113 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{omesh},@var{nodelist},@var{elementlist}]} = MSH2Msubmesh(@var{imesh},@var{intrfc},@var{sdl}) -## -## Gives as output a specified submesh, and the lists of nodes and element that are mantained. -## -## Input: -## @itemize @minus -## @item @var{imesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{intrfc}: row vector containing the number of the internal interface sides (numbering referred to mesh.e(5,:)) that should be mantained. -## Could be empty. -## @item @var{sdl} (subdomain list): row vector containing the list of the subdomain that are going to be extracted. -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{omesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{nodelist}: list of the node of the original mesh that are present in the restricted one. -## @item @var{elementlist}: list of the element of the original mesh that are present in the restricted one. -## @end itemize -## -## @seealso{MSH2Mstructmesh,MSH2Mjoinstructm,MSH2Mgmsh} -## @end deftypefn - -function [omesh,nodelist,elementlist] = MSH2Msubmesh(imesh,intrfc,sdl) - - nsd = length(sdl); - - ## set list of output triangles - elementlist=[]; - for isd=1:nsd - elementlist = [elementlist find(imesh.t(4,:) == sdl(isd))]; - endfor - - omesh.t = imesh.t(:,elementlist); - - ## set list of output nodes - nodelist = unique(reshape(imesh.t(1:3,elementlist),1,[])); - omesh.p = imesh.p(:,nodelist); - - ## use new node numbering in connectivity matrix - indx(nodelist) = [1:length(nodelist)]; - iel = [1:length(elementlist)]; - omesh.t(1:3,iel) = indx(omesh.t(1:3,iel)); - - ## set list of output edges - omesh.e =[]; - for isd=1:nsd - omesh.e = [omesh.e imesh.e(:,imesh.e(7,:)==sdl(isd))]; - omesh.e = [omesh.e imesh.e(:,imesh.e(6,:)==sdl(isd))]; - endfor - omesh.e=unique(omesh.e',"rows")'; - - ## use new node numbering in boundary segment list - ied = [1:size(omesh.e,2)]; - omesh.e(1:2,ied) = indx(omesh.e(1:2,ied)); - -endfunction - -%!test -%! [mesh1] = MSH2Mstructmesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); -%! [mesh2] = MSH2Mstructmesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); -%! [mesh] = MSH2Mjoinstructm(mesh1,mesh2,2,4); -%! [omesh,nodelist,elementlist] = MSH2Msubmesh(mesh,[],2); -%! p = [1.00000 1.00000 1.00000 1.50000 1.50000 1.50000 2.00000 2.00000 2.00000 -%! 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000]; -%! e = [1 1 2 3 4 6 7 8 -%! 2 4 3 6 7 9 8 9 -%! 0 0 0 0 0 0 0 0 -%! 0 0 0 0 0 0 0 0 -%! 2 5 2 7 5 7 6 6 -%! 2 0 2 0 0 0 0 0 -%! 1 2 1 2 2 2 2 2]; -%! t = [1 2 4 5 2 3 5 6 -%! 4 5 7 8 4 5 7 8 -%! 2 3 5 6 5 6 8 9 -%! 2 2 2 2 2 2 2 2]; -%! nl = [7 8 9 10 11 12 13 14 15]; -%! el = [9 10 11 12 13 14 15 16]; -%! toll = 1e-4; -%! assert(omesh.p,p,toll); -%! assert(omesh.e,e); -%! assert(omesh.t,t); -%! assert(nodelist,nl); -%! assert(elementlist,el);
--- a/extra/msh/inst/MSH2Mtopprop.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,296 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{varargout}]} = MSH2Mtopprop(@var{mesh}, [@var{string1}, @var{string2},...]) -## -## Computes some topological properties of the mesh @var{mesh}. -## -## Input: -## @itemize @minus -## @item @var{mesh}: mesh structure with (at least) fields @code{p}, -## @code{e} and @code{t} a detailed in the help for @code{MSH2Mstructmesh}. -## @item @var{string1}, @var{string2},...: properties to compute, -## available choices are: -## -## @itemize @bullet -## @item "n" (neighbouring triangles): a matrix @code{M} of size 3 by -## the number of elements. Element @code{M(i,j)} will be the index of -## the mesh element sharing with elemnt number @code{j} its -## @code{i}-th edge if no such element exists (i.e. for boundary -## edges) a value of @code{NaN} is set. -## @item "sides" (global edge matrix): a matrix @code{M} of size 2 by -## total number of edges. Element @code{M(i,j)} will be the index of -## the vertex at the @code{i}-th end of the @code{j}-th edge. -## @item "ts" (triangle sides matrix): a matrix @code{M} of size 3 by -## number of elements. Element @code{M(i,j)} will be the index of the -## @code{i}-th side of the @code{j}-th mesh element. -## @item "tws" (trg with sides matrix):a matrix @code{M} of size 2 by -## the total number of edges. -## @code{M(:,j)} will the set of indeces of the mesh -## elements to whose boundary the @code{j}-th mesh edge belongs. -## For an edge @code{j} belonging to one triangle onlly (an external -## boundary edge) @code{M(2,j) = NaN}. -## @item "coinc" (coincident circumcenter matrix): a matrix @code{M} -## with 2 rows. -## Each column will contain the indices of two triangles sharing the -## same circumcenter. -## @item "boundary" (boundary edge matrix): a matrix @code{M} of size 2 -## times the total number of boundary edges. -## @code{M(1,:)} will be the set of mesh elements (possibly repeated) -## with at least one edge belonging to the external boundary. -## @code{0 < M(2,j) < 4} will be the (local) index of an edge of the -## @code{j}-th mesh element that lies on the external boundary. -## @end itemize -## @end itemize -## -## The output will contain the geometrical properties requested in the input in the same order specified in the function call -## @seealso{MSHM2geomprop} -## @end deftypefn - -function [varargout] = MSH2Mtopprop(mesh,varargin) - - p = mesh.p; e = mesh.e; t = mesh.t; - ## Number of elements in the mesh - nelem = columns(t); - [n,ts,tws,sides] = neigh(t,nelem); - - for nn = 1:length(varargin) - request = varargin{nn}; - switch request - - case "n" - if isfield(mesh,'n') - varargout{nn} = mesh.n; - else - varargout{nn} = n; - endif - - case "sides" - if isfield(mesh,'sides') - varargout{nn} = mesh.sides; - else - varargout{nn} = sides; - endif - - case "ts" - if isfield(mesh,'ts') - varargout{nn} = mesh.ts; - else - varargout{nn} = ts; - endif - - case "tws" - if isfield(mesh,'tws') - varargout{nn} = mesh.tws; - else - varargout{nn} = tws; - endif - - case "coinc" - if isfield(mesh,'coinc') - varargout{nn} = mesh.coinc; - else - ##Compute the matrix listing the trgs that shares the circum-centre - if isfield(mesh,'cdist') - d = mesh.cdist; - else - [d] = MSH2Mgeomprop(mesh,'cdist'); - endif - [b] = coinc(n,d); - varargout{nn} = b; - clear b - endif - - case "boundary" - if isfield(mesh,'boundary') - varargout{nn} = mesh.boundary; - else - [b] = borderline(e,t); - varargout{nn} = b; - clear b - endif - - otherwise - warning("Unexpected value in passed string. Empty vector passed as output.") - varargout{nn} = []; - endswitch - - endfor - -endfunction - -function [n,ts,triwside,sides] = neigh(t,nelem) - - n = nan*ones(3,nelem); - t = t(1:3,:); - - s3 = sort(t(1:2,:),1); - s1 = sort(t(2:3,:),1); - s2 = sort(t([3,1],:),1); - - allsides = [s1 s2 s3]'; - [sides, ii, jj] = unique( allsides,'rows'); - sides = sides'; - - ts = reshape(jj,[],3)'; - - triwside = zeros(2,columns(sides)); - for kk =1:3 - triwside(1,ts(kk,1:end)) = 1:nelem; - triwside(2,ts(4-kk,end:-1:1)) = nelem:-1:1; - endfor - - triwside(2,triwside(1,:)==triwside(2,:)) = NaN; - - n(1,:) = triwside(1,ts(1,:)); - n(1,n(1,:)==1:nelem) = triwside(2,ts(1,:))(n(1,:)==1:nelem); - n(2,:) = triwside(1,ts(2,:)); - n(2,n(2,:)==1:nelem) = triwside(2,ts(2,:))(n(2,:)==1:nelem); - n(3,:) = triwside(1,ts(3,:)); - n(3,n(3,:)==1:nelem) = triwside(2,ts(3,:))(n(3,:)==1:nelem); - -endfunction - -function [output] = coinc(n,d); - ##Compute the matrix listing the trgs that shares the circum-centre - - ##Tolerance value for considering two point to be coincident - toll = 1e-10; - ##Check the presence of more than two trgs sharing the same circum centre - degen = d < toll; res = sum(degen); - [check] = find(res > 1); - ##Index of the sharing pairs - [ii, jj] = find(degen >= 1); - if isempty(jj) == 0 - temp = zeros(2,length(jj)); - temp(1,:) = jj'; - temp(2,:) = diag(n(ii,jj))'; - temp = sort(temp); - temp = temp'; - [output] = unique(temp,'rows'); - output = output'; - if isempty(check) == 0 - warning("More than two trgs sharing the same circum-centre.") - ##FIXME if more than two trgs shares the same circen ---> construct a cell array - endif - else - output = []; - endif -endfunction - -function [output] = borderline(e,t) - ##Compute a vector which contains the element number for the corresponding edge - - nelem = columns(e); - t = t(1:3,:); - output = zeros(4,nelem); - for ii = 1:nelem - - point = ( e(1,ii) == t ); - point += ( e(2,ii) == t ); - - [jj1] = find( sum(point(2:3,:)) == 2); - [jj2] = find( sum(point([3 1],:)) == 2); - [jj3] = find( sum(point(1:2,:)) == 2); - - assert( (length(jj1) + length(jj2) + length(jj3)) <= 2 ); - - numtrg = 0; - for jj=1:length(jj1) - output(2*numtrg+1,ii) = jj1(jj); - output(2*numtrg+2,ii) = 1; - numtrg += 1; - endfor - for jj=1:length(jj2) - output(2*numtrg+1,ii) = jj2(jj); - output(2*numtrg+2,ii) = 2; - numtrg += 1; - endfor - for jj=1:length(jj3) - output(2*numtrg+1,ii) = jj3(jj); - output(2*numtrg+2,ii) = 3; - numtrg += 1; - endfor - - endfor -endfunction - -%!test -%! [mesh] = MSH2Mstructmesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); -%! [mesh.n,mesh.sides,mesh.ts,mesh.tws,mesh.coinc,mesh.boundary] = MSH2Mtopprop(mesh,'n','sides','ts','tws','coinc','boundary'); -%! n = [5 6 7 8 3 4 NaN NaN -%! NaN NaN 5 6 2 NaN 4 NaN -%! NaN 5 NaN 7 1 2 3 4]; -%! sides = [1 1 2 2 2 3 3 4 4 5 5 5 6 6 7 8 -%! 2 4 3 4 5 5 6 5 7 6 7 8 8 9 8 9]; -%! ts = [4 6 11 13 8 10 15 16 -%! 1 3 8 10 5 7 12 14 -%! 2 5 9 12 4 6 11 13]; -%! tws = [ 1 1 2 5 2 6 6 3 3 4 7 4 8 8 7 8 -%! NaN NaN NaN 1 5 2 NaN 5 NaN 6 3 7 4 NaN NaN NaN]; -%! coinc = [1 2 3 4 -%! 5 6 7 8]; -%! boundary =[ 1 3 7 8 6 8 1 2 -%! 3 3 1 1 2 2 2 2 -%! 0 0 0 0 0 0 0 0 -%! 0 0 0 0 0 0 0 0]; -%! assert(mesh.n,n); -%! assert(mesh.sides,sides); -%! assert(mesh.ts,ts); -%! assert(mesh.tws,tws); -%! assert(mesh.coinc,coinc); -%! assert(mesh.boundary,boundary); - -%!test -%! mesh.p = []; mesh.e = []; -%! mesh.t = [3 9 10 1 6 9 10 9 8 9 -%! 9 3 1 10 10 10 7 5 9 8 -%! 6 5 7 8 2 6 2 4 4 10 -%! 6 6 6 6 6 6 6 6 6 6]; -%! [mesh.n] = MSH2Mtopprop(mesh,'n'); -%! n = [6 NaN NaN 10 7 5 NaN NaN 8 4 -%! NaN 8 7 NaN NaN 1 5 9 NaN 6 -%! 2 1 4 3 6 10 3 2 10 9]; -%! assert(mesh.n,n); - - -%!test -%! mesh.p = []; mesh.e = []; -%! mesh.t =[ -%! 10 3 6 11 10 3 6 11 1 7 5 9 2 5 11 9 13 6 -%! 14 7 10 15 15 8 11 16 5 11 9 13 6 6 12 10 14 7 -%! 15 8 11 16 11 4 7 12 2 8 6 10 3 2 8 6 10 3 -%! 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; -%! [mesh.n] = MSH2Mtopprop(mesh,'n'); -%! n =[ -%! NaN 10 5 NaN 4 NaN 10 NaN 14 15 16 17 18 13 NaN 3 1 2 -%! 5 6 7 8 3 NaN 18 15 NaN 2 14 16 NaN 9 10 11 12 13 -%! 17 18 16 5 1 2 3 4 NaN 7 NaN NaN 14 11 8 12 NaN 7]; -%! assert(mesh.n,n);
--- a/extra/msh/inst/MSH3Mgeomprop.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,177 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{varargout}]} = MSH3Mgeomprop(@var{mesh},[@var{string1},@var{string2},...]) -## -## Computes geometrical properties of the specified mesh -## -## Input: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{string1}, @var{string2},...: identifier of the property to compute. Possible choices are listed below. -## @itemize @bullet -## @item "wjacdet" : weighted determinant of jacobian trasformation to -## reference tetrahedron -## @item "shg": gradient of the P1 shape functions for BIM method -## @item "shp": value of the P1 shape functions for BIM method -## @end itemize -## @end itemize -## -## The output will contain the geometrical properties requested in the input in the same order specified in the function call -## @end deftypefn - - -function [varargout] = MSH3Mgeomprop(imesh,varargin) - - ## check if varargin is composed of strings as expected - if iscellstr(varargin) == 0 - warning("Unexpected input. See help for more information."); - print_usage; - endif - - ## extract tetrahedra node coordinates - x1 = imesh.p(1,imesh.t(1,:)); - y1 = imesh.p(2,imesh.t(1,:)); - z1 = imesh.p(3,imesh.t(1,:)); - x2 = imesh.p(1,imesh.t(2,:)); - y2 = imesh.p(2,imesh.t(2,:)); - z2 = imesh.p(3,imesh.t(2,:)); - x3 = imesh.p(1,imesh.t(3,:)); - y3 = imesh.p(2,imesh.t(3,:)); - z3 = imesh.p(3,imesh.t(3,:)); - x4 = imesh.p(1,imesh.t(4,:)); - y4 = imesh.p(2,imesh.t(4,:)); - z4 = imesh.p(3,imesh.t(4,:)); - - for nn = 1:length(varargin) - - request = varargin{nn}; - switch request - case "wjacdet" - b = wjacdet(x1,y1,z1,\ - x2,y2,z2,\ - x3,y3,z3,\ - x4,y4,z4); - varargout{nn} = b; - clear b - case "area" - tmp = wjacdet(x1,y1,z1,\ - x2,y2,z2,\ - x3,y3,z3,\ - x4,y4,z4); - b = sum(tmp,1); - varargout{nn} = b; - case "shg" - b = shg(x1,y1,z1,\ - x2,y2,z2,\ - x3,y3,z3,\ - x4,y4,z4); - varargout{nn} = b; - clear b - case "shp" - varargout{nn} = eye(4); - otherwise - warning("Unexpected value in passed string. Empty vector passed as output.") - varargout{nn} = []; - endswitch - - endfor - -endfunction - -function [b] = wjacdet(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) - ## COMPUTE WEIGHTED YACOBIAN DETERMINANT - - ## weight - weight = [1/4 1/4 1/4 1/4]'; - - Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); - Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); - Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); - - ## Determinant of the Jacobian of the - ## transformation from the base tetrahedron - ## to the tetrahedron K - detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; - ## Volume of the reference tetrahedron - Kkvolume = 1/6; - - b(:,:) = Kkvolume * weight * detJ; - -endfunction - -function [b] = shg(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) - ## COMPUTE GRADIENT OF SHAPE FUNCTIONS - - - Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); - Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); - Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); - - ## Determinant of the Jacobian of the - ## transformation from the base tetrahedron - ## to the tetrahedron K - detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; - - ## Shape function gradients follow - ## First index represents space direction - ## Second index represents the shape function - ## Third index represents the tetrahedron number - b(1,1,:) = (y2.*(z4-z3) + y3.*(z2-z4) + y4.*(z3-z2))./ detJ; - b(2,1,:) = (x2.*(z3-z4) + x3.*(z4-z2) + x4.*(z2-z3))./ detJ; - b(3,1,:) = (x2.*(y4-y3) + x3.*(y2-y4) + x4.*(y3-y2))./ detJ; - - b(1,2,:) = ( Nb2 ) ./ detJ; - b(2,2,:) = (x1.*(z4-z3) + x3.*(z1-z4) + x4.*(z3-z1)) ./ detJ; - b(3,2,:) = (x1.*(y3-y4) + x3.*(y4-y1) + x4.*(y1-y3)) ./ detJ; - - b(1,3,:) = ( Nb3 ) ./ detJ; - b(2,3,:) = (x1.*(z2-z4) + x2.*(z4-z1) + x4.*(z1-z2)) ./ detJ; - b(3,3,:) = (x1.*(y4-y2) + x2.*(y1-y4) + x4.*(y2-y1)) ./ detJ; - - b(1,4,:) = ( Nb4) ./ detJ; - b(2,4,:) = (x1.*(z3-z2) + x2.*(z1-z3) + x3.*(z2-z1)) ./ detJ; - b(3,4,:) = (x1.*(y2-y3) + x2.*(y3-y1) + x3.*(y1-y2)) ./ detJ; -endfunction - -%!shared mesh,wjacdet,shg,shp -% x = y = z = linspace(0,1,2); -% [mesh] = MSH3Mstructmesh(x,y,z,1,1:6) -% [wjacdet] = MSH3Mgeomprop(mesh,"wjacdet") -% [shg] = MSH3Mgeomprop(mesh,"shg") -% [shp] = MSH3Mgeomprop(mesh,"shp") -%!test -% assert(columns(mesh.t),columns(wjacdet)) -%!test -% assert(size(shg),[3 4 6]) -%!test -% assert(shp,eye(4)) -%!test -% fail(MSH3Mgeomprop(mesh,"samanafattababbudoiu"),"warning","Unexpected value in passed string. Empty vector passed as output.") \ No newline at end of file
--- a/extra/msh/inst/MSH3Mgmsh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,146 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH3Mgmsh(@var{geometry},@var{option},@var{value},...) -## -## Construct an unstructured 3D mesh making use of the free software gmsh. Give as output the PDE-tool like mesh structure. -## -## Input: -## @itemize @minus -## @item @var{geometry}: name of the ".geo" file describing the 2D geometry. Required by gmsh to start the meshing operation. -## @item @var{option}: option to be used by gmsh -## @item @var{value}: value of the option -## @end itemize -## For more information refer to gmsh manual, or gmsh site: -## -## http://www.geuz.org/gmsh/ -## -## Output: mesh basic structure, composed of the following fields -## @itemize @minus -## @item @var{p}: matrix with size 3 times number of mesh point. -## @itemize @bullet -## @item 1st row: x-coordinates of the points. -## @item 2nd row: y-coordinates of the points. -## @item 3rd row: z-coordinates of the points. -## @end itemize -## @item @var{e}: matrix with size 10 times number of mesh border edges. -## @itemize @bullet -## @item 1st row: p-matrix column number of the first edge-vertex. -## @item 2nd row: p-matrix column number of the second edge-vertex. -## @item 3rd row: p-matrix column number of the third edge-vertex. -## @item 4th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 5th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 6th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 7th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 8th row: number of the region to the right of the referred mesh edge. -## @item 9th row: number of the region to the left of the referred mesh edge. -## @item 10th row: number of the geometrical border upon which the referred mesh edge is lying on. -## @end itemize -## @item @var{t}: -## @itemize @bullet -## @item 1st row: p-matrix column number of the first tetrahedra vertex. -## @item 2nd row: p-matrix column number of the second tetrahedra vertex. -## @item 3rd row: p-matrix column number of the third tetrahedra vertex. -## @item 4th row: p-matrix column number of the fourth tetrahedra vertex. -## @item 5th row: number of the region upon which the referred trg is lying on. -## @end itemize -## @end itemize -## -## @seealso{MSH2Mgmsh,MSH2Mjoinstructm,MSH2Msubmesh} -## @end deftypefn - -function [mesh] = MSH3Mgmsh(geometry,varargin) - - ## Number of passed options - noptions = (nargin - 1) / 2; - - optstring = ""; - for ii = 1:noptions - option = varargin{2*(ii)-1}; - value = varargin{2*ii}; - if !ischar(value) - value = num2str(value); - endif - optstring = [optstring," -",option," ",value]; - endfor - - ## Generate mesh using Gmsh - printf("\n"); - printf("Generating mesh...\n"); - system(["gmsh -format msh -3 -o " geometry ".msh" optstring " " geometry ".geo"]); - - printf("Processing gmsh data...\n"); - ## Points - com_p = "awk '/\\$Nodes/,/\\$EndNodes/ {print $2, $3, $4 > ""msh_p.txt""}' "; - ## Surface edges - com_e = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""2"") print $7, $8, $9, $5 > ""msh_e.txt""}' "; - ## Tetrahedra - com_t = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""4"") print $7, $8, $9, $10, $5 > ""msh_t.txt""}' "; - ## Side edges - com_s = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""1"") print $7, $8, $5 > ""msh_s.txt""}' "; - - command = [com_p,geometry,".msh ; "]; - command = [command,com_e,geometry,".msh ; "]; - command = [command,com_t,geometry,".msh ; "]; - command = [command,com_s,geometry,".msh"]; - - system(command); - - ## Create PDE-tool like structure - printf("Creating PDE-tool like mesh...\n"); - ## Mesh-points - p = load("msh_p.txt")'; - ## Mesh side-edges - s = load("msh_s.txt")'; - ## Mesh surface-edges - tmp = load("msh_e.txt")'; - be = zeros(10,columns(tmp)); - be([1,2,3,10],:) = tmp; - ## Mesh tetrahedra - t = load("msh_t.txt")'; - - - ## Remove hanging nodes - printf("Check for hanging nodes...\n"); - nnodes = columns(p); - in_msh = intersect( 1:nnodes , t(1:4,:) ); - if length(in_msh) != nnodes - new_num(in_msh) = [1:length(in_msh)]; - t(1:4,:) = new_num(t(1:4,:)); - be(1:3,:) = new_num(be(1:3,:)); - p = p(:,in_msh); - endif - - mesh = struct("p",p,"s",s,"e",be,"t",t); - - printf("Deleting temporary files...\n"); - system(["rm -f msh_p.txt msh_e.txt msh_t.txt msh_s.txt *.msh"]); - -endfunction \ No newline at end of file
--- a/extra/msh/inst/MSH3Mjoinstructm.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,157 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH3Mjoinstructm(@var{mesh1},@var{mesh2},@var{s1},@var{s2}) -## -## Join two structured meshes (created by MSH3Mstructmesh) into one -## mesh structure variable. -## -## Input: -## @itemize @minus -## @item @var{mesh1}, @var{mesh2}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{s1}, @var{s2}: number of the corresponding geometrical border edge for respectively mesh1 and mesh2. -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @end itemize -## -## WARNING: on the common edge the two meshes must share the same vertexes. -## -## @seealso{MSH3Mstructmesh,MSH3Mgmsh,MSH3Msubmesh} -## @end deftypefn - -function mesh = MSH3Mjoinstructm(mesh1,mesh2,s1,s2) - - ## outside world must always be on the same side of the - ## boundary of mesh1 - [mesh1.e(8:9,:),I] = sort(mesh1.e(8:9,:)); - - ## IF THE REGIONS ARE INVERTED THE VERTEX ORDER SHOULD ALSO BE - ## INVERTED!! - - ## get interface nodes - intfcnodes1 = MSH3Mnodesonfaces(mesh1,s1)'; - intfcnodes2 = MSH3Mnodesonfaces(mesh2,s2)'; - - ## sort interface nodes by position - [tmp,I] = sort(mesh1.p(1,intfcnodes1)); - intfcnodes1 = intfcnodes1(I); - [tmp,I] = sort(mesh1.p(2,intfcnodes1)); - intfcnodes1 = intfcnodes1(I); - [tmp,I] = sort(mesh1.p(3,intfcnodes1)); - intfcnodes1 = intfcnodes1(I); - - [tmp,I] = sort(mesh2.p(1,intfcnodes2)); - intfcnodes2 = intfcnodes2(I); - [tmp,I] = sort(mesh2.p(2,intfcnodes2)); - intfcnodes2 = intfcnodes2(I); - [tmp,I] = sort(mesh2.p(3,intfcnodes2)); - intfcnodes2 = intfcnodes2(I); - - ## delete redundant boundary faces - ## but first remeber what region - ## they were connected to - for is = 1:length(s2) - ii = find( mesh2.e(10,:)==s2(is) ); - adreg(is,:) = unique(mesh2.e(9,ii)); - endfor - - for is = 1:length(s2) - mesh2.e(:,find( mesh2.e(10,:)==s2(is) )) = []; - endfor - - ## change face numbers - idx = []; - consecutives = []; - idx = unique(mesh2.e(10,:)); - consecutives (idx) = [1:length(idx)] + max(mesh1.e(10,:)); - mesh2.e(10,:) = consecutives(mesh2.e(10,:)); - - ## change node indices in connectivity matrix - ## and edge list - idx = []; - consecutives = []; - idx = 1:size(mesh2.p,2); - offint = setdiff(idx,intfcnodes2); - consecutives (offint) = [1:length(offint)]+size(mesh1.p,2); - - consecutives (intfcnodes2) = intfcnodes1; - mesh2.e(1:3,:) = consecutives(mesh2.e(1:3,:)); - mesh2.t(1:4,:) = consecutives(mesh2.t(1:4,:)); - - ## delete redundant points - mesh2.p(:,intfcnodes2) = []; - - ## set region numbers - regions = unique(mesh1.t(5,:)); - newregions(regions) = 1:length(regions); - mesh1.t(5,:) = newregions(mesh1.t(5,:)); - - ## set region numbers - regions = unique(mesh2.t(5,:)); - newregions(regions) = [1:length(regions)]+max(mesh1.t(5,:)); - mesh2.t(5,:) = newregions(mesh2.t(5,:)); - - ## set adjacent region numbers in face structure 2 - [i,j] = find(mesh2.e(8:9,:)); - i += 7; - - mesh2.e(i,j) = newregions(mesh2.e(i,j)); - - ## set adjacent region numbers in edge structure 1 - for is = 1:length(s1) - ii = find( mesh1.e(10,:)==s1(is) ); - mesh1.e(8,ii) = newregions(regions(adreg(is,:))); - endfor - - ## build nbew mesh structure - mesh.p = [mesh1.p mesh2.p]; - mesh.e = [mesh1.e mesh2.e]; - mesh.t = [mesh1.t mesh2.t]; - -endfunction - -%!shared mesh1,mesh2,jmesh -% x = y = z = linspace(0,1,2); -% x2 = linspace(1,2,2); -% [mesh1] = MSH3Mstructmesh(x,y,z,1,1:6); -% [mesh2] = MSH3Mstructmesh(x2,y,z,3,1:6); -% [jmesh] = MSH3Mjoinstructm(mesh1,mesh2,2,1); -%!test -% assert(columns(jmesh.p),12) -%!test -% tmp = sort(unique(jmesh.e(10,:))); -% assert(tmp,1:11) -%!test -% assert(columns(jmesh.t),columns(mesh1.t)+columns(mesh2.t)) -%!test -% assert(unique(jmesh.e(8:9,:)),0:2)
--- a/extra/msh/inst/MSH3Mnodesonfaces.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{nodelist}]} =@ -## MSH3Mnodesonfaces(@var{mesh}, @var{facelist}) -## -## Returns a list of the nodes lying on the faces @var{facelist} of -## the mesh @var{mesh}. -## -## Input: -## @itemize @minus -## @item @var{mesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{facelist}: row vector containing the number of the faces (numbering referred to mesh.e(10,:)). -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{nodelist}: list of the nodes that lies on the specified faces. -## @end itemize -## -## @seealso{MSH3Mgeomprop,MSH3Mtopprop} -## @end deftypefn - -function [nodelist] = MSH3Mnodesonfaces(mesh,facelist); - - facefaces = []; - - for ii=1:length(facelist) - facefaces = [facefaces,find(mesh.e(10,:)==facelist(ii))]; - endfor - - facenodes = mesh.e(1:3,facefaces); - nodelist = unique(facenodes(:)); - -endfunction - -%!shared x,y,z,mesh -% x = y = z = linspace(0,1,2); -% [mesh] = MSH3Mstructmesh(x,y,z,1,1:6); -%!test -% nodelist = MSH3Mnodesonfaces(mesh,1); -% assert(nodelist,[1 2 5 6]') -%!test -% nodelist = MSH3Mnodesonfaces(mesh,2); -% assert(nodelist,[3 4 7 8]') -%!test -% nodelist = MSH3Mnodesonfaces(mesh,3); -% assert(nodelist,[1 3 5 7]') -%!test -% nodelist = MSH3Mnodesonfaces(mesh,[1 2 3]); -% assert(nodelist,[1:8]') \ No newline at end of file
--- a/extra/msh/inst/MSH3Mstructmesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,237 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{mesh}]} = MSH3Mstructmesh(@var{x},@var{y},@var{z},@var{region},@var{sides}) -## -## Constructs a structured tetrahedral 3D mesh on a parallelepipedal domain, -## and returns a PDEtool-like mesh structure. -## -## Input: -## @itemize @minus -## @item @var{x}: vector representing the 1D meshing of the side -## parallel to x axis. -## @item @var{y}: vector representing the 1D meshing of the side -## parallel to y axis. -## @item @var{z}: vector representing the 1D meshing of the side -## parallel to z axis. -## @item @var{region}: number assigned to the meshed region. -## @item @var{sides}: row vector containing the six numbers assigned to the geometrical surface-edges. -## @end itemize -## -## Output: mesh basic structure, composed of the following fields -## @itemize @minus -## @item @var{p}: matrix with size 3 times number of mesh point. -## @itemize @bullet -## @item 1st row: x-coordinates of the points. -## @item 2nd row: y-coordinates of the points. -## @item 3rd row: z-coordinates of the points. -## @end itemize -## @item @var{e}: matrix with size 10 times number of mesh border edges. -## @itemize @bullet -## @item 1st row: p-matrix column number of the first edge-vertex. -## @item 2nd row: p-matrix column number of the second edge-vertex. -## @item 3rd row: p-matrix column number of the third edge-vertex. -## @item 4th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 5th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 6th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 7th row: not initialized, only for compatibility with standard PDE-tool like mesh. -## @item 8th row: number of the region to the right of the referred mesh -## edge. -## @item 9th row: number of the region to the left of the referred mesh edge. -## @item 10th row: number of the geometrical border upon which the referred mesh edge is lying on. -## @end itemize -## @item @var{t}: -## @itemize @bullet -## @item 1st row: p-matrix column number of the first tetrahedra vertex. -## @item 2nd row: p-matrix column number of the second tetrahedra vertex. -## @item 3rd row: p-matrix column number of the third tetrahedra vertex. -## @item 4th row: p-matrix column number of the fourth tetrahedra vertex. -## @item 5th row: number of the region upon which the referred trg is lying on. -## @end itemize -## @end itemize -## -## @seealso{MSH2Mgmsh,MSH2Mjoinstructm,MSH2Msubmesh} -## @end deftypefn - -function [mesh] = MSH3Mstructmesh(x,y,z,region,sides) - - ## check for correct input - scalar = ( isscalar(x) || isscalar(y) || isscalar(z) ); - matrix = ( min(size(x)) + min(size(y)) + min(size(z)) != 3 ); - if scalar - warning("x, y, z cannot be scalar numbers!") - print_usage; - endif - if matrix - warning("x, y, z cannot be matrices!") - print_usage; - endif - - ## sort point coordinates - x = sort(x); - y = sort(y); - z = sort(z); - ## compute # of points in each direction - nx = length(x); - ny = length(y); - nz = length(z); - - ## generate verticeces - [XX,YY,ZZ] = meshgrid(x,y,z); - p = [XX(:),YY(:),ZZ(:)]'; - - iiv (ny,nx,nz)=0; - iiv(:)=1:nx*ny*nz; - iiv(end,:,:)=[]; - iiv(:,end,:)=[]; - iiv(:,:,end)=[]; - iiv=iiv(:)'; - - ## generate connections: - - ## bottom faces - n1 = iiv; - n2 = iiv + 1; - n3 = iiv + ny; - n4 = iiv + ny + 1; - - ## top faces - N1 = iiv + nx * ny; - N2 = N1 + 1; - N3 = N1 + ny; - N4 = N3 + 1; - - t = [... - [n1; n3; n2; N2],... - [N1; N2; N3; n3],... - [N1; N2; n3; n1],... - [N2; n3; n2; n4],... - [N3; n3; N2; N4],... - [N4; n3; N2; n4],... - ]; - - ## generate boundary face list: - - ## left - T = t; - T(:) = p(1,t)'==x(1); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e1(:,jj) = t(order(:,ii(jj)),ii(jj)); - endfor - e1(10,:) = sides(1); - - ## right - T(:) = p(1,t)'==x(end); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e2(:,jj) = t(order(:,ii(jj)),ii(jj)); - end - e2(10,:) = sides(2); - - ## front - T(:) = p(2,t)'==y(1); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e3(:,jj) = t(order(:,ii(jj)),ii(jj)); - endfor - e3(10,:) = sides(3); - - ## back - T(:) = p(2,t)'==y(end); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e4(:,jj) = t(order(:,ii(jj)),ii(jj)); - endfor - e4(10,:) = sides(4); - - ## bottom - T = t; - T(:) = p(3,t)'==z(1); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e5(:,jj) = t(order(:,ii(jj)),ii(jj)); - endfor - e5(10,:) = sides(5); - - ## top - T = t; - T(:) = p(3,t)'==z(end); - [ignore,order] = sort(T,1); - ii = (find(sum(T,1)==3)); - order(1,:) = []; - for jj=1:length(ii) - e6(:,jj) = t(order(:,ii(jj)),ii(jj)); - endfor - e6(10,:) = sides(6); - - - mesh.e = [e1,e2,e3,e4,e5,e6]; - mesh.t = t; - mesh.e (9,:) = region; - mesh.t (5,:) = region; - mesh.p = p; - -endfunction - -%!test -% x = y = z = linspace(0,1,2) -% [mesh] = MSH3Mstructmesh(x,y,z,1,1:6) -% assert = (columns(mesh.p),8) -% assert = (columns(mesh.t),6) -% assert = (columns(mesh.e),12) -%!test -% x = y = z = linspace(0,1,3) -% [mesh] = MSH3Mstructmesh(x,y,z,1,1:6) -% assert = (columns(mesh.p),27) -% assert = (columns(mesh.t),48) -% assert = (columns(mesh.e),48) -%!test -% x = y = z = linspace(0,1,4) -% [mesh] = MSH3Mstructmesh(x,y,z,1,1:6) -% assert = (columns(mesh.p),64) -% assert = (columns(mesh.t),162) -% assert = (columns(mesh.e),108) -%!test -% x = y = z = linspace(0,1,1) -% fail([mesh] = MSH3Mstructmesh(x,y,z,1,1:6),"warning","x, y, z cannot be scalar numbers!") -%!test -% x = y = z = eye(2) -% fail([mesh] = MSH3Mstructmesh(x,y,z,1,1:6),"warning","x, y, z cannot be matrices!") \ No newline at end of file
--- a/extra/msh/inst/MSH3Msubmesh.m Mon Feb 01 03:57:21 2010 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,95 +0,0 @@ -## Copyright (C) 2007,2008 Carlo de Falco, Massimiliano Culpo -## -## This file is part of -## -## MSH - Meshing Software Package for Octave -## -## 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, see <http://www.gnu.org/licenses/>. -## -## -## AUTHORS: -## Carlo de Falco <cdf _AT_ users.sourceforge.net> -## -## Culpo Massimiliano -## Bergische Universitaet Wuppertal -## Fachbereich C - Mathematik und Naturwissenschaften -## Arbeitsgruppe fuer Angewandte MathematD-42119 Wuppertal Gaussstr. 20 -## D-42119 Wuppertal, Germany - -## -*- texinfo -*- -## @deftypefn {Function File} {[@var{omesh},@var{nodelist},@var{elementlist}]} = MSH3Msubmesh(@var{imesh},@var{intrfc},@var{sdl}) -## -## Gives as output a specified submesh, and the lists of nodes and element that are mantained. -## -## Input: -## @itemize @minus -## @item @var{imesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{intrfc}: row vector containing the number of the internal -## interface sides (numbering referred to mesh.e(10,:)) that should be -## mantained. Could be empty. -## @item @var{sdl} (subdomain list): row vector containing the list of the subdomain that are going to be extracted. -## @end itemize -## -## Output: -## @itemize @minus -## @item @var{omesh}: standard PDEtool-like mesh, with field "p", "e", "t". -## @item @var{nodelist}: list of the node of the original mesh that are present in the restricted one. -## @item @var{elementlist}: list of the element of the original mesh that are present in the restricted one. -## @end itemize -## -## @seealso{MSH3Mstructmesh,MSH3Mjoinstructm,MSH3Mgmsh} -## @end deftypefn - -function [omesh,nodelist,elementlist] = MSH3Msubmesh(imesh,intrfc,sdl) - - ## build element list - elementlist=[]; - for ir = 1:length(sdl) - elementlist = [ elementlist find(imesh.t(5,:)==sdl(ir)) ]; - endfor - - ## build nodelist - nodelist = reshape(imesh.t(1:4,elementlist),1,[]); - nodelist = unique(nodelist); - - ## extract submesh - omesh.p = imesh.p (:,nodelist); - indx(nodelist) = 1:length (nodelist); - omesh.t = imesh.t (:,elementlist); - omesh.t(1:4,:) = indx(omesh.t(1:4,:)); - - omesh.e = []; - for ifac = 1:size(imesh.e,2) - if (length(intersect(imesh.e(1:3,ifac),nodelist) )== 3) - omesh.e = [omesh.e imesh.e(:,ifac)]; - endif - endfor - - omesh.e(1:3,:) = indx(omesh.e(1:3,:)); - -endfunction - -%!shared mesh1,mesh2,jmesh,exmesh,nodelist,elemlist -% x = y = z = linspace(0,1,2); -% x2 = linspace(1,2,2); -% [mesh1] = MSH3Mstructmesh(x,y,z,1,1:6); -% [mesh2] = MSH3Mstructmesh(x2,y,z,1,1:6); -% [jmesh] = MSH3Mjoinstructm(mesh1,mesh2,2,1); -% [exmesh,nodelist,elemlist] = MSH3Msubmesh(jmesh,2,1); -%!test -% assert(size(exmesh.p),size(mesh1.p)) -%!test -% assert(size(exmesh.t),size(mesh1.t)) -%!test -% assert(size(exmesh.e),size(mesh1.e)) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_displacement_smoothing.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,155 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{Ax},@var{Ay}]} = @ +## msh2m_displacement_smoothing(@var{msh},@var{k}) +## +## Displace the boundary of a 2D mesh setting a spring with force/length +## constant @var{k} along each edge and enforcing 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 to make the system +## solvable, e.g.: +## +## @example +## msh = msh2m_structured_mesh(linspace(0,1,10),@ +## linspace(0,1,10),@ +## 1,1:4,"left"); +## +## dnodes = msh2m_nodes_on_sides(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); +## dxtot = dytot = -.5*sin(xd.*yd*pi/2); +## Nsteps = 10; +## +## for ii = 1:Nsteps +## dx(dnodes) = dxtot; +## dy(dnodes) = dytot; +## [Ax,Ay] = msh2m_displacement_smoothing(msh,1); +## dx(varnodes) = Ax(varnodes,varnodes) \ ... +## (-Ax(varnodes,dnodes)*dx(dnodes)); +## dy(varnodes) = Ay(varnodes,varnodes) \ ... +## (-Ay(varnodes,dnodes)*dy(dnodes)); +## msh.p += [ dx'/Nsteps; dy'/Nsteps ] ; +## triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); +## pause(.01) +## endfor +## @end example +## +## @seealso{msh2m_jiggle_mesh} +## +## @end deftypefn + +function [Ax,Ay] = msh2m_displacement_smoothing(msh, k) + + ## Check input + if nargin != 2 # Number of input parameters + error("msh2m_displacement_smoothing: wrong number of input parameters."); + elseif !(isstruct(msh) && isfield(msh,"p") && + isfield(msh,"t") && isfield(msh,"e")) + error("msh2m_displacement_smoothing: first input is not a valid mesh structure."); + elseif !isscalar(k) + error("msh2m_displacement_smoothing: k must be a valid scalar"); + endif + + ## Construct matrices + x = msh.p(1,:); + y = msh.p(2,:); + + dx2 = (x(msh.t([1 2 3],:))-x(msh.t([2 3 1],:))).^2; + dy2 = (y(msh.t([1 2 3],:))-y(msh.t([2 3 1],:))).^2; + + l2 = dx2 + dy2; + + Ax = spalloc(length(x),length(x),1); + Ay = spalloc(length(x),length(x),1); + + ax = zeros(3,3,columns(msh.t)); + ay = zeros(3,3,columns(msh.t)); + + for inode=1:3 + for jnode=1:3 + ginode(inode,jnode,:)=msh.t(inode,:); + gjnode(inode,jnode,:)=msh.t(jnode,:); + endfor + endfor + + for ii=1:3 + for jj=ii+1:3 + + ax(ii,jj,:) = ax(jj,ii,:) = reshape(-k * dx2(ii,:)./l2(ii,:),1,1,[]); + ay(ii,jj,:) = ay(jj,ii,:) = reshape(-k * dy2(ii,:)./l2(ii,:),1,1,[]); + + 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(:)); + +endfunction + +%!demo +%! msh = msh2m_structured_mesh(linspace(0,1,10), +%! linspace(0,1,10), +%! 1,1:4,"left"); +%! dnodes = msh2m_nodes_on_sides(msh,1:4); +%! varnodes = setdiff([1:columns(msh.p)],dnodes); +%! +%! xd = msh.p(1,dnodes)'; +%! yd = msh.p(2,dnodes)'; +%! +%! dy = zeros(columns(msh.p),1); +%! dx = dy; +%! +%! dxtot = -.5*sin(xd.*yd*pi/2); +%! dytot = -.5*sin(xd.*yd*pi/2); +%! +%! Nsteps = 5; +%! for ii=1:Nsteps +%! +%! dx(dnodes) = dxtot; +%! dy(dnodes) = dytot; +%! +%! [Ax,Ay] = msh2m_displacement_smoothing(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; +%! +%! if mod(ii,2)==0 +%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); +%! pause(.01) +%! endif +%! endfor
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_equalize_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,132 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh2m_equalize_mesh(@var{mesh}) +## +## Apply a baricentric regularization to equalize the size of triangle +## edges, i.e. move each node to the center of mass of the patch of +## triangles to which it belongs. +## +## May be useful when distorting a mesh. +## Type @code{demo msh2m_equalize_mesh} to see some examples. +## +## @seealso{msh2m_displacement_smoothing} +## +## @end deftypefn + +function [msh] = msh2m_equalize_mesh(msh) + + ## Check input + if nargin != 1 # Number of input parameters + error("msh2m_equalize_mesh: wrong number of input parameters."); + elseif !(isstruct(msh) && isfield(msh,"p") && + isfield(msh,"t") && isfield(msh,"e")) + error("msh2m_equalize_mesh: first input is not a valid mesh structure."); + endif + + ## Apply regularization + nel= columns(msh.t); + + x = msh.p(1,:)'; + y = msh.p(2,:)'; + + dnodes = unique(msh.e(1:2,:)(:)); + varnodes = setdiff([1:columns(msh.p)],dnodes); + + Ax = spalloc(length(x),length(x),1); + Ay = spalloc(length(x),length(x),1); + + ax = zeros(3,3,nel); + ay = zeros(3,3,nel); + + 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,:); + endfor + endfor + + for ii=1:3 + 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(:)); + + x(varnodes) = Ax(varnodes,varnodes) \ (-Ax(varnodes,dnodes)*x(dnodes)); + y(varnodes) = Ay(varnodes,varnodes) \ (-Ay(varnodes,dnodes)*y(dnodes)); + msh.p(1,:) = x'; + msh.p(2,:) = y'; + +endfunction + +%!demo +%! ### equalize a structured mesh without moving boundary nodes +%! msh = msh2m_structured_mesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"random"); +%! dnodes = msh2m_nodes_on_sides(msh,1:4); +%! varnodes = setdiff([1:columns(msh.p)],dnodes); +%! x = msh.p(1,:)'; +%! y = msh.p(2,:)'; +%! msh = msh2m_equalize_mesh(msh); +%! 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 = msh2m_structured_mesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"random"); +%! dnodes = msh2m_nodes_on_sides(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 = 10; +%! for ii=1:Nsteps +%! dx(dnodes) = dxtot; +%! dy(dnodes) = dytot; +%! [Ax,Ay] = msh2m_displacement_smoothing(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,:)'; +%! msh = msh2m_equalize_mesh(msh); +%! hold on;triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)');hold off +%! pause(.5) +%! endfor \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_geometrical_properties.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,450 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{varargout}]} = @ +## msh2m_geometrical_properties(@var{mesh},[@var{string1},@var{string2},...]) +## +## Compute @var{mesh} geometrical properties identified by input strings. +## +## Valid properties are: +## @itemize @bullet +## @item @code{"bar"}: return a matrix with size 2 times the number of mesh +## elements containing the center of mass coordinates. +## @item @code{"cir"}: return a matrix with size 2 times the number of +## mesh elements containing the circumcenter coordinates. +## @item @code{"emidp"}: return a matrix with size 2 times the number of +## side edges containing their midpoint coordinates. +## @item @code{"slength"}: return a matrix with size 3 times the number +## of mesh elements containing the length of each element side. +## @item @code{"cdist"}: return a matrix of size 3 times the number of +## mesh elements containing the distance among circumcenters of +## neighbouring elements. If the corresponding side lies on the edge, +## the distance between circumcenter and border edge is returned in the +## matrix. +## @item @code{"wjacdet"}: return the weigthed Jacobian determinant used +## for the numerical integration with trapezoidal rule over an element. +## @item @code{"shg"}: return a matrix of size 3 times the number of +## elements matrix containing the gradient of P1 shape functions. +## @item @code{"area"}: return a row vector containing the area of every +## element. +## @item @code{"midedge"}: return a multi-dimensional array with size 2 +## times 3 times the number of elements containing the coordinates of +## the midpoint of every edge. +## @end itemize +## +## The output will contain the geometrical properties requested in the +## input in the same order specified in the function call. +## +## If an unexpected string is given as input, an empty vector is +## returned in output. +## +## @seealso{msh2m_topological_properties, msh3m_geometrical_properties} +## @end deftypefn + +function [varargout] = msh2m_geometrical_properties(mesh,varargin) + + ## Check input + if nargin < 2 # Number of input parameters + error("msh2m_geometrical_properties: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield(mesh,"t") && isfield(mesh,"e")) + error("msh2m_geometrical_properties: first input is not a valid mesh structure."); + elseif !iscellstr(varargin) + error("msh2m_geometrical_properties: only string value admitted for properties."); + endif + + ## Compute properties + p = mesh.p; + e = mesh.e; + t = mesh.t; + nelem = columns(t); # Number of elements in the mesh + + [k,j,w] = coeflines(p,t,nelem); # Edges coefficients + + for nn = 1:length(varargin) + + request = varargin{nn}; + switch request + + case "bar" # Center of mass coordinates + if isfield(mesh,"bar") + varargout{nn} = mesh.bar; + else + [b] = coordinates(p,t,nelem,j,w,k,"bar"); + varargout{nn} = b; + clear b; + endif + + case "cir" # Circum center coordinates + if isfield(mesh,"cir") + varargout{nn} = mesh.cir; + else + [b] = coordinates(p,t,nelem,j,w,k,"cir"); + varargout{nn} = b; + clear b; + endif + + case "emidp" # Boundary edges midpoint coordinates + if isfield(mesh,"emidp") + varargout{nn} = mesh.emidp; + else + [b] = midpoint(p,e); + varargout{nn} = b; + clear b; + endif + + case "slength" # Length of every side + if isfield(mesh,"slength") + varargout{nn} = mesh.slength; + else + [b] = sidelength(p,t,nelem); + varargout{nn} = b; + clear b; + endif + + case "cdist" # Distance among circumcenters of neighbouring elements + if isfield(mesh,"cdist") + varargout{nn} = mesh.cdist; + else + + if isfield(mesh,"cir") + cir = mesh.cir; + else + [cir] = coordinates(p,t,nelem,j,w,k,"cir"); + endif + + if isfield(mesh,"n") + n = mesh.n; + else + [n] = msh2m_topological_properties(mesh,"n"); + endif + + [b] = distance(cir,n,nelem); + [semib] = semidistance(cir,nelem,j,w,k); + border = isnan(n); + [index1] = find( border(1,:) ); + [index2] = find( border(2,:) ); + [index3] = find( border(3,:) ); + b(1,index1) = semib(1,index1); + b(2,index2) = semib(2,index2); + b(3,index3) = semib(3,index3); + varargout{nn} = b; + clear b semib index1 index2 index3 border; + endif + + case "wjacdet" # Weighted Jacobian determinant + if isfield(mesh,"wjacdet") + varargout{nn} = mesh.wjacdet; + else + [b] = computearea(p,e,t,"wjac"); + varargout{nn} = b; + clear b + endif + + case "area" # Area of the elements + if isfield(mesh,"area") + varargout{nn} = mesh.area; + else + [b] = computearea(p,e,t,"area"); + varargout{nn} = b; + clear b + endif + + case "shg" # Gradient of hat functions + if isfield(mesh,"shg") + varargout{nn} = mesh.shg; + else + [b] = shapegrad(p,t); + varargout{nn} = b; + clear b + endif + + case "midedge" # Mid-edge coordinates + if isfield(mesh,"midedge") + varargout{nn} = mesh.midedge; + else + [b] = midedge(p,t,nelem); + varargout{nn} = b; + clear b; + endif + + otherwise + warning("msh2m_geometrical_properties: unexpected value in property string. Empty vector passed as output.") + varargout{nn} = []; + endswitch + + endfor + +endfunction + +function [k,j,w] = coeflines(p,t,nelem) + + ## Edges are described by the analytical expression: + ## + ## k*x + j*y + w = 0 + ## + ## Coefficients k,j,w are stored in matrixes + + ## i-th edge list, i =1,2,3 + s1 = sort(t(2:3,:),1); + s2 = sort(t([3,1],:),1); + s3 = sort(t(1:2,:),1); + ## Initialization of the matrix data-structure + k = ones(3,nelem); + j = ones(3,nelem); + w = ones(3,nelem); + ## Searching for lines parallel to x axis + [i1] = find( (p(2,s1(2,:)) - p(2,s1(1,:))) != 0 ); + noti1 = setdiff([1:nelem], i1); + [i2] = find( (p(2,s2(2,:)) - p(2,s2(1,:))) != 0 ); + noti2 = setdiff([1:nelem], i2); + [i3] = find( (p(2,s3(2,:)) - p(2,s3(1,:))) != 0 ); + noti3 = setdiff([1:nelem], i3); + ## Computation of the coefficients + ## Edge 1 + j(1,i1) = (p(1,s1(1,i1)) - p(1,s1(2,i1)))./(p(2,s1(2,i1)) - p(2,s1(1,i1))); + w(1,i1) = -(p(1,s1(1,i1)) + p(2,s1(1,i1)).*j(1,i1)); + k(1,noti1) = 0; + j(1,noti1) = 1; + w(1,noti1) = - p(2,s1(1,noti1)); + ## Edge 2 + j(2,i2) = (p(1,s2(1,i2)) - p(1,s2(2,i2)))./(p(2,s2(2,i2)) - p(2,s2(1,i2))); + w(2,i2) = -(p(1,s2(1,i2)) + p(2,s2(1,i2)).*j(2,i2)); + k(2,noti2) = 0; + j(2,noti2) = 1; + w(2,noti2) = - p(2,s2(1,noti2)); + ## Edge 3 + j(3,i3) = (p(1,s3(1,i3)) - p(1,s3(2,i3)))./(p(2,s3(2,i3)) - p(2,s3(1,i3))); + w(3,i3) = -(p(1,s3(1,i3)) + p(2,s3(1,i3)).*j(3,i3)); + k(3,noti3) = 0; + j(3,noti3) = 1; + w(3,noti3) = - p(2,s3(1,noti3)); +endfunction + +function [b] = coordinates(p,t,nelem,j,w,k,string) + + ## Compute the coordinates of the geometrical entity specified by string + + ## Initialization of the output vectors + b = zeros(2, nelem); + switch string + case "bar" + b(1,:) = ( p(1,t(1,:)) + p(1,t(2,:)) + p(1,t(3,:)) )/3; + b(2,:) = ( p(2,t(1,:)) + p(2,t(2,:)) + p(2,t(3,:)) )/3; + case "cir" + ## Computation of the midpoint of the first two edges + mid1 = zeros(2,nelem); + mid2 = zeros(2,nelem); + ## X coordinate + mid1(1,:) = ( p(1,t(2,:)) + p(1,t(3,:)) )/2; + mid2(1,:) = ( p(1,t(3,:)) + p(1,t(1,:)) )/2; + ## Y coordinate + mid1(2,:) = ( p(2,t(2,:)) + p(2,t(3,:)) )/2; + mid2(2,:) = ( p(2,t(3,:)) + p(2,t(1,:)) )/2; + ## Computation of the intersect between axis 1 and axis 2 + ## Searching for element with edge 1 parallel to x-axes + [parx] = find( j(1,:) == 0); + [notparx] = setdiff(1:nelem, parx); + coefy = zeros(1,nelem); + ## If it is not parallel + coefy(notparx) = ((j(2,notparx)./j(1,notparx)).*k(1,notparx) - k(2,notparx)).^(-1); + b(2,notparx) = coefy(1,notparx).*( j(2,notparx).*mid2(1,notparx) - k(2,notparx).*mid2(2,notparx) + k(1,notparx)./j(1,notparx).*j(2,notparx).*mid1(2,notparx) - j(2,notparx).*mid1(1,notparx) ); + b(1,notparx) =( k(1,notparx).*b(2,notparx) + j(1,notparx).*mid1(1,notparx) - k(1,notparx).*mid1(2,notparx) )./j(1,notparx); + ## If it is parallel + b(2,parx) = mid1(2,parx); + b(1,parx) = k(2,parx)./j(2,parx).*( b(2,parx) - mid2(2,parx) ) + mid2(1,parx); + endswitch +endfunction + +function [b] = midpoint(p,e) + + ## Compute the coordinates of the midpoint on the boundary edges + + b = zeros(2,columns(e)); + b(1,:) = ( p(1,e(1,:)) + p(1,e(2,:)) )./2; + b(2,:) = ( p(2,e(1,:)) + p(2,e(2,:)) )./2; +endfunction + +function [l] = sidelength(p,t,nelem) + + ## Compute the length of every side + + l = zeros(3, nelem); + ## i-th edge list, i =1,2,3 + s1 = sort(t(2:3,:),1); + s2 = sort(t([3,1],:),1); + s3 = sort(t(1:2,:),1); + ## First side length + l(1,:) = sqrt( (p(1,s1(1,:)) - p(1,s1(2,:))).^2 + (p(2,s1(1,:)) - p(2,s1(2,:))).^2 ); + ## Second side length + l(2,:) = sqrt( (p(1,s2(1,:)) - p(1,s2(2,:))).^2 + (p(2,s2(1,:)) - p(2,s2(2,:))).^2 ); + ## Third side length + l(3,:) = sqrt( (p(1,s3(1,:)) - p(1,s3(2,:))).^2 + (p(2,s3(1,:)) - p(2,s3(2,:))).^2 ); +endfunction + +function [d] = semidistance(b,nelem,j,w,k) + + ## Compute the distance to the sides of the nodes with coordinates b + ## The edges are described by the analytical expression: + ## + ## k*x + j*y + w = 0 + ## + ## The coefficients k,j,w are stored in matrixes + + ## Initialization of the distance output vector + d = zeros(3, nelem); + ## Computation of the distances from the geometrical entity to the edges + d(1,:) = abs( k(1,:).*b(1,:) + j(1,:).*b(2,:) + w(1,:))./(sqrt( k(1,:).^2 + j(1,:).^2)); + d(2,:) = abs( k(2,:).*b(1,:) + j(2,:).*b(2,:) + w(2,:))./(sqrt( k(2,:).^2 + j(2,:).^2)); + d(3,:) = abs( k(3,:).*b(1,:) + j(3,:).*b(2,:) + w(3,:))./(sqrt( k(3,:).^2 + j(3,:).^2)); +endfunction + +function [d] = distance(b,n,nelem) + + ## Compute the distance between two neighbouring entities + + ## Initialization of the distance output vector + d = NaN(3, nelem); + ## Trg not on the geometrical border + border = isnan(n); + [index1] = find( border(1,:) == 0 ); + [index2] = find( border(2,:) == 0 ); + [index3] = find( border(3,:) == 0 ); + ## Computation of the distances between two neighboring geometrical entities + d(1,index1) = sqrt( ( b(1,index1) - b(1,n(1,index1)) ).^2 + ( b(2,index1) - b(2,n(1,index1)) ).^2 ); + d(2,index2) = sqrt( ( b(1,index2) - b(1,n(2,index2)) ).^2 + ( b(2,index2) - b(2,n(2,index2)) ).^2 ); + d(3,index3) = sqrt( ( b(1,index3) - b(1,n(3,index3)) ).^2 + ( b(2,index3) - b(2,n(3,index3)) ).^2 ); +endfunction + +function [b] = computearea(p,e,t,string) + + ## Compute the area of every element in the mesh + + weight = [1/3 1/3 1/3]; + areakk = 1/2; + Nelements = columns(t); + + jac([1,2],:) = [p(1,t(2,:))-p(1,t(1,:)); + p(1,t(3,:))-p(1,t(1,:))]; + jac([3,4],:) = [p(2,t(2,:))-p(2,t(1,:)); + p(2,t(3,:))-p(2,t(1,:))]; + jacdet = jac(1,:).*jac(4,:)-jac(2,:).*jac(3,:); + + degen=find(jacdet <= 0); + if ~isempty(degen) + ## XXX FIXME: there should be a -verbose option to allow to see this + ## fprintf(1,"invalid mesh element: %d fixing...\n",degen); + t(1:3,degen) = t([2,1,3],degen); + jac([1,2],degen) = [p(1,t(2,degen))-p(1,t(1,degen)); + p(1,t(3,degen))-p(1,t(1,degen))]; + jac([3,4],degen) = [p(2,t(2,degen))-p(2,t(1,degen)); + p(2,t(3,degen))-p(2,t(1,degen))]; + jacdet(degen) = jac(1,degen).*jac(4,degen)-jac(2,degen).*jac(3,degen); + endif + + for inode = 1:3 + wjacdet(inode,:) = areakk .* jacdet .* weight(inode); + endfor + + if string == "wjac" + b = wjacdet(); + elseif string == "area" + b = sum(wjacdet)'; + endif + +endfunction + +function [d] = midedge(p,t,nelem) + ## Compute the midpoint coordinates for every edge + s1 = t(2:3,:); s2 = t([3,1],:); s3 = t(1:2,:); + edge = cell(3,1); + edge(1) = s1; edge(2) = s2; edge(3) = s3; + d = zeros(2,3,nelem); #Lati * Coordinate * Elementi + for jj = 1:3 + tempx = (p(1,edge{jj}(1,:)) + p(1,edge{jj}(2,:)))/2; + tempy = (p(2,edge{jj}(1,:)) + p(2,edge{jj}(2,:)))/2; + temp = [tempx; tempy]; + d(:,jj,:) = temp; + endfor +endfunction + +function [shg] = shapegrad(p,t) + + ## Compute the gradient of the hat functions + + x0 = p(1,t(1,:)); + y0 = p(2,t(1,:)); + x1 = p(1,t(2,:)); + y1 = p(2,t(2,:)); + x2 = p(1,t(3,:)); + y2 = p(2,t(3,:)); + + denom = (-(x1.*y0) + x2.*y0 + x0.*y1 - x2.*y1 - x0.*y2 + x1.*y2); + shg(1,1,:) = (y1 - y2)./denom; + shg(2,1,:) = -(x1 - x2)./denom; + shg(1,2,:) = -(y0 - y2)./denom; + shg(2,2,:) = (x0 - x2)./denom; + shg(1,3,:) = (y0 - y1)./denom; + shg(2,3,:) = -(x0 - x1)./denom; +endfunction + +%!test +%! [mesh] = msh2m_structured_mesh(0:.5:1, 0:.5:1, 1, 1:4, "left"); +%! [mesh.bar, mesh.cir, mesh.emidp, mesh.slength, mesh.cdist, mesh.area,mesh.midedge] = msh2m_geometrical_properties(mesh,"bar","cir","emidp","slength","cdist","area","midedge"); +%! bar = [0.16667 0.16667 0.66667 0.66667 0.33333 0.33333 0.83333 0.83333 +%! 0.16667 0.66667 0.16667 0.66667 0.33333 0.83333 0.33333 0.83333]; +%! cir = [0.25000 0.25000 0.75000 0.75000 0.25000 0.25000 0.75000 0.75000 +%! 0.25000 0.75000 0.25000 0.75000 0.25000 0.75000 0.25000 0.75000]; +%! emidp =[0.25000 0.75000 1.00000 1.00000 0.25000 0.75000 0.00000 0.00000 +%! 0.00000 0.00000 0.25000 0.75000 1.00000 1.00000 0.25000 0.75000]; +%! slength =[0.70711 0.70711 0.70711 0.70711 0.50000 0.50000 0.50000 0.50000 +%! 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 0.50000 +%! 0.50000 0.50000 0.50000 0.50000 0.70711 0.70711 0.70711 0.70711]; +%! cdist = [0.00000 0.00000 0.00000 0.00000 0.50000 0.50000 0.25000 0.25000 +%! 0.25000 0.25000 0.50000 0.50000 0.50000 0.25000 0.50000 0.25000 +%! 0.25000 0.50000 0.25000 0.50000 0.00000 0.00000 0.00000 0.00000]; +%! area = [ 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500 ; 0.12500]; +%! midedge = zeros(2,3,8); +%! midedge(:,:,1) = [0.25000 0.00000 0.25000 +%! 0.25000 0.25000 0.00000]; +%! midedge(:,:,2) = [0.25000 0.00000 0.25000 +%! 0.75000 0.75000 0.50000]; +%! midedge(:,:,3) = [0.75000 0.50000 0.75000 +%! 0.25000 0.25000 0.00000]; +%! midedge(:,:,4) = [0.75000 0.50000 0.75000 +%! 0.75000 0.75000 0.50000]; +%! midedge(:,:,5) = [0.50000 0.25000 0.25000 +%! 0.25000 0.50000 0.25000]; +%! midedge(:,:,6) = [0.50000 0.25000 0.25000 +%! 0.75000 1.00000 0.75000]; +%! midedge(:,:,7) = [1.00000 0.75000 0.75000 +%! 0.25000 0.50000 0.25000]; +%! midedge(:,:,8) = [1.00000 0.75000 0.75000 +%! 0.75000 1.00000 0.75000]; +%! toll = 1e-4; +%! assert(mesh.bar,bar,toll); +%! assert(mesh.cir,cir,toll); +%! assert(mesh.emidp,emidp,toll); +%! assert(mesh.slength,slength,toll); +%! assert(mesh.cdist,cdist,toll); +%! assert(mesh.area,area,toll); +%! assert(mesh.midedge,midedge,toll);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_gmsh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,150 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh2m_gmsh(@var{geometry},@var{option},@var{value},...) +## +## Construct an unstructured triangular 2D mesh making use of the free +## software gmsh. +## +## The compulsory argument @var{geometry} is the basename of the +## @code{*.geo} file to be meshed. +## +## The optional arguments @var{option} and @var{value} identify +## respectively a gmsh option and its value. For more information +## regarding the possible option to pass, refer to gmsh manual or gmsh +## site @url{http://www.geuz.org/gmsh/}. +## +## The returned value @var{mesh} is a PDE-tool like mesh structure. +## +## @seealso{msh2m_structured_mesh, msh3m_gmsh, msh2m_mesh_along_spline} +## @end deftypefn + +function [mesh] = msh2m_gmsh(geometry,varargin) + + ## Check input + if !mod(nargin,2) # Number of input parameters + error("msh2m_gmsh: wrong number of input parameters."); + endif + ## FIXME: add input type check? + + ## Build mesh + noptions = (nargin - 1) / 2; # Number of passed options + + ## Construct system command string + verbose = 1; + optstring = ""; + for ii = 1:noptions + option = varargin{2*(ii)-1}; + value = varargin{2*ii}; + ## Check for verbose option + if strcmp(option,"v") + verbose = value; + endif + if !ischar(value) + value = num2str(value); + endif + optstring = [optstring," -",option," ",value]; + endfor + + ## Invoke gmsh + if (verbose) + printf("\n"); + printf("Generating mesh...\n"); + endif + system(["gmsh -format msh -2 -o " geometry ".msh" optstring " " geometry ".geo"]); + + ## Build structure fields + if (verbose) + printf("Processing gmsh data...\n"); + endif + ## Points + com_p = "awk '/\\$Nodes/,/\\$EndNodes/ {print $2, $3 > ""msh_p.txt""}' "; + ## Side edges + com_e = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""1"") print $7, $8,$5 > ""msh_e.txt""}' "; + ## Triangles + com_t = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""2"") print $7, $8, $9, $5 > ""msh_t.txt""}' "; + + command = [com_p,geometry,".msh ; "]; + command = [command,com_e,geometry,".msh ; "]; + command = [command,com_t,geometry,".msh"]; + + system(command); + + ## Create PDE-tool like structure + if (verbose) + printf("Creating PDE-tool like mesh...\n"); + endif + p = load("msh_p.txt")'; # Mesh-points + tmp = load("msh_e.txt")'; # Mesh surface-edges + be = zeros(7,columns(tmp)); + be([1,2,5],:) = tmp; + t = load("msh_t.txt")'; # Mesh tetrahedra + + ## Remove hanging nodes + if (verbose) + printf("Check for hanging nodes...\n"); + endif + nnodes = columns(p); + in_msh = intersect( 1:nnodes , t(1:3,:) ); + if length(in_msh) != nnodes + new_num(in_msh) = [1:length(in_msh)]; + t(1:3,:) = new_num(t(1:3,:)); + be(1:2,:) = new_num(be(1:2,:)); + p = p(:,in_msh); + endif + + ## Set region numbers in edge structure + if (verbose) + printf("Setting region number in edge structure...\n"); + endif + msh = struct("p",p,"t",t,"e",be); + msh.be = msh2m_topological_properties(msh, "boundary"); + msh.e(6,:) = msh.t(4,msh.be(1,:)); + jj = find (sum(msh.be>0)==4); + msh.e(7,jj) = msh.t(4,msh.be(3,jj)); + + mesh = struct("p",p,"e",be,"t",t); + + ## Delete temporary files + if (verbose) + printf("Deleting temporary files...\n"); + endif + system(["rm -f msh_p.txt msh_e.txt msh_t.txt msh_s.txt *.msh"]); + +endfunction + +%!test +%! fid = fopen("circle.geo","w"); +%! fprintf(fid,"Point(1) = {0, 0, 0, 1};\n"); +%! fprintf(fid,"Point(2) = {1, 0, 0, 1};\n"); +%! fprintf(fid,"Point(3) = {-1, 0, 0, 1};\n"); +%! fprintf(fid,"Circle(1) = {3, 1, 2};\n"); +%! fprintf(fid,"Circle(2) = {2, 1, 3};\n"); +%! fprintf(fid,"Line Loop(4) = {2, 1};\n"); +%! fprintf(fid,"Plane Surface(4) = {4};"); +%! fclose(fid); +%! mesh = msh2m_gmsh("circle","v",0); +%! system("rm circle.geo"); +%! nnodest = length(unique(mesh.t)); +%! nnodesp = columns(mesh.p); +%! assert(nnodest,nnodesp); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_jiggle_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,118 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{newmsh}]} = @ +## msh2m_jiggle_mesh(@var{msh},@var{steps}) +## +## Equalize the size of triangle edges setting a spring of rest length +## @var{factor}*@var{area} along each edge of the mesh and solving for +## static equilibrium. +## +## The non-linear eqautions of the system obtained are solved via a +## non-linear Gauss-Seidel method. @var{step} is the number of steps of +## the method to be applied. +## +## May be useful when distorting a mesh, type @code{demo +## msh2m_jiggle_mesh} to see some examples. +## +## @seealso{msh2m_displacement_smoothing, msh2m_equalize_mesh} +## +## @end deftypefn + +function [msh] = msh2m_jiggle_mesh(msh,steps) + + ## Check input + if nargin != 2 # Number of input parameters + error("msh2m_jiggle_mesh: wrong number of input parameters."); + elseif !(isstruct(msh) && isfield(msh,"p") && + isfield(msh,"t") && isfield(msh,"e")) + error("msh2m_jiggle_mesh: first input is not a valid mesh structure."); + elseif !isscalar(steps) + error("msh2m_jiggle_mesh: second argument is not a valid scalar"); + endif + + ## Solve for static equilibrium + nel= columns(msh.t); + nnodes = columns(msh.p); + + x = msh.p(1,:)'; + y = msh.p(2,:)'; + + dnodes = unique(msh.e(1:2,:)(:)); + vnodes = setdiff(1:nnodes,dnodes); + + ## Find node neighbours + ## FIXME: should this go into msh2m_topological_properties ? + sides = msh2m_topological_properties(msh,"sides"); + for inode = 1:nnodes + neig{inode} = (sides(:, sides(1,:) == inode | sides(2,:) == inode))(:); + neig{inode} (neig{inode} == inode) = []; + endfor + + for istep = 1:steps + for inode =vnodes + xx = x(neig{inode}) * ones(size(neig{inode}))'; + lx = abs ( xx - xx' )(:); + mx = ( xx + xx' )(:)/2; + x(inode) = sum(mx.*lx)/sum(lx); + + yy = y(neig{inode}) * ones(size(neig{inode}))'; + ly = abs ( yy - yy' )(:); + my = (yy + yy')(:)/2; + y(inode) = sum(my.*ly)/sum(ly); + endfor + endfor + + msh.p = [x';y']; + +endfunction + +%!demo +%! ### distort a mesh on a square equalizing at each step +%! msh = msh2m_structured_mesh(linspace(0,1,10),linspace(0,1,10),1,1:4,"right"); +%! dnodes = msh2m_nodes_on_sides(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 = -.4*sin(x(dnodes).*y(dnodes)*pi/2); +%! Nsteps = 30; +%! for ii=1:Nsteps +%! dx(dnodes) = dxtot; +%! dy(dnodes) = dytot; +%! [Ax,Ay] = msh2m_displacement_smoothing(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,:)'; +%! msh = msh2m_jiggle_mesh(msh,10); +%! hold on; +%! triplot(msh.t(1:3,:)',msh.p(1,:)',msh.p(2,:)'); +%! hold off; +%! pause(.5) +%! endfor \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_join_structured_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,163 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh2m_join_structured_mesh(@var{mesh1},@var{mesh2},@var{s1},@var{s2}) +## +## Join the two structured meshes @var{mesh1} and @var{mesh2} into one +## single mesh. +## +## The two meshes must share a common edge identified by @var{s1} and +## @var{s2}. +## +## @strong{WARNING}: the two meshes must share the same vertexes on the +## common edge. +## +## @seealso{msh2m_structured_mesh, msh2m_gmsh, msh2m_submesh, +## msh3m_join_structured_mesh} +## @end deftypefn + +function [mesh] = msh2m_join_structured_mesh(mesh1,mesh2,s1,s2) + + ## Check input + if nargin != 4 # Number of input parameters + error("msh2m_join_structured_mesh: wrong number of input parameters."); + elseif !(isstruct(mesh1) && isfield(mesh1,"p") && + isfield (mesh1,"e") && isfield(mesh1,"t") && + isstruct(mesh2) && isfield(mesh2,"p") && + isfield (mesh2,"e") && isfield(mesh2,"t") ) + error("msh2m_join_structured_mesh: invalid mesh structure passed as input."); + elseif !(isvector(s1) && isvector(s2)) + error("msh2m_join_structured_mesh: shared geometrical sides are not vectors."); + elseif (length(s1) != length(s2)) + error("msh2m_join_structured_mesh: vectors containing shared geometrical sides are not of the same length."); + endif + + ## Join meshes + + ## Make sure that the outside world is always on the same side of the + ## boundary of mesh1 + [mesh1.e(6:7,:),I] = sort(mesh1.e(6:7,:)); + for ic=1:size(mesh1.e,2) + mesh1.e(1:2,ic) = mesh1.e(I(:,ic),ic); + endfor + + ## FIXME: here a check could be added to see whether + ## the coordinate points of the two meshes coincide on the + ## side edges + intnodes1=[]; + intnodes2=[]; + + ## FIXME: Can the following cycle be replaced by + ## msh2m_nodes_on_sides? + j1=[];j2=[]; + for is=1:length(s1) + side1 = s1(is); + side2 = s2(is); + [i,j] = find(mesh1.e(5,:)==side1); + j1=[j1 j]; + [i,j] = find(mesh2.e(5,:)==side2); + oldregion(side1) = max(max(mesh2.e(6:7,j))); + j2=[j2 j]; + endfor + + intnodes1=[mesh1.e(1,j1),mesh1.e(2,j1)]; + intnodes2=[mesh2.e(1,j2),mesh2.e(2,j2)]; + + intnodes1 = unique(intnodes1); + [tmp,I] = sort(mesh1.p(1,intnodes1)); + intnodes1 = intnodes1(I); + [tmp,I] = sort(mesh1.p(2,intnodes1)); + intnodes1 = intnodes1(I); + + intnodes2 = unique(intnodes2); + [tmp,I] = sort(mesh2.p(1,intnodes2)); + intnodes2 = intnodes2(I); + [tmp,I] = sort(mesh2.p(2,intnodes2)); + intnodes2 = intnodes2(I); + + ## Delete redundant edges + mesh2.e(:,j2) = []; + + ## Change edge numbers + indici=[]; + consecutivi=[]; + indici = unique(mesh2.e(5,:)); + consecutivi (indici) = [1:length(indici)]+max(mesh1.e(5,:)); + mesh2.e(5,:)=consecutivi(mesh2.e(5,:)); + + ## Change node indices in connectivity matrix and edge list + indici=[]; consecutivi=[]; + indici = 1:size(mesh2.p,2); + offint = setdiff(indici,intnodes2); + consecutivi (offint) = [1:length(offint)]+size(mesh1.p,2); + consecutivi (intnodes2) = intnodes1; + mesh2.e(1:2,:)=consecutivi(mesh2.e(1:2,:)); + mesh2.t(1:3,:)=consecutivi(mesh2.t(1:3,:)); + + ## Delete redundant points + mesh2.p(:,intnodes2) = []; + + ## Set region numbers + regions = unique(mesh1.t(4,:)); # Mesh 1 + newregions(regions) = 1:length(regions); + mesh1.t(4,:) = newregions(mesh1.t(4,:)); + + regions = unique(mesh2.t(4,:)); # Mesh 2 + newregions(regions) = [1:length(regions)]+max(mesh1.t(4,:)); + mesh2.t(4,:) = newregions(mesh2.t(4,:)); + + ## Set adjacent region numbers in edge structure 2 + [i,j] = find(mesh2.e(6:7,:)); + i = i+5; + mesh2.e(i,j) = newregions(mesh2.e(i,j)); + ## Set adjacent region numbers in edge structure 1 + mesh1.e(6,j1) = newregions(oldregion(mesh1.e(5,j1))); + + ## Make the new p structure + mesh.p = [mesh1.p mesh2.p]; + mesh.e = [mesh1.e mesh2.e]; + mesh.t = [mesh1.t mesh2.t]; + +endfunction + +%!test +%! [mesh1] = msh2m_structured_mesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); +%! [mesh2] = msh2m_structured_mesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); +%! [mesh] = msh2m_join_structured_mesh(mesh1,mesh2,2,4); +%! p = [0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 1.00000 1.00000 1.00000 1.50000 1.50000 1.50000 2.00000 2.00000 2.00000 +%! 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000]; +%! e = [1 4 7 8 3 6 1 2 7 10 13 14 9 12 +%! 4 7 8 9 6 9 2 3 10 13 14 15 12 15 +%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +%! 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +%! 1 1 2 2 3 3 4 4 5 5 6 6 7 7 +%! 0 0 2 2 0 0 0 0 0 0 0 0 0 0 +%! 1 1 1 1 1 1 1 1 2 2 2 2 2 2]; +%! t = [1 2 4 5 2 3 5 6 7 8 10 11 8 9 11 12 +%! 4 5 7 8 4 5 7 8 10 11 13 14 10 11 13 14 +%! 2 3 5 6 5 6 8 9 8 9 11 12 11 12 14 15 +%! 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]; +%! toll = 1e-4; +%! assert(mesh.p,p,toll); +%! assert(mesh.e,e,toll); +%! assert(mesh.t,t,toll); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_mesh_along_spline.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,93 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh2m_mesh_along_spline(@var{xc},@var{yc},@var{Nnx},@var{Nny},@var{sigma}) +## +## Generate a structured mesh in a thin layer of size @var{sigma} +## sitting on a natural Catmull-Rom type cubic spline with control +## points @var{xc}, @var{yc}. +## +## If @var{Nnx} and @var{Nny} are scalars, the mesh has @var{Nnx} nodes +## in the direction along the spline and @var{Nny} in the normal +## direction. +## +## If @var{Nnx} and @var{Nny} are vectors they indicate the curvilinear +## coordinates of the mesh nodes. +## +## The returned value @var{mesh} is a PDE-tool like mesh structure. +## +## Be aware that if @var{sigma} is not much smaller than the curvature +## of the line the resulting mesh may be invalid. +## +## @seealso{msh2m_structured_mesh, msh2m_gmsh, msh3m_structured_mesh} +## @end deftypefn + + +function msh2 = msh2m_mesh_along_spline(xc,yc,Nnx,Nny,sigma) + + ## Check input + ## FIXME: input type not checked for the first 4 arguments + if (nargin != 5) # Number of input parameters + error("msh2m_mesh_along_spline: wrong number of input parameters."); + elseif (!isscalar(sigma)) + error("msh2m_mesh_along_spline: sigma must be a valid scalar value."); + endif + + ## Construct mesh + s = [0:length(xc)-1]; + xsPP = catmullrom ( s, xc ); + ysPP = catmullrom ( s, yc ); + + if (length(Nnx)>1) + ss = Nnx(:).'; + else + ss = linspace(0,s(end),Nnx); + endif + xs = ppval(xsPP,ss); + ys = ppval(ysPP,ss); + + dxsPP = fnder(xsPP,1); + dysPP = fnder(ysPP,1); + + + nx = -ppval(dysPP,ss)'; + ny = ppval(dxsPP,ss)'; + + nx = nx ./ sqrt(nx.^2+ny.^2); + ny = ny ./ sqrt(nx.^2+ny.^2); + + if (length(Nny)>1) + ssy = Nny(:).'; + else + ssy = linspace(0,1,Nny); + endif + + msh2 = msh2m_structured_mesh([1:length(ss)], ssy, 1, 1:4); + + jj = (msh2.p(1,:)); + p(1,:) = xs(jj) + sigma*nx(jj)' .* msh2.p(2,:); + p(2,:) = ys(jj) + sigma*ny(jj)' .* msh2.p(2,:); + + msh2.p = p; + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_nodes_on_sides.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,65 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{nodelist}]} = @ +## msh2m_nodes_on_sides(@var{mesh},@var{sidelist}) +## +## Return a list of @var{mesh} nodes lying on the sides specified in +## @var{sidelist}. +## +## @seealso{msh2m_geometrical_properties, msh2m_topological_properties, +## msh3m_nodes_on_faces} +## @end deftypefn + +function [nodelist] = msh2m_nodes_on_sides(mesh,sidelist) + + ## Check input + if nargin != 2 # Number of input parameters + error("msh2m_nodes_on_sides: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield(mesh,"t") && isfield(mesh,"e")) + error("msh2m_nodes_on_sides: first input is not a valid mesh structure."); + elseif !isnumeric(sidelist) + error("msh2m_nodes_on_sides: only numeric value admitted as sidelist."); + endif + + ## Search nodes + + edgelist = []; + + for ii = 1:length(sidelist) + edgelist=[edgelist,find(mesh.e(5,:)==sidelist(ii))]; + endfor + + nodelist = mesh.e(1:2,edgelist); + nodelist = [nodelist(1,:) nodelist(2,:)]; + nodelist = unique(nodelist); + +endfunction + +%!test +%! [mesh1] = msh2m_structured_mesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); +%! [mesh2] = msh2m_structured_mesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); +%! [mesh] = msh2m_join_structured_mesh(mesh1,mesh2,2,4); +%! [nodelist] = msh2m_nodes_on_sides(mesh,[1 2]); +%! reallist = [1 4 7 8 9]; +%! assert(nodelist,reallist);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_structured_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,249 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh2m_structured_mesh(@var{x},@var{y},@var{region},@var{sides},@var{string}) +## +## Construct a structured triangular 2D mesh on a rectangular domain. +## +## @itemize @bullet +## @item @var{x} and @var{y} are the one dimensional mesh vector of the +## corresponding Cartesian axis. +## @item @var{region} is a number identifying the geometrical surface +## region, while @var{sides} is a 4 components vector containing the +## numbers used to identify the geometrical side edges. +## @item @var{string} is an optional value specifying the orientation of +## the diagonal edge of the structured mesh. It may take the value +## @code{"right"} (default), @code{"left"}, @code{"random"}. +## @end itemize +## +## The returned value @var{mesh} is a PDE-tool like mesh structure +## composed of the following fields: +## @itemize @minus +## @item @var{p}: matrix with size 2 times number of mesh points. +## @itemize @bullet +## @item 1st row: x-coordinates of the points. +## @item 2nd row: y-coordinates of the points. +## @end itemize +## @item @var{e}: matrix with size 7 times number of mesh side edges. +## @itemize @bullet +## @item 1st row: number of the first vertex of the side edge. +## @item 2nd row: number of the second vertex of the side edge. +## @item 3rd row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 4th row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 5th row: number of the geometrical border containing the side +## edge. +## @item 6th row: number of the geometrical surface to the right of +## side edge. +## @item 7th row: number of the geometrical surface to the left of the +## side edge. +## @end itemize +## @item @var{t}: matrix with size 4 times number of mesh elements. +## @itemize @bullet +## @item 1st row: number of the first vertex of the element. +## @item 2nd row: number of the second vertex of the element. +## @item 3rd row: number of the third vertex of the element. +## @item 4th row: number of the geometrical surface containing the element. +## @end itemize +## @end itemize +## +## @seealso{msh3m_structured_mesh, msh2m_gmsh, msh2m_mesh_along_spline, +## msh2m_join_structured_mesh, msh2m_submesh} +## @end deftypefn + +function [mesh] = msh2m_structured_mesh(x,y,region,sides,varargin) + + ## Check input + if ((nargin < 4) || (nargin > 5)) # Number of input parameters + error("msh2m_structured_mesh: wrong number of input parameters."); + elseif !(isvector(x) && isnumeric(x) && isvector(y) && isnumeric(y)) + error("msh2m_structured_mesh: X and Y must be valid numeric vectors."); + elseif !isscalar(region) + error("msh2m_structured_mesh: REGION must be a valid scalar."); + elseif !(isvector(sides) && (length(sides) == 4)) + error("msh2m_structured_mesh: SIDES must be a 4 components vector."); + endif + + ## Build mesh + default = "right"; + + ## Check if any orientation is given + if length(varargin)==0 + string = default; + else + string = varargin{1}; + endif + + ## Construct mesh + switch string + case "right" + [mesh] = Ustructmesh_right(x, y, region, sides); + case "left" + [mesh] = Ustructmesh_left(x, y, region, sides); + case "random" + [mesh] = Ustructmesh_random(x, y, region, sides); + otherwise + error("msh2m_structured_mesh: STRING has not a valid value."); + endswitch + +endfunction + +## Right diagonal structured mesh +function [mesh]=Ustructmesh_right(x,y,region,sides) + + x = sort(x); + y = sort(y); + + nx = length(x); + ny = length(y); + [XX,YY] = meshgrid(x,y); + p = [XX(:),YY(:)]'; + iiv (ny,nx)=0; + iiv(:)=1:nx*ny; + iiv(end,:)=[]; + iiv(:,end)=[]; + iiv=iiv(:)'; + t = [[iiv;iiv+ny;iiv+ny+1],[iiv;iiv+ny+1;iiv+1] ]; + t (4,:)=region; + + l1 = 1+ny*([1:nx]-1); + l4 = 1:ny; + l2 = ny*(nx-1)+1:nx*ny; + l3 = ny + l1 -1; + + e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) + l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region + ]; + mesh.p = p; + mesh.e = e; + mesh.t = t; + +endfunction + +## Left diagonal structured mesh +function [mesh]=Ustructmesh_left(x,y,region,sides) + + x = sort(x); + y = sort(y); + + nx = length(x); + ny = length(y); + [XX,YY] = meshgrid(x,y); + p = [XX(:),YY(:)]'; + iiv (ny,nx)=0; + iiv(:)=1:nx*ny; + iiv(end,:)=[]; + iiv(:,end)=[]; + iiv=iiv(:)'; + t = [[iiv;iiv+ny;iiv+1],[iiv+1;iiv+ny;iiv+ny+1] ]; + t (4,:)=region; + + l1 = 1+ny*([1:nx]-1); + l4 = 1:ny; + l2 = ny*(nx-1)+1:nx*ny; + l3 = ny + l1 -1; + + e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) + l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region + ]; + mesh.p = p; + mesh.e = e; + mesh.t = t; + +endfunction + +## Random diagonal structured mesh +function [mesh]=Ustructmesh_random(x,y,region,sides) + + x = sort(x); + y = sort(y); + + nx = length(x); + ny = length(y); + [XX,YY] = meshgrid(x,y); + p = [XX(:),YY(:)]'; + iiv (ny,nx)=0; + iiv(:)=1:nx*ny; + iiv(end,:)=[]; + iiv(:,end)=[]; + iiv=iiv(:)'; + + niiv = length(iiv); + theperm = iiv(randperm(niiv)); + first = theperm(1:floor(niiv/2)); + second = theperm(floor(niiv/2)+1:end); + + t = [[first;first+ny;first+ny+1],[first;first+ny+1;first+1] ]; + t = [t,[second;second+ny;second+1],[second+ny;second+ny+1;second+1] ]; + + t (4,:)=region; + + l1 = 1+ny*([1:nx]-1); + l4 = 1:ny; + l2 = ny*(nx-1)+1:nx*ny; + l3 = ny + l1 -1; + + e = [ l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1]) + l1([2:end]) l2([2:end]) l3([2:end]) l4([2:end]) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + l1([1:end-1])*0+sides(1) l2([1:end-1])*0+sides(2) l3([1:end-1])*0+sides(3) l4([1:end-1])*0+sides(4) + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0 + [l1([1:end-1]) l2([1:end-1]) l3([1:end-1]) l4([1:end-1])]*0+region + ]; + mesh.p = p; + mesh.e = e; + mesh.t = t; + +endfunction + +%!test +%! x = y = linspace(0,1,3); +%! msh = msh2m_structured_mesh(x,y,1,[1:4]); +%! p = [0.00000 0.00000 0.00000 0.50000 0.50000 0.50000 \ +%! 1.00000 1.00000 1.00000 +%! 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 \ +%! 0.00000 0.50000 1.00000]; +%! assert(msh.p,p) +%! e = [1 4 7 8 3 6 1 2 +%! 4 7 8 9 6 9 2 3 +%! 0 0 0 0 0 0 0 0 +%! 0 0 0 0 0 0 0 0 +%! 1 1 2 2 3 3 4 4 +%! 0 0 0 0 0 0 0 0 +%! 1 1 1 1 1 1 1 1]; +%! assert(msh.e,e) +%! t =[1 2 4 5 1 2 4 5 +%! 4 5 7 8 5 6 8 9 +%! 5 6 8 9 2 3 5 6 +%! 1 1 1 1 1 1 1 1]; +%! assert(msh.t,t) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_submesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,111 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{omesh},@var{nodelist},@var{elementlist}]} = @ +## msh2m_submesh(@var{imesh},@var{intrfc},@var{sdl}) +## +## Extract the subdomain(s) in @var{sdl} from @var{imesh}. +## +## The row vector @var{intrfc} contains the internal interface sides to +## be maintained (field @code{mesh.e(5,:)}). It can be empty. +## +## Return the vectors @var{nodelist} and @var{elementlist} containing +## respectively the list of nodes and elements of the original mesh that +## are part of the selected subdomain(s). +## +## @seealso{msh2m_join_structured_mesh, msh3m_submesh, +## msh3e_surface_mesh} +## @end deftypefn + +function [omesh,nodelist,elementlist] = msh2m_submesh(imesh,intrfc,sdl) + + ## Check input + if nargin != 3 + error("msh2m_submesh: wrong number of input parameters."); + endif + if !isstruct(imesh) + error("msh2m_submesh: first input is not a valid mesh structure."); + endif + if !isvector(sdl) + error("msh2m_submesh: third input is not a valid vector."); + endif + + ## Extract sub-mesh + nsd = length(sdl); # number of subdomains + + ## Set list of output triangles + elementlist=[]; + for isd=1:nsd + elementlist = [elementlist find(imesh.t(4,:) == sdl(isd))]; + endfor + + omesh.t = imesh.t(:,elementlist); + + ## Set list of output nodes + nodelist = unique(reshape(imesh.t(1:3,elementlist),1,[])); + omesh.p = imesh.p(:,nodelist); + + ## Use new node numbering in connectivity matrix + indx(nodelist) = [1:length(nodelist)]; + iel = [1:length(elementlist)]; + omesh.t(1:3,iel) = indx(omesh.t(1:3,iel)); + + ## Set list of output edges + omesh.e =[]; + for isd=1:nsd + omesh.e = [omesh.e imesh.e(:,imesh.e(7,:)==sdl(isd))]; + omesh.e = [omesh.e imesh.e(:,imesh.e(6,:)==sdl(isd))]; + endfor + omesh.e=unique(omesh.e',"rows")'; + + ## Use new node numbering in boundary segment list + ied = [1:size(omesh.e,2)]; + omesh.e(1:2,ied) = indx(omesh.e(1:2,ied)); + +endfunction + +%!test +%! [mesh1] = msh2m_structured_mesh(0:.5:1, 0:.5:1, 1, 1:4, 'left'); +%! [mesh2] = msh2m_structured_mesh(1:.5:2, 0:.5:1, 1, 1:4, 'left'); +%! [mesh] = msh2m_join_structured_mesh(mesh1,mesh2,2,4); +%! [omesh,nodelist,elementlist] = msh2m_submesh(mesh,[],2); +%! p = [1.00000 1.00000 1.00000 1.50000 1.50000 1.50000 2.00000 2.00000 2.00000 +%! 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000 0.00000 0.50000 1.00000]; +%! e = [1 1 2 3 4 6 7 8 +%! 2 4 3 6 7 9 8 9 +%! 0 0 0 0 0 0 0 0 +%! 0 0 0 0 0 0 0 0 +%! 2 5 2 7 5 7 6 6 +%! 2 0 2 0 0 0 0 0 +%! 1 2 1 2 2 2 2 2]; +%! t = [1 2 4 5 2 3 5 6 +%! 4 5 7 8 4 5 7 8 +%! 2 3 5 6 5 6 8 9 +%! 2 2 2 2 2 2 2 2]; +%! nl = [7 8 9 10 11 12 13 14 15]; +%! el = [9 10 11 12 13 14 15 16]; +%! toll = 1e-4; +%! assert(omesh.p,p,toll); +%! assert(omesh.e,e); +%! assert(omesh.t,t); +%! assert(nodelist,nl); +%! assert(elementlist,el);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2m_topological_properties.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,288 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{varargout}]} = @ +## msh2m_topological_properties(@var{mesh},[@var{string1},@var{string2},...]) +## +## Compute @var{mesh} topological properties identified by input strings. +## +## Valid properties are: +## @itemize @bullet +## @item @code{"n"}: return a matrix with size 3 times the number of +## mesh elements containing the list of its neighbours. The entry +## @code{M(i,j)} in this matrix is the mesh element sharing the side +## @code{i} of triangle @code{j}. If no such element exists (i.e. for +## boundary edges) a value of @code{NaN} is set. +## @item @code{"sides"}: return a matrix with size 2 times number of +## sides.The entry @code{M(i,j)} is the index of the i-th vertex of j-th +## side. +## @item @code{"ts"}: return a matrix with size 3 times the number of +## mesh elements containing the sides associated with each element. +## @item @code{"tws"}:return a matrix with size 2 times the number of +## mesh sides containing the elements associated with each side. For a +## side belonging to one triangle only a value of @code{NaN} is set. +## @item @code{"coinc"}: return a matrix with 2 rows. Each column +## contains the indices of two triangles sharing the same circumcenter. +## @item @code{"boundary"}: return a matrix with size 2 times the number +## of side edges. The first row contains the mesh element to which the +## side belongs, the second row is the local index of this edge. +## @end itemize +## +## The output will contain the geometrical properties requested in the +## input in the same order specified in the function call. +## +## If an unexpected string is given as input, an empty vector is +## returned in output. +## +## @seealso{mshm2m_geometrical_properties, msh3m_geometrical_properties} +## @end deftypefn + +function [varargout] = msh2m_topological_properties(mesh,varargin) + + ## Check input + if nargin < 2 # Number of input parameters + error("msh2m_topological_properties: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield(mesh,"t") && isfield(mesh,"e")) + error("msh2m_topological_properties: first input is not a valid mesh structure."); + elseif !iscellstr(varargin) + error("msh2m_topological_properties: only string value admitted for properties."); + endif + + ## Compute properties + p = mesh.p; + e = mesh.e; + t = mesh.t; + + nelem = columns(t); # Number of elements in the mesh + [n,ts,tws,sides] = neigh(t,nelem); + + for nn = 1:length(varargin) + request = varargin{nn}; + switch request + + case "n" # Neighbouring triangles + if isfield(mesh,"n") + varargout{nn} = mesh.n; + else + varargout{nn} = n; + endif + + case "sides" # Global edge matrix + if isfield(mesh,"sides") + varargout{nn} = mesh.sides; + else + varargout{nn} = sides; + endif + + case "ts" # Triangle sides matrix + if isfield(mesh,"ts") + varargout{nn} = mesh.ts; + else + varargout{nn} = ts; + endif + + case "tws" # Trg with sides matrix + if isfield(mesh,"tws") + varargout{nn} = mesh.tws; + else + varargout{nn} = tws; + endif + + case "coinc" # Coincident circumcenter matrix + if isfield(mesh,"coinc") + varargout{nn} = mesh.coinc; + else + if isfield(mesh,"cdist") + d = mesh.cdist; + else + [d] = msh2m_geometrical_properties(mesh,"cdist"); + endif + [b] = coinc(n,d); + varargout{nn} = b; + clear b + endif + + case "boundary" # Boundary edge matrix + if isfield(mesh,"boundary") + varargout{nn} = mesh.boundary; + else + [b] = borderline(e,t); + varargout{nn} = b; + clear b + endif + + otherwise + warning("msh2m_topological_properties: unexpected value in property string. Empty vector passed as output.") + varargout{nn} = []; + endswitch + + endfor + +endfunction + +function [n,ts,triwside,sides] = neigh(t,nelem) + + n = nan*ones(3,nelem); + t = t(1:3,:); + + s3 = sort(t(1:2,:),1); + s1 = sort(t(2:3,:),1); + s2 = sort(t([3,1],:),1); + + allsides = [s1 s2 s3]'; + [sides, ii, jj] = unique( allsides,"rows"); + sides = sides'; + + ts = reshape(jj,[],3)'; + + triwside = zeros(2,columns(sides)); + for kk =1:3 + triwside(1,ts(kk,1:end)) = 1:nelem; + triwside(2,ts(4-kk,end:-1:1)) = nelem:-1:1; + endfor + + triwside(2,triwside(1,:)==triwside(2,:)) = NaN; + + n(1,:) = triwside(1,ts(1,:)); + n(1,n(1,:)==1:nelem) = triwside(2,ts(1,:))(n(1,:)==1:nelem); + n(2,:) = triwside(1,ts(2,:)); + n(2,n(2,:)==1:nelem) = triwside(2,ts(2,:))(n(2,:)==1:nelem); + n(3,:) = triwside(1,ts(3,:)); + n(3,n(3,:)==1:nelem) = triwside(2,ts(3,:))(n(3,:)==1:nelem); + +endfunction + +function [output] = coinc(n,d); + + ## Tolerance value for considering two point to be coincident + toll = 1e-10; + ## Check the presence of more than two trgs sharing the same circum centre + degen = d < toll; res = sum(degen); + [check] = find(res > 1); + ## Index of the sharing pairs + [ii, jj] = find(degen >= 1); + if isempty(jj) == 0 + temp = zeros(2,length(jj)); + temp(1,:) = jj'; + temp(2,:) = diag(n(ii,jj))'; + temp = sort(temp); + temp = temp'; + [output] = unique(temp,"rows"); + output = output'; + if isempty(check) == 0 + warning("More than two trgs sharing the same circum-centre.") + ## FIXME if more than two trgs shares the same circen ---> construct a cell array + endif + else + output = []; + endif +endfunction + +function [output] = borderline(e,t) + + nelem = columns(e); + t = t(1:3,:); + output = zeros(4,nelem); + for ii = 1:nelem + + point = ( e(1,ii) == t ); + point += ( e(2,ii) == t ); + + [jj1] = find( sum(point(2:3,:)) == 2); + [jj2] = find( sum(point([3 1],:)) == 2); + [jj3] = find( sum(point(1:2,:)) == 2); + + assert( (length(jj1) + length(jj2) + length(jj3)) <= 2 ); + + numtrg = 0; + for jj=1:length(jj1) + output(2*numtrg+1,ii) = jj1(jj); + output(2*numtrg+2,ii) = 1; + numtrg += 1; + endfor + for jj=1:length(jj2) + output(2*numtrg+1,ii) = jj2(jj); + output(2*numtrg+2,ii) = 2; + numtrg += 1; + endfor + for jj=1:length(jj3) + output(2*numtrg+1,ii) = jj3(jj); + output(2*numtrg+2,ii) = 3; + numtrg += 1; + endfor + + endfor +endfunction + +%!test +%! [mesh] = msh2m_structured_mesh(0:.5:1, 0:.5:1, 1, 1:4, "left"); +%! [mesh.n,mesh.sides,mesh.ts,mesh.tws,mesh.coinc,mesh.boundary] = msh2m_topological_properties(mesh,"n","sides","ts","tws","coinc","boundary"); +%! n = [5 6 7 8 3 4 NaN NaN +%! NaN NaN 5 6 2 NaN 4 NaN +%! NaN 5 NaN 7 1 2 3 4]; +%! sides = [1 1 2 2 2 3 3 4 4 5 5 5 6 6 7 8 +%! 2 4 3 4 5 5 6 5 7 6 7 8 8 9 8 9]; +%! ts = [4 6 11 13 8 10 15 16 +%! 1 3 8 10 5 7 12 14 +%! 2 5 9 12 4 6 11 13]; +%! tws = [ 1 1 2 5 2 6 6 3 3 4 7 4 8 8 7 8 +%! NaN NaN NaN 1 5 2 NaN 5 NaN 6 3 7 4 NaN NaN NaN]; +%! coinc = [1 2 3 4 +%! 5 6 7 8]; +%! boundary =[ 1 3 7 8 6 8 1 2 +%! 3 3 1 1 2 2 2 2 +%! 0 0 0 0 0 0 0 0 +%! 0 0 0 0 0 0 0 0]; +%! assert(mesh.n,n); +%! assert(mesh.sides,sides); +%! assert(mesh.ts,ts); +%! assert(mesh.tws,tws); +%! assert(mesh.coinc,coinc); +%! assert(mesh.boundary,boundary); + +%!test +%! mesh.p = []; mesh.e = []; +%! mesh.t = [3 9 10 1 6 9 10 9 8 9 +%! 9 3 1 10 10 10 7 5 9 8 +%! 6 5 7 8 2 6 2 4 4 10 +%! 6 6 6 6 6 6 6 6 6 6]; +%! [mesh.n] = msh2m_topological_properties(mesh,"n"); +%! n = [6 NaN NaN 10 7 5 NaN NaN 8 4 +%! NaN 8 7 NaN NaN 1 5 9 NaN 6 +%! 2 1 4 3 6 10 3 2 10 9]; +%! assert(mesh.n,n); + + +%!test +%! mesh.p = []; mesh.e = []; +%! mesh.t =[ +%! 10 3 6 11 10 3 6 11 1 7 5 9 2 5 11 9 13 6 +%! 14 7 10 15 15 8 11 16 5 11 9 13 6 6 12 10 14 7 +%! 15 8 11 16 11 4 7 12 2 8 6 10 3 2 8 6 10 3 +%! 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]; +%! [mesh.n] = msh2m_topological_properties(mesh,"n"); +%! n =[ +%! NaN 10 5 NaN 4 NaN 10 NaN 14 15 16 17 18 13 NaN 3 1 2 +%! 5 6 7 8 3 NaN 18 15 NaN 2 14 16 NaN 9 10 11 12 13 +%! 17 18 16 5 1 2 3 4 NaN 7 NaN NaN 14 11 8 12 NaN 7]; +%! assert(mesh.n,n);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh2p_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,52 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} msh2p_mesh(@var{mesh}, @var{linespec}) +## +## Plot @var{mesh} with the line specification in @var{linespec} using +## @code{triplot}. +## +## @seealso{triplot} +## +## @end deftypefn + +function msh2p_mesh(mesh,linespec) + + ## Check input + if nargin > 2 + error("msh2p_mesh: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield (mesh,"t") && isfield(mesh,"e")) + error("msh2p_mesh: first input is not a valid mesh structure."); + endif + + tri = mesh.t(1:3,:)'; + x = mesh.p(1,:)'; + y = mesh.p(2,:)'; + + if ~exist("linespec") + linespec = "r"; + endif + + triplot(tri,x,y,linespec); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3e_surface_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,102 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{emesh},@var{snodes},@var{ssides},@var{striangles}]} = @ +## msh3e_surface_mesh(@var{mesh},@var{nsrf},@var{nsides}) +## +## Extract the plane surface @var{nsrf} delimited by @var{nsides} from +## @var{mesh}. +## +## Return the vector @var{snodes} containing the references to input +## mesh nodes (field @code{mesh.p}), the vector @var{ssides} containing +## the references to input mesh side (field @code{mesh.s}) and the +## vector @var{striangles} containing the references to input mesh side +## edges (field @code{mesh.e}). +## +## @strong{WARNING}: the suface MUST be ortogonal to either X, Y or Z +## axis. This should be changed to account for generic 2D surface. +## +## @end deftypefn + +function [emesh,snodes,ssides,striangles] = msh3e_surface_mesh(mesh,nsrf,nsides) + + ## Check input + if nargin != 3 + error("msh3e_surface_mesh: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield (mesh,"t") && isfield(mesh,"e")) + error("msh3e_surface_mesh: first input is not a valid mesh structure."); + elseif !isscalar(nsrf) + error("msh3e_surface_mesh: second input is not a valid scalar."); + elseif !(isvector(nsides) && isnumeric(nsides)) + error("msh3e_surface_mesh: third input is not a valid numeric vector."); + endif + + ## Surface extraction + + ## Extraction of 2D surface elements + striangles = find( mesh.e(10,:) == nsrf ); + t = mesh.e(1:3,striangles); + tmp = reshape(t,[],1); + ## Renumbering + [snodes,ii,jj] = unique(tmp); + nds = 1:length(snodes); + emesh.t = reshape(nds(jj),3,[]); + + ## Extraction of 2D mesh points + points = mesh.p(:,snodes); + + ## Test for normals + ## FIXME: this should disappear as soon as 2D mesh are not supposed to + ## lie on a plane. + if length(unique(points(1,:))) == 1 + xyz = [2,3]; # normal to X coordinate + elseif length(unique(points(2,:))) == 1 + xyz = [1,3]; # normal to Y coordinate + else + xyz = [1,2]; # normal to Z coordinate + endif + + emesh.p = points(xyz,:); + + ## Extraction of 1D side edges + ssides = []; + + for ll = nsides + tmp = find ( mesh.s(3,:) == ll ); + ssides = [ssides,tmp]; + endfor + + nedges = length(ssides); + emesh.e = zeros(7,nedges); + emesh.e(5,:) = mesh.s(3,ssides); + + tmp = reshape(mesh.s(1:2,ssides),[],1); + [enodes,nn,mm] = unique(tmp); + [tmp1, nds, tmp2] = intersect(snodes,enodes); + emesh.e(1:2,:) = reshape(nds(mm),2,[]); + + ## Compute mesh properties + ## FIXME: this has to be removed. MSH should not depend on BIM. + emesh = bim2c_mesh_properties(emesh); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_geometrical_properties.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,181 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{varargout}]} = @ +## msh3m_geometrical_properties(@var{mesh},[@var{string1},@var{string2},...]) +## +## +## Compute @var{mesh} geometrical properties identified by input strings. +## +## Valid properties are: +## @itemize @bullet +## @item @code{"wjacdet"}: return the weigthed Jacobian determinant used +## for the numerical integration with trapezoidal rule over an element. +## @item @code{"shg"}: return a matrix of size 3 times the number of +## elements matrix containing the gradient of P1 shape functions. +## @item @code{"shp"}: return a matrix containing the the value of P1 shape +## functions. +## @end itemize +## +## The output will contain the geometrical properties requested in the +## input in the same order specified in the function call. +## +## If an unexpected string is given as input, an empty vector is +## returned in output. +## +## @seealso{msh2m_topological_properties, msh2m_geometrical_properties} +## @end deftypefn + + +function [varargout] = msh3m_geometrical_properties(imesh,varargin) + + ## Check input + if nargin < 2 # Number of input parameters + error("msh3m_geometrical_properties: wrong number of input parameters."); + elseif !(isstruct(imesh) && isfield(imesh,"p") && + isfield(imesh,"t") && isfield(imesh,"e")) + error("msh3m_geometrical_properties: first input is not a valid mesh structure."); + elseif !iscellstr(varargin) + error("msh3m_geometrical_properties: only string value admitted for properties."); + endif + + ## Compute properties + + ## Extract tetrahedra node coordinates + x1 = imesh.p(1,imesh.t(1,:)); + y1 = imesh.p(2,imesh.t(1,:)); + z1 = imesh.p(3,imesh.t(1,:)); + x2 = imesh.p(1,imesh.t(2,:)); + y2 = imesh.p(2,imesh.t(2,:)); + z2 = imesh.p(3,imesh.t(2,:)); + x3 = imesh.p(1,imesh.t(3,:)); + y3 = imesh.p(2,imesh.t(3,:)); + z3 = imesh.p(3,imesh.t(3,:)); + x4 = imesh.p(1,imesh.t(4,:)); + y4 = imesh.p(2,imesh.t(4,:)); + z4 = imesh.p(3,imesh.t(4,:)); + + for nn = 1:length(varargin) + + request = varargin{nn}; + switch request + case "wjacdet" # Weighted Jacobian determinant + b = wjacdet(x1,y1,z1,\ + x2,y2,z2,\ + x3,y3,z3,\ + x4,y4,z4); + varargout{nn} = b; + clear b + case "area" # Element area + tmp = wjacdet(x1,y1,z1,\ + x2,y2,z2,\ + x3,y3,z3,\ + x4,y4,z4); + b = sum(tmp,1); + varargout{nn} = b; + case "shg" # Gradient of shape functions + b = shg(x1,y1,z1,\ + x2,y2,z2,\ + x3,y3,z3,\ + x4,y4,z4); + varargout{nn} = b; + clear b + case "shp" # Value of shape functions + varargout{nn} = eye(4); + otherwise + warning("msh3m_geometrical_properties: unexpected value in property string. Empty vector passed as output.") + varargout{nn} = []; + endswitch + + endfor + +endfunction + +function [b] = wjacdet(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) + + ## Compute weighted yacobian determinant + + weight = [1/4 1/4 1/4 1/4]'; + + Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); + Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); + Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); + + ## Determinant of the Jacobian of the + ## transformation from the base tetrahedron + ## to the tetrahedron K + detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; + ## Volume of the reference tetrahedron + Kkvolume = 1/6; + + b(:,:) = Kkvolume * weight * detJ; + +endfunction + +function [b] = shg(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4) + + ## Compute gradient of shape functions + + Nb2 = y1.*(z3-z4) + y3.*(z4-z1) + y4.*(z1-z3); + Nb3 = y1.*(z4-z2) + y2.*(z1-z4) + y4.*(z2-z1); + Nb4 = y1.*(z2-z3) + y2.*(z3-z1) + y3.*(z1-z2); + + ## Determinant of the Jacobian of the + ## transformation from the base tetrahedron + ## to the tetrahedron K + detJ = (x2-x1).*Nb2 +(x3-x1).*Nb3 +(x4-x1).*Nb4; + + ## Shape function gradients follow + ## First index represents space direction + ## Second index represents the shape function + ## Third index represents the tetrahedron number + b(1,1,:) = (y2.*(z4-z3) + y3.*(z2-z4) + y4.*(z3-z2))./ detJ; + b(2,1,:) = (x2.*(z3-z4) + x3.*(z4-z2) + x4.*(z2-z3))./ detJ; + b(3,1,:) = (x2.*(y4-y3) + x3.*(y2-y4) + x4.*(y3-y2))./ detJ; + + b(1,2,:) = ( Nb2 ) ./ detJ; + b(2,2,:) = (x1.*(z4-z3) + x3.*(z1-z4) + x4.*(z3-z1)) ./ detJ; + b(3,2,:) = (x1.*(y3-y4) + x3.*(y4-y1) + x4.*(y1-y3)) ./ detJ; + + b(1,3,:) = ( Nb3 ) ./ detJ; + b(2,3,:) = (x1.*(z2-z4) + x2.*(z4-z1) + x4.*(z1-z2)) ./ detJ; + b(3,3,:) = (x1.*(y4-y2) + x2.*(y1-y4) + x4.*(y2-y1)) ./ detJ; + + b(1,4,:) = ( Nb4) ./ detJ; + b(2,4,:) = (x1.*(z3-z2) + x2.*(z1-z3) + x3.*(z2-z1)) ./ detJ; + b(3,4,:) = (x1.*(y2-y3) + x2.*(y3-y1) + x3.*(y1-y2)) ./ detJ; +endfunction + +%!shared mesh,wjacdet,shg,shp +% x = y = z = linspace(0,1,2); +% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6) +% [wjacdet] = msh3m_geometrical_properties(mesh,"wjacdet") +% [shg] = msh3m_geometrical_properties(mesh,"shg") +% [shp] = msh3m_geometrical_properties(mesh,"shp") +%!test +% assert(columns(mesh.t),columns(wjacdet)) +%!test +% assert(size(shg),[3 4 6]) +%!test +% assert(shp,eye(4)) +%!test +% fail(msh3m_geometrical_properties(mesh,"samanafattababbudoiu"),"warning","Unexpected value in passed string. Empty vector passed as output.") \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_gmsh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,133 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh3m_gmsh(@var{geometry},@var{option},@var{value},...) +## +## Construct an unstructured tetrahedral 3D mesh making use of the free +## software gmsh. +## +## The compulsory argument @var{geometry} is the basename of the +## @code{*.geo} file to be meshed. +## +## The optional arguments @var{option} and @var{value} identify +## respectively a gmsh option and its value. For more information +## regarding the possible option to pass, refer to gmsh manual or gmsh +## site @url{http://www.geuz.org/gmsh/}. +## +## The returned value @var{mesh} is a PDE-tool like mesh structure. +## +## @seealso{msh3m_structured_mesh, msh2m_gmsh, msh2m_mesh_along_spline} +## @end deftypefn + +function [mesh] = msh3m_gmsh(geometry,varargin) + + ## Check input + ## Number of input + if !mod(nargin,2) + warning("WRONG NUMBER OF INPUT."); + print_usage; + endif + ## FIXME: add input type check? + + ## Build mesh + noptions = (nargin - 1) / 2; # Number of passed options + + ## Construct system command string + verbose = 1; + optstring = ""; + for ii = 1:noptions + option = varargin{2*(ii)-1}; + value = varargin{2*ii}; + ## Check for verbose option + if strcmp(option,"v") + verbose = value; + endif + if !ischar(value) + value = num2str(value); + endif + optstring = [optstring," -",option," ",value]; + endfor + + ## Generate mesh using Gmsh + if (verbose) + printf("\n"); + printf("Generating mesh...\n"); + endif + system(["gmsh -format msh -3 -o " geometry ".msh" optstring " " geometry ".geo"]); + + if (verbose) + printf("Processing gmsh data...\n"); + endif + ## Points + com_p = "awk '/\\$Nodes/,/\\$EndNodes/ {print $2, $3, $4 > ""msh_p.txt""}' "; + ## Surface edges + com_e = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""2"") print $7, $8, $9, $5 > ""msh_e.txt""}' "; + ## Tetrahedra + com_t = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""4"") print $7, $8, $9, $10, $5 > ""msh_t.txt""}' "; + ## Side edges + com_s = "awk '/\\$Elements/,/\\$EndElements/ {if ($2 == ""1"") print $7, $8, $5 > ""msh_s.txt""}' "; + + command = [com_p,geometry,".msh ; "]; + command = [command,com_e,geometry,".msh ; "]; + command = [command,com_t,geometry,".msh ; "]; + command = [command,com_s,geometry,".msh"]; + + system(command); + + ## Create PDE-tool like structure + if (verbose) + printf("Creating PDE-tool like mesh...\n"); + endif + ## Mesh-points + p = load("msh_p.txt")'; + ## Mesh side-edges + s = load("msh_s.txt")'; + ## Mesh surface-edges + tmp = load("msh_e.txt")'; + be = zeros(10,columns(tmp)); + be([1,2,3,10],:) = tmp; + ## Mesh tetrahedra + t = load("msh_t.txt")'; + + + ## Remove hanging nodes + if (verbose) + printf("Check for hanging nodes...\n"); + endif + nnodes = columns(p); + in_msh = intersect( 1:nnodes , t(1:4,:) ); + if length(in_msh) != nnodes + new_num(in_msh) = [1:length(in_msh)]; + t(1:4,:) = new_num(t(1:4,:)); + be(1:3,:) = new_num(be(1:3,:)); + p = p(:,in_msh); + endif + + mesh = struct("p",p,"s",s,"e",be,"t",t); + + if (verbose) + printf("Deleting temporary files...\n"); + endif + system(["rm -f msh_p.txt msh_e.txt msh_t.txt msh_s.txt *.msh"]); + +endfunction \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_join_structured_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,161 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh3m_join_structured_mesh(@var{mesh1},@var{mesh2},@var{s1},@var{s2}) +## +## Join the two structured meshes @var{mesh1} and @var{mesh2} into one +## single mesh. +## +## The two meshes must share a common face identified by @var{s1} and +## @var{s2}. +## +## @strong{WARNING}: the two meshes must share the same vertexes on the +## common face. +## +## @seealso{msh3m_structured_mesh, msh3m_gmsh, msh3m_submesh, +## msh2m_join_structured_mesh} +## @end deftypefn + +function mesh = msh3m_join_structured_mesh(mesh1,mesh2,s1,s2) + + ## Check input + if nargin != 4 # Number of input parameters + error("msh3m_join_structured_mesh: wrong number of input parameters."); + elseif !(isstruct(mesh1) && isfield(mesh1,"p") && + isfield (mesh1,"e") && isfield(mesh1,"t") && + isstruct(mesh2) && isfield(mesh2,"p") && + isfield (mesh2,"e") && isfield(mesh2,"t") ) + error("msh3m_join_structured_mesh: invalid mesh structure passed as input."); + elseif !(isvector(s1) && isvector(s2)) + error("msh3m_join_structured_mesh: shared geometrical sides are not vectors."); + elseif (length(s1) != length(s2)) + error("msh3m_join_structured_mesh: vectors containing shared geometrical sides are not of the same length."); + endif + + ## Join meshes + + ## Make sure that the outside world is always on the same side of the + ## boundary of mesh1 + [mesh1.e(8:9,:),I] = sort(mesh1.e(8:9,:)); + + ## IF THE REGIONS ARE INVERTED THE VERTEX ORDER SHOULD ALSO BE + ## INVERTED!! + + ## FIXME: here a check could be added to see whether + ## the coordinate points of the two meshes coincide on the + ## side edges + + ## Get interface nodes + intfcnodes1 = msh3m_nodes_on_faces(mesh1,s1)'; + intfcnodes2 = msh3m_nodes_on_faces(mesh2,s2)'; + + ## Sort interface nodes by position + [tmp,I] = sort(mesh1.p(1,intfcnodes1)); + intfcnodes1 = intfcnodes1(I); + [tmp,I] = sort(mesh1.p(2,intfcnodes1)); + intfcnodes1 = intfcnodes1(I); + [tmp,I] = sort(mesh1.p(3,intfcnodes1)); + intfcnodes1 = intfcnodes1(I); + + [tmp,I] = sort(mesh2.p(1,intfcnodes2)); + intfcnodes2 = intfcnodes2(I); + [tmp,I] = sort(mesh2.p(2,intfcnodes2)); + intfcnodes2 = intfcnodes2(I); + [tmp,I] = sort(mesh2.p(3,intfcnodes2)); + intfcnodes2 = intfcnodes2(I); + + ## Delete redundant boundary faces but first remeber what region they + ## were connected to + for is = 1:length(s2) + ii = find( mesh2.e(10,:)==s2(is) ); + adreg(is,:) = unique(mesh2.e(9,ii)); + endfor + + for is = 1:length(s2) + mesh2.e(:,find( mesh2.e(10,:)==s2(is) )) = []; + endfor + + ## Change face numbers + idx = []; + consecutives = []; + idx = unique(mesh2.e(10,:)); + consecutives (idx) = [1:length(idx)] + max(mesh1.e(10,:)); + mesh2.e(10,:) = consecutives(mesh2.e(10,:)); + + ## Change node indices in connectivity matrix and edge list + idx = []; + consecutives = []; + idx = 1:size(mesh2.p,2); + offint = setdiff(idx,intfcnodes2); + consecutives (offint) = [1:length(offint)]+size(mesh1.p,2); + + consecutives (intfcnodes2) = intfcnodes1; + mesh2.e(1:3,:) = consecutives(mesh2.e(1:3,:)); + mesh2.t(1:4,:) = consecutives(mesh2.t(1:4,:)); + + ## Delete redundant points + mesh2.p(:,intfcnodes2) = []; + + ## Set region numbers + regions = unique(mesh1.t(5,:));# Mesh 1 + newregions(regions) = 1:length(regions); + mesh1.t(5,:) = newregions(mesh1.t(5,:)); + + regions = unique(mesh2.t(5,:));# Mesh 2 + newregions(regions) = [1:length(regions)]+max(mesh1.t(5,:)); + mesh2.t(5,:) = newregions(mesh2.t(5,:)); + + ## Set adjacent region numbers in face structure 2 + [i,j] = find(mesh2.e(8:9,:)); + i += 7; + + mesh2.e(i,j) = newregions(mesh2.e(i,j)); + + ## Set adjacent region numbers in edge structure 1 + for is = 1:length(s1) + ii = find( mesh1.e(10,:)==s1(is) ); + mesh1.e(8,ii) = newregions(regions(adreg(is,:))); + endfor + + ## Build new mesh structure + mesh.p = [mesh1.p mesh2.p]; + mesh.e = [mesh1.e mesh2.e]; + mesh.t = [mesh1.t mesh2.t]; + +endfunction + +%!shared mesh1,mesh2,jmesh +% x = y = z = linspace(0,1,2); +% x2 = linspace(1,2,2); +% [mesh1] = msh3m_structured_mesh(x,y,z,1,1:6); +% [mesh2] = msh3m_structured_mesh(x2,y,z,3,1:6); +% [jmesh] = msh3m_join_structured_mesh(mesh1,mesh2,2,1); +%!test +% assert(columns(jmesh.p),12) +%!test +% tmp = sort(unique(jmesh.e(10,:))); +% assert(tmp,1:11) +%!test +% assert(columns(jmesh.t),columns(mesh1.t)+columns(mesh2.t)) +%!test +% assert(unique(jmesh.e(8:9,:)),0:2)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_nodes_on_faces.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,70 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{nodelist}]} = @ +## msh3m_nodes_on_faces(@var{mesh},@var{facelist}) +## +## Return a list of @var{mesh} nodes lying on the faces specified in +## @var{facelist}. +## +## @seealso{msh3m_geometrical_properties, msh2m_nodes_on_faces} +## @end deftypefn + +function [nodelist] = msh3m_nodes_on_faces(mesh,facelist); + + ## Check input + if nargin != 2 # Number of input parameters + error("msh3m_nodes_on_faces: wrong number of input parameters."); + elseif !(isstruct(mesh) && isfield(mesh,"p") && + isfield(mesh,"t") && isfield(mesh,"e")) + error("msh3m_nodes_on_faces: first input is not a valid mesh structure."); + elseif !isnumeric(facelist) + error("msh3m_nodes_on_faces: only numeric value admitted as facelist."); + endif + + ## Search nodes + facefaces = []; + + for ii=1:length(facelist) + facefaces = [facefaces,find(mesh.e(10,:)==facelist(ii))]; + endfor + + facenodes = mesh.e(1:3,facefaces); + nodelist = unique(facenodes(:)); + +endfunction + +%!shared x,y,z,mesh +% x = y = z = linspace(0,1,2); +% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6); +%!test +% nodelist = msh3m_nodes_on_faces(mesh,1); +% assert(nodelist,[1 2 5 6]') +%!test +% nodelist = msh3m_nodes_on_faces(mesh,2); +% assert(nodelist,[3 4 7 8]') +%!test +% nodelist = msh3m_nodes_on_faces(mesh,3); +% assert(nodelist,[1 3 5 7]') +%!test +% nodelist = msh3m_nodes_on_faces(mesh,[1 2 3]); +% assert(nodelist,[1:8]') \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_structured_mesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,230 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{mesh}]} = @ +## msh3m_structured_mesh(@var{x},@var{y},@var{z},@var{region},@var{sides}) +## +## Construct a structured tetrahedral 3D mesh on a parallelepipedal +## domain. +## +## @itemize @bullet +## @item @var{x}, @var{y} and @var{z} are the one dimensional mesh +## vector of the corresponding Cartesian axis. +## @item @var{region} is a number identifying the geometrical volume, +## while @var{sides} is a 6 components vector containing the numbers +## used to identify the geometrical face edges. +## @end itemize +## +## The returned value @var{mesh} is a PDE-tool like mesh structure +## composed of the following fields: +## @itemize @minus +## @item @var{p}: matrix with size 3 times number of mesh points. +## @itemize @bullet +## @item 1st row: x-coordinates of the points. +## @item 2nd row: y-coordinates of the points. +## @item 3rd row: z-coordinates of the points. +## @end itemize +## @item @var{e}: matrix with size 10 times number of mesh face edges. +## @itemize @bullet +## @item 1st row: number of the first vertex of the face edge. +## @item 2nd row: number of the second vertex of the face edge. +## @item 3rd row: number of the third vertex of the face edge. +## @item 4th row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 5th row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 6th row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 7th row: set to 0, present for compatibility with MatLab PDE-tool. +## @item 8th row: number of the geometrical volume to the right of the +## face edge. +## @item 9th row: number of the geometrical volume to the left of the +## face edge. +## @item 10th row: number of the geometrical border containing the face +## edge. +## @end itemize +## @item @var{t}: matrix with size 5 times number of mesh elements. +## @itemize @bullet +## @item 1st row: number of the first vertex of the element. +## @item 2nd row: number of the second vertex of the element. +## @item 3rd row: number of the third vertex of the element. +## @item 4th row: number of the fourth vertex of the element. +## @item 5th row: number of the geometrical volume containing the element. +## @end itemize +## @end itemize +## +## @seealso{msh2m_structured_mesh, msh3m_gmsh, msh2m_mesh_along_spline, +## msh3m_join_structured_mesh, msh3m_submesh} +## @end deftypefn + +function [mesh] = msh3m_structured_mesh(x,y,z,region,sides) + + ## Check input + if (nargin != 5) # Number of input parameters + error("msh3m_structured_mesh: wrong number of input parameters."); + elseif !(isvector(x) && isnumeric(x) && + isvector(y) && isnumeric(y) && + isvector(z) && isnumeric(z) ) + error("msh3m_structured_mesh: X, Y, Z must be valid numeric vectors."); + elseif !isscalar(region) + error("msh3m_structured_mesh: REGION must be a valid scalar."); + elseif !(isvector(sides) && (length(sides) == 4)) + error("msh3m_structured_mesh: SIDES must be a 4 components vector."); + endif + + ## Build mesh + ## Sort point coordinates + x = sort(x); + y = sort(y); + z = sort(z); + ## Compute # of points in each direction + nx = length(x); + ny = length(y); + nz = length(z); + + ## Generate vertices + [XX,YY,ZZ] = meshgrid(x,y,z); + p = [XX(:),YY(:),ZZ(:)]'; + + iiv (ny,nx,nz)=0; + iiv(:)=1:nx*ny*nz; + iiv(end,:,:)=[]; + iiv(:,end,:)=[]; + iiv(:,:,end)=[]; + iiv=iiv(:)'; + + ## Generate connections: + + n1 = iiv; # bottom faces + n2 = iiv + 1; + n3 = iiv + ny; + n4 = iiv + ny + 1; + + N1 = iiv + nx * ny; # top faces + N2 = N1 + 1; + N3 = N1 + ny; + N4 = N3 + 1; + + t = [... + [n1; n3; n2; N2],... + [N1; N2; N3; n3],... + [N1; N2; n3; n1],... + [N2; n3; n2; n4],... + [N3; n3; N2; N4],... + [N4; n3; N2; n4],... + ]; + + ## Generate boundary face list + + ## left + T = t; + T(:) = p(1,t)'==x(1); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e1(:,jj) = t(order(:,ii(jj)),ii(jj)); + endfor + e1(10,:) = sides(1); + + ## right + T(:) = p(1,t)'==x(end); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e2(:,jj) = t(order(:,ii(jj)),ii(jj)); + end + e2(10,:) = sides(2); + + ## front + T(:) = p(2,t)'==y(1); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e3(:,jj) = t(order(:,ii(jj)),ii(jj)); + endfor + e3(10,:) = sides(3); + + ## back + T(:) = p(2,t)'==y(end); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e4(:,jj) = t(order(:,ii(jj)),ii(jj)); + endfor + e4(10,:) = sides(4); + + ## bottom + T = t; + T(:) = p(3,t)'==z(1); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e5(:,jj) = t(order(:,ii(jj)),ii(jj)); + endfor + e5(10,:) = sides(5); + + ## top + T = t; + T(:) = p(3,t)'==z(end); + [ignore,order] = sort(T,1); + ii = (find(sum(T,1)==3)); + order(1,:) = []; + for jj=1:length(ii) + e6(:,jj) = t(order(:,ii(jj)),ii(jj)); + endfor + e6(10,:) = sides(6); + + ## Assemble structure + mesh.e = [e1,e2,e3,e4,e5,e6]; + mesh.t = t; + mesh.e (9,:) = region; + mesh.t (5,:) = region; + mesh.p = p; + +endfunction + +%!test +% x = y = z = linspace(0,1,2) +% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6) +% assert = (columns(mesh.p),8) +% assert = (columns(mesh.t),6) +% assert = (columns(mesh.e),12) +%!test +% x = y = z = linspace(0,1,3) +% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6) +% assert = (columns(mesh.p),27) +% assert = (columns(mesh.t),48) +% assert = (columns(mesh.e),48) +%!test +% x = y = z = linspace(0,1,4) +% [mesh] = msh3m_structured_mesh(x,y,z,1,1:6) +% assert = (columns(mesh.p),64) +% assert = (columns(mesh.t),162) +% assert = (columns(mesh.e),108) +%!test +% x = y = z = linspace(0,1,1) +% fail([mesh] = msh3m_structured_mesh(x,y,z,1,1:6),"warning","x, y, z cannot be scalar numbers!") +%!test +% x = y = z = eye(2) +% fail([mesh] = msh3m_structured_mesh(x,y,z,1,1:6),"warning","x, y, z cannot be matrices!") \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/msh/inst/msh3m_submesh.m Mon Feb 01 07:38:09 2010 +0000 @@ -0,0 +1,92 @@ +## Copyright (C) 2006,2007,2008,2009,2010 Carlo de Falco, Massimiliano Culpo +## +## This file is part of: +## MSH - Meshing Software Package for Octave +## +## 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, see <http://www.gnu.org/licenses/>. +## +## author: Carlo de Falco <cdf _AT_ users.sourceforge.net> +## author: Massimiliano Culpo <culpo _AT_ users.sourceforge.net> + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{omesh},@var{nodelist},@var{elementlist}]} = @ +## msh3m_submesh(@var{imesh},@var{intrfc},@var{sdl}) +## +## Extract the subdomain(s) in @var{sdl} from @var{imesh}. +## +## The row vector @var{intrfc} contains the internal interface sides to +## be maintained (field @code{mesh.e(5,:)}). It can be empty. +## +## Return the vectors @var{nodelist} and @var{elementlist} containing +## respectively the list of nodes and elements of the original mesh that +## are part of the selected subdomain(s). +## +## @seealso{msh3m_join_structured_mesh, msh2m_join_structured_mesh, +## msh3m_submesh} +## @end deftypefn + +function [omesh,nodelist,elementlist] = msh3m_submesh(imesh,intrfc,sdl) + + ## Check input + if nargin != 3 + error("msh3m_submesh: wrong number of input parameters."); + elseif !(isstruct(imesh) && isfield(imesh,"p") && + isfield (imesh,"t") && isfield(imesh,"e")) + error("msh3m_submesh: first input is not a valid mesh structure."); + elseif !isvector(sdl) + error("msh3m_submesh: third input is not a valid vector."); + endif + + ## Extract sub-mesh + + ## Build element list + elementlist=[]; + for ir = 1:length(sdl) + elementlist = [ elementlist find(imesh.t(5,:)==sdl(ir)) ]; + endfor + + ## Build nodelist + nodelist = reshape(imesh.t(1:4,elementlist),1,[]); + nodelist = unique(nodelist); + + ## Extract submesh + omesh.p = imesh.p (:,nodelist); + indx(nodelist) = 1:length (nodelist); + omesh.t = imesh.t (:,elementlist); + omesh.t(1:4,:) = indx(omesh.t(1:4,:)); + + omesh.e = []; + for ifac = 1:size(imesh.e,2) + if (length(intersect(imesh.e(1:3,ifac),nodelist) )== 3) + omesh.e = [omesh.e imesh.e(:,ifac)]; + endif + endfor + + omesh.e(1:3,:) = indx(omesh.e(1:3,:)); + +endfunction + +%!shared mesh1,mesh2,jmesh,exmesh,nodelist,elemlist +% x = y = z = linspace(0,1,2); +% x2 = linspace(1,2,2); +% [mesh1] = msh3m_structured_mesh(x,y,z,1,1:6); +% [mesh2] = msh3m_structured_mesh(x2,y,z,1,1:6); +% [jmesh] = msh3m_join_structured_mesh(mesh1,mesh2,2,1); +% [exmesh,nodelist,elemlist] = msh3m_submesh(jmesh,2,1); +%!test +% assert(size(exmesh.p),size(mesh1.p)) +%!test +% assert(size(exmesh.t),size(mesh1.t)) +%!test +% assert(size(exmesh.e),size(mesh1.e)) \ No newline at end of file