Index: /issm/trunk-jpl/src/m/classes/model.js
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.js	(revision 22718)
+++ /issm/trunk-jpl/src/m/classes/model.js	(revision 22719)
@@ -360,4 +360,126 @@
 			}
 		} //}}}
+		this.extract = function(md,area) { //{{{
+			//extract - 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 
+			//   extract2d, add '~' to the name of the domain file (ex: '~HO.exp');
+			//   an empty string '' will be considered as an empty domain
+			//   a string 'all' will be considered as the entire domain
+			//
+			//   Usage:
+			//      md2=extract(md,area);
+			//
+			//   Examples:
+			//      md2=extract(md,'Domain.exp');
+			//
+			//   See also: EXTRUDE, COLLAPSE
+
+			//some checks on list of arguments
+			var argc = arguments.length;
+
+			if (!((argc == 2) | (argc == 1))) {
+				throw "extract error message: bad usage";
+			}
+			
+			//get elements that are inside area
+			flag_elem=FlagElements(this,area);
+
+			if (!ArrayAnyEqual(flag_elem,1))throw "extracted model is empty!"
+
+			/*kick out all elements with 3 dirichlets: actually, not so fast, not good for javscript usually.
+			spc_elem=find(~flag_elem);
+			spc_node=sort(unique(md1.mesh.elements(spc_elem,:)));
+			flag=ones(md1.mesh.numberofvertices,1);
+			flag(spc_node)=0;
+			pos=find(sum(flag(md1.mesh.elements),2)==0);
+			flag_elem(pos)=0;*/
+
+			//extracted elements and nodes lists
+			var pos_elem = ArrayFind(flag_elem,1);
+			var dup_nodes= new Array(pos_elem.length*3);
+			for(var i=0;i<pos_elem.length;i++){
+				dup_nodes[3*i]=md.mesh.elements[pos_elem[i]][0]-1;
+				dup_nodes[3*i+1]=md.mesh.elements[pos_elem[i]][1]-1;
+				dup_nodes[3*i+2]=md.mesh.elements[pos_elem[i]][2]-1;
+			}
+			pos_node=ArrayUnique(dup_nodes); pos_node=ArraySort(pos_node);
+
+			//keep track of some fields
+			var numberofvertices1=md.mesh.numberofvertices;
+			var numberofelements1=md.mesh.numberofelements;
+			var numberofvertices2=pos_node.length;
+			var numberofelements2=pos_elem.length;
+			var flag_node=NewArrayFill(numberofvertices1,0);
+			for (var i=0;i<pos_node.length;i++)flag_node[pos_node[i]]=1;
+
+			//Create Pelem and Pnode (transform old nodes in new nodes and same thing for the elements)
+			Pelem=NewArrayFill(numberofelements1,0);
+			for (var i=0;i<numberofelements2;i++) Pelem[pos_elem[i]]=i;
+			Pnode=NewArrayFill(numberofvertices1,0);
+			for (var i=0;i<numberofvertices2;i++) Pnode[pos_node[i]]=i;
+			//renumber the elements (some nodes won't exist anymore)
+			var elements_2=NewArrayFill2D(numberofelements2,3,0);
+			for (var i=0;i<numberofelements2;i++){
+				for (var j=0;j<3;j++){
+					elements_2[i][j]=Pnode[md.mesh.elements[i][j]-1]+1;
+				}
+			}
+			
+			/*if isa(md.mesh,'mesh3dprisms'),
+				elements_2(:,4)=Pnode(elements_2(:,4));
+				elements_2(:,5)=Pnode(elements_2(:,5));
+				elements_2(:,6)=Pnode(elements_2(:,6));
+			end
+			*/
+			
+			//OK, now create the new model!
+
+			//take every field from model
+			var md2=md.deepcopy(md);
+			//var md2=new model(); md2.mesh=new mesh3dsurface();
+			
+			//deal with mesh: {{{
+			md2.mesh.numberofvertices=numberofvertices2;
+			md2.mesh.numberofelements=numberofelements2;
+			md2.mesh.elements=elements_2;
+
+			md2.mesh.x=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.x[i]=md.mesh.x[pos_node[i]];
+			md2.mesh.y=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.y[i]=md.mesh.y[pos_node[i]];
+			md2.mesh.z=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.z[i]=md.mesh.z[pos_node[i]];
+			md2.mesh.lat=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.lat[i]=md.mesh.lat[pos_node[i]];
+			md2.mesh.long=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.long[i]=md.mesh.long[pos_node[i]];
+			md2.mesh.r=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.mesh.r[i]=md.mesh.r[pos_node[i]];
+			//}}}
+			//deal with geometry: {{{
+			md2.geometry.base=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.base[i]=md.geometry.base[pos_node[i]];
+			md2.geometry.thickness=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.thickness[i]=md.geometry.thickness[pos_node[i]];
+			md2.geometry.surface=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.surface[i]=md.geometry.surface[pos_node[i]];
+			md2.geometry.bed=new Array(numberofvertices2); for (var i=0;i<numberofvertices2;i++)md2.geometry.bed[i]=md.geometry.bed[pos_node[i]];
+			//}}}
+
+			//Keep track of pos_node and pos_elem
+			for (var i=0;i<md2.mesh.numberofvertices;i++)pos_node[i]=pos_node[i]+1;
+			for (var i=0;i<md2.mesh.numberofelements;i++)pos_elem[i]=pos_elem[i]+1;
+			md2.mesh.extractedvertices=pos_node;
+			md2.mesh.extractedelements=pos_elem;
+
+			return md2;
+
+
+
+			//automatically modify fields
+
+			//loop over model fields
+			//model_fields=fields(md);
+
+
+
+
+
+		} //}}}
 		this.collapse = function(md) { //{{{
 			/*