Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[6818] trunk/octave-forge/extra/msh
|
From: <cu...@us...> - 2010-02-01 07:38:17
|
Revision: 6818
http://octave.svn.sourceforge.net/octave/?rev=6818&view=rev
Author: culpo
Date: 2010-02-01 07:38:09 +0000 (Mon, 01 Feb 2010)
Log Message:
-----------
Updated OF package to new naming convention. Deleted old functions. Updated INDEX and DESCRIPTION.
Modified Paths:
--------------
trunk/octave-forge/extra/msh/DESCRIPTION
trunk/octave-forge/extra/msh/INDEX
Added Paths:
-----------
trunk/octave-forge/extra/msh/inst/msh2m_displacement_smoothing.m
trunk/octave-forge/extra/msh/inst/msh2m_equalize_mesh.m
trunk/octave-forge/extra/msh/inst/msh2m_geometrical_properties.m
trunk/octave-forge/extra/msh/inst/msh2m_gmsh.m
trunk/octave-forge/extra/msh/inst/msh2m_jiggle_mesh.m
trunk/octave-forge/extra/msh/inst/msh2m_join_structured_mesh.m
trunk/octave-forge/extra/msh/inst/msh2m_mesh_along_spline.m
trunk/octave-forge/extra/msh/inst/msh2m_nodes_on_sides.m
trunk/octave-forge/extra/msh/inst/msh2m_structured_mesh.m
trunk/octave-forge/extra/msh/inst/msh2m_submesh.m
trunk/octave-forge/extra/msh/inst/msh2m_topological_properties.m
trunk/octave-forge/extra/msh/inst/msh2p_mesh.m
trunk/octave-forge/extra/msh/inst/msh3e_surface_mesh.m
trunk/octave-forge/extra/msh/inst/msh3m_geometrical_properties.m
trunk/octave-forge/extra/msh/inst/msh3m_gmsh.m
trunk/octave-forge/extra/msh/inst/msh3m_join_structured_mesh.m
trunk/octave-forge/extra/msh/inst/msh3m_nodes_on_faces.m
trunk/octave-forge/extra/msh/inst/msh3m_structured_mesh.m
trunk/octave-forge/extra/msh/inst/msh3m_submesh.m
Removed Paths:
-------------
trunk/octave-forge/extra/msh/inst/MSH2Mdisplacementsmoothing.m
trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m
trunk/octave-forge/extra/msh/inst/MSH2Mgeomprop.m
trunk/octave-forge/extra/msh/inst/MSH2Mgmsh.m
trunk/octave-forge/extra/msh/inst/MSH2Mjigglemesh.m
trunk/octave-forge/extra/msh/inst/MSH2Mjoinstructm.m
trunk/octave-forge/extra/msh/inst/MSH2Mmeshalongspline.m
trunk/octave-forge/extra/msh/inst/MSH2Mnodesonsides.m
trunk/octave-forge/extra/msh/inst/MSH2Mstructmesh.m
trunk/octave-forge/extra/msh/inst/MSH2Msubmesh.m
trunk/octave-forge/extra/msh/inst/MSH2Mtopprop.m
trunk/octave-forge/extra/msh/inst/MSH3Mgeomprop.m
trunk/octave-forge/extra/msh/inst/MSH3Mgmsh.m
trunk/octave-forge/extra/msh/inst/MSH3Mjoinstructm.m
trunk/octave-forge/extra/msh/inst/MSH3Mnodesonfaces.m
trunk/octave-forge/extra/msh/inst/MSH3Mstructmesh.m
trunk/octave-forge/extra/msh/inst/MSH3Msubmesh.m
Modified: trunk/octave-forge/extra/msh/DESCRIPTION
===================================================================
--- trunk/octave-forge/extra/msh/DESCRIPTION 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/DESCRIPTION 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Modified: trunk/octave-forge/extra/msh/INDEX
===================================================================
--- trunk/octave-forge/extra/msh/INDEX 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/INDEX 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mdisplacementsmoothing.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mdisplacementsmoothing.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mdisplacementsmoothing.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mequalizemesh.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mgeomprop.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mgeomprop.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mgeomprop.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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);
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mgmsh.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mgmsh.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mgmsh.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mjigglemesh.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mjigglemesh.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mjigglemesh.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mjoinstructm.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mjoinstructm.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mjoinstructm.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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);
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mmeshalongspline.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mmeshalongspline.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mmeshalongspline.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mnodesonsides.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mnodesonsides.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mnodesonsides.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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);
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Mstructmesh.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Mstructmesh.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Mstructmesh.m 2010-02-01 07:38:09 UTC (rev 6818)
@@ -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
Deleted: trunk/octave-forge/extra/msh/inst/MSH2Msubmesh.m
===================================================================
--- trunk/octave-forge/extra/msh/inst/MSH2Msubmesh.m 2010-02-01 03:57:21 UTC (rev 6817)
+++ trunk/octave-forge/extra/msh/inst/MSH2Msubmesh.m 2010-02-01 07:38:09 UTC (r...
[truncated message content] |
