Index: /issm/trunk/src/m/classes/public/modelextract2.m
===================================================================
--- /issm/trunk/src/m/classes/public/modelextract2.m	(revision 1348)
+++ /issm/trunk/src/m/classes/public/modelextract2.m	(revision 1348)
@@ -0,0 +1,243 @@
+function md2=modelextract2(md1,area,varargin)
+%modelextract2 - extract a model according to an Argus contour or flag list
+%
+%   This routine extracts a submodel from a bigger model with respect to a given contour
+%   md must be followed by the corresponding exp file or flags list
+%   It can either be a domain file (argus type, .exp extension), or an array of element flags. 
+%   If user wants every element outside the domain to be 
+%   modelextract2d, add '~' to the name of the domain file (ex: '~Pattyn.exp');
+%   an empty string '' will be considered as an empty domain
+%   a string 'all' will be considered as the entire domain
+%   add an argument 0 if you do not want the elements to be checked (faster)
+%
+%   Usage:
+%      md2=modelextract2(md1,area,varargin);
+%
+%   Examples:
+%      md2=modelextract2(md,'Domain.exp');
+%      md2=modelextract2(md,'Domain.exp',1);
+%      md2=modelextract2(md,md.elementoniceshelf);
+%
+%   See also: EXTRUDE, COLLAPSE
+
+%some checks
+if ((nargin~=2 & nargin~=3) | (nargout~=1)),
+	help modelextract2
+	error('modelextract2 error message: bad usage');
+end
+if strcmpi(md1.type,'3d'),
+	error('modelextract2 error message: only 2d model supported yet. Use BasinConstrain instead');
+end
+
+%get check option
+if (nargin==3 & varargin{1}==0),
+	checkoutline=0;
+else
+	checkoutline=1;
+end
+
+%first
+if ischar(area),
+	if isempty(area),
+		flag_elem=zeros(md1.numberofelements,1);
+		invert=0;
+	elseif strcmpi(area,'all')
+		flag_elem=ones(md1.numberofelements,1);
+		invert=0;
+	else
+		%make sure that we actually don't want the elements outside the domain outline!
+		if strcmpi(area(1),'~'),
+			area=area(2:length(area));
+			invert=1;
+		else
+			invert=0;
+		end
+
+		%ok, flag_elem elements
+		A=expread(area,1);
+		flag_elem=ContourToMesh(md1.elements(:,1:3),md1.x,md1.y,A,'element',1);
+	end
+	if invert, flag_elem=~flag_elem; end
+
+elseif isfloat(area),
+	if size(area,1)~=md1.numberofelements,
+		setelementstypeusage();
+		error('Flags must be of same size as number of elements in model');
+	end
+	flag_elem=area;
+else
+	error('Invalide option');
+end
+
+%element and grid position
+pos_elem=find(flag_elem);
+pos_grid=sort(unique(md1.elements(pos_elem,:)));
+
+%keep track of some fields
+numberofgrids1=md1.numberofgrids;
+numberofelements1=md1.numberofelements;
+numberofgrids2=length(pos_grid);
+numberofelements2=length(pos_elem);
+
+%Create Pelem and Pgrid (transform old grids in new grids and same thing for the elements)
+Pelem=zeros(numberofelements1,1);
+Pelem(pos_elem)=[1:numberofelements2]';
+Pgrid=zeros(numberofgrids1,1);
+Pgrid(pos_grid)=[1:numberofgrids2]';
+
+%renumber the elements (some grid won't exist anymore)
+elements_1=md1.elements;
+elements_2=elements_1(pos_elem,:);
+elements_2(:,1)=Pgrid(elements_2(:,1));
+elements_2(:,2)=Pgrid(elements_2(:,2));
+elements_2(:,3)=Pgrid(elements_2(:,3));
+if strcmpi(md1.type,'3d'),
+	elements_2(:,4)=Pgrid(elements_2(:,4));
+	elements_2(:,5)=Pgrid(elements_2(:,5));
+	elements_2(:,6)=Pgrid(elements_2(:,6));
+end
+
+%OK, now create the new model !
+
+	%take every fields from model
+	md2=md1;
+
+	%modify some specific fields
+
+	%Mesh
+	md2.numberofelements=numberofelements2;
+	md2.numberofgrids=numberofgrids2;
+	md2.elements=elements_2;
+
+	%Initial 2d mesh 
+	if strcmpi(md1.type,'3d')
+		flag_elem_2d=flag_elem(1:md1.numberofelements2d);
+		pos_elem_2d=find(flag_elem_2d);
+		flag_grid_2d=flag_grid(1:md1.numberofgrids2d);
+		pos_grid_2d=find(flag_grid_2d);
+
+		md2.numberofelements2d=length(pos_elem_2d);
+		md2.numberofgrids2d=length(pos_grid_2d);
+		md2.elements2d=md1.elements2d(pos_elem_2d,:);
+		md2.elements2d(:,1)=Pgrid(md2.elements2d(:,1));
+		md2.elements2d(:,2)=Pgrid(md2.elements2d(:,2));
+		md2.elements2d(:,3)=Pgrid(md2.elements2d(:,3));
+
+		if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d,:); end;
+		md2.x2d=md1.x(pos_grid_2d);
+		md2.y2d=md1.y(pos_grid_2d);
+		md2.z2d=md1.z(pos_grid_2d);
+	end
+
+	%Penalties
+	if ~isnan(md2.penalties),
+		for i=1:size(md1.penalties,1);
+			md2.penalties(i,:)=Pgrid(md1.penalties(i,:));
+		end
+		md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
+	end
+
+	%recreate segments
+	md2.segments=findsegments(md2);
+	md2.gridonboundary=zeros(numberofgrids2,1); md1.gridonboundary(md2.segments(:,1:2))=1;
+
+	if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
+
+	%Dirichlets Diag
+	%Catch the elements that have not been extracted
+	orphans_elem=find(~flag_elem);
+	orphans_grid=unique(md1.elements(orphans_elem,:))';
+	%Figure out which grid are on the boundary between md2 and md1
+	gridstoflag1=intersect(orphans_grid,pos_grid);
+	gridstoflag2=Pgrid(gridstoflag1);
+	if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
+		md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
+		md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
+	else
+		md2.dirichletvalues_diag=zeros(numberofgrids2,2);
+		disp(' ')
+		disp('!! modelextract warning: dirichlet values should be checked !!')
+		disp(' ')
+	end
+
+	%Diagnostic
+	if ~isnan(md2.segmentonneumann_diag)
+		md2.segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1)); 
+		md2.segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2)); 
+		md2.segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end)); 
+		if strcmpi(md1.type,'3d')
+			md2.segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3)); 
+			md2.segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4)); 
+		end
+		md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);
+	end
+	md2.neumannvalues_diag=NaN;
+
+	%Transient
+	if ~isnan(md2.segmentonneumann_prog)
+		md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1)); 
+		md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2)); 
+		md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end)); 
+		md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:);
+	end
+	md2.neumannvalues_prog=NaN;
+	if ~isnan(md2.segmentonneumann_prog2)
+		md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1)); 
+		md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2)); 
+		md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end)); 
+		md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:);
+	end
+
+	%Results fields
+	if isstruct(md1.results),
+		solutionfields=fields(md1.results);
+		for i=1:length(solutionfields),
+			%get subfields
+			solutionsubfields=fields(md1.results.(solutionfields(i)));
+			for j=1:length(solutionsubfields),
+				field=md1.results.(solutionfields(i)).(solutionsubfields(j));
+				if length(field)==numberofgrids1,
+					md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_grid);
+				elseif length(field)==numberofelements1,
+					md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_elem);
+				end
+			end
+		end
+	end
+
+	%reinitialize output parameters
+	md2.viscousheating=NaN;
+	md2.pressure_elem=NaN;
+	md2.stress=NaN;
+	md2.stress_surface=NaN;
+	md2.stress_bed=NaN;
+	md2.deviatoricstress=NaN;
+	md2.strainrate=NaN;
+
+	%other fields
+
+	%loop over model fields
+	model_fields=fields(md1);
+	for i=1:length(model_fields),
+
+		%get field
+		field=md1.(model_fields(i));
+
+		%size = number of grids * 1
+		if (size(field,1)==numberofgrids1 & size(field,2)==1),
+			md2.(model_fields(i))=field(pos_grid);
+
+		%size = number of grids * 2
+		elseif (size(field,1)==numberofgrids1 & size(field,2)==2), 
+			md2.(model_fields(i))=[field(pos_grid,1) field(pos_grid,2)];
+
+		%size = number of elements * 1
+		elseif (size(field,1)==numberofelements1 & size(field,2)==1),
+			md2.(model_fields(i))=field(pos_elem);
+
+		end
+	end
+
+%Keep track of pos_grid and pos_elem
+md2.extractedgrids=pos_grid;
+md2.extractedelements=pos_elem;