Per i futuri lettori, ecco una funzione python (usando numpy) che crea una lista degli indici di punti sul confine.
Sto iniziando con una rete generale di punti collegati (non necessariamente triangolati), quindi i miei dati sono nella forma di (xy, NL, KL) con NL che è la lista dei vicini e KL che è la lista di connettività (vedi sotto) .
def extract_boundary_from_NL(xy,NL,KL):
'''
Extract the boundary of a 2D network (xy,NL,KL).
Parameters
----------
xy : #pts x 2 float array
point set in 2D
NL : #pts x max # nearest neighbors int array
Neighbor list. The ith row contains the indices of xy that are the bonded pts to the ith pt.
Nonexistent bonds are replaced by zero.
KL : #pts x max # nearest neighbors int array
Connectivity list. The jth column of the ith row ==1 if pt i is bonded to pt NL[i,j].
The jth column of the ith row ==0 if pt i is not bonded to point NL[i,j].
Returns
----------
boundary : #points on boundary x 1 int array
indices of points living on boundary of the network
'''
# Initialize the list of boundary indices to be larger than necessary
bb = np.zeros(len(xy), dtype=int)
# Start with the rightmost point, which is guaranteed to be
# at the convex hull and thus also at the outer edge.
# Then take the first step to be along the minimum angle bond
rightIND = np.where(xy[:,0]== max(xy[:,0]))[0]
# If there are more than one rightmost point, choose one
if rightIND.size >1:
rightIND = rightIND[0]
# Grab the true neighbors of this starting point
neighbors = NL[rightIND,np.argwhere(KL[rightIND]).ravel()]
# Compute the angles of the neighbor bonds
angles = np.mod(np.arctan2(xy[neighbors,1]-xy[rightIND,1],xy[neighbors,0]-xy[rightIND,0]).ravel(), 2*np.pi)
# Take the second particle to be the one with the lowest bond angle (will be >= pi/2)
nextIND = neighbors[angles==min(angles)][0]
bb[0] = rightIND
dmyi = 1
# as long as we haven't completed the full outer edge/boundary, add nextIND
while nextIND != rightIND:
bb[dmyi] = nextIND
n_tmp = NL[nextIND,np.argwhere(KL[nextIND]).ravel()]
# Exclude previous boundary particle from the neighbors array
# since its angle with itself is zero!
neighbors = np.delete(n_tmp, np.where(n_tmp == bb[dmyi-1])[0])
angles = np.mod( np.arctan2(xy[neighbors,1]-xy[nextIND,1],xy[neighbors,0]-xy[nextIND,0]).ravel() \
- np.arctan2( xy[bb[dmyi-1],1]-xy[nextIND,1], xy[bb[dmyi-1],0]-xy[nextIND,0] ).ravel(), 2*np.pi)
nextIND = neighbors[angles==min(angles)][0]
dmyi += 1
# Truncate the list of boundary indices
boundary = bb[0:dmyi]
return boundary
Per inserire una triangolazione nel formato (NL, KL) per utilizzare quanto sopra, puoi fare:
BL = TRI2BL(TRI)
NL,KL = BL2NLandKL(BL,NN=10)
def BL2NLandKL(BL, NP='auto', NN=6):
'''Convert bond list (#bonds x 2) to neighbor list (NL) and connectivity list (KL) for lattice of bonded points.
Returns KL as ones where there is a bond and zero where there is not.
(Even if you just want NL from BL, you have to compute KL anyway.)
Note that this makes no attempt to preserve any previous version of NL, which in the philosophy of these simulations should remain constant during a simulation.
If NL is known, use BL2KL instead, which creates KL according to the existing NL.
Parameters
----------
BL : array of dimension #bonds x 2
Each row is a bond and contains indices of connected points
NP : int
number of points (defines the length of NL and KL)
NN : int
maximum number of neighbors (defines the width of NL and KL)
Returns
----------
NL : array of dimension #pts x max(#neighbors)
The ith row contains indices for the neighbors for the ith point.
KL : array of dimension #pts x (max number of neighbors)
Spring constant list, where 1 corresponds to a true connection while 0 signifies that there is not a connection.
'''
if NP=='auto':
if BL.size>0:
NL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
KL = np.zeros((max(BL.ravel())+1,NN),dtype=np.intc)
else:
raise RuntimeError('ERROR: there is no BL to use to define NL and KL, so cannot run BL2NLandKL()')
else:
NL = np.zeros((NP,NN),dtype=np.intc)
KL = np.zeros((NP,NN),dtype=np.intc)
if BL.size > 0:
for row in BL:
col = np.where(KL[row[0],:]==0)[0][0]
NL[row[0],col] = row[1]
KL[row[0],col] = 1
col = np.where(KL[row[1],:]==0)[0][0]
NL[row[1],col] = row[0]
KL[row[1],col] = 1
return NL, KL
def TRI2BL(TRI):
'''
Convert triangulation index array (Ntris x 3) to Bond List (Nbonds x 2) array.
Parameters
----------
TRI : Ntris x 3 int array
Triangulation of a point set. Each row gives indices of vertices of single triangle.
Returns
----------
BL : Nbonds x 2 int array
Bond list
'''
# each edge is shared by 2 triangles unless at the boundary.
# each row contains 3 edges.
# An upper bound on the number bonds is 3*len(TRI)
BL = np.zeros((3*len(TRI),2),dtype = int)
dmyi = 0
for row in TRI:
BL[dmyi] = [row[0], row[1]]
BL[dmyi+1] = [row[1], row[2]]
BL[dmyi+2] = [row[0], row[2]]
dmyi += 3
# Sort each row to be ascending
BL_sort = np.sort(BL, axis=1)
BLtrim = unique_rows(BL_sort)
return BLtrim
