mesh

PURPOSE ^

MESH - create model mesh

SYNOPSIS ^

function md=mesh(md,domainname,varargin)

DESCRIPTION ^

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=mesh(md,domainname,resolution)
   or md=mesh(md,domainname,riftname, resolution)

   Examples:
      md=mesh(md,'DOmainOutline.exp',1000);
      md=mesh(md,'DOmainOutline.exp','Rifts.exp',1500);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function md=mesh(md,domainname,varargin)
0002 %MESH - create model mesh
0003 %
0004 %   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
0005 %   where md is a @model object, domainname is the name of an Argus domain outline file,
0006 %   and resolution is a characteristic length for the mesh (same unit as the domain outline
0007 %   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
0008 %
0009 %   Usage:
0010 %      md=mesh(md,domainname,resolution)
0011 %   or md=mesh(md,domainname,riftname, resolution)
0012 %
0013 %   Examples:
0014 %      md=mesh(md,'DOmainOutline.exp',1000);
0015 %      md=mesh(md,'DOmainOutline.exp','Rifts.exp',1500);
0016 
0017 %Figure out a characteristic area. Resolution is a grid oriented concept (ex a 1000m  resolution grid would
0018 %be made of 1000*1000 area squares).
0019 if (nargin==3),
0020     resolution=varargin{1};
0021     riftname='';
0022 end
0023 if (nargin==4),
0024     riftname=varargin{1};
0025     resolution=varargin{2};
0026 end
0027 
0028 %Check that mesh was not already run, and warn user:
0029 if subsref(md,struct('type','.','subs','counter'))>=1,
0030     choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s');
0031     if ~strcmp(choice,'y')
0032         error('no meshing done ... exiting');
0033     end
0034 end
0035 
0036 area=resolution^2;
0037 
0038 %Mesh using TriMesh
0039 if strcmp(riftname,''),
0040     [md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes');
0041 else
0042     [elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes');
0043 
0044     %check that all the created grids belong to at least one element
0045     orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
0046     for i=1:length(orphan),
0047         %get rid of the orphan grid i
0048         %update x and y
0049         x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)];
0050         y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)];
0051         %update elements
0052         pos=find(elements>orphan(i)-(i-1));
0053         elements(pos)=elements(pos)-1;
0054         %update segments
0055         pos1=find(segments(:,1)>orphan(i)-(i-1));
0056         pos2=find(segments(:,2)>orphan(i)-(i-1));
0057         segments(pos1,1)=segments(pos1,1)-1;
0058         segments(pos2,2)=segments(pos2,2)-1;
0059     end
0060 
0061     %plug into md
0062     md.x=x;
0063     md.y=y;
0064     md.elements=elements;
0065     md.segments=segments;
0066     md.segmentmarkers=segmentmarkers;
0067 end
0068 
0069 %Fill in rest of fields:
0070 md.numberofelements=length(md.elements);
0071 md.numberofgrids=length(md.x);
0072 md.z=zeros(md.numberofgrids,1);
0073 md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
0074 md.gridonbed=ones(md.numberofgrids,1);
0075 md.gridonsurface=ones(md.numberofgrids,1);
0076 md.elementonbed=ones(md.numberofelements,1);
0077 md.elementonsurface=ones(md.numberofelements,1);
0078 
0079 %outline names
0080 md.domainoutline=readfile(domainname);
0081 if strcmp(riftname,''),
0082     md.riftoutline='';
0083 else
0084     md.riftoutline=readfile(riftname);
0085 end
0086 
0087 %type of model
0088 md.type='2d';
0089     
0090 %augment counter  keeping track of what has been done to this model
0091 md.counter=1;

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003