Index: /issm/trunk/src/m/classes/public/modelextract.m
===================================================================
--- /issm/trunk/src/m/classes/public/modelextract.m	(revision 1493)
+++ /issm/trunk/src/m/classes/public/modelextract.m	(revision 1494)
@@ -1,4 +1,4 @@
-function md2=modelextract(md1,area,varargin)
-%MODELEXTRACT - extract a model according to an Argus contour or flag list
+function md2=modelextract2(md1,area)
+%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
@@ -6,5 +6,5 @@
 %   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 
-%   modelextractd, add '~' to the name of the domain file (ex: '~Pattyn.exp');
+%   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
@@ -12,20 +12,16 @@
 %
 %   Usage:
-%      md2=modelextract(md1,area,varargin);
+%      md2=modelextract2(md1,area);
 %
 %   Examples:
-%      md2=modelextract(md,'Domain.exp');
-%      md2=modelextract(md,'Domain.exp',1);
-%      md2=modelextract(md,md.elementoniceshelf);
+%      md2=modelextract2(md,'Domain.exp');
+%      md2=modelextract2(md,md.elementoniceshelf);
 %
 %   See also: EXTRUDE, COLLAPSE
 
 %some checks
-if ((nargin~=2 & nargin~=3) | (nargout~=1)),
-	help modelextract
-	error('modelextract error message: bad usage');
-end
-if strcmpi(md1.type,'3d'),
-	error('modelextract error message: only 2d model supported yet. Use BasinConstrain instead');
+if ((nargin~=2) | (nargout~=1)),
+	help modelextract2
+	error('modelextract2 error message: bad usage');
 end
 
@@ -46,5 +42,4 @@
 		invert=0;
 	else
-
 		%make sure that we actually don't want the elements outside the domain outline!
 		if strcmpi(area(1),'~'),
@@ -58,21 +53,7 @@
 		A=expread(area,1);
 		flag_elem=ContourToMesh(md1.elements(:,1:3),md1.x,md1.y,A,'element',1);
-
-		%check that the elements do not cross the contour -> flag as 0
-		if checkoutline,
-			pos=find(flag_elem);
-			for contour=1:length(A)
-				for i=1: length(pos)
-					test=cross_extract(md1.x(md1.elements(pos(i),:)),md1.y(md1.elements(pos(i),:)),A(contour).x,A(contour).y);
-					if test
-						flag_elem(pos(i))=0;
-					end
-				end
-			end
-		end
-	end
-	if invert,
-		flag_elem=~flag_elem;
-	end
+	end
+	if invert, flag_elem=~flag_elem; end
+
 elseif isfloat(area),
 	if size(area,1)~=md1.numberofelements,
@@ -84,9 +65,16 @@
 	error('Invalide option');
 end
+
+%kick out all elements with 3 dirichlets
+spc_elem=find(~flag_elem);
+spc_grid=sort(unique(md1.elements(spc_elem,:)));
+flag=ones(md1.numberofgrids,1);
+flag(spc_grid)=0;
+pos=find(sum(flag(md1.elements),2)==0);
+flag_elem(pos)=0;
+
+%extracted elements and nodes lists
 pos_elem=find(flag_elem);
-
-%get the flags of grids
-flag_grid=ismember([1:md1.numberofgrids]',sort(unique(md1.elements(pos_elem,:)))');
-pos_grid=find(flag_grid);
+pos_grid=sort(unique(md1.elements(pos_elem,:)));
 
 %keep track of some fields
@@ -95,4 +83,6 @@
 numberofgrids2=length(pos_grid);
 numberofelements2=length(pos_elem);
+flag_grid=zeros(numberofgrids1,1);
+flag_grid(pos_grid)=1;
 
 %Create Pelem and Pgrid (transform old grids in new grids and same thing for the elements)
@@ -105,5 +95,4 @@
 elements_1=md1.elements;
 elements_2=elements_1(pos_elem,:);
-
 elements_2(:,1)=Pgrid(elements_2(:,1));
 elements_2(:,2)=Pgrid(elements_2(:,2));
@@ -120,16 +109,41 @@
 	md2=md1;
 
+	%automatically modify fields
+
+	%loop over model fields
+	model_fields=fields(md1);
+	for i=1:length(model_fields),
+
+		%get field
+		field=md1.(model_fields(i));
+		fieldsize=size(field);
+
+		%size = number of grids * n
+		if fieldsize(1)==numberofgrids1
+			md2.(model_fields(i))=field(pos_grid,:);
+
+		%size = number of elements * n
+		elseif fieldsize(1)==numberofelements1
+			md2.(model_fields(i))=field(pos_elem,:);
+		end
+	end
+
+	%modify some specific fields
+
 	%Mesh
 	md2.numberofelements=numberofelements2;
 	md2.numberofgrids=numberofgrids2;
 	md2.elements=elements_2;
-	if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
-	md2.x=md1.x(pos_grid);
-	md2.y=md1.y(pos_grid);
-	md2.z=md1.z(pos_grid);
-	if ~isnan(md2.bed_slopex), md2.bed_slopex=md1.bed_slopex(pos_grid); end;
-	if ~isnan(md2.bed_slopey), md2.bed_slopex=md1.bed_slopey(pos_grid); end;
-	if ~isnan(md2.surface_slopex), md2.surface_slopex=md1.surface_slopex(pos_grid); end;
-	if ~isnan(md2.surface_slopey), md2.surface_slopex=md1.surface_slopey(pos_grid); end;
+
+	%uppernodes lowernodes
+	if strcmpi(md1.type,'3d')
+		md2.uppergrids=md1.uppergrids(pos_grid);
+		pos=find(~isnan(md2.uppergrids));
+		md2.uppergrids(pos)=Pgrid(md2.uppergrids(pos));
+
+		md2.lowergrids=md1.lowergrids(pos_grid);
+		pos=find(~isnan(md2.lowergrids));
+		md2.lowergrids(pos)=Pgrid(md2.lowergrids(pos));
+	end
 
 	%Initial 2d mesh 
@@ -152,17 +166,5 @@
 		md2.z2d=md1.z(pos_grid_2d);
 	end
-		
-	%Elements
-	if ~isnan(md2.elementonhutter), md2.elementonhutter=md1.elementonhutter(pos_elem); end;
-	if ~isnan(md2.elementonmacayeal), md2.elementonmacayeal=md1.elementonmacayeal(pos_elem); end;
-	if ~isnan(md2.elementonpattyn), md2.elementonpattyn=md1.elementonpattyn(pos_elem); end;
-	if ~isnan(md2.elementonstokes), md2.elementonstokes=md1.elementonstokes(pos_elem); end;
-
-	%Nodes
-	if ~isnan(md2.gridonhutter), md2.gridonhutter=md1.gridonhutter(pos_grid); end;
-	if ~isnan(md2.gridonmacayeal), md2.gridonmacayeal=md1.gridonmacayeal(pos_grid); end;
-	if ~isnan(md2.gridonpattyn), md2.gridonpattyn=md1.gridonpattyn(pos_grid); end;
-	if ~isnan(md2.gridonstokes), md2.gridonstokes=md1.gridonstokes(pos_grid); end;
-	
+
 	%Penalties
 	if ~isnan(md2.penalties),
@@ -172,39 +174,51 @@
 		md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
 	end
-	md2.rifts=NaN; %????????????????????????????
-	md2.numrifts=0;%????????????????????????????
-
-	%Projections
-	if ~isnan(md2.uppergrids), md2.uppergrids=Pgrid(md1.uppergrids(pos_grid)); end;
-	if ~isnan(md2.lowergrids), md2.lowergrids=Pgrid(md1.lowergrids(pos_grid)); end;
-	if ~isnan(md2.deadgrids), md2.deadgrids=md1.deadgrids(pos_grid);end;
-	
-	%Extrusion
-	md2.elementonbed=md1.elementonbed(pos_elem);
-	md2.elementonsurface=md1.elementonsurface(pos_elem);
-	md2.gridonbed=md1.gridonbed(pos_grid);
-	md2.gridonsurface=md1.gridonsurface(pos_grid);
-	if ~isnan(md2.firn_layer), md2.firn_layer=md1.firn_layer(pos_grid); end;
-	
-	%Physical parameters
-	if ~isnan(md2.drag), md2.drag=md1.drag(pos_grid);end;
-	if ~isnan(md2.p), md2.p=md1.p(pos_elem);end;
-	if ~isnan(md2.q), md2.q=md1.q(pos_elem);end;
-	if ~isnan(md2.n), md2.n=md1.n(pos_elem);end;
-	if ~isnan(md2.B), md2.B=md1.B(pos_grid);end;
-
-	%Geometrical parameters
-	if ~isnan(md2.elementoniceshelf), md2.elementoniceshelf=md1.elementoniceshelf(pos_elem);end;
-	if ~isnan(md2.elementonicesheet), md2.elementonicesheet=md1.elementonicesheet(pos_elem);end;
-	if ~isnan(md2.elementonwater), md2.elementonwater=md1.elementonwater(pos_elem);end;
-	if ~isnan(md2.gridoniceshelf), md2.gridoniceshelf=md1.gridoniceshelf(pos_grid);end;
-	if ~isnan(md2.gridonwater), md2.gridonwater=md1.gridonwater(pos_grid);end;
-	if ~isnan(md2.gridonicesheet), md2.gridonicesheet=md1.gridonicesheet(pos_grid);end;
-	if ~isnan(md2.surface), md2.surface=md1.surface(pos_grid); end;
-	if ~isnan(md2.thickness), md2.thickness=md1.thickness(pos_grid); end;
-	if ~isnan(md2.bed), md2.bed=md1.bed(pos_grid); end;
-
-	%Boundary conditions
-	md2.gridonboundary=md1.gridonboundary(pos_grid);
+
+	%recreate segments
+	if strcmpi(md1.type,'2d')
+		md2.segments=findsegments(md2);
+		md2.gridonboundary=zeros(numberofgrids2,1); md1.gridonboundary(md2.segments(:,1:2))=1;
+	end
+
+	%Boundary conditions: Dirichlets on new boundary
+	%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.gridondirichlet_diag),
+		md2.gridondirichlet_diag(gridstoflag2)=1;
+		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
+	end
+	if ~isnan(md1.gridondirichlet_thermal),
+		md2.gridondirichlet_thermal(gridstoflag2)=1;
+		if ~isnan(md1.observed_temperature)
+			md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2); 
+		else
+			md2.dirichletvalues_thermal=zeros(numberofgrids2,2);
+			disp(' ')
+			disp('!! modelextract warning: dirichlet values should be checked !!')
+			disp(' ')
+		end
+	end
+
+	%Recreate border stokes
+	if md1.isstokes,
+		pos=find(md2.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
+		md2.gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront
+		stokes_elements=find(md2.elements_type(:,2)==StokesEnum()); %find the elements on the stokes domain
+		borderflags=zeros(numberofgrids2,1); 
+		borderflags(md2.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
+		md2.borderstokes=borderflags-md2.gridonstokes;   %remove stokes grids from this list
+	end
 
 	%Diagnostic
@@ -219,16 +233,15 @@
 		md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);
 	end
-	md2.neumannvalues_diag=NaN; %???????????????????????????
-	if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag=md1.gridondirichlet_diag(pos_grid); end;
-	if ~isnan(md2.dirichletvalues_diag),
-		md2.dirichletvalues_diag=[];
-		for i=1:size(md1.dirichletvalues_diag,2)
-			md2.dirichletvalues_diag(:,i)=md1.dirichletvalues_diag(pos_grid,i);
-		end
-	end
-
-	%Thermal
-	if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal=md1.gridondirichlet_thermal(pos_grid); end;
-	if ~isnan(md2.dirichletvalues_thermal), md2.dirichletvalues_thermal=md1.dirichletvalues_thermal(pos_grid); end;
+	if ~isnan(md2.segmentonneumann_diag_stokes)
+		md2.segmentonneumann_diag_stokes(:,1)=Pgrid(md1.segmentonneumann_diag_stokes(:,1)); 
+		md2.segmentonneumann_diag_stokes(:,2)=Pgrid(md1.segmentonneumann_diag_stokes(:,2)); 
+		md2.segmentonneumann_diag_stokes(:,end)=Pelem(md1.segmentonneumann_diag_stokes(:,end)); 
+		if strcmpi(md1.type,'3d')
+			md2.segmentonneumann_diag_stokes(:,3)=Pgrid(md1.segmentonneumann_diag_stokes(:,3)); 
+			md2.segmentonneumann_diag_stokes(:,4)=Pgrid(md1.segmentonneumann_diag_stokes(:,4)); 
+		end
+		md2.segmentonneumann_diag_stokes=md2.segmentonneumann_diag_stokes(find(md2.segmentonneumann_diag_stokes(:,1) & md2.segmentonneumann_diag_stokes(:,2)),:);
+	end
+	md2.neumannvalues_diag=NaN;
 
 	%Transient
@@ -239,5 +252,5 @@
 		md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:);
 	end
-	md2.neumannvalues_prog=NaN; %???????????????????????????
+	md2.neumannvalues_prog=NaN;
 	if ~isnan(md2.segmentonneumann_prog2)
 		md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1)); 
@@ -246,24 +259,4 @@
 		md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:);
 	end
-	md2.neumannvalues_prog=NaN; %???????????????????????????
-	if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog=md1.gridondirichlet_prog(pos_grid); end;
-	if ~isnan(md2.dirichletvalues_prog), md2.dirichletvalues_prog=md1.dirichletvalues_prog(pos_grid); end;
-
-	%Observations
-	if ~isnan(md2.vx_obs), md2.vx_obs=md1.vx_obs(pos_grid); end;
-	if ~isnan(md2.vy_obs), md2.vy_obs=md1.vy_obs(pos_grid); end;
-	if ~isnan(md2.vel_obs), md2.vel_obs=md1.vel_obs(pos_grid); end;
-	if ~isnan(md2.accumulation), md2.accumulation=md1.accumulation(pos_grid); end;
-	if ~isnan(md2.geothermalflux), md2.geothermalflux=md1.geothermalflux(pos_grid); end;
-	if ~isnan(md2.observed_temperature), md2.observed_temperature=md1.observed_temperature(pos_grid); end;
-	
-	%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;
 
 	%Results fields
@@ -283,105 +276,15 @@
 		end
 	end
-	if ~isnan(md2.vx), md2.vx=md1.vx(pos_grid); end;
-	if ~isnan(md2.vy), md2.vy=md1.vy(pos_grid); end;
-	if ~isnan(md2.vz), md2.vz=md1.vz(pos_grid); end;
-	if ~isnan(md2.vel), md2.vel=md1.vel(pos_grid); end;
-	if ~isnan(md2.temperature), md2.temperature=md1.temperature(pos_grid); end;
-	if ~isnan(md2.melting), md2.melting=md1.melting(pos_grid); end;
-	if ~isnan(md2.pressure), md2.pressure=md1.pressure(pos_grid); end;
-
-%Now take care of the boudaries (Dirichlet).
-
-%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);
-
-%Same thing for the elements
-flaglist=ismember(md1.elements,gridstoflag1);
-elementstoflag1=find((flaglist(:,1) | flaglist(:,2) | flaglist(:,3)) & ~flag_elem);
-flaglist=ismember(md2.elements,gridstoflag2);
-elementstoflag2=find(flaglist(:,1) | flaglist(:,2) | flaglist(:,3));
-
-%These grids are "on boundary" and hence "on dirichlet"
-test=0;
-if ~isnan(md2.gridonboundary), md2.gridonboundary(gridstoflag2)=1; end;
-if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag(gridstoflag2)=1; end;
-if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog(gridstoflag2)=1; end;
-if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal(gridstoflag2)=1; end;
-if ~isnan(md2.dirichletvalues_diag), 
-	test=1; 
-	if strcmpi(md1.type,'3d')
-		if ~isnan(md2.vx) & ~isnan(md2.vy)
-			md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx(gridstoflag2); 
-			md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy(gridstoflag2);
-		end
-	else
-		if ~isnan(md2.vx_obs) & ~isnan(md2.vy_obs)
-			md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
-			md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
-		end
-	end
-end
-if ~isnan(md2.dirichletvalues_prog), test=1; end;
-if ~isnan(md2.dirichletvalues_thermal), test=1; end;
-
-if test
-	disp(' ')
-	disp('!! modelextract warning: dirichlet values should be checked !!')
-	disp(' ')
-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;
 
 %Keep track of pos_grid and pos_elem
 md2.extractedgrids=pos_grid;
 md2.extractedelements=pos_elem;
-
-%get segments
-elementsonboundary=find(sum(md2.gridonboundary(md2.elements),2)>=2);
-md2.segments=[];
-for i=1:length(elementsonboundary),
-	el=elementsonboundary(i);
-	for j=1:3,
-		grid1=md2.elements(el,j);
-		if (j==3),
-			grid2=md2.elements(el,1);
-		else
-			grid2=md2.elements(el,j+1);
-		end
-		%ok, [grid1,grid2] are a segment only if they belong to one element
-		pos=find(    ((md2.elements(:,1)==grid1) | (md2.elements(:,2)==grid1) | (md2.elements(:,3)==grid1)) & ...
-		            ((md2.elements(:,1)==grid2) | (md2.elements(:,2)==grid2) | (md2.elements(:,3)==grid2))...
-				);
-	    if length(pos)==1,
-			md2.segments=[md2.segments; [grid1 grid2 pos]];
-		end
-	end
-end
-md2.segmentmarkers=ones(size(md2.segments,1));
-
-
-end %function
-
-function bool=cross_extract(x1,y1,x2,y2);
-%return 1 if the polygon defined by (x1,y1) crosses the polygon defined by (x2,y2)
-	x1 = x1(:); y1 = y1(:);
-	x2 = x2(:); y2 = y2(:);
-	l1 = length(x1);
-	l2 = length(x2);
-
-	% Indexing
-	[i11,i12] = meshgrid(1:l1,1:l2);
-	[i21,i22] = meshgrid([2:l1 1],[2:l2 1]);
-	i11 = i11(:); i12 = i12(:);
-	i21 = i21(:); i22 = i22(:);
-
-	% Calculate matrix of intersections
-	S = iscross([x1(i11) x1(i21)]',[y1(i11) y1(i21)]',...
-		[x2(i12) x2(i22)]',[y2(i12) y2(i22)]')';
-	S(find(S==0.5))=0;
-
-	bool=any(any(S));
-end
Index: sm/trunk/src/m/classes/public/modelextract2.m
===================================================================
--- /issm/trunk/src/m/classes/public/modelextract2.m	(revision 1493)
+++ 	(revision )
@@ -1,290 +1,0 @@
-function md2=modelextract2(md1,area)
-%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);
-%
-%   Examples:
-%      md2=modelextract2(md,'Domain.exp');
-%      md2=modelextract2(md,md.elementoniceshelf);
-%
-%   See also: EXTRUDE, COLLAPSE
-
-%some checks
-if ((nargin~=2) | (nargout~=1)),
-	help modelextract2
-	error('modelextract2 error message: bad usage');
-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
-
-%kick out all elements with 3 dirichlets
-spc_elem=find(~flag_elem);
-spc_grid=sort(unique(md1.elements(spc_elem,:)));
-flag=ones(md1.numberofgrids,1);
-flag(spc_grid)=0;
-pos=find(sum(flag(md1.elements),2)==0);
-flag_elem(pos)=0;
-
-%extracted elements and nodes lists
-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);
-flag_grid=zeros(numberofgrids1,1);
-flag_grid(pos_grid)=1;
-
-%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;
-
-	%automatically modify fields
-
-	%loop over model fields
-	model_fields=fields(md1);
-	for i=1:length(model_fields),
-
-		%get field
-		field=md1.(model_fields(i));
-		fieldsize=size(field);
-
-		%size = number of grids * n
-		if fieldsize(1)==numberofgrids1
-			md2.(model_fields(i))=field(pos_grid,:);
-
-		%size = number of elements * n
-		elseif fieldsize(1)==numberofelements1
-			md2.(model_fields(i))=field(pos_elem,:);
-		end
-	end
-
-	%modify some specific fields
-
-	%Mesh
-	md2.numberofelements=numberofelements2;
-	md2.numberofgrids=numberofgrids2;
-	md2.elements=elements_2;
-
-	%uppernodes lowernodes
-	if strcmpi(md1.type,'3d')
-		md2.uppergrids=md1.uppergrids(pos_grid);
-		pos=find(~isnan(md2.uppergrids));
-		md2.uppergrids(pos)=Pgrid(md2.uppergrids(pos));
-
-		md2.lowergrids=md1.lowergrids(pos_grid);
-		pos=find(~isnan(md2.lowergrids));
-		md2.lowergrids(pos)=Pgrid(md2.lowergrids(pos));
-	end
-
-	%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
-	if strcmpi(md1.type,'2d')
-		md2.segments=findsegments(md2);
-		md2.gridonboundary=zeros(numberofgrids2,1); md1.gridonboundary(md2.segments(:,1:2))=1;
-	end
-
-	%Boundary conditions: Dirichlets on new boundary
-	%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.gridondirichlet_diag),
-		md2.gridondirichlet_diag(gridstoflag2)=1;
-		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
-	end
-	if ~isnan(md1.gridondirichlet_thermal),
-		md2.gridondirichlet_thermal(gridstoflag2)=1;
-		if ~isnan(md1.observed_temperature)
-			md2.dirichletvalues_thermal(gridstoflag2,1)=md2.observed_temperature(gridstoflag2); 
-		else
-			md2.dirichletvalues_thermal=zeros(numberofgrids2,2);
-			disp(' ')
-			disp('!! modelextract warning: dirichlet values should be checked !!')
-			disp(' ')
-		end
-	end
-
-	%Recreate border stokes
-	if md1.isstokes,
-		pos=find(md2.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
-		md2.gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront
-		stokes_elements=find(md2.elements_type(:,2)==StokesEnum()); %find the elements on the stokes domain
-		borderflags=zeros(numberofgrids2,1); 
-		borderflags(md2.elements(stokes_elements,:))=1; %find all the grids of the elements on stokes domain, ie stokes grids and borderstokes
-		md2.borderstokes=borderflags-md2.gridonstokes;   %remove stokes grids from this list
-	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
-	if ~isnan(md2.segmentonneumann_diag_stokes)
-		md2.segmentonneumann_diag_stokes(:,1)=Pgrid(md1.segmentonneumann_diag_stokes(:,1)); 
-		md2.segmentonneumann_diag_stokes(:,2)=Pgrid(md1.segmentonneumann_diag_stokes(:,2)); 
-		md2.segmentonneumann_diag_stokes(:,end)=Pelem(md1.segmentonneumann_diag_stokes(:,end)); 
-		if strcmpi(md1.type,'3d')
-			md2.segmentonneumann_diag_stokes(:,3)=Pgrid(md1.segmentonneumann_diag_stokes(:,3)); 
-			md2.segmentonneumann_diag_stokes(:,4)=Pgrid(md1.segmentonneumann_diag_stokes(:,4)); 
-		end
-		md2.segmentonneumann_diag_stokes=md2.segmentonneumann_diag_stokes(find(md2.segmentonneumann_diag_stokes(:,1) & md2.segmentonneumann_diag_stokes(:,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;
-
-%Keep track of pos_grid and pos_elem
-md2.extractedgrids=pos_grid;
-md2.extractedelements=pos_elem;