import numpy as np
from ElementsFromEdge import ElementsFromEdge
import MatlabFuncs as m
[docs]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,md.mesh.y,domainoutline,'node',0)
tips=rift.tips
outsidetips=tips[np.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=np.nonzero(rift.segments[:,0]==tip)[0]
if tipindex:
tipindex=tipindex[0]
node_connected_to_tip=rift.segments[tipindex,1]
else:
tipindex=np.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=np.empty(0,int)
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=m.ismember(edgeelements,elements)
nextelement=edgeelements(np.nonzero(np.logical_not(already_detected))[0])
#add new detected element to the list of elements we are looking for.
elements=np.concatenate((elements,nextelement))
#new B:
B=md.mesh.elements[nextelement-1,np.nonzero(np.logical_not(m.ismember(md.mesh.elements[nextelement-1,:],np.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=np.size(md.mesh.x)+1
md.mesh.x=np.concatenate((md.mesh.x,md.mesh.x[tip]))
md.mesh.y=np.concatenate((md.mesh.y,md.mesh.y[tip]))
md.mesh.numberofvertices=num
#replace tip in elements
newelements=md.mesh.elements[elements-1,:]
pos=np.nonzero(newelements==tip)
newelements[pos]=num
md.mesh.elements[elements-1,:]=newelements
rift.tips=np.concatenate((rift.tips,num))
#deal with segments
tipsegments=np.nonzero(np.logical_or(md.mesh.segments[:,0]==tip,md.mesh.segments[:,1]==tip))[0]
for segment_index in tipsegments:
pos=np.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=np.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=np.size(md.mesh.elements,axis=0)
md.mesh.numberofvertices=np.size(md.mesh.x)
md.mesh.vertexonboundary=np.zeros(np.size(md.mesh.x),bool)
md.mesh.vertexonboundary[md.mesh.segments[:,0:2]-1]=True
md.rifts.numrifts=length(md.rifts.riftstruct)
return md
[docs]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
# }}}