Index: /issm/trunk-jpl/src/m/mesh/bamg.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.m (revision 13494) +++ /issm/trunk-jpl/src/m/mesh/bamg.m (revision 13495) @@ -137,5 +137,5 @@ %Get position of the two nodes of the edge in domain i1=j; - i2=mod(j+1,domain(1).nods); + i2=j+1; %rift is crossing edge [i1 i2] of the domain @@ -152,5 +152,5 @@ disp('moving tip-domain intersection point'); - %Get position of the closest point + %Get position of the closer point if tipdis(1)>tipdis(2), pos=i2; @@ -165,5 +165,5 @@ bamg_geometry.Edges=[bamg_geometry.Edges;... pos count+1 (1+i);... - [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; + [transpose(count+1:count+nods-1) transpose(count+2:count+nods) (1+i)*ones(nods-1,1)]]; count=count+nods; @@ -175,10 +175,10 @@ count=count+1; - %Decompose the crossing edge in 2 subedges + %Decompose the crossing edge into 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 + if isempty(pos) error('bamg error message: a problem occurred...'); 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) count bamg_geometry.Edges(pos,3);... + count bamg_geometry.Edges(pos,2) bamg_geometry.Edges(pos,3);... bamg_geometry.Edges(pos+1:end,:)]; @@ -188,5 +188,5 @@ bamg_geometry.Edges=[bamg_geometry.Edges;... count count+1 2 ;... - [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; + [transpose(count+1:count+nods-1) transpose(count+2:count+nods) (1+i)*ones(nods-1,1)]]; count=count+nods; @@ -198,5 +198,5 @@ nods=rift(i).nods-1; bamg_geometry.Vertices=[bamg_geometry.Vertices; [rift(i).x(:) rift(i).y(:) ones(nods+1,1)]]; - bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose([count+2:count+nods+1]) (1+i)*ones(nods,1)]]; + bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods) transpose(count+2:count+nods+1) (1+i)*ones(nods,1)]]; count=count+nods+1; end @@ -224,5 +224,5 @@ nods=size(track,1); bamg_geometry.Vertices=[bamg_geometry.Vertices; track]; - bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) 3.*ones(nods-1,1)]]; + bamg_geometry.Edges=[bamg_geometry.Edges; [transpose(count+1:count+nods-1) transpose(count+2:count+nods) 3.*ones(nods-1,1)]]; %update counter Index: /issm/trunk-jpl/src/m/mesh/bamg.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13494) +++ /issm/trunk-jpl/src/m/mesh/bamg.py (revision 13495) @@ -7,4 +7,5 @@ from expread import * from expwrite import * +from SegIntersect import * from MatlabFuncs import * from BamgMesher import * @@ -121,91 +122,91 @@ elif numpy.any(numpy.logical_not(flags)): - raise RuntimeError("bamg.py for rifts is not complete.") #We LOTS of work to do print("Rift tip outside of or on the domain has been detected and is being processed...") - """ - - %check that only one point is outside (for now) - if sum(~flags)~=1, - error('bamg error message: only one point outside of the domain is supported yet'); - end - - %Move tip outside to the first position - if flags(1)==0, - %OK, first point is outside (do nothing), - elseif (flags(end)==0), - rift(i).x=flipud(rift(i).x); - rift(i).y=flipud(rift(i).y); - else - error('bamg error message: only a rift tip can be outside of the domain'); - end - - %Get cordinate of intersection point - x1=rift(i).x(1); y1=rift(i).y(1); - x2=rift(i).x(2); y2=rift(i).y(2); - for j=1:length(domain[0]['x'])-1; - if SegIntersect([x1 y1; x2 y2],[domain[0]['x'](j) domain[0]['y'](j); domain[0]['x'](j+1) domain[0]['y'](j+1)]), - - %Get position of the two nodes of the edge in domain - i1=j; - i2=mod(j+1,domain[0]['nods']); - - %rift is crossing edge [i1 i2] of the domain - %Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) - x3=domain[0]['x'](i1); y3=domain[0]['y'](i1); - x4=domain[0]['x'](i2); y4=domain[0]['y'](i2); - x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); - 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]); - - 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, now 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;... - pos count+1 (1+i);... - [transpose(count+1:count+nods-1) transpose([count+2:count+nods]) (1+i)*ones(nods-1,1)]]; + + #check that only one point is outside (for now) + if numpy.sum(numpy.logical_not(flags).astype(int))!=1: + raise RuntimeError("bamg error message: only one point outside of the domain is supported yet") + + #Move tip outside to the first position + if not flags[0]: + #OK, first point is outside (do nothing), + pass + elif not flags[-1]: + rifti['x']=numpy.flipud(rifti['x']) + rifti['y']=numpy.flipud(rifti['y']) + else: + raise RuntimeError("bamg error message: only a rift tip can be outside of the domain") + + #Get cordinate of intersection point + x1=rifti['x'][0] + y1=rifti['y'][0] + x2=rifti['x'][1] + y2=rifti['y'][1] + for j in xrange(0,numpy.size(domain[0]['x'])-1): + if SegIntersect(numpy.array([[x1,y1],[x2,y2]]),numpy.array([[domain[0]['x'][j],domain[0]['y'][j]],[domain[0]['x'][j+1],domain[0]['y'][j+1]]])): + + #Get position of the two nodes of the edge in domain + i1=j + i2=j+1 + + #rift is crossing edge [i1 i2] of the domain + #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) + x3=domain[0]['x'][i1] + y3=domain[0]['y'][i1] + x4=domain[0]['x'][i2] + y4=domain[0]['y'][i2] +# x=det([det([x1 y1; x2 y2]) x1-x2;det([x3 y3; x4 y4]) x3-x4])/det([x1-x2 y1-y2;x3-x4 y3-y4]); +# 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]); + x=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])); + y=numpy.linalg.det(numpy.array([[numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4]]))/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]])); + + segdis= sqrt((x4-x3)**2+(y4-y3)**2) + tipdis=numpy.array([sqrt((x-x3)**2+(y-y3)**2),sqrt((x-x4)**2+(y-y4)**2)]) + + if numpy.min(tipdis)/segdis < options.getfieldvalue('toltip',0): + print("moving tip-domain intersection point") + + #Get position of the closer point + if tipdis[0]>tipdis[1]: + pos=i2 + else: + pos=i1 + + #This point is only in Vertices (number pos). + #OK, now we can add our own rift + nods=rifti['nods']-1 + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) + bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ + numpy.array([[pos,count+1,(1+i)]]),\ + numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) + count+=nods + + break + + else: + #Add intersection point to Vertices + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.array([[x,y,1]]))) + count+=1 + + #Decompose the crossing edge into 2 subedges + pos=numpy.nonzero(numpy.logical_and(bamg_geometry.Edges[:,0]==i1,bamg_geometry.Edges[:,1]==i2))[0] + if not pos: + raise RuntimeError("bamg error message: a problem occurred...") + bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges[0:pos-1,:],\ + numpy.array([[bamg_geometry.Edges[pos,0],count ,bamg_geometry.Edges[pos,2]]]),\ + numpy.array([[count ,bamg_geometry.Edges[pos,1],bamg_geometry.Edges[pos,2]]]),\ + bamg_geometry.Edges[pos+1:,:])) + + #OK, now we can add our own rift + nods=rifti['nods']-1 + bamg_geometry.Vertices=numpy.vstack((bamg_geometry.Vertices,numpy.hstack((rifti['x'][1:].reshape(-1,1),rifti['y'][1:].reshape(-1,1),numpy.ones((nods,1)))))) + bamg_geometry.Edges=numpy.vstack((bamg_geometry.Edges,\ + numpy.array([[count,count+1,2]]),\ + numpy.hstack((numpy.arange(count+1,count+nods).reshape(-1,1),numpy.arange(count+2,count+nods+1).reshape(-1,1),(1+i)*numpy.ones((nods-1,1)))))) count=count+nods; - 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, now 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]) (1+i)*ones(nods-1,1)]]; - count=count+nods; - - break; - end - end - end - """ + break + else: nods=rifti['nods']-1 @@ -384,6 +385,6 @@ #Get coordinate of intersection point (http://mathworld.wolfram.com/Line-LineIntersection.html) - x=det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) - y=det(numpy.array([det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) + x=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),x1-x2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),x3-x4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) + y=numpy.linalg.det(numpy.array([numpy.linalg.det(numpy.array([[x1,y1],[x2,y2]])),y1-y2],[numpy.linalg.det(numpy.array([[x3,y3],[x4,y4]])),y3-y4])/numpy.linalg.det(numpy.array([[x1-x2,y1-y2],[x3-x4,y3-y4]]))) #Add vertex to the list of vertices