Index: /issm/trunk/src/m/classes/public/bamg.m
===================================================================
--- /issm/trunk/src/m/classes/public/bamg.m	(revision 3304)
+++ /issm/trunk/src/m/classes/public/bamg.m	(revision 3305)
@@ -4,5 +4,5 @@
 %   Available options (for more details see ISSM website http://issm.jpl.nasa.gov/):
 %
-%   - domain: followed by an ARGUS file that prescribes the domain outline (rifts and holes)
+%   - domain: followed by an ARGUS file that prescribes the domain outline
 %   - hmin  : minimum edge length (default is 10^-100)
 %   - hmax  : maximum esge length (default is 10^100)
@@ -37,4 +37,7 @@
 %   - verbose     : level of verbosity (default is 1)
 %
+%   - rifts : followed by an ARGUS file that prescribes the rifts
+%   - toltip: tolerance to move tip on an existing point of the domain outline
+%
 %   Examples:
 %      md=bamg(md,'domain','DomainOutline.exp','hmax',3000);
@@ -139,25 +142,52 @@
 						y=det([det([x1 y1; x2 y2])  y1-y2;det([x3 y3; x4 y4])  y3-y4])/det([x1-x2 y1-y2;x3-x4 y3-y4]);
 
-						%Add intersection point to Vertices
-						bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1];
-						count=count+1;
-
-						%Decompose the crossing edge in 2 subedges
-						pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2);
-						if isempty(pos) error('bamg error message: a problem occured...'); end
-						bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);...
-							bamg_geometry.Edges(pos,1) count                           bamg_geometry.Edges(pos,3);...
-							count                      bamg_geometry.Edges(pos,2)   bamg_geometry.Edges(pos,3);...
-							bamg_geometry.Edges(pos+1:end,:)];
-
-						%OK, no we can add our own rift
-						nods=rift(i).nods-1;
-						bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
-						bamg_geometry.Edges=[bamg_geometry.Edges;...
-							count  count+1  2 ;...
-							[transpose(count+1:count+nods-1) transpose([count+2:count+nods])  2*ones(nods-1,1)]];
-						count=count+nods;
-
-						break;
+						segdis= sqrt((x4-x3)^2+(y4-y3)^2);
+						tipdis=[sqrt((x-x3)^2+(y-y3)^2)  sqrt((x-x4)^2+(y-y4)^2)];
+
+						if (min(tipdis)/segdis) < getfieldvalue(options,'toltip',0),
+							disp('moving tip-domain intersection point');
+
+							%Get position of the closest point
+							if tipdis(1)>tipdis(2),
+								pos=i2;
+							else
+								pos=i1;
+							end
+
+							%This point is only in Vertices (number pos).
+							%OK, no we can add our own rift
+							nods=rift(i).nods-1;
+							bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
+							bamg_geometry.Edges=[bamg_geometry.Edges;...
+								i1 count+1  2 ;...
+								[transpose(count+1:count+nods-1) transpose([count+2:count+nods])  2*ones(nods-1,1)]];
+							count=count+nods;
+							bamg_geometry.Edges
+
+							break;
+
+						else
+							%Add intersection point to Vertices
+							bamg_geometry.Vertices=[bamg_geometry.Vertices; x y 1];
+							count=count+1;
+
+							%Decompose the crossing edge in 2 subedges
+							pos=find(bamg_geometry.Edges(:,1)==i1 & bamg_geometry.Edges(:,2)==i2);
+							if isempty(pos) error('bamg error message: a problem occured...'); end
+							bamg_geometry.Edges=[bamg_geometry.Edges(1:pos-1,:);...
+								bamg_geometry.Edges(pos,1) count                           bamg_geometry.Edges(pos,3);...
+								count                      bamg_geometry.Edges(pos,2)   bamg_geometry.Edges(pos,3);...
+								bamg_geometry.Edges(pos+1:end,:)];
+
+							%OK, no we can add our own rift
+							nods=rift(i).nods-1;
+							bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(2:end) rift(i).y(2:end) ones(nods,1)]];
+							bamg_geometry.Edges=[bamg_geometry.Edges;...
+								count  count+1  2 ;...
+								[transpose(count+1:count+nods-1) transpose([count+2:count+nods])  2*ones(nods-1,1)]];
+							count=count+nods;
+
+							break;
+						end
 					end
 				end