parameters: (v, t)
return_statement: return A, B
docstring: computeAB(v,t) computes the two sparse symmetric matrices representing
             the Laplace Beltrami Operator for a given triangle mesh using
             the linear finite element method (assuming a closed mesh or 
             the Neumann boundary condition).
  
      Inputs:   v - vertices : list of lists of 3 floats
                t - triangles: list of lists of 3 int of indices (>=0) into v array
  
      Outputs:  A - sparse sym. (n x n) positive semi definite numpy matrix 
                B - sparse sym. (n x n) positive definite numpy matrix (inner product)
  
      Can be used to solve sparse generalized Eigenvalue problem: A x = lambda B x
      or to solve Poisson equation: A x = B f (where f is function on mesh vertices)
      or to solve Laplace equation: A x = 0
      or to model the operator's action on a vector x:   y = B\(Ax)
docstring_summary: computeAB(v,t) computes the two sparse symmetric matrices representing
             the Laplace Beltrami Operator for a given triangle mesh using
             the linear finite element method (assuming a closed mesh or 
             the Neumann boundary condition).
docstring_tokens: computeAB ( v t ) computes the two sparse symmetric matrices representing the Laplace Beltrami Operator for a given triangle mesh using the linear finite element method ( assuming a closed mesh or the Neumann boundary condition ) .
function: def computeAB(v, t):
      """
      computeAB(v,t) computes the two sparse symmetric matrices representing
             the Laplace Beltrami Operator for a given triangle mesh using
             the linear finite element method (assuming a closed mesh or 
             the Neumann boundary condition).
  
      Inputs:   v - vertices : list of lists of 3 floats
                t - triangles: list of lists of 3 int of indices (>=0) into v array
  
      Outputs:  A - sparse sym. (n x n) positive semi definite numpy matrix 
                B - sparse sym. (n x n) positive definite numpy matrix (inner product)
  
      Can be used to solve sparse generalized Eigenvalue problem: A x = lambda B x
      or to solve Poisson equation: A x = B f (where f is function on mesh vertices)
      or to solve Laplace equation: A x = 0
      or to model the operator's action on a vector x:   y = B\(Ax) 
      """  
      v = np.array(v)
      t = np.array(t)
      tnum = t.shape[0]
  
  
      # Linear local matrices on unit triangle:
      tB = (np.ones((3,3)) + np.eye(3)) / 24.0
  
      tA00 = np.array([[ 0.5,-0.5, 0.0],
                       [-0.5, 0.5, 0.0],
                       [ 0.0, 0.0, 0.0]])
  
      tA11 = np.array([[ 0.5, 0.0,-0.5],
                       [ 0.0, 0.0, 0.0],
                       [-0.5, 0.0, 0.5]])
  
      tA0110 = np.array([[ 1.0,-0.5,-0.5],
                         [-0.5, 0.0, 0.5],
                         [-0.5, 0.5, 0.0]])
  
      # Compute vertex coordinates and a difference vector for each triangle:
      v1 = v[t[:, 0], :]
      v2 = v[t[:, 1], :]
      v3 = v[t[:, 2], :]
      v2mv1 = v2 - v1
      v3mv1 = v3 - v1
      
      # Compute length^2 of v3mv1 for each triangle:
      a0 = np.sum(v3mv1 * v3mv1, axis=1)
  
      # Compute length^2 of v2mv1 for each triangle:
      a1 = np.sum(v2mv1 * v2mv1, axis=1)
  
      # Compute dot product (v2mv1*v3mv1) for each triangle:
      a0110 = np.sum(v2mv1 * v3mv1, axis=1)
  
      # Compute cross product and 2*vol for each triangle:
      cr  = np.cross(v2mv1,v3mv1)
      vol = np.sqrt(np.sum(cr*cr, axis=1))
      # zero vol will cause division by zero below, so set to small value:
      vol_mean = 0.001*np.mean(vol)
      vol = [vol_mean if x == 0 else x for x in vol]
  
      # Construct all local A and B matrices (guess: for each triangle):
      localB = np.array([x*tB for x in vol])
      localA = np.array([(1.0/xv) * (xa0*tA00 + xa1*tA11 - xa0110*tA0110) for xv, xa0, xa1, xa0110 in zip(vol,a0,a1,a0110)])
  
      # Construct row and col indices.
      I = np.array([np.tile(x, (3,1)) for x in t])
      J = np.array([np.transpose(np.tile(x, (3,1))) for x in t])
  
      # Flatten arrays and swap I and J:
      I = I.flatten()
      J = J.flatten()
      localA = localA.flatten()
      localB = localB.flatten()
  
      # Construct sparse matrix:
      A = sparse.csr_matrix((localA, (I, J)))
      B = sparse.csr_matrix((localB, (I, J)))
  
      return A, B
function_tokens: def computeAB ( v , t ) : v = np . array ( v ) t = np . array ( t ) tnum = t . shape [ 0 ] # Linear local matrices on unit triangle: tB = ( np . ones ( ( 3 , 3 ) ) + np . eye ( 3 ) ) / 24.0 tA00 = np . array ( [ [ 0.5 , - 0.5 , 0.0 ] , [ - 0.5 , 0.5 , 0.0 ] , [ 0.0 , 0.0 , 0.0 ] ] ) tA11 = np . array ( [ [ 0.5 , 0.0 , - 0.5 ] , [ 0.0 , 0.0 , 0.0 ] , [ - 0.5 , 0.0 , 0.5 ] ] ) tA0110 = np . array ( [ [ 1.0 , - 0.5 , - 0.5 ] , [ - 0.5 , 0.0 , 0.5 ] , [ - 0.5 , 0.5 , 0.0 ] ] ) # Compute vertex coordinates and a difference vector for each triangle: v1 = v [ t [ : , 0 ] , : ] v2 = v [ t [ : , 1 ] , : ] v3 = v [ t [ : , 2 ] , : ] v2mv1 = v2 - v1 v3mv1 = v3 - v1 # Compute length^2 of v3mv1 for each triangle: a0 = np . sum ( v3mv1 * v3mv1 , axis = 1 ) # Compute length^2 of v2mv1 for each triangle: a1 = np . sum ( v2mv1 * v2mv1 , axis = 1 ) # Compute dot product (v2mv1*v3mv1) for each triangle: a0110 = np . sum ( v2mv1 * v3mv1 , axis = 1 ) # Compute cross product and 2*vol for each triangle: cr = np . cross ( v2mv1 , v3mv1 ) vol = np . sqrt ( np . sum ( cr * cr , axis = 1 ) ) # zero vol will cause division by zero below, so set to small value: vol_mean = 0.001 * np . mean ( vol ) vol = [ vol_mean if x == 0 else x for x in vol ] # Construct all local A and B matrices (guess: for each triangle): localB = np . array ( [ x * tB for x in vol ] ) localA = np . array ( [ ( 1.0 / xv ) * ( xa0 * tA00 + xa1 * tA11 - xa0110 * tA0110 ) for xv , xa0 , xa1 , xa0110 in zip ( vol , a0 , a1 , a0110 ) ] ) # Construct row and col indices. I = np . array ( [ np . tile ( x , ( 3 , 1 ) ) for x in t ] ) J = np . array ( [ np . transpose ( np . tile ( x , ( 3 , 1 ) ) ) for x in t ] ) # Flatten arrays and swap I and J: I = I . flatten ( ) J = J . flatten ( ) localA = localA . flatten ( ) localB = localB . flatten ( ) # Construct sparse matrix: A = sparse . csr_matrix ( ( localA , ( I , J ) ) ) B = sparse . csr_matrix ( ( localB , ( I , J ) ) ) return A , B