<s> import numpy as np <EOL> from scipy import sparse as sp <EOL> from matutils import mkvc , ndgrid , sub2ind , sdiag <EOL> def volTetra ( xyz , A , B , C , D ) : <EOL> """<STR_LIT>""" <EOL> AD = xyz [ A , : ] - xyz [ D , : ] <EOL> BD = xyz [ B , : ] - xyz [ D , : ] <EOL> CD = xyz [ C , : ] - xyz [ D , : ] <EOL> V = ( BD [ : , <NUM_LIT:0> ] * CD [ : , <NUM_LIT:1> ] - BD [ : , <NUM_LIT:1> ] * CD [ : , <NUM_LIT:0> ] ) * AD [ : , <NUM_LIT:2> ] - ( BD [ : , <NUM_LIT:0> ] * CD [ : , <NUM_LIT:2> ] - BD [ : , <NUM_LIT:2> ] * CD [ : , <NUM_LIT:0> ] ) * AD [ : , <NUM_LIT:1> ] + ( BD [ : , <NUM_LIT:1> ] * CD [ : , <NUM_LIT:2> ] - BD [ : , <NUM_LIT:2> ] * CD [ : , <NUM_LIT:1> ] ) * AD [ : , <NUM_LIT:0> ] <EOL> return V / <NUM_LIT:6> <EOL> def indexCube ( nodes , gridSize , n = None ) : <EOL> """<STR_LIT>""" <EOL> assert type ( nodes ) == str , "<STR_LIT>" <EOL> assert isinstance ( gridSize , np . ndarray ) , "<STR_LIT>" <EOL> nodes = nodes . upper ( ) <EOL> possibleNodes = '<STR_LIT>' if gridSize . size == <NUM_LIT:2> else '<STR_LIT>' <EOL> for node in nodes : <EOL> assert node in possibleNodes , "<STR_LIT>" % possibleNodes <EOL> dim = gridSize . size <EOL> if n is None : <EOL> n = gridSize - <NUM_LIT:1> <EOL> if dim == <NUM_LIT:2> : <EOL> ij = ndgrid ( np . arange ( n [ <NUM_LIT:0> ] ) , np . arange ( n [ <NUM_LIT:1> ] ) ) <EOL> i , j = ij [ : , <NUM_LIT:0> ] , ij [ : , <NUM_LIT:1> ] <EOL> elif dim == <NUM_LIT:3> : <EOL> ijk = ndgrid ( np . arange ( n [ <NUM_LIT:0> ] ) , np . arange ( n [ <NUM_LIT:1> ] ) , np . arange ( n [ <NUM_LIT:2> ] ) ) <EOL> i , j , k = ijk [ : , <NUM_LIT:0> ] , ijk [ : , <NUM_LIT:1> ] , ijk [ : , <NUM_LIT:2> ] <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> nodeMap = { '<STR_LIT:A>' : [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:0> ] , '<STR_LIT:B>' : [ <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:0> ] , '<STR_LIT:C>' : [ <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:0> ] , '<STR_LIT:D>' : [ <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:0> ] , <EOL> '<STR_LIT:E>' : [ <NUM_LIT:0> , <NUM_LIT:0> , <NUM_LIT:1> ] , '<STR_LIT:F>' : [ <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT:1> ] , '<STR_LIT>' : [ <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:1> ] , '<STR_LIT:H>' : [ <NUM_LIT:1> , <NUM_LIT:0> , <NUM_LIT:1> ] } <EOL> out = ( ) <EOL> for node in nodes : <EOL> shift = nodeMap [ node ] <EOL> if dim == <NUM_LIT:2> : <EOL> out += ( sub2ind ( gridSize , np . c_ [ i + shift [ <NUM_LIT:0> ] , j + shift [ <NUM_LIT:1> ] ] ) . flatten ( ) , ) <EOL> elif dim == <NUM_LIT:3> : <EOL> out += ( sub2ind ( gridSize , np . c_ [ i + shift [ <NUM_LIT:0> ] , j + shift [ <NUM_LIT:1> ] , k + shift [ <NUM_LIT:2> ] ] ) . flatten ( ) , ) <EOL> return out <EOL> def faceInfo ( xyz , A , B , C , D , average = True , normalizeNormals = True ) : <EOL> """<STR_LIT>""" <EOL> assert type ( average ) is bool , '<STR_LIT>' <EOL> assert type ( normalizeNormals ) is bool , '<STR_LIT>' <EOL> AB = xyz [ B , : ] - xyz [ A , : ] <EOL> BC = xyz [ C , : ] - xyz [ B , : ] <EOL> CD = xyz [ D , : ] - xyz [ C , : ] <EOL> DA = xyz [ A , : ] - xyz [ D , : ] <EOL> def cross ( X , Y ) : <EOL> return np . c_ [ X [ : , <NUM_LIT:1> ] * Y [ : , <NUM_LIT:2> ] - X [ : , <NUM_LIT:2> ] * Y [ : , <NUM_LIT:1> ] , <EOL> X [ : , <NUM_LIT:2> ] * Y [ : , <NUM_LIT:0> ] - X [ : , <NUM_LIT:0> ] * Y [ : , <NUM_LIT:2> ] , <EOL> X [ : , <NUM_LIT:0> ] * Y [ : , <NUM_LIT:1> ] - X [ : , <NUM_LIT:1> ] * Y [ : , <NUM_LIT:0> ] ] <EOL> nA = cross ( AB , DA ) <EOL> nB = cross ( BC , AB ) <EOL> nC = cross ( CD , BC ) <EOL> nD = cross ( DA , CD ) <EOL> length = lambda x : np . sqrt ( x [ : , <NUM_LIT:0> ] ** <NUM_LIT:2> + x [ : , <NUM_LIT:1> ] ** <NUM_LIT:2> + x [ : , <NUM_LIT:2> ] ** <NUM_LIT:2> ) <EOL> normalize = lambda x : x / np . kron ( np . ones ( ( <NUM_LIT:1> , x . shape [ <NUM_LIT:1> ] ) ) , mkvc ( length ( x ) , <NUM_LIT:2> ) ) <EOL> if average : <EOL> N = ( nA + nB + nC + nD ) / <NUM_LIT:4> <EOL> N = normalize ( N ) <EOL> else : <EOL> if normalizeNormals : <EOL> N = [ normalize ( nA ) , normalize ( nB ) , normalize ( nC ) , normalize ( nD ) ] <EOL> else : <EOL> N = [ nA , nB , nC , nD ] <EOL> area = ( length ( nA ) + length ( nB ) + length ( nC ) + length ( nD ) ) / <NUM_LIT:4> <EOL> return N , area </s>
