[docs]defgetNodeNeighbors(NodeCoords,NodeConn,ElemType='auto'):""" Determines the adjacent nodes for each node in the mesh Parameters ---------- NodeCoords : array_like List of nodal coordinates. NodeConn : array_like Nodal connectivity list. ElemType : str, optional Type of element contained in the mesh, by default 'auto'. See converter.solid2edges() for details. 'auto' is suitable for most element types and mixed-element meshes, 4-node elements are assumed to be tets, not quads, unless ElemType is set to 'quad' or 'surf'. Returns ------- NodeNeighbors : list List of neighboring nodes for each node in NodeCoords. """Edges,EdgeElem=converter.solid2edges(NodeCoords,NodeConn,return_EdgeElem=True,ElemType=ElemType)UEdges,idx,inv=converter.edges2unique(Edges,return_idx=True,return_inv=True)NotInMesh=set(range(len(NodeCoords))).difference(np.unique(UEdges))Neighbors=np.append(UEdges.flatten(order='F'),np.repeat(-1,len(NotInMesh)))Idx=np.append(np.fliplr(UEdges).flatten(order='F'),list(NotInMesh))arg=Idx.argsort()key_func=lambdax:x[0]NodeNeighbors=[[zfory,zinx[1]ifz!=-1]forxinitertools.groupby(zip(Idx[arg],Neighbors[arg]),key_func)]returnNodeNeighbors
[docs]defgetElemConnectivity(NodeCoords,NodeConn):""" Determines the elements connected to each node in the mesh Parameters ---------- NodeCoords : array_like List of nodal coordinates. NodeConn : array_like Nodal connectivity list. Returns ------- ElemConn : list List of elements connected to each node. """nodes,elems=zip(*[(n,i)fori,eleminenumerate(NodeConn)forninelem])NotInMesh=set(range(len(NodeCoords))).difference(nodes)nodes+=tuple(NotInMesh)elems+=tuple(itertools.repeat(-1,len(nodes)))ElemConn=[list(set(elemfor_,elemingroupifelem!=-1))fornode,groupinitertools.groupby(sorted(zip(nodes,elems),key=lambdax:x[0]),key=lambdax:x[0])]returnElemConn
[docs]defgetNodeNeighborhood(NodeCoords,NodeConn,nRings):""" Gives the connected nodes in an n ring neighborhood for each node in the mesh Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list Nodal connectivity list. nRings : int Number of rings to include. Returns ------- NodeNeighborhoods : list List of neighboring nodes in an n ring neighborhood around each node in NodeCoords. """NodeNeighbors=getNodeNeighbors(NodeCoords,NodeConn)NodeNeighborhoods=[[jforjinNodeNeighbors[i]]foriinrange(len(NodeNeighbors))]ifnRings==1:returnNodeNeighborhoodselse:forninrange(nRings-1):# For each ring, loop through and add the neighbors of the nodes in the neighborhood to the neighborhoodforiinrange(len(NodeNeighborhoods)):temp=[jforjinNodeNeighborhoods[i]]forjintemp:forkinrange(len(NodeNeighbors[j])):if(NodeNeighbors[j][k]notinNodeNeighborhoods[i])and(NodeNeighbors[j][k]!=i):NodeNeighborhoods[i].append(NodeNeighbors[j][k])returnNodeNeighborhoods
[docs]defgetNodeNeighborhoodByRadius(NodeCoords,NodeConn,Radius):""" Gives the connected nodes in a neighborhood with a specified radius for each node in the mesh. Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list Nodal connectivity list. radius : float Radius of around each node. Returns ------- NodeNeighborhoods : list List of neighboring nodes in an neighborhood around each node in NodeCoords with the neighborhoods specified by a radius. """NodeNeighbors=getNodeNeighbors(NodeCoords,NodeConn)NodeNeighborhoods=[[]foriinrange(len(NodeNeighbors))]foriinrange(len(NodeNeighborhoods)):thisNode=NodeCoords[i]thinking=TrueNodeNeighborhoods[i]=[jforjinNodeNeighbors[i]if \
np.sqrt((thisNode[0]-NodeCoords[j][0])**2+ \
(thisNode[1]-NodeCoords[j][1])**2+ \
(thisNode[2]-NodeCoords[j][2])**2)<=Radius]whilethinking:thinking=Falsetemp=[jforjinNodeNeighborhoods[i]]forjintemp:forkinrange(len(NodeNeighbors[j])):if(NodeNeighbors[j][k]notinNodeNeighborhoods[i])and(NodeNeighbors[j][k]!=i):otherNode=NodeCoords[NodeNeighbors[j][k]]ifnp.sqrt((thisNode[0]-otherNode[0])**2+(thisNode[1]-otherNode[1])**2+(thisNode[2]-otherNode[2])**2)<=Radius:thinking=TrueNodeNeighborhoods[i].append(NodeNeighbors[j][k])returnNodeNeighborhoods
[docs]defgetElemNeighbors(NodeCoords,NodeConn,mode='face',ElemConn=None):""" Get list of neighboring elements for each element in the mesh. Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list List of nodal connectivities. mode : str, optional Neighbor mode, will determine what type of connectivity constitutes an element neighbor, by default 'face'. 'node' : Any elements that share at least one node are considered neighbors. TODO: Not currently implemented. 'edge' : Any elements that share an edge are considered neighbors. 'face' : Any elements that share a face are considered neighbors. NOTE that in surface meshes, no elements share faces. ElemConn : list, optional Node-Element connectivity of the mesh as obtained by getNodeNeighbors. If supplied, won't require an additional call to getNodeNeighbors. Only relevant if mode == 'node', by default None. Returns ------- ElemNeighbors : list List of element neighbors. For each element, there is a list of the indices of the neighboring elements. """# Get Element neighbors ifmode=='node':ElemNeighbors=[set()foriinrange(len(NodeConn))]ElemConn=getElemConnectivity(NodeCoords,NodeConn)fori,eleminenumerate(NodeConn):forninelem:ElemNeighbors[i].update(ElemConn[n])ElemNeighbors=[list(s)forsinElemNeighbors]elifmode=='edge':ElemNeighbors=[set()foriinrange(len(NodeConn))]Edges,EdgeConn,EdgeElem=converter.solid2edges(NodeCoords,NodeConn,return_EdgeElem=True,return_EdgeConn=True)UEdges,idx,inv=converter.edges2unique(Edges,return_idx=True,return_inv=True)inv=np.append(inv,-1)UEdgeConn=inv[PadRagged(EdgeConn)]UEdgeElem=EdgeElem[idx]EdgeElemConn=np.nan*(np.ones((len(UEdges),2)))# Elements attached to each edger=np.repeat(np.arange(len(UEdgeConn))[:,None],UEdgeConn.shape[1],axis=1)EECidx=(UEdgeElem[UEdgeConn]==r).astype(int)EdgeElemConn[UEdgeConn,EECidx]=rforiinrange(len(EdgeElemConn)):ifnotany(np.isnan(EdgeElemConn[i])):ElemNeighbors[int(EdgeElemConn[i][0])].add(int(EdgeElemConn[i][1]))ElemNeighbors[int(EdgeElemConn[i][1])].add(int(EdgeElemConn[i][0]))ElemNeighbors=[list(s)forsinElemNeighbors]elifmode=='face':#TODO: This needs to updated, can be made faster, should use converter.faces2uniquefaces,faceconn,faceelem=converter.solid2faces(NodeCoords,NodeConn,return_FaceConn=True,return_FaceElem=True)####### if v == 1:sortface=[tuple(sorted(face))forfaceinfaces]FaceElem=collections.defaultdict(set)fori,facekeyinenumerate(sortface):FaceElem[facekey].add(faceelem[i])ElemNeighborDict=dict()fori,fsinenumerate(faceconn):neighbors={elemforfinfsforeleminFaceElem[sortface[f]]}neighbors.discard(i)ElemNeighborDict[i]=neighborsElemNeighbors=[list(ElemNeighborDict[i])foriinrange(len(faceconn))]###### elif v == 0:# ElemNeighbors = [set() for i in range(len(NodeConn))]# # Pad Ragged arrays in case of mixed-element meshes# Rfaces = PadRagged(faces)# Rfaceconn = PadRagged(faceconn)# # Get all unique element faces (accounting for flipped versions of faces)# _,idx,inv = np.unique(np.sort(Rfaces,axis=1),axis=0,return_index=True,return_inverse=True)# RFaces = Rfaces[idx]# FaceElem = faceelem[idx]# RFaces = np.append(RFaces, np.repeat(-1,RFaces.shape[1])[None,:],axis=0)# inv = np.append(inv,-1)# RFaceConn = inv[Rfaceconn] # Faces attached to each element# # Face-Element Connectivity# FaceElemConn = np.nan*(np.ones((len(RFaces),2)))# FECidx = (FaceElem[RFaceConn] == np.repeat(np.arange(len(NodeConn))[:,None],RFaceConn.shape[1],axis=1)).astype(int)# FaceElemConn[RFaceConn,FECidx] = np.repeat(np.arange(len(NodeConn))[:,None],RFaceConn.shape[1],axis=1)# FaceElemConn = [[int(x) if not np.isnan(x) else x for x in y] for y in FaceElemConn[:-1]]# for i in range(len(FaceElemConn)):# if np.any(np.isnan(FaceElemConn[i])): continue# ElemNeighbors[FaceElemConn[i][0]].add(FaceElemConn[i][1])# ElemNeighbors[FaceElemConn[i][1]].add(FaceElemConn[i][0])# ElemNeighbors = [list(s) for s in ElemNeighbors] else:raiseException('Invalid mode. Must be "edge" or "face".')returnElemNeighbors
[docs]defgetConnectedNodes(NodeCoords,NodeConn,NodeNeighbors=None,BarrierNodes=set()):""" Identifies groups of connected nodes. For a fully connected mesh, a single region will be identified Parameters ---------- NodeCoords : list of lists List of nodal coordinates. NodeConn : list of lists Nodal connectivity list. NodeNeighbors : list, optional List of neighboring nodes for each node in NodeCoords. The default is None. If no value is provided, it will be computed with getNodeNeighbors BarrierNodes : set, optional Set of nodes that can separate regions. Returns ------- NodeRegions : list of sets Each set in the list contains a region of connected nodes. Sorted by size of region such that the region with the most nodes is first in the list. """NodeRegions=[]ifnotNodeNeighbors:NodeNeighbors=getNodeNeighbors(NodeCoords,NodeConn)iflen(BarrierNodes)>0:NodeNeighbors=[[]ifiinBarrierNodeselsenfori,ninenumerate(NodeNeighbors)]NeighborSets=[set(n)forninNodeNeighbors]AllNodes=set(range(len(NodeCoords)))DetachedNodes=AllNodes.difference(set(np.unique([nforeleminNodeConnforninelem])))todo=AllNodes.difference(DetachedNodes).difference(BarrierNodes)whilelen(todo)>0:seed=todo.pop()region={seed}new={seed}nOld=0nCurrent=len(region)k=0whilenOld!=nCurrent:k+=1nOld=nCurrentold=copy.copy(new)new=set()foriinold:new.update(NeighborSets[i])new.difference_update(region)region.update(new)nCurrent=len(region)todo.difference_update(region)NodeRegions.append(region)NodeRegions=[NodeRegions[i]foriinnp.argsort([len(region)forregioninNodeRegions])[::-1]]returnNodeRegions
[docs]defgetConnectedElements(NodeCoords,NodeConn,ElemNeighbors=None,mode='edge',BarrierElems=set()):""" Identifies groups of connected nodes. For a fully connected mesh, a single region will be identified Parameters ---------- NodeCoords : list of lists List of nodal coordinates. NodeConn : list of lists Nodal connectivity list. ElemNeighbors : list, optional List of neighboring elements for each element in NodeConn. The default is None. If no value is provided, it will be computed with getNodeNeighbors mode : str, optional Connectivity method to be used for getElemNeighbors. The default is 'edge'. BarrierElems : set, optional Set of barrier elements that the connected region cannot move past. They can be included in a region, but will not connect to their neighbors Returns ------- ElemRegions : list of sets Each set in the list contains a region of connected nodes. Sorted by size of region such that the region with the most nodes is first in the list. """ElemRegions=[]ifnotElemNeighbors:ElemNeighbors=getElemNeighbors(NodeCoords,NodeConn,mode=mode)iflen(BarrierElems)>0:ElemNeighbors=[[]ifiinBarrierElemselseefori,einenumerate(ElemNeighbors)]NeighborSets=[set(n)forninElemNeighbors]todo=set(range(len(NodeConn))).difference(BarrierElems)whilelen(todo)>0:seed=todo.pop()region={seed}new={seed}nOld=0nCurrent=len(region)k=0whilenOld!=nCurrent:k+=1nOld=nCurrentold=copy.copy(new)new=set()foriinold:new.update(NeighborSets[i])new.difference_update(region)region.update(new)nCurrent=len(region)todo.difference_update(region)ElemRegions.append(region)ElemRegions=[ElemRegions[i]foriinnp.argsort([len(region)forregioninElemRegions])[::-1]]returnElemRegions
[docs]defCentroids(NodeCoords,NodeConn):""" Calculate element centroids. Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list List of nodal connectivities. Returns ------- centroids : list list of element centroids. """iflen(NodeConn)==0:return[]ArrayCoords=np.vstack([NodeCoords,[np.nan,np.nan,np.nan]])R=PadRagged(NodeConn,fillval=-1)Points=ArrayCoords[R]centroids=np.nanmean(Points,axis=1)returncentroids
[docs]defCalcFaceNormal(NodeCoords,SurfConn):""" Calculates normal vectors on the faces of a triangular surface mesh. Assumes triangles are in counter-clockwise when viewed from the outside Parameters ---------- NodeCoords : list List of nodal coordinates. SurfConn : list Nodal connectivity list of a triangular surface mesh. Returns ------- ElemNormals list List of element surface normals . """ArrayCoords=np.asarray(NodeCoords)_,TriConn,inv=converter.surf2tris(NodeCoords,SurfConn,return_inv=True)points=ArrayCoords[TriConn]TriNormals=_tri_normals(points)ElemNormals=np.zeros((len(SurfConn),3))np.add.at(ElemNormals,inv,TriNormals)ElemNormals/=np.bincount(inv)[:,None]returnElemNormals
[docs]defFace2NodeNormal(NodeCoords,NodeConn,ElemConn,ElemNormals,method='Angle'):""" Calculate node normal vectors based on the element face normals. Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list Nodal connectivity list. ElemConn : list List of elements connected to each node. ElemNormals : list List of element normal vectors. method : str, optional Method used to determine node normals. The default is 'Angle'. - Angle - performs an angle weighted average of connected element normals :cite:p:`Thurrner1998` - Average - performs a simple averaging of connected element normals - MostVisible - determines the most visible normal :cite:p:`Aubry2008a` - MostVisible_Loop - non-vectorized version of MostVisible, slower but more readable - MostVisible_Iter - iterative method for determining the most visible normal :cite:p:`Aubry2008a` MostVisible_Loop and MostVisible_Iter are included for completeness, but in general, MostVisible should be used instead. Returns ------- NodeNormals : list Unit normal vectors for each node. """if(method.lower()=='angle'):# Based on: Grit Thürrner & Charles A. Wüthrich (1998)# Perform angle weighted average to compute vertex normals# Calculate the angles to use as weight# Cast ElemConn into a rectangular matrix# Warning: This code is very vectorized - it might be difficult to debugNodeSet=np.unique(PadRagged(NodeConn,fillval=-1))if-1inNodeSet:NodeSet=np.delete(NodeSet,0)ArrayCoords=np.vstack([NodeCoords,[np.nan,np.nan,np.nan]])R=PadRagged(ElemConn,fillval=-1)[NodeSet]Mask0=(R>=0).astype(int)Masknan=Mask0.astype(float)Masknan[Mask0==0]=np.nanNs=np.vstack([ElemNormals,[np.nan,np.nan,np.nan]])[R]RNodeConn=PadRagged(NodeConn,fillval=-1)ArrayConn=np.vstack([RNodeConn,-1*np.ones((1,RNodeConn.shape[1]),dtype=int)])IncidentNodes=ArrayConn[R]x=(ArrayCoords[IncidentNodes]-ArrayCoords[NodeSet,None,None,:])x[np.all(x==[0,0,0],axis=-1)]=np.nan# For each node and for each incident element on that node, dot product of the two edges of the element that meet at the nodedots=np.sum((np.nanprod(x,axis=2)*Mask0[:,:,None]),axis=2)# For each node and for each incident element on that node, the product of the norms of the two edges of the element that meet at the nodenorms=np.nanprod(np.linalg.norm(x,axis=3),axis=2)# cos(alpha) = dot(u,v)/(norm(u)*norm(v))cosAlpha=dots/normsalpha=np.arccos(cosAlpha)*MasknansumAlphaN=np.nansum(alpha[:,:,None]*Ns,axis=1)NodeNormals=np.nan*np.ones_like(NodeCoords)NodeNormals[NodeSet]=sumAlphaN/np.linalg.norm(sumAlphaN,axis=1)[:,None]elif(method.lower()=='average')or(method=='none')or(method==None):# Cast ElemConn into a rectangular matrixNodeSet=np.unique(PadRagged(NodeConn,fillval=-1))R=PadRagged(ElemConn,fillval=-1)[NodeSet]Ns=np.array(np.append(ElemNormals,[[np.nan,np.nan,np.nan]],axis=0))[R]NodeNormals=np.nan*np.ones_like(NodeCoords)NodeNormals[NodeSet]=np.nanmean(Ns,axis=1)NodeNormals[NodeSet]=NodeNormals[NodeSet]/np.linalg.norm(NodeNormals[NodeSet],axis=1)[:,None]elifmethod=='MostVisible':# Note: this code uses dot(Ni,Nj) as a surrogate for radius; since Ni,Nj are both unit vectors# cos(theta) = dot(Ni,Nj) -> theta = arccos(dot(Ni,Nj)). Since arccos is a monotonically # decreasing function, if dot(Ni,Nj) < dot(Ni,Nk), then rij > rikeps=-1e-8NodeSet=np.unique(PadRagged(NodeConn,fillval=-1))if-1inNodeSet:NodeSet=np.delete(NodeSet,0)R=PadRagged(ElemConn,fillval=-1)ifnp.shape(R)[1]<3:# Handling for case where no nodes have more than 2 connected elementstempR=-1*np.ones((len(ElemConn),3),dtype='int')tempR[:,:np.shape(R)[1]]=RR=tempRNs=np.vstack([ElemNormals,[np.nan,np.nan,np.nan]])[R]# 2 Point Circlesscalmin=-1Combos2=Ns[NodeSet][:,np.array(list(itertools.combinations(range(Ns.shape[1]),2)))]Nb=np.sum(Combos2,axis=2)Nb=Nb/np.linalg.norm(Nb,axis=2)[:,:,None]scal2=np.sum(Nb*Combos2[:,:,0,:],axis=2)# 3 Point CirclesCombos3=Ns[NodeSet][:,np.array(list(itertools.combinations(range(Ns.shape[1]),3)))]Ni=Combos3[:,:,0,:]Nj=Combos3[:,:,1,:]Nk=Combos3[:,:,2,:]denom=2*np.linalg.norm(np.cross(Ni-Nk,Nj-Nk),axis=2)**2withnp.errstate(divide='ignore',invalid='ignore'):Nc=np.cross(((np.linalg.norm(Ni-Nk,axis=2)**2)[:,:,None]*(Nj-Nk))-(np.linalg.norm(Nj-Nk,axis=2)**2)[:,:,None]*(Ni-Nk),np.cross(Ni-Nk,Nj-Nk))/denom[:,:,None]+NkNc=Nc/np.linalg.norm(Nc,axis=2)[:,:,None]scal3=np.sum(Nc*Ni,axis=2)Nc[scal3<0]=-Nc[scal3<0]scal3[scal3<0]=-scal3[scal3<0]scal23=np.hstack([scal2,scal3])Nbc=np.hstack([Nb,Nc])check=np.any((np.einsum('lij,ljk->lik',Nbc,np.swapaxes(Ns[NodeSet],1,2))-scal23[:,:,None])<eps,axis=2)scal23[check]=scalmin# Indices of the smallest radius that contains all pointsIdx=scal23==np.nanmax(scal23,axis=1)[:,None]# In case of duplicates, only taking the first onenewIdx=np.zeros_like(Idx)newIdx[np.arange(len(Idx)),Idx.argmax(axis=1)]=Idx[np.arange(len(Idx)),Idx.argmax(axis=1)]NodeNormals=np.nan*np.ones_like(NodeCoords)NodeNormals[NodeSet]=Nbc[newIdx]else:raiseValueError(f'Invalid method: {method:s}')# NodeNormals = [[] for i in range(len(NodeCoords))] # Normal vectors for each vertex# NodeSet = {n for elem in NodeConn for n in elem}# for i in range(len(NodeCoords)):# if i not in NodeSet:# NodeNormals[i] = [np.nan,np.nan,np.nan]# continue# angles = [0 for j in range(len(ElemConn[i]))]# elemnormals = [np.array(ElemNormals[elem]) for elem in ElemConn[i]]# if method == 'MostVisible_Loop':# # This is kept for readability; 'MostVisible' is a vectorized equivalent that performs significantly faster# # Note: this code uses dot(Ni,Nj) as a surrogate for radius; since Ni,Nj are both unit vectors# # cos(theta) = dot(Ni,Nj) -> theta = arccos(dot(Ni,Nj)). Since arccos is a monotonically # # decreasing function, if dot(Ni,Nj) < dot(Ni,Nk), then rij > rik# eps = -1e-8# scalmin = -1# C = [np.nan,np.nan,np.nan]# for ii in range(len(elemnormals)-1):# # Check the 2 point circles# Ni = np.array(elemnormals[ii])# for j in range(ii+1,len(elemnormals)):# Nj = np.array(elemnormals[j])# Nb = Ni+Nj# Nb = Nb/np.linalg.norm(Nb)# scal = np.dot(Nb,Ni)# if scal < scalmin: # pass# elif any((np.dot(Nl,Nb) - scal) < eps for Nl in elemnormals):# pass# else:# C = Nb.tolist()# scalmin = scal# for ii in range(len(elemnormals)-2): # # Check the 3 point circles# Ni = elemnormals[ii]# for j in range(ii+1,len(elemnormals)-1):# Nj = elemnormals[j]# for k in range(j+1,len(elemnormals)):# Nk = elemnormals[k]# denom = (2*np.linalg.norm(np.cross(Ni-Nk,Nj-Nk))**2) # if denom == 0:# continue# Nc = np.cross(np.linalg.norm(Ni-Nk)**2 * (Nj-Nk) - np.linalg.norm(Nj-Nk)**2 * (Ni-Nk), np.cross(Ni-Nk,Nj-Nk))/denom + Nk# nNc = np.linalg.norm(Nc)# if nNc == 0:# continue# Nc = Nc/nNc# scal = np.dot(Nc, Ni)# if scal < 0:# Nc = [-1*n for n in Nc]# scal = -scal# if scal < scalmin:# pass# elif any((np.dot(Nl,Nc) - scal) < eps for Nl in elemnormals):# pass # else:# C = Nc# scalmin = scal# NodeNormals[i] = C # if np.any(np.isnan(C)):# print(i)# elif method == 'MostVisible_Iter':# conv = 1e-3# beta = 0.5# # Initial weights# ws = [1/len(elemnormals) for i in range(len(elemnormals))]# # Compute initial guess normal# Sp = sum([w*n for w,n in zip(ws,elemnormals)])# Np = Sp/np.linalg.norm(Sp)# k = 0# thinking = True# while thinking:# k+=1# alphas = [np.arccos(np.clip(np.dot(Np,Ni),-1,1)) for Ni in elemnormals]# Salpha = sum(alphas)# if Salpha == 0:# thinking = False# else:# ws = [w*alpha/Salpha for w,alpha in zip(ws,alphas)]# Sw = sum(ws)# ws = [w/Sw for w in ws]# Spnew = sum([w*n for w,n in zip(ws,elemnormals)])# if np.linalg.norm(Spnew) == 0:# print('merp3')# Npnew = Spnew/np.linalg.norm(Spnew)# # Relax# Nprel = beta*Npnew + (1-beta)*Np# if np.linalg.norm(Np-Nprel) < conv or k > 100:# thinking = False# Np = Nprel# if any(np.isnan(Np)) and len(elemnormals)>0:# merp = 2# NodeNormals[i] = Np.tolist()returnNodeNormals
[docs]@try_njitdefBaryTri(Nodes,Pt):""" Returns the bary centric coordinates of a point (Pt) relative to a triangle (Nodes) Parameters ---------- Nodes : np.ndarray List of coordinates of the triangle vertices. Pt : np.ndarray Coordinates of the point. Returns ------- alpha : float First barycentric coordinate. beta : float Second barycentric coordinate. gamma : float Third barycentric coordinate. """Nodes=np.asarray(Nodes,dtype=np.float64)Pt=np.asarray(Pt,dtype=np.float64)A=Nodes[0]B=Nodes[1]C=Nodes[2]BA=np.subtract(B,A)# CB = np.subtract(C,B)# AC = np.subtract(A,C)CA=np.subtract(C,A)BABA=np.dot(BA,BA)BACA=np.dot(BA,CA)CACA=np.dot(CA,CA)PABA=np.dot(np.subtract(Pt,A),BA)PACA=np.dot(np.subtract(Pt,A),CA)d=(BABA*CACA-BACA*BACA)denom=1/dbeta=(CACA*PABA-BACA*PACA)*denomgamma=(BABA*PACA-BACA*PABA)*denomalpha=1-gamma-betareturnalpha,beta,gamma
[docs]defBaryTris(Tris,Pt):""" Returns the barycentric coordinates of a point or points relative to a triangle. This can either compare a set of n triangles to a single point, or pairwise comparison between n triangles and n points. Parameters ---------- Tris : array_like nx3x3 coordinates of the vertices. The array should be formatted as if obtained by indexing NodeCoords[NodeConn] for a purely triangular mesh. Pt : array_like Coordinates of the point or points. For a single point, this should have the a shape = (3,), for a set of points, this should have a shape = (n,3) where n is equal to the number of triangles. Returns ------- alpha : float First barycentric coordinate. beta : float Second barycentric coordinate. gamma : float Third barycentric coordinate. """A=Tris[:,0]B=Tris[:,1]C=Tris[:,2]BA=np.subtract(B,A)# CB = np.subtract(C,B)# AC = np.subtract(A,C)CA=np.subtract(C,A)BABA=np.sum(BA*BA,axis=1)BACA=np.sum(BA*CA,axis=1)CACA=np.sum(CA*CA,axis=1)PABA=np.sum(np.subtract(Pt,A)*BA,axis=1)PACA=np.sum(np.subtract(Pt,A)*CA,axis=1)withnp.errstate(divide='ignore',invalid='ignore'):denom=1/(BABA*CACA-BACA*BACA)beta=(CACA*PABA-BACA*PACA)*denomgamma=(BABA*PACA-BACA*PABA)*denomalpha=1-gamma-beta;returnalpha,beta,gamma
[docs]@try_njitdefBaryTet(Nodes,Pt):""" Returns the bary centric coordinates of a point (Pt) relative to a tetrahedron (Nodes) Parameters ---------- Nodes : list List of coordinates of the tetrahedral vertices. Pt : list Coordinates of the point. Returns ------- alpha : float First barycentric coordinate. beta : float Second barycentric coordinate. gamma : float Third barycentric coordinate. delta : float Fourth barycentric coordinate. """A=Nodes[0]B=Nodes[1]C=Nodes[2]D=Nodes[3]T=np.array([[A[0]-D[0],B[0]-D[0],C[0]-D[0]],[A[1]-D[1],B[1]-D[1],C[1]-D[1]],[A[2]-D[2],B[2]-D[2],C[2]-D[2]]])alpha,beta,gamma=np.linalg.solve(T,np.subtract(Pt,D))delta=1-(alpha+beta+gamma)returnalpha,beta,gamma,delta
[docs]defProject2Surface(Points,Normals,NodeCoords,SurfConn,tol=np.inf,Octree='generate'):""" Projects a node along its normal vector onto a surface. Returns the index of the element (elemID) that contains the projected node and the barycentric coordinates (alpha, beta, gamma) of that projection within that element. Parameters ---------- Point : list or np.ndarray Coordinates of the point to be projected on to the surface. Normal : list or np.ndarray Vector along which the point will be projected. NodeCoords : list or np.ndarray Node coordinates list of the mesh that the point is being projected to. SurfConn : list or np.ndarray Nodal connectivity of the surface mesh that the point is being projected to. tol : float, optional Tolerance value, if the projection distance is greater than tol, the projection will be exculded, default is np.inf Octree : str (or octree.OctreeNode), optional octree options. An octree representation of the surface can significantly improve mapping speeds, by default 'generate'. 'generate' - Will generate an octree for use in surface mapping. 'none' or None - Won't generate an octree and will use a brute force approach. octree.OctreeNode - Provide a precompute octree structure corresponding to the surface mesh. Should be created by octree.Surface2Octree(NodeCoords,SurfConn) Returns ------- MappingMatrix : np.ndarray nx4 array consisting of the element ID (in SurfConn) and three barycentric coordinates (alpha, beta, gamma) for each point in Points """iftype(NodeCoords)islist:NodeCoords=np.array(NodeCoords)iftype(SurfConn)islist:SurfConn=np.array(SurfConn)intersections,distances,intersectionPts=rays.RaysSurfIntersection(Points,Normals,NodeCoords,SurfConn,Octree=Octree)argmindist=[np.argmin(np.abs(x))iflen(x)>0else-1forxindistances]mindist=np.array([x[argmindist[i]]iflen(x)>0elsenp.inffori,xinenumerate(distances)])elemID=np.array([intersections[i][argmindist[i]]iflen(x)>0else-1fori,xinenumerate(distances)])ps=np.array([intersectionPts[i][argmindist[i]]iflen(x)>0else[np.nan,np.nan,np.nan]fori,xinenumerate(distances)])mappedbool=(elemID>=0)&(mindist<=tol)alphas,betas,gammas=BaryTris(NodeCoords[SurfConn[elemID[mappedbool]]],ps[mappedbool,:])alpha=-1*np.ones(len(Points))beta=-1*np.ones(len(Points))gamma=-1*np.ones(len(Points))alpha[mappedbool]=alphasbeta[mappedbool]=betasgamma[mappedbool]=gammasMappingMatrix=np.column_stack([elemID,alpha,beta,gamma])returnMappingMatrix
[docs]defSurfMapping(NodeCoords1,SurfConn1,NodeCoords2,SurfConn2,tol=np.inf,verbose=False,Octree='generate',return_octree=False,npts=np.inf):""" Generate a mapping matrix from to map data from one surface to another using barycentric interpolation. Each row of the mapping matrix contains an element ID followed by barycentric coordinates alpha, beta, gamma that define the position of the nodes of surface 1 (NodeCoords1) relative to the specified surface element of surface 2 (SurfConn2). An element ID of -1 indicates a failed mapping. NOTE: Only triangular surface meshes are supported. Parameters ---------- NodeCoords1 : list List of nodal coordinates. SurfConn1 : list List of nodal connectivities. NodeCoords2 : list List of nodal coordinates. SurfConn2 : list List of nodal connectivities. tol : float, optional Tolerance value, if the projection distance is greater than tol, the projection will be exculded, default is np.inf verbose : bool, optional If true, will print mapping statistics, by default False. Octree : str (or octree.OctreeNode), optional octree options. An octree representation of surface 2 can significantly improve mapping speeds, by default 'generate'. 'generate' - Will generate an octree for use in surface mapping. 'none' or None - Won't generate an octree and will use a brute force approach. octree.OctreeNode - Provide a precompute octree structure corresponding to surface 2. Should be created by octree.Surface2Octree(NodeCoords2,SurfConn2) return_octree : bool, optional If true, will return the generated or provided octree, by default False. npts : int, optional Number of points to map. A random sample of min(npts, len(NodeCoords1)) from Surface 1 will be mapped , by default np.inf (all points). Returns ------- MappingMatrix : list min(npts, len(NodeCoords1))x4 matrix of of barycentric coordinates, defining NodeCoords1 in terms of the triangular surface elements of Surface 2. Octree : octree.OctreeNode, optional The generated or provided octree structure corresponding to Surface 2. """iftype(NodeCoords1)islist:NodeCoords1=np.array(NodeCoords1)iftype(NodeCoords2)islist:NodeCoords2=np.array(NodeCoords2)iftype(SurfConn1)islist:SurfConn1=np.array(SurfConn1)iftype(SurfConn2)islist:SurfConn2=np.array(SurfConn2)Surf1Nodes=np.unique(SurfConn1.flatten())ifnpts>=len(Surf1Nodes):N=len(NodeCoords1)NodeIds=Surf1Nodeselse:N=nptsidx=np.random.choice(range(len(Surf1Nodes)),size=N,replace=False)NodeIds=Surf1Nodes[idx]assertSurfConn1.shape[1]==SurfConn2.shape[1]==3,'Currently only triangular surfaces are supported.'ElemConn1=getElemConnectivity(NodeCoords1,SurfConn1)ElemNormals1=CalcFaceNormal(NodeCoords1,SurfConn1)NodeNormals1=Face2NodeNormal(NodeCoords1,SurfConn1,ElemConn1,ElemNormals1,method='angle')ifOctree=='generate':Octree=octree.Surface2Octree(NodeCoords2,SurfConn2)MappingMatrix=-1*np.ones((len(NodeCoords1),4))MappingMatrix[NodeIds,:]=Project2Surface(NodeCoords1[NodeIds,:],NodeNormals1[NodeIds,:],NodeCoords2,SurfConn2,tol=tol,Octree=Octree)ifverbose:failcount=np.sum(MappingMatrix[NodeIds,0]==-1)print('{:.3f}% of nodes mapped'.format((len(NodeIds)-failcount)/len(NodeIds)*100))ifreturn_octree:returnMappingMatrix,OctreereturnMappingMatrix
[docs]defValueMapping(NodeCoords1,SurfConn1,NodeVals1,NodeCoords2,SurfConn2,tol=np.inf,Octree='generate',MappingMatrix=None,verbose=False,return_MappingMatrix=False,return_octree=False,npts=np.inf):""" Maps nodal values one surface to another. This currently only supports triangluar surface meshes TODO: Multi-value mapping may produce errors - need to better verify. Parameters ---------- NodeCoords1 : List of lists Contains coordinates for each node in surface 1. Ex. [[x1,y1,z1],...] SurfConn1 : List of lists Contains the nodal connectivity defining the surface elements. NodeVals1 : List or List of lists Scalar nodal values associated with surface 1. For multiple values: [[x1,x2,x3,...],[y1,y2,y3,...],[z1,z2,z3,...],...] NodeCoords2 : List of lists Contains coordinates for each node in surface 2. Ex. [[x1,y1,z1],...]. SurfConn2 : List of lists Contains the nodal connectivity defining the surface elements. tol : float, optional Tolerance value, if the projection distance is greater than tol, the projection will be exculded, default is np.inf Octree : str (or octree.OctreeNode), optional octree options. An octree representation of surface 1 can significantly improve mapping speeds, by default 'generate'. 'generate' - Will generate an octree for use in surface mapping. 'none' or None - Won't generate an octree and will use a brute force approach. octree.OctreeNode - Provide a precompute octree structure corresponding to surface 1. Should be created by octree.Surface2Octree(NodeCoords1,SurfConn1) MappingMatrix : list len(NodeCoords2)x4 matrix of of barycentric coordinates, defining NodeCoords2 in terms of the triangular surface elements of Surface 1. verbose : bool, optional If true, will print mapping statistics, by default False. return_MappingMatrix : bool, optional If true, will return MappingMatrix, by default False. return_octree : bool, optional If true, will return generated or provided octree, by defualt False. NOTE if MappingMatrix is provided, the octree structure won't be generated. In this cases, if Octree='generate' and return_octree=True, the returned value for octree will simply be the string 'generate'. npts : int, optional Number of points to map. Values from Surface 1 will be mapped to random sample of min(npts, len(NodeCoords2) in Surface 2, by default np.inf (all points). Returns ------- NodeVals2 : List Scalar nodal values associated with surface 2, mapped from surface 1. """# if type(NodeVals1[0]) is list or type(NodeVals1[0]) is np.ndarray:# singleVal = False# # NodeVals2 = [[0 for j in range(len(NodeCoords2))] for i in range(len(NodeVals1))]# else:# singleVal = True# NodeVals2 = [0 for i in range(len(NodeCoords2))]# Map the coordinates from surface 2 to surface 1ifMappingMatrixisNone:MappingMatrix,Octree=SurfMapping(NodeCoords2,SurfConn2,NodeCoords1,SurfConn1,Octree=Octree,tol=tol,verbose=verbose,return_octree=True,npts=npts)# if singleVal:iflen(np.shape(NodeVals1))==1:# 1D data_NodeVals1=np.append(NodeVals1,np.nan)alpha=MappingMatrix[:,1]beta=MappingMatrix[:,2]gamma=MappingMatrix[:,3]else:# ND data_NodeVals1=np.append(NodeVals1,[np.repeat(np.nan,np.shape(NodeVals1)[1])],axis=0)alpha=MappingMatrix[:,1][:,None]beta=MappingMatrix[:,2][:,None]gamma=MappingMatrix[:,3][:,None]# NodeVals2 = np.nan*np.ones(np.shape(NodeVals1))elemID=MappingMatrix[:,0].astype(int)ArrayConn=np.append(SurfConn1,[[-1,-1,-1]],axis=0)NodeVals2=alpha*_NodeVals1[ArrayConn[elemID][:,0]]+ \
beta*_NodeVals1[ArrayConn[elemID][:,1]]+ \
gamma*_NodeVals1[ArrayConn[elemID][:,2]]ifreturn_MappingMatrixandreturn_octree:returnNodeVals2,MappingMatrix,Octreeelifreturn_MappingMatrix:returnNodeVals2,MappingMatrixelifreturn_octree:returnNodeVals2,OctreereturnNodeVals2
[docs]defDeleteDuplicateNodes(NodeCoords,NodeConn,tol=1e-12,return_idx=False,return_inv=False):""" Remove nodes that are duplicated in the mesh, either at exactly the same location as another node or a distance < tol apart. Nodes are renumbered and elements reconnected such that the geometry and structure of the mesh remains unchanged. Parameters ---------- NodeCoords : list Contains coordinates for each node. Ex. [[x1,y1,z1],...] NodeConn : list Nodal connectivity list. tol : float, optional Tolerance value to be used when determining if two nodes are the same. The default is 1e-14. return_idx : bool, optional Returns the indices of each row of NodeCoords in the order that they're sorted place into the new array, by default False. return_inv : bool, optional Returns the indices that reverse the operation, by default False. Returns ------- NewCoords : list Updated node coordinates without duplicates. NewConn : list Updated node connectivity without duplicate nodes. idx : np.ndarray Array of indices that convert from the original node coordinates to the new node coordinates (NewCoords = [NodeCoords[i] for i in idx]) inv : np.ndarray Array of indices that can reverse the operation to convert from the new node coordinates to old node coordinates (NodeCoords = [NewCoords[i] for i in inv]). Examples -------- >>> NodeCoords = [[0.,0.,0.], [0.,1.,0.], [1.,1.,0.], [0.,0.,0.], [1.,1.,0.], [1.,0.,0.]] >>> NodeConn = [[0,1,2],[3,4,5]] >>> NewCoords, NewConn, idx, inv = utils.DeleteDuplicateNodes(NodeCoords,NodeConn, return_idx=True,return_inv=True) >>> NewConn [[0, 1, 3], [0, 3, 2]] >>> NewCoords == [NodeCoords[i] for i in idx] True >>> NodeCoords == [NewCoords[i] for i in inv] True """iflen(NodeCoords)==0:ifreturn_idxandreturn_inv:returnNodeCoords,NodeConn,np.array([]),np.array([])elifreturn_idxorreturn_inv:returnNodeCoords,NodeConn,np.array([])returnNodeCoords,NodeConniftol>0:arrayCoords=np.round(np.array(NodeCoords)/tol)*tolelse:arrayCoords=np.array(NodeCoords)unq,idx,inv=np.unique(arrayCoords,return_index=True,return_inverse=True,axis=0)iftype(NodeCoords)islist:NewCoords=np.asarray(NodeCoords)[idx].tolist()else:NewCoords=np.asarray(NodeCoords)[idx]iflen(NodeConn)>0:tempIds=np.append(inv,-1)R=PadRagged(NodeConn,fillval=-1)NewConn=ExtractRagged(tempIds[R],delval=-1)else:NewConn=NodeConnreturns=(True,True,return_idx,return_inv)returntuple(outputfori,outputinenumerate((NewCoords,NewConn,idx,inv))ifreturns[i])
defRemoveNodes(NodeCoords,NodeConn):""" Removes nodes that aren't held by any element Parameters ---------- NodeCoords : array_like Node coordinates NodeConn : array_like Node connectivity Returns ------- NewNodeCoords : array_like New set of node coordinates where unused nodes have been removed NewNodeConn : array_like Renumbered set of node connectivities to be consistent with NewNodeCoords OriginalIds : np.ndarray The indices the original IDs of the nodes still in the mesh. This can be used to remove entries in associated node data (ex. new_data = old_data[OriginalIds]). """iftype(NodeConn)isnp.ndarray:F=NodeConn.flatten()else:F=np.array([nforeleminNodeConnforninelem])node_mask=np.zeros(len(NodeCoords),dtype=np.uint8)node_mask[F]=1OriginalIds=np.where(node_mask)[0]replace=np.zeros(len(NodeCoords),dtype=int)replace[OriginalIds]=np.arange(len(OriginalIds))NewNodeCoords=np.asarray(NodeCoords)[OriginalIds]New=replace[F]iftype(NodeConn)isnp.ndarray:NewNodeConn=New.reshape(NodeConn.shape)else:Newiter=iter(New)NewNodeConn=[list(itertools.islice(Newiter,len(elem)))foreleminNodeConn]returnNewNodeCoords,NewNodeConn,OriginalIds
[docs]defRelabelNodes(NodeCoords,NodeConn,newIds,faces=None):""" Relabel the nodes in the mesh according to the newIds list Parameters ---------- NodeCoords : array_like List of nodal coordinates. NodeConn : array_like Nodal connectivity list. newIds : list list of node ids where the new index is located at the old index faces : list, optional list of face elements, that will also be relabel, by default None Returns ------- NewCoords : array_like Relabeled of nodal coordinates. NewConn : array_like Relabeled nodal connectivity list. """NewConn=ExtractRagged(np.append(newIds,[-1])[PadRagged(NodeConn)],dtype=int)iffaces!=None:iflen(faces)>0:NewFaces=ExtractRagged(np.append(newIds,[-1])[PadRagged(faces)],dtype=int)else:NewFaces=facesNewCoords=np.nan*np.ones(np.shape(NodeCoords))NewCoords[newIds.astype(int)]=np.array(NodeCoords)iffaces!=None:returnNewCoords,NewConn,NewFaceselse:returnNewCoords,NewConn
[docs]defDeleteDegenerateElements(NodeCoords,NodeConn,tol=1e-12,angletol=1e-3,strict=False):""" Deletes degenerate elements from a mesh. TODO: Currently only valid for triangles. Parameters ---------- NodeCoords : list List of nodal coordinates. NodeConn : list List of nodal connectivities. angletol : float, optional Tolerance value for determining what constitutes a degenerate element, by default 1e-3. Degenerate elements will be those who have an angle greater than or equal to 180-180*angletol (default 179.82 degrees) Returns ------- NodeCoords : list List of nodal coordinates. NodeConn : list List of nodal connectivities. """# Remove elements that have a collapsed edge - i.e. two collinear edgesiflen(NodeConn)==0:returnNodeCoords,NodeConnifstrict:NewCoords=NodeCoordsNewConn=[elemforeleminNodeConniflen(elem)==len(set(elem))]else:NewCoords,NewConn=DeleteDegenerateElements(NodeCoords,NodeConn,strict=True)iflen(NewConn)==0:returnNewCoords,NewConnifangletol==0:warnings.warn("Change to strict=True")thetal=np.pi-np.pi*angletol# Maximum angle threshold defdo_split(NewCoords,NewConn,EdgeSort,ConnSort,AngleSort,i):elem0=NewConn[ConnSort[i,0]]elem1=NewConn[ConnSort[i,1]]NotShared0=set(elem0).difference(EdgeSort[i]).pop()NotShared1=set(elem1).difference(EdgeSort[i]).pop()# Get the node not belonging to the edgeif(AngleSort[i,0]>=thetalandConnSort[i,0]>=0andtype(NewConn[ConnSort[i,0]][0])!=list)and(AngleSort[i,1]>=thetalandConnSort[i,1]>=0andtype(NewConn[ConnSort[i,1]][0])!=list):# Both connected elements are degenerateNewNode=NewCoords[NotShared0]elif(AngleSort[i,0]>=thetalandConnSort[i,0]>=0andtype(NewConn[ConnSort[i,0]][0])!=list):NewNode=NewCoords[NotShared0]elif(AngleSort[i,1]>=thetalandConnSort[i,1]>=0andtype(NewConn[ConnSort[i,1]][0])!=list):NewNode=NewCoords[NotShared1]else:returnNewCoords,NewConnNewId=len(NewCoords)NewCoords=np.vstack([NewCoords,NewNode])ifConnSort[i,0]>=0:whileelem0[0]!=NotShared0:elem0=[elem0[-1]]+elem0[0:-1]# cycle the element definition so that it starts with the non-shared node (Might be unnecessarily slow)NewConn[ConnSort[i,0]]=[[elem0[0],elem0[1],NewId],[elem0[0],NewId,elem0[2]]]ifConnSort[i,1]>=0:whileelem1[0]!=NotShared1:elem1=[elem1[-1]]+elem1[0:-1]NewConn[ConnSort[i,1]]=[[elem1[0],elem1[1],NewId],[elem1[0],NewId,elem1[2]]]returnNewCoords,NewConniftype(NewConn)isnp.ndarray:NewConn=NewConn.tolist()Thinking=Truek=0;maxiter=3whileThinkingandk<3:k+=1NewCoords=np.array(NewCoords)Edges,EdgeConn,EdgeElem=converter.solid2edges(NewCoords,NewConn,return_EdgeConn=True,return_EdgeElem=True)UEdges,UIdx,UInv=converter.edges2unique(Edges,return_idx=True,return_inv=True)UEdgeElem=np.asarray(EdgeElem)[UIdx]UEdgeConn=UInv[PadRagged(EdgeConn)]EECidx=(UEdgeElem[UEdgeConn]==np.repeat(np.arange(len(EdgeConn))[:,None],UEdgeConn.shape[1],axis=1)).astype(int)EdgeElemConn=-1*(np.ones((len(UEdges),2),dtype=int))EdgeElemConn[UEdgeConn,EECidx]=np.repeat(np.arange(len(EdgeConn))[:,None],UEdgeConn.shape[1],axis=1)Edges=np.asarray(Edges);EdgeConn=np.asarray(EdgeConn)EdgeVectors=NewCoords[Edges[:,1]]-NewCoords[Edges[:,0]]EdgeLengths=np.linalg.norm(EdgeVectors,axis=1)ElemVectors=EdgeVectors[EdgeConn]ElemLengths=EdgeLengths[EdgeConn]OppositeAngles=-1*np.ones(ElemLengths.shape)withnp.errstate(divide='ignore',invalid='ignore'):OppositeAngles[:,0]=np.clip(np.sum(ElemVectors[:,2]*-ElemVectors[:,1],axis=1)/(ElemLengths[:,1]*ElemLengths[:,2]),-1,1)OppositeAngles[:,1]=np.clip(np.sum(ElemVectors[:,0]*-ElemVectors[:,2],axis=1)/(ElemLengths[:,0]*ElemLengths[:,2]),-1,1)OppositeAngles[:,2]=np.clip(np.sum(ElemVectors[:,1]*-ElemVectors[:,0],axis=1)/(ElemLengths[:,1]*ElemLengths[:,0]),-1,1)OppositeAngles=np.arccos(OppositeAngles)EdgeOppositeAngles=-1*np.ones((len(UEdges),2))EdgeOppositeAngles[UEdgeConn,EECidx]=OppositeAnglessortkey=np.argsort(EdgeLengths[UIdx])[::-1]LengthSort=EdgeLengths[UIdx][sortkey]AngleSort=EdgeOppositeAngles[sortkey]EdgeSort=np.asarray(UEdges)[sortkey]ConnSort=np.array(EdgeElemConn)[sortkey]AbsLargeAngle=np.any(AngleSort>=thetal,axis=1)todo=np.where(AbsLargeAngle)[0]# Splitsrepeat=Falseforiintodo:iftype(NewConn[ConnSort[i,0]][0])islistortype(NewConn[ConnSort[i,1]][0])islist:repeat=TruecontinueNewCoords,NewConn=do_split(NewCoords,NewConn,EdgeSort,ConnSort,AngleSort,i)NewConn=[elemif(type(elem[0])!=list)elseelem[0]foreleminNewConn]+[elem[1]foreleminNewConnif(type(elem[0])==list)]NewCoords=NewCoords.tolist()ifrepeat:Thinking=Trueelse:Thinking=FalseNewCoords,NewConn=DeleteDuplicateNodes(NewCoords,NewConn,tol=tol)NewCoords,NewConn=DeleteDegenerateElements(NewCoords,NewConn,strict=True)returnNewCoords,NewConn
[docs]defCleanupDegenerateElements(NodeCoords,NodeConn,Type='auto',return_idx=False):""" Checks for elements with degenerate edges and either changes the element type or removes the element depending on how degenerate it is. Elements with less than 3 (for surface meshes) or 4 (for volume meshes) unique nodes will be deleted, others will be reduced (ex. a quad with 3 unique nodes will be converted to a triangle). The ordering of nodes will be kept. This function only changes the mesh of an element in NodeConn has the same node number more than once. For meshes that have two differently numbered nodes at the same location, first use utils.DeleteDuplicateNodes. Parameters ---------- NodeCoords : array_like Node coordinates NodeConn : list, array_like Node connectivity Type : str, optional Specifies whether the mesh contains surface elements (tris, quads) or volume elements (tets, hexs, etc.). Must be either "auto", "surf" or "vol". If "auto", Type will be inferred using :func:`identify_type`. By default "auto". Returns ------- NodeCoords : array_like Node coordinates (these are simply passed through from the input) NewConn : list, array_like Updated node connectivity idx : np.ndarray Array of indices that convert from the original list of elements IDs to the new list of element IDs """defrowunique(NodeConn,min_node):# based on unutbu's answer to https://stackoverflow.com/questions/26958233/numpy-row-wise-unique-elementsPadConn=PadRagged(NodeConn)weight=1j*np.linspace(0,PadConn.shape[1],PadConn.shape[0],endpoint=False)uConn=PadConn+weight[:,np.newaxis]u,ind=np.unique(uConn,return_index=True)uConn=-1*np.ones_like(PadConn)np.put(uConn,ind,PadConn.flat[ind])to_delete=np.sum(uConn!=-1,axis=1)<min_nodeifPadConn.shape[1]>=6:# Special attention need for degenerate wedge elementswedge_rows=np.sum(PadConn!=-1,axis=1)==6wedge2tet=np.where((wedge_rows)&(np.sum(uConn!=-1,axis=1)==4))[0]wedge2pyr=np.where((wedge_rows)&(np.sum(uConn!=-1,axis=1)==5))[0]tetints=np.sum((uConn[wedge2tet,:6]==-1)*2**np.arange(0,6)[::-1],axis=1)# Note that the number of possible cases is much less than the maximum 6 digit binary (63) # since unique always keeps the first occurence of a duplicate, and there can only be two "1"s# Cases where a quad face has collapsed make the pyramid degenerate plane, should be removedto_delete[wedge2tet[np.isin(tetints,(9,18))]]=True# Reordering for proper tetsuConn[wedge2tet[tetints==10],:6]=uConn[wedge2tet[tetints==10]][:,[0,3,1,5,2,4]]# node 2, 4 removeduConn[wedge2tet[tetints==12],:6]=uConn[wedge2tet[tetints==12]][:,[0,4,1,5,2,3]]# node 2, 3 removed# Okay cases: 3, 5, 6, 17, 24ifnp.any(~np.isin(tetints,(3,5,6,9,10,12,17,18,24))):warnings.warn(f'Unaccounted for wedge-to-tet case(s) in CleanupDegenerateElements: {str(np.unique(tetints[~np.isin(tetints,(3,5,6,9,10,12,17,18,24))])):s}. This is a bug, please report.')pyrints=np.sum((uConn[wedge2pyr,:6]==-1)*2**np.arange(0,6)[::-1],axis=1)# Pyramids need to be reordered (note case 32 where node 0 is removed never occurs since unique always keeps the first occurence of a duplicate)uConn[wedge2pyr[pyrints==1],:6]=uConn[wedge2pyr[pyrints==1]][:,[0,3,4,1,2,5]]# node 5 removeduConn[wedge2pyr[pyrints==2],:6]=uConn[wedge2pyr[pyrints==2]][:,[0,2,5,3,1,4]]# node 4 removeduConn[wedge2pyr[pyrints==4],:6]=uConn[wedge2pyr[pyrints==4]][:,[1,4,5,2,0,3]]# node 3 removeduConn[wedge2pyr[pyrints==8],:6]=uConn[wedge2pyr[pyrints==8]][:,[0,3,4,1,5,2]]# node 2 removeduConn[wedge2pyr[pyrints==16],:6]=uConn[wedge2pyr[pyrints==16]][:,[0,2,5,3,4,1]]# node 1 removedifnp.any(~np.isin(pyrints,[1,2,4,8,16])):warnings.warn(f'Unaccounted for wedge-to-pyr case(s) in CleanupDegenerateElements: {str(np.unique(pyrints[~np.isin(pyrints,[1,2,4,8,16])])):s}. This is a bug, please report.')ifPadConn.shape[1]>=8:# Special attention need for degenerate hex elementshex_rows=np.sum(PadConn!=-1,axis=1)==8hex2tet=np.where((hex_rows)&(np.sum(uConn!=-1,axis=1)==4))[0]hex2pyr=np.where((hex_rows)&(np.sum(uConn!=-1,axis=1)==5))[0]hex2wdg=np.where((hex_rows)&(np.sum(uConn!=-1,axis=1)==6))[0]tetints=np.sum((uConn[hex2tet,:8]==-1)*2**np.arange(0,8)[::-1],axis=1)pyrints=np.sum((uConn[hex2pyr,:8]==-1)*2**np.arange(0,8)[::-1],axis=1)wdgints=np.sum((uConn[hex2wdg,:8]==-1)*2**np.arange(0,8)[::-1],axis=1)# Wedge cases: TODO: Not all cases accounted for# Case 3 : Face 1 vertical collapse (2==6, 3==7)uConn[hex2wdg[wdgints==3],:8]=uConn[hex2wdg[wdgints==3]][:,[0,3,4,1,2,5,6,7]]# Case 12 : Face 1 vertical collapse (0==4, 1==5)uConn[hex2wdg[wdgints==12],:8]=uConn[hex2wdg[wdgints==12]][:,[0,3,7,1,2,6,4,5]]# Pyramid cases:# Case 112 : Face 0 collapse (0==1==2==3)uConn[hex2pyr[pyrints==112],:8]=uConn[hex2pyr[pyrints==112]][:,[7,6,5,4,0,1,2,3]]# Case 76 : Face 1 collapse (0==1==4==5)uConn[hex2pyr[pyrints==76],:8]=uConn[hex2pyr[pyrints==76]][:,[2,6,7,3,0,1,4,5]]# Case 38 : Face 2 collapse (1==2==5==6)uConn[hex2pyr[pyrints==38],:8]=uConn[hex2pyr[pyrints==38]][:,[0,3,7,4,1,2,5,6]]# Case 25 : Face 4 collapse (0==3==4==7)uConn[hex2pyr[pyrints==25],:8]=uConn[hex2pyr[pyrints==25]][:,[1,5,6,2,0,3,4,7]]# Case 19 : Face 3 collapse (2==3==6==7)uConn[hex2pyr[pyrints==19],:8]=uConn[hex2pyr[pyrints==19]][:,[0,4,5,1,2,3,6,7]]# Case 7 : Face 5 collapse (4==5==6==7)uConn[hex2pyr[pyrints==7],:8]=uConn[hex2pyr[pyrints==7]][:,[0,1,2,3,4,5,6,7]]ifnp.any((wdgints!=3)&(wdgints!=12)):warnings.warn(f'Unaccounted for hex-to-wedge case(s) in CleanupDegenerateElements. This is a bug, please report.')ifnp.any((pyrints!=7)&(pyrints!=19)&(pyrints!=25)&(pyrints!=38)&(pyrints!=76)&(pyrints!=112)):warnings.warn(f'Unaccounted for hex-to-pyr case(s) in CleanupDegenerateElements. This is a bug, please report.')# if len(tetints) > 0:# warnings.warn(f'Unaccounted for hex-to-tet case(s) in CleanupDegenerateElements. This is a bug, please report.')uConn=uConn[~to_delete]NewConn=ExtractRagged(uConn)returnNewConn,np.where(~to_delete)[0]ifType.lower()=='auto':Type=identify_type(NodeCoords,NodeConn)ifType.lower()=='surf':NewConn,idx=rowunique(NodeConn,3)elifType.lower()=='vol':NewConn,idx=rowunique(NodeConn,4)else:raiseValueError(f'Type must be "surf" or "vol", not {Type:s}.')ifreturn_idx:returnNodeCoords,NewConn,idxreturnNodeCoords,NewConn
[docs]defMirrorMesh(NodeCoords,NodeConn,x=None,y=None,z=None):""" Creates a mirrored copy of a mesh by mirroring about the planes defined by X=x, Y=y, and Z=z Parameters ---------- NodeCoords : list Nodal Coordinates. NodeConn : list Nodal Connectivity. x : float, optional YZ plane at X = x. The default is None. y : float, optional XZ plane at Y = y. The default is None. z : float, optional XY plane at Z = z. The default is None. Returns ------- MirroredCoords : list Mirrored Nodal Coordinates. MirroredConn : list Nodal Connectivity of Mirrored Elements. """ifxisNoneandyisNoneandzisNone:warnings.warn('No mirror plane was specified, specify at least one of x, y, or z.')MirroredCoords=np.copy(NodeCoords)ifx!=None:MirroredCoords[:,0]=-(MirroredCoords[:,0]-x)+xify!=None:MirroredCoords[:,1]=-(MirroredCoords[:,1]-y)+yifz!=None:MirroredCoords[:,2]=-(MirroredCoords[:,2]-z)+zreturnMirroredCoords,NodeConn
[docs]defMergeMesh(NodeCoords1,NodeConn1,NodeCoords2,NodeConn2,NodeVals1=[],NodeVals2=[],cleanup=True):""" Merge two meshes together Parameters ---------- NodeCoords1 : list List of nodal coordinates for mesh 1. NodeConn1 : list List of nodal connectivities for mesh 1. NodeCoords2 : list List of nodal coordinates for mesh 2. NodeConn2 : list List of nodal connectivities for mesh 2. NodeVals1 : list, optional List of node data associated with mesh 1, by default [] NodeVals2 : list, optional List of node data associated with mesh 2, by default [] cleanup : bool, optional If true, duplicate nodes will be deleted and renumbered accordingly, by default True. Returns ------- MergedCoords : list List of nodal coordinates of the merged mesh. Nodes from mesh 1 appear first, followed by those of mesh 2. MergedConn : list List of nodal connectivities of the merged mesh. MergedVals : list, optional If provided, merged list of NodeVals. """iftype(NodeCoords1)==np.ndarray:NodeCoords1=NodeCoords1.tolist()iftype(NodeConn1)==np.ndarray:NodeConn1=NodeConn1.tolist()iftype(NodeCoords2)==np.ndarray:NodeCoords2=NodeCoords2.tolist()iftype(NodeConn2)==np.ndarray:NodeConn2=NodeConn2.tolist()iftype(NodeVals1)==np.ndarray:NodeVals1=NodeVals1.tolist()iftype(NodeVals2)==np.ndarray:NodeVals2=NodeVals2.tolist()MergeCoords=NodeCoords1+NodeCoords2MergeConn=NodeConn1+[[node+len(NodeCoords1)fornodeinelem]foreleminNodeConn2]iflen(NodeVals1)>0:assertlen(NodeVals1)==len(NodeCoords1),'NodeVals lists must contain the number of entries as nodes.'assertlen(NodeVals2)==len(NodeCoords2),'NodeVals lists must contain the number of entries as nodes.'MergeVals=[[]foriinrange(len(NodeVals1))]foriinrange(len(NodeVals1)):MergeVals[i]=NodeVals1[i]+NodeVals2[i]ifcleanup:MergeCoords,MergeConn,inv=DeleteDuplicateNodes(MergeCoords,MergeConn,return_inv=True)foriinrange(len(MergeVals)):MergeVals[i]=[MergeVals[i][j]forjinidx]returnMergeCoords,MergeConn,MergeValselifcleanup:MergeCoords,MergeConn=DeleteDuplicateNodes(MergeCoords,MergeConn)returnMergeCoords,MergeConn
[docs]defDetectFeatures(NodeCoords,SurfConn,angle=25):""" Classifies nodes as edges or corners if the angle between adjacent surface elements is less than or equal to `angle`. Parameters ---------- NodeCoords : list List of nodal coordinates. SurfConn : list List of nodal connectivities of a surface mesh. angle : float, optional Dihedral angle threshold (in degrees) used to determine whether an edge exists between two adjacent faces, by default 25. Returns ------- edges : list list of nodes identified to lie on an edge of the geometry. corners : list list of nodes identified to lie on a corner of the geometry. Examples -------- .. plot:: background = primitives.Grid([0,1,0,1,0,1], .02, ElemType='tet') S = implicit.TetMesh(implicit.thickenf(implicit.gyroid,1), [0,1,0,1,0,1], .03, background=background) edges, corners = utils.DetectFeatures(S.NodeCoords, S.SurfConn) features = np.zeros(S.NNode) features[edges] = 1 features[corners] = 2 S.plot(scalars=features, color='coolwarm', bgcolor='w') """ElemNormals=np.asarray(CalcFaceNormal(NodeCoords,SurfConn))Edges,EdgeConn,EdgeElem=converter.solid2edges(NodeCoords,SurfConn,return_EdgeConn=True,return_EdgeElem=True)UEdges,UIdx,UInv=converter.edges2unique(Edges,return_idx=True,return_inv=True)UEdgeElem=np.asarray(EdgeElem)[UIdx]UEdgeConn=UInv[PadRagged(EdgeConn)]EECidx=(UEdgeElem[UEdgeConn]==np.repeat(np.arange(len(EdgeConn))[:,None],UEdgeConn.shape[1],axis=1)).astype(int)EdgeElemConn=-1*(np.ones((len(UEdges),2),dtype=int))EdgeElemConn[UEdgeConn,EECidx]=np.repeat(np.arange(len(EdgeConn))[:,None],UEdgeConn.shape[1],axis=1)ConnectedNormals=ElemNormals[EdgeElemConn]angles=quality.dihedralAngles(ConnectedNormals[:,0],ConnectedNormals[:,1],Abs=False)FeatureEdges=np.where(angles>angle*np.pi/180)[0]FeatureNodes=[nforedgeinFeatureEdgesforninUEdges[edge]]unq,counts=np.unique(FeatureNodes,return_counts=True)corners=unq[counts>2].tolist()edges=unq[counts<=2].tolist()returnedges,corners
[docs]defmakePyramidLayer(VoxelCoords,VoxelConn,PyramidHeight=None):""" Generate a set of pyramid elements that cover the surface of the voxel mesh. To merge the pyramid layer with the voxel mesh, use :func:`MergeMesh`. Parameters ---------- VoxelCoords : list Contains coordinates for each node in a voxel mesh. Ex. [[x0,y0,z0],...]. VoxelConn : List Nodal connectivity list. The voxel mesh is assumed to consist of a set of uniform cubic hexahedral elements. PyramidHeight : float (or None), optional Height of the pyramids. The default is None. If no height as assigned, it will default to 1/2 of the voxel size Returns ------- PyramidCoords : list List of nodal coordinates for the pyramid elements. PyramidConn : list List of nodal connectivities for the pyramid elements. """ifPyramidHeight==None:PyramidHeight=abs(VoxelCoords[VoxelConn[0][0]][0]-VoxelCoords[VoxelConn[0][1]][0])/2SurfConn=converter.solid2surface(VoxelCoords,VoxelConn)SurfCoords,SurfConn,_=RemoveNodes(VoxelCoords,SurfConn)FaceNormals=CalcFaceNormal(SurfCoords,SurfConn)ArrayCoords=np.array(SurfCoords)PyramidConn=[[]foriinrange(len(SurfConn))]PyramidCoords=SurfCoordsfori,faceinenumerate(SurfConn):nodes=ArrayCoords[face]centroid=np.mean(nodes,axis=0)tipCoord=centroid+PyramidHeight*np.array(FaceNormals[i])PyramidConn[i]=face+[len(PyramidCoords)]PyramidCoords.append(tipCoord.tolist())returnPyramidCoords,PyramidConn
[docs]defErodeVoxel(NodeCoords,NodeConn,nLayers=1):""" Removes the specified number of layers from a hexahedral mesh Parameters ---------- NodeCoords : list of lists Contains coordinates for each node in a voxel mesh. Ex. [[x1,y1,z1],...]. The mesh is assumed to consist of only hexahedral elements. NodeConn : List of lists Nodal connectivity list. nLayers : int, optional Number of layers to peel. The default is 1. Returns ------- PeeledCoords : List Node coordinates for each node in the peeled mesh. PeeledConn : list Nodal connectivity for each element in the peeled mesh. PeelCoords : list Node coordinates for each node in the layers of the mesh that have been removed. PeelConn : list Nodal connectivity for each element in the layers of the mesh that have been removed. """NewCoords=copy.copy(NodeCoords)NewConn=copy.copy(NodeConn)PeelConn=[]foriinrange(nLayers):HexSurfConn=converter.solid2surface(NewCoords,NewConn)SurfNodes=np.unique(HexSurfConn)SurfNodeSet=set(SurfNodes)PeelConn+=[NewConn[i]foriinrange(len(NewConn))if(set(NewConn[i])&SurfNodeSet)]NewConn=[NewConn[i]foriinrange(len(NewConn))ifnot(set(NewConn[i])&SurfNodeSet)]PeelCoords,PeelConn,_=RemoveNodes(NewCoords,PeelConn)PeeledCoords,PeeledConn,_=RemoveNodes(NewCoords,NewConn)returnPeeledCoords,PeeledConn,PeelCoords,PeelConn
[docs]defDilateVoxel(VoxelCoords,VoxelConn):""" For a given voxel mesh, will generate a layer of voxels that wrap around the current voxel mesh. NOTE: This has the potential to create overlapping voxels Parameters ---------- VoxelCoords : list of lists Contains coordinates for each node in a voxel mesh. Ex. [[x1,y1,z1],...]. The voxel mesh is assumed to consist of a set of uniform cubic hexahedral elements. VoxelConn : List of lists Nodal connectivity list. Returns ------- LayerCoords : list New node coordinates. LayerConn : TYPE New node connectivity. """VoxelSize=abs(VoxelCoords[VoxelConn[0][0]][0]-VoxelCoords[VoxelConn[0][1]][0])SurfConn=converter.solid2surface(VoxelCoords,VoxelConn)SurfCoords,SurfConn,_=RemoveNodes(VoxelCoords,SurfConn)FaceNormals=CalcFaceNormal(SurfCoords,SurfConn)ArrayCoords=np.array(SurfCoords)LayerConn=[[]foriinrange(len(SurfConn))]LayerCoords=SurfCoordsfori,faceinenumerate(SurfConn):nodes=ArrayCoords[face]coord0=nodes[0]+VoxelSize*np.array(FaceNormals[i])coord1=nodes[1]+VoxelSize*np.array(FaceNormals[i])coord2=nodes[2]+VoxelSize*np.array(FaceNormals[i])coord3=nodes[3]+VoxelSize*np.array(FaceNormals[i])LayerConn[i]=face+[len(LayerCoords),len(LayerCoords)+1,len(LayerCoords)+2,len(LayerCoords)+3]LayerCoords.append(coord0.tolist())LayerCoords.append(coord1.tolist())LayerCoords.append(coord2.tolist())LayerCoords.append(coord3.tolist())returnLayerCoords,LayerConn
[docs]defTriSurfVol(NodeCoords,SurfConn):""" Calculates the volume contained within a surface mesh. Based on 'Efficient feature extraction for 2D/3D objects in mesh representation.' - Zhang, C. and Chen, T., 2001 Parameters ---------- NodeCoords : list of lists Contains coordinates for each node. Ex. [[x1,y1,z1],...]. SurfConn : List of lists Nodal connectivity list for a triangular surface mesh. Returns ------- V : float Volume contained within the surface mesh. """defTriSignedVolume(nodes):return1/6*(-nodes[2][0]*nodes[1][1]*nodes[0][2]+nodes[1][0]*nodes[2][1]*nodes[0][2]+nodes[2][0]*nodes[0][1]*nodes[1][2]-nodes[0][0]*nodes[2][1]*nodes[1][2]-nodes[1][0]*nodes[0][1]*nodes[2][2]+nodes[0][0]*nodes[1][1]*nodes[2][2])V=sum([TriSignedVolume([NodeCoords[node]fornodeinelem])foreleminSurfConn])returnV
[docs]defTetMeshVol(NodeCoords,NodeConn):""" Calculates the volume contained within a tetrahedral mesh Parameters ---------- NodeCoords : list of lists Contains coordinates for each node. Ex. [[x1,y1,z1],...]. NodeConn : List of lists Nodal connectivity list for a tetrahedral mesh. Returns ------- V : float Volume contained within the tetrahedral mesh. """vs=quality.Volume(NodeCoords,NodeConn)V=np.sum(vs)returnV
[docs]defMVBB(Points,return_matrix=False):""" Calculate the minimum volume bounding box of the set of points Parameters ---------- Points : array_like nx3 point coordinates. return_matrix : bool, optional option to return the rotation matrix that aligns the input Points with the local coordinate system of the MVBB, by default False. Returns ------- mvbb : np.ndarray Coordinates of the corners of the MVBB mat : np.ndarray, optional Rotation matrix that aligns the input Points with the local coordinate system of the MVBB """hull=scipy.spatial.ConvexHull(Points)hull_points,hull_facets,_=RemoveNodes(Points,hull.simplices)hull_points=np.asarray(hull_points)# Calculate rotation matrices to align each hull facet with [0,0,-1] (so that it's rotated to the minimal z plane)normals=np.asarray(CalcFaceNormal(hull_points,hull_facets))rot_axes=np.cross(normals,[0,0,-1])rot_axes[np.all(normals==[0,0,-1],axis=1)]=[0,0,-1]rot_axes[np.all(normals==[0,0,1],axis=1)]=[1,0,0]rot_axes=rot_axes/np.linalg.norm(rot_axes,axis=1)[:,None]thetas=np.arccos(np.sum(normals*[0,0,-1],axis=1))thetas[np.all(normals==[0,0,1],axis=1)]=np.piouter_prod=rot_axes[:,np.newaxis,:]*rot_axes[:,:,np.newaxis]cross_prod_matrices=np.zeros((len(hull_facets),3,3))cross_prod_matrices[:,0,1]=-rot_axes[:,2]cross_prod_matrices[:,1,0]=rot_axes[:,2]cross_prod_matrices[:,0,2]=rot_axes[:,1]cross_prod_matrices[:,2,0]=-rot_axes[:,1]cross_prod_matrices[:,1,2]=-rot_axes[:,0]cross_prod_matrices[:,2,1]=rot_axes[:,0]rot_matrices=np.cos(thetas)[:,None,None]*np.repeat([np.eye(3)],len(hull_facets),axis=0)+np.sin(thetas)[:,None,None]*cross_prod_matrices+(1-np.cos(thetas))[:,None,None]*outer_prod# NOTE: might be able to reduce memory usage by not explicitly obtaining the rotation matrices (see https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle)# For each possible rotation, rotate all of the pointsrotated_points=rot_matrices@hull_points.T[None,:,:]# Get the local coordinate system axis-aligned bounding boxes for each rotationmins=np.min(rotated_points,axis=2)maxs=np.max(rotated_points,axis=2)# Calculate the box volumes and find the smallestside_lengths=maxs-minsvolumes=np.prod(side_lengths,axis=1)min_idx=np.argmin(volumes)# Get local coordinates of the MVBBrotated_bb=np.array([[mins[min_idx,0],mins[min_idx,1],mins[min_idx,2]],[maxs[min_idx,0],mins[min_idx,1],mins[min_idx,2]],[maxs[min_idx,0],maxs[min_idx,1],mins[min_idx,2]],[mins[min_idx,0],maxs[min_idx,1],mins[min_idx,2]],[mins[min_idx,0],mins[min_idx,1],maxs[min_idx,2]],[maxs[min_idx,0],mins[min_idx,1],maxs[min_idx,2]],[maxs[min_idx,0],maxs[min_idx,1],maxs[min_idx,2]],[mins[min_idx,0],maxs[min_idx,1],maxs[min_idx,2]],])# Return the MVBB to the original coordinate systemmat=rot_matrices[min_idx]mvbb=(np.linalg.inv(mat)@rotated_bb.T).Tifreturn_matrix:returnmvbb,matreturnmvbb
[docs]defAABB(Points):""" Calculate the axis-aligned bounding box of a set of points Parameters ---------- Points : array_like nx3 point coordinates. Returns ------- aabb : np.ndarray Coordinates of the corners of the AABB """mins=np.min(Points,axis=0)maxs=np.max(Points,axis=0)aabb=np.array([[mins[0],mins[1],mins[2]],[maxs[0],mins[1],mins[2]],[maxs[0],maxs[1],mins[2]],[mins[0],maxs[1],mins[2]],[mins[0],mins[1],maxs[2]],[maxs[0],mins[1],maxs[2]],[maxs[0],maxs[1],maxs[2]],[mins[0],maxs[1],maxs[2]],])returnaabb
[docs]defSortRaggedByLength(In,return_idx=False,return_inv=False,return_separators=False):""" Sorted a ragged list of lists by the length of each sublist Parameters ---------- In : list List of lists to be sorted return_idx : bool, optional Returns the indices of each row of In in the order that they're sorted into Out, by default False. return_inv : bool, optional Returns the indices that reverse the sorting operation, by default False. return_separators : bool, optional Returns the indices that separate sections of the list by length. Determining these separators requires a small amount of additional work, by default False. Returns ------- Out : list List of lists sorted by row length idx : np.ndarray, optional Indices used to reorder In to Out. These are the indices of each row of In in the order that they're sorted into Out. Returned if return_idx is True. inv : np.ndarray, optional Indices to recover the original List (in) from the output (Out). Return if return_inv us True. separators : np.ndarray, optional Indices of Out that separate sections of the list by length. These separators will always include 0 as the first separator and len(Out) as the last separator. With the exception of the last separator, each separator is the start of a new section and are set such that the sublists of equal-length lists can be accessed by slices with two adjacent separators. Examples -------- >>> In = [[0,1], [2, 3, 4, 5], [6, 7], [8, 9, 10]] >>> Out, idx, inv = utils.SortRaggedByLength(In, return_idx=True, return_inv=True) >>> Out >>> [In[i] for i in idx] == Out >>> [Out[i] for i in inv] == In """lengths=np.array(list(map(len,In)))idx=lengths.argsort()Out=[In[i]foriinidx]ifreturn_separators:separators=np.concatenate([[0],np.where(np.diff(lengths[idx])!=0)[0]+1,[len(Out)]])else:separators=Noneifreturn_inv:inv=np.zeros(len(idx),dtype=int)inv[idx]=np.arange(len(idx),dtype=int)else:inv=Nonereturns=(True,return_idx,return_inv,return_separators)ifsum(returns)>1:returntuple(outputfori,outputinenumerate((Out,idx,inv,separators))ifreturns[i])returnOut
[docs]defSplitRaggedByLength(In,return_idx=False,return_inv=False):""" Split a ragged list of lists into a list of array_like groupings of the original list in which all rows are equal length. The returned list will be the length of the number of unique row lengths of the original list of lists. Parameters ---------- In : list List of lists to be sorted return_idx : bool, optional Returns the indices of each row of In in the order that they're sorted into Out, by default False. Returns ------- Out : list List of array_like groupings of the original list in which all rows are equal length. idx : np.ndarray, optional Indices used to reorder In to Out. These are the indices of each row of In in the order that they're sorted into Out. Returned if return_idx is True. inv : np.ndarray, optional Indices to recover the original List (in) from the output (Out). Return if return_inv us True. Examples -------- >>> In = [[0,1], [2, 3, 4, 5], [6, 7], [8, 9, 10]] >>> Out = utils.SplitRaggedByLength(In) >>> Out """out=SortRaggedByLength(In,return_idx=return_idx,return_inv=return_inv,return_separators=True)out=list(out)[::-1]In_sorted,In_idx,In_inv,separators=[out.pop()ifbelseNoneforbin(True,return_idx,return_inv,True)]Out=[In_sorted[separators[i]:separators[i+1]]foriinrange(len(separators)-1)]ifreturn_idx:idx=[In_idx[separators[i]:separators[i+1]]foriinrange(len(separators)-1)]else:idx=Noneifreturn_inv:inv=[In_inv[separators[i]:separators[i+1]]foriinrange(len(separators)-1)]else:inv=Nonereturns=(True,return_idx,return_inv)ifsum(returns)>1:returntuple(outputfori,outputinenumerate((Out,idx,inv))ifreturns[i])returnOut
[docs]defPadRagged(In,fillval=-1):""" Pads a 2d list of lists with variable length into a rectangular numpy array with specified fill value. Parameters ---------- In : list Input list of lists to be padded. fillval : int (or other), optional Value used to pad the ragged array, by default -1 Returns ------- Out : np.ndarray Padded array. """# Out = np.array(list(itertools.zip_longest(*In,fillvalue=fillval))).TmaxL=max(len(row)forrowinIn)Out=np.full((len(In),maxL),fillval)fori,rowinenumerate(In):Out[i,:len(row)]=rowreturnOut
[docs]defExtractRagged(In,delval=-1,dtype=None):""" Convert a padded numpy array to a ragged list of list by removing entries that match the specified value. Parameters ---------- In : np.ndarray Input array delval : int, optional Value to remove from the input array, by default -1 dtype : type, optional Data type to cast the array to, by default the data type is unchanged. Returns ------- Out : list Output list of lists with the specified value removed. """ifdtype:iftype(In)islist:In=np.array(In)In=In.astype(dtype)delval=np.array([delval]).astype(dtype)[0]ifnp.isnan(delval):delval=np.nanmax(In)+1In=np.copy(In)In[np.isnan(In)]=delvalwhere=In!=delvalifnotnp.all(where):iflen(In.shape)==2:Out=np.split(In[where],np.cumsum(np.sum(where,axis=1)))[:-1]eliflen(In.shape)==3:Out=[[[xforxinyifx!=delval]foryinzifall([x!=delvalforxiny])]forzinIn]else:raiseException('Currently only supported for 2- or 3D matrices')else:Out=In.tolist()returnOut
[docs]defidentify_type(NodeCoords,NodeConn):""" Classify the mesh as either a surface or volume. A mesh is classified as a volume mesh (``vol``) if any elements are unambiguous volume elements - pyramid (5 nodes), wedge (6), hexahedron (8), or if any of a random sample of 10 elements (or all elements if NElem < 10) has a volume less than machine precision (``np.finfo(float).eps``). Alternatively, a surface mesh (``surf``) is identified if any of the elements is a triangle (3 nodes). In the case of a mesh containing both triangular and volume elements, the mesh will be classified arbitrarily by whichever appears first in NodeConn. A ``line`` mesh is identified if any line (2 node) elements are present. This approach has a chance of mislabeling the mesh in some rare or non-standard scenarios, but attempts to be as efficient as possible to minimize overhead when creating a mesh. Potentially ambiguous meshes are: - Meshes containing both triangular elements and volume elements (this should generally be avoided as most functions aren't set up to handle that case). - Tetrahedral meshes with many degenerate elements with abs(vol) < machine precision. - Quadrilateral meshes with non-planar elements. In such cases, Type should be specified explicitly when creating the mesh object. Returns ------- Type : str Mesh type, either 'line', 'surf', 'vol', or 'empty'. """# Check if the mesh is emptyNNode=len(NodeCoords)NElem=len(NodeConn)ifNNode==0orNElem==0:Type='empty'returnType# Check node dimensions, if it's 2D, then it must be a surface or lineiflen(NodeCoords[0])==2:iflen(NodeConn[0])==2:Type='line'else:Type='surf'returnType# Check element lengthslengths=(len(e)foreinNodeConn)vol_elem_set={5,6,8}# Set of unambiguous volume element lengthsforlinlengths:# NOTE: any mesh containing triangle and volume elements will be # arbitrarily classified by whichever comes first.ifl==2:# The presence of any edge elements triggers 'line'Type='line'returnTypeelifl==3:# The presence of any triangular elements triggers 'surf'Type='surf'returnTypeeliflinvol_elem_set:# The presence of any unambiguous volume elements triggers 'vol'Type='vol'returnType# If still ambiguous, check volumes of a random selection of elementsifNElem>10:TempConn=[NodeConn[i]foriinnp.random.randint(NElem,size=(10))]else:TempConn=NodeConnv=quality.Volume(NodeCoords,TempConn,verbose=False)ifnp.max(np.abs(v))>np.finfo(float).eps:# If any of the sampled element volumes exceed machine precision,# mesh is assumed to be a volume mesh, otherwise a surfaceType='vol'else:Type='surf'returnType
