Index: /issm/trunk-jpl/src/m/classes/model.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/model.m	(revision 26160)
+++ /issm/trunk-jpl/src/m/classes/model.m	(revision 26161)
@@ -818,4 +818,119 @@
 			md2.mesh.extractedvertices=pos_node;
 			md2.mesh.extractedelements=pos_elem;
+		end % }}}
+		function md2 = refine(md) % {{{
+			%refine - split all triangles into 3 to refine the mesh everywhere
+			%
+			%   This function only works for 2d triangle meshes
+			%
+			%   Usage:
+			%      md2=refine(md);
+			%
+			%   See also: EXTRUDE, COLLAPSE, EXTRACT
+
+			%Check incoming 
+			if ~strcmp(elementtype(md.mesh),'Tria')
+				error('not supported for 3d meshes');
+			end
+
+			%copy model
+			md2=md;
+
+			disp('Getting edges');
+			%initialization of some variables
+			nbe   = md.mesh.numberofelements;
+			nbv   = md.mesh.numberofvertices;
+			index = md.mesh.elements;
+			elementslist=1:nbe;
+			%1: list of edges
+			edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
+			%2: find unique edges
+			[edges,I,J]=unique(sort(edges,2),'rows');
+			%3: unique edge numbers
+			vec=J;
+			%4: unique edges numbers in each triangle (2 triangles sharing the same edge will have the same edge number)
+			edges_tria=[vec(elementslist+nbe) vec(elementslist+2*nbe) vec(elementslist)];
+
+			% We divide each element as follows
+			%
+			%                   e2
+			%    n1 ------------+------------ n3
+			%       \          / \          /
+			%        \    1   /   \   3    /
+			%         \      /     \      /
+			%          \    /   2   \    /
+			%           \  /         \  /
+			%         e3 +____________\/ e1
+			%             \           /
+			%              \         /
+			%               \   4   /
+			%                \     /
+			%                 \   /
+			%                   n2
+
+			%Create new coordinates
+			disp('Remeshing...');
+			x_edges = 0.5*(md.mesh.x(edges(:,1)) + md.mesh.x(edges(:,2)));
+			y_edges = 0.5*(md.mesh.y(edges(:,1)) + md.mesh.y(edges(:,2)));
+			md2.mesh.x = [md2.mesh.x;x_edges];
+			md2.mesh.y = [md2.mesh.y;y_edges];
+			md2.mesh.elements = [...
+				index(:,1)          nbv+edges_tria(:,3) nbv+edges_tria(:,2);...
+				nbv+edges_tria(:,2) nbv+edges_tria(:,3) nbv+edges_tria(:,1);...
+				nbv+edges_tria(:,2) nbv+edges_tria(:,1) index(:,3);...
+				nbv+edges_tria(:,3) index(:,2)          nbv+edges_tria(:,1)];
+			md2.mesh.numberofelements = 4*nbe;
+			md2.mesh.numberofvertices = nbv + size(edges,1);
+			disp(['   Old number of elements: ' num2str(nbe)]);
+			disp(['   New number of elements: ' num2str(4*nbe)]);
+
+			disp('Interpolate all fields');
+			numberofvertices1 = md.mesh.numberofvertices;
+			numberofelements1 = md.mesh.numberofelements;
+			nbv2 = md2.mesh.numberofvertices;
+			%Create transformation vectors
+			nbedges = size(edges,1);
+			Pelem = sparse(1:4*nbe,repmat([1:nbe],1,4),ones(4*nbe,1),4*nbe,nbe);
+			Pnode = sparse([1:nbv,repmat([nbv+1:nbv+nbedges],1,2)],[1:nbv edges(:)'],[ones(nbv,1);1/2*ones(2*nbedges,1)],md2.mesh.numberofvertices,nbv);
+			%loop over model fields
+			model_fields=setxor(fields(md),{'mesh'});
+			%remove mesh from this field
+			for i=1:length(model_fields),
+				%get field
+				field=md.(model_fields{i});
+				fieldsize=size(field);
+				if isobject(field), %recursive call
+					object_fields=fields(md.(model_fields{i}));
+					for j=1:length(object_fields),
+						%get field
+						field=md.(model_fields{i}).(object_fields{j});
+						fieldsize=size(field);
+						%size = number of nodes * n
+						if fieldsize(1)==numberofvertices1
+							md2.(model_fields{i}).(object_fields{j})=Pnode*field;
+						elseif (fieldsize(1)==numberofvertices1+1)
+							md2.(model_fields{i}).(object_fields{j})=[Pnode*field(1:end-1,:); field(end,:)];
+							%size = number of elements * n
+						elseif fieldsize(1)==numberofelements1
+							md2.(model_fields{i}).(object_fields{j})=Pelem*field;
+						elseif (fieldsize(1)==numberofelements1+1)
+							md2.(model_fields{i}).(object_fields{j})=[Pelem*field(1:end-1,:); field(end,:)];
+						end
+					end
+				else
+					%size = number of nodes * n
+					if fieldsize(1)==numberofvertices1
+						md2.(model_fields{i})=Pnode*field;
+					elseif (fieldsize(1)==numberofvertices1+1)
+						md2.(model_fields{i})=[Pnode*field(1:end-1,:); field(end,:)];
+						%size = number of elements * n
+					elseif fieldsize(1)==numberofelements1
+						md2.(model_fields{i})=Pelem*field;
+					elseif (fieldsize(1)==numberofelements1+1)
+						md2.(model_fields{i})=[Pelem*field(1:end-1,:); field(end,:)];
+					end
+				end
+			end
+
 		end % }}}
 		function md = extrude(md,varargin) % {{{