Index: sm/trunk/src/m/model/mesh/mesh.m =================================================================== --- /issm/trunk/src/m/model/mesh/mesh.m (revision 9703) +++ (revision ) @@ -1,85 +1,0 @@ -function md=setmesh(md,domainname,varargin) -%MESH - create model mesh -% -% This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution -% where md is a @model object, domainname is the name of an Argus domain outline file, -% and resolution is a characteristic length for the mesh (same unit as the domain outline -% unit). Riftname is an optional argument (Argus domain outline) describing rifts. -% -% Usage: -% md=setmesh(md,domainname,resolution) -% or md=setmesh(md,domainname,riftname, resolution) -% -% Examples: -% md=setmesh(md,'DomainOutline.exp',1000); -% md=setmesh(md,'DomainOutline.exp','Rifts.exp',1500); - -%Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m resolution node would -%be made of 1000*1000 area squares). -if (nargin==3), - resolution=varargin{1}; - riftname=''; -end -if (nargin==4), - riftname=varargin{1}; - resolution=varargin{2}; -end - -%Check that mesh was not already run, and warn user: -if md.numberofelements~=0, - choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s'); - if ~strcmp(choice,'y') - disp('no meshing done ... exiting'); - return - end -end - -area=resolution^2; - -%Mesh using TriMesh -if strcmp(riftname,''), - [md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes'); -else - [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes'); - - %check that all the created nodes belong to at least one element - orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); - for i=1:length(orphan), - %get rid of the orphan node i - %update x and y - x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)]; - y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)]; - %update elements - pos=find(elements>orphan(i)-(i-1)); - elements(pos)=elements(pos)-1; - %update segments - pos1=find(segments(:,1)>orphan(i)-(i-1)); - pos2=find(segments(:,2)>orphan(i)-(i-1)); - segments(pos1,1)=segments(pos1,1)-1; - segments(pos2,2)=segments(pos2,2)-1; - end - - %plug into md - md.x=x; - md.y=y; - md.elements=elements; - md.segments=segments; - md.segmentmarkers=segmentmarkers; -end - -%Fill in rest of fields: -md.numberofelements=length(md.elements); -md.numberofnodes=length(md.x); -md.z=zeros(md.numberofnodes,1); -md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1; -md.nodeonbed=ones(md.numberofnodes,1); -md.nodeonsurface=ones(md.numberofnodes,1); -md.elementonbed=ones(md.numberofelements,1); -md.elementonsurface=ones(md.numberofelements,1); - -%Now, build the connectivity tables for this mesh. -md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes); -md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity); - -%type of model -md.dim=2; Index: /issm/trunk/src/m/model/mesh/setmesh.m =================================================================== --- /issm/trunk/src/m/model/mesh/setmesh.m (revision 9704) +++ /issm/trunk/src/m/model/mesh/setmesh.m (revision 9704) @@ -0,0 +1,85 @@ +function md=setmesh(md,domainname,varargin) +%MESH - create model mesh +% +% This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution +% where md is a @model object, domainname is the name of an Argus domain outline file, +% and resolution is a characteristic length for the mesh (same unit as the domain outline +% unit). Riftname is an optional argument (Argus domain outline) describing rifts. +% +% Usage: +% md=setmesh(md,domainname,resolution) +% or md=setmesh(md,domainname,riftname, resolution) +% +% Examples: +% md=setmesh(md,'DomainOutline.exp',1000); +% md=setmesh(md,'DomainOutline.exp','Rifts.exp',1500); + +%Figure out a characteristic area. Resolution is a node oriented concept (ex a 1000m resolution node would +%be made of 1000*1000 area squares). +if (nargin==3), + resolution=varargin{1}; + riftname=''; +end +if (nargin==4), + riftname=varargin{1}; + resolution=varargin{2}; +end + +%Check that mesh was not already run, and warn user: +if md.numberofelements~=0, + choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s'); + if ~strcmp(choice,'y') + disp('no meshing done ... exiting'); + return + end +end + +area=resolution^2; + +%Mesh using TriMesh +if strcmp(riftname,''), + [md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes'); +else + [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes'); + + %check that all the created nodes belong to at least one element + orphan=find(~ismember([1:length(x)],sort(unique(elements(:))))); + for i=1:length(orphan), + %get rid of the orphan node i + %update x and y + x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)]; + y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)]; + %update elements + pos=find(elements>orphan(i)-(i-1)); + elements(pos)=elements(pos)-1; + %update segments + pos1=find(segments(:,1)>orphan(i)-(i-1)); + pos2=find(segments(:,2)>orphan(i)-(i-1)); + segments(pos1,1)=segments(pos1,1)-1; + segments(pos2,2)=segments(pos2,2)-1; + end + + %plug into md + md.x=x; + md.y=y; + md.elements=elements; + md.segments=segments; + md.segmentmarkers=segmentmarkers; +end + +%Fill in rest of fields: +md.numberofelements=length(md.elements); +md.numberofnodes=length(md.x); +md.z=zeros(md.numberofnodes,1); +md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1; +md.nodeonbed=ones(md.numberofnodes,1); +md.nodeonsurface=ones(md.numberofnodes,1); +md.elementonbed=ones(md.numberofelements,1); +md.elementonsurface=ones(md.numberofelements,1); + +%Now, build the connectivity tables for this mesh. +md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes); +md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity); + +%type of model +md.dim=2;