Index: /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.m (revision 13523) +++ /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.m (revision 13524) @@ -1,15 +1,15 @@ -function edgeelements=ElementsFromEdge(elements,A,B) +function edgeelements=ElementsFromEdge(elements,A,B) %ELEMENTSFROMEDGE: find elements connected to one edge defined by nodes A and B % -% Usage: edgeelements=ElementsFromEdge(elements,A,B) +% Usage: edgeelements=ElementsFromEdge(elements,A,B) % -% Eg: edgeelements=ElementsFromEdge(md.mesh.elements,tip1,tip2) +% Eg: edgeelements=ElementsFromEdge(md.mesh.elements,tip1,tip2) % % edgeelements=find(... -(elements(:,1)==A & elements(:,2)==B )| ... -(elements(:,1)==A & elements(:,3)==B )| ... -(elements(:,2)==A & elements(:,3)==B )| ... -(elements(:,2)==A & elements(:,1)==B )| ... -(elements(:,3)==A & elements(:,1)==B )| ... -(elements(:,3)==A & elements(:,2)==B )); + (elements(:,1)==A & elements(:,2)==B )| ... + (elements(:,1)==A & elements(:,3)==B )| ... + (elements(:,2)==A & elements(:,3)==B )| ... + (elements(:,2)==A & elements(:,1)==B )| ... + (elements(:,3)==A & elements(:,1)==B )| ... + (elements(:,3)==A & elements(:,2)==B )); Index: /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py (revision 13524) +++ /issm/trunk-jpl/src/m/mesh/ElementsFromEdge.py (revision 13524) @@ -0,0 +1,27 @@ +import numpy + +def ElementsFromEdge(elements,A,B): + """ + ELEMENTSFROMEDGE: find elements connected to one edge defined by nodes A and B + + Usage: edgeelements=ElementsFromEdge(elements,A,B) + + Eg: edgeelements=ElementsFromEdge(md.mesh.elements,tip1,tip2) + + """ + + edgeelements=numpy.nonzero(\ + numpy.logical_or( \ + numpy.logical_or( \ + numpy.logical_or(numpy.logical_and(elements[:,0]==A,elements[:,1]==B), \ + numpy.logical_and(elements[:,0]==A,elements[:,2]==B)) \ + , \ + numpy.logical_or(numpy.logical_and(elements[:,1]==A,elements[:,2]==B), \ + numpy.logical_and(elements[:,1]==A,elements[:,0]==B)) \ + ), \ + numpy.logical_or(numpy.logical_and(elements[:,2]==A,elements[:,0]==B), \ + numpy.logical_and(elements[:,2]==A,elements[:,1]==B)) \ + ))[0]+1 + + return edgeelements + Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13523) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.m (revision 13524) @@ -41,5 +41,5 @@ elements=[]; - while flags(B), %as long as B does not belong to the domain outline, keep looking. + while flags(B), %as long as B does not belong to the domain outline, keep looking. %detect elements on edge A,B: edgeelements=ElementsFromEdge(md.mesh.elements,A,B); @@ -69,5 +69,5 @@ %deal with segments tipsegments=find((md.mesh.segments(:,1)==tip) | (md.mesh.segments(:,2)==tip)); - for k=1:length(tipsegments), + for k=1:length(tipsegments), segment_index=tipsegments(k); pos=find(md.mesh.segments(segment_index,1:2)~=tip); Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 13524) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessoutsiderifts.py (revision 13524) @@ -0,0 +1,110 @@ +import numpy +from ElementsFromEdge import * +from MatlabFuncs import * + +def meshprocessoutsiderifts(md,domainoutline): + """ + MESHPROCESSOUTSIDERIFTS - process rifts when they touch the domain outline + + Usage: + md=meshprocessoutsiderifts(md,domain) + + """ + + #go through rifts, and figure out which ones touch the domain outline + for rift in md.rifts.riftstruct: + + #first, flag nodes that belong to the domain outline + flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) + + tips=rift.tips + outsidetips=tips[numpy.nonzero(flags[rift.tips-1])[0]] + + #we have found outsidetips, tips that touch the domain outline. go through them + for tip in outsidetips: + + #find tip in the segments, take first segment (there should be 2) that holds tip, + #and node_connected_to_tip is the other node on this segment: + tipindex=numpy.nonzero(rift.segments[:,0]==tip)[0] + if tipindex: + tipindex=tipindex[0] + node_connected_to_tip=rift.segments[tipindex,1] + else: + tipindex=numpy.nonzero(rift.segments[:,1]==tip)[0] + tipindex=tipindex[0] + node_connected_to_tip=rift.segments[tipindex,1] + + #ok, we have the tip node, and the first node connected to it, on the rift. Now, + #identify all the elements that are connected to the tip, and that are on the same + #side of the rift. + A=tip + B=node_connected_to_tip + + elements=numpy.empty(0) + + while flags(B): #as long as B does not belong to the domain outline, keep looking. + #detect elements on edge A,B: + edgeelements=ElementsFromEdge(md.mesh.elements,A,B) + #rule out those we already detected + already_detected=ismember(edgeelements,elements) + nextelement=edgeelements(numpy.nonzero(numpy.logical_not(already_detected))[0]) + #add new detected element to the list of elements we are looking for. + elements=numpy.concatenate((elements,nextelement)) + #new B: + B=md.mesh.elements[nextelement-1,numpy.nonzero(numpy.logical_not(ismember(md.mesh.elements[nextelement-1,:],numpy.array([A,B]))))] + + #take the list of elements on one side of the rift that connect to the tip, + #and duplicate the tip on them, so as to open the rift to the outside. + num=numpy.size(md.mesh.x)+1 + md.mesh.x=numpy.concatenate((md.mesh.x,md.mesh.x[tip])) + md.mesh.y=numpy.concatenate((md.mesh.y,md.mesh.y[tip])) + md.mesh.numberofvertices=num + + #replace tip in elements + newelements=md.mesh.elements[elements.astype(int)-1,:] + pos=numpy.nonzero(newelements==tip) + newelements[pos]=num + md.mesh.elements[elements.astype(int)-1,:]=newelements + rift.tips=numpy.concatenate((rift.tips,num)) + + #deal with segments + tipsegments=numpy.nonzero(numpy.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0] + for segment_index in tipsegments: + pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]!=tip)[0] + other_node=md.mesh.segments[segment_index,pos] + if not isconnected(md.mesh.elements,other_node,tip): + pos=numpy.nonzero(md.mesh.segments[segment_index,0:2]==tip)[0] + md.mesh.segments[segment_index,pos]=num + + #Fill in rest of fields: + md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=numpy.size(md.mesh.x) + md.mesh.z=numpy.zeros(md.mesh.numberofvertices) + md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x)) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.rifts.numrifts=length(md.rifts.riftstruct) + md.flowequation.element_equation=3.*numpy.ones(md.mesh.numberofelements) + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) + + return md + +def isconnected(elements,A,B): # {{{ + """ + ISCONNECTED: are two nodes connected by a triangulation? + + Usage: flag=isconnected(elements,A,B) + + """ + + elements=ElementsFromEdge(elements,A,B) + if not elements: + flag=0 + else: + flag=1 + + return flag + # }}} + Index: /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py =================================================================== --- /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 13524) +++ /issm/trunk-jpl/src/m/mesh/rifts/meshprocessrifts.py (revision 13524) @@ -0,0 +1,64 @@ +import numpy +#from TriMeshProcessRifts import * +from ContourToMesh import * +from meshprocessoutsiderifts import * +from GetAreas import * + +def meshprocessrifts(md,domainoutline): + """ + MESHPROCESSRIFTS - process mesh when rifts are present + + split rifts inside mesh (rifts are defined by presence of + segments inside the domain outline) + if domain outline is provided, check for rifts that could touch it, and open them up. + + Usage: + md=meshprocessrifts(md,domainoutline) + + Ex: + md=meshprocessrifts(md,'DomainOutline.exp'); + + """ + + #Call MEX file + [md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.mesh.x,md.mesh.y,md.mesh.segments,md.mesh.segmentmarkers) + if not isinstance(md.rifts.riftstruct,'list') or not md.rifts.riftstruct: + raise RuntimeError("TriMeshProcessRifts did not find any rift") + + #Fill in rest of fields: + numrifts=len(md.rifts.riftstruct) + md.mesh.numberofelements=numpy.size(md.mesh.elements,axis=0) + md.mesh.numberofvertices=numpy.size(md.mesh.x) + md.mesh.z=numpy.zeros(md.mesh.numberofvertices) + md.mesh.vertexonboundary=numpy.zeros(numpy.size(md.mesh.x)) + md.mesh.vertexonboundary[md.mesh.segments[:,0:2].astype(int)-1]=1 + md.mesh.vertexonbed=numpy.ones(md.mesh.numberofvertices) + md.mesh.vertexonsurface=numpy.ones(md.mesh.numberofvertices) + md.mesh.elementonbed=numpy.ones(md.mesh.numberofelements) + md.mesh.elementonsurface=numpy.ones(md.mesh.numberofelements) + + #get coordinates of rift tips + for rift in md.rifts.riftstruct: + rift.tip1coordinates=numpy.hstack((md.mesh.x[rift.tips[0].astype(int)-1].reshape(-1,1),md.mesh.y[rift.tips[0].astype(int)-1].reshape(-1,1))) + rift.tip2coordinates=numpy.hstack((md.mesh.x[rift.tips[1].astype(int)-1].reshape(-1,1),md.mesh.y[rift.tips[1].astype(int)-1].reshape(-1,1))) + + #In case we have rifts that open up the domain outline, we need to open them: + flags=ContourToMesh(md.mesh.elements,md.mesh.x.reshape(-1,1),md.mesh.y.reshape(-1,1),domainoutline,'node',0) + found=0 + for rift in md.rifts.riftstruct: + if flags[rift.tips[0].astype(int)-1]==0: + found=1 + break + if flags[rift.tips[1].astype(int)-1]==0: + found=1 + break + if found: + md=meshprocessoutsiderifts(md,domainoutline) + + #get elements that are not correctly oriented in the correct direction: + aires=GetAreas(md.mesh.elements,md.mesh.x,md.mesh.y) + pos=numpy.nonzero(aires<0)[0] + md.mesh.elements[pos,:]=numpy.array([[md.mesh.elements[pos,1],md.mesh.elements[pos,0],md.mesh.elements[pos,2]]]) + + return md +