Computes the dual of the mesh.  This is basically an optimized version of buildIndex for the situation where from_cells is just the list of vertices
186926
function dual(cells, vertex_count) {
  if(!vertex_count) {
    return incidence(unique(skeleton(cells, 0)), cells, 0)
  }
  var res = new Array(vertex_count)
  for(var i=0; i<vertex_count; ++i) {
    res[i] = []
  }
  for(var i=0, len=cells.length; i<len; ++i) {
    var c = cells[i]
    for(var j=0, cl=c.length; j<cl; ++j) {
      res[c[j]].push(i)
    }
  }
  return res
}
156700
def get_vert_connectivity(mesh_v, mesh_f):
    """"""

    vpv = sp.csc_matrix((len(mesh_v),len(mesh_v)))

    # for each column in the faces...
    for i in range(3):
        IS = mesh_f[:,i]
        JS = mesh_f[:,(i+1)%3]
        data = np.ones(len(IS))
        ij = np.vstack((row(IS.flatten()), row(JS.flatten())))
        mtx = sp.csc_matrix((data, ij), shape=vpv.shape)
        vpv = vpv + mtx + mtx.T

    return vpv
180462
function (options) {
  options = options || {}
  if (('points' in options) !== ('faces' in options)) {
    throw new Error("polyhedron needs 'points' and 'faces' arrays")
  }
  let vertices = parseOptionAs3DVectorList(options, 'points', [
            [1, 1, 0],
            [1, -1, 0],
            [-1, -1, 0],
            [-1, 1, 0],
            [0, 0, 1]
  ])
        .map(function (pt) {
          return new Vertex3(pt)
        })
  let faces = parseOption(options, 'faces', [
            [0, 1, 4],
            [1, 2, 4],
            [2, 3, 4],
            [3, 0, 4],
            [1, 0, 3],
            [2, 1, 3]
  ])
    // Openscad convention defines inward normals - so we have to invert here
  faces.forEach(function (face) {
    face.reverse()
  })
  let polygons = faces.map(function (face) {
    return new Polygon3(face.map(function (idx) {
      return vertices[idx]
    }))
  })

    // TODO: facecenters as connectors? probably overkill. Maybe centroid
    // the re-tesselation here happens because it's so easy for a user to
    // create parametrized polyhedrons that end up with 1-2 dimensional polygons.
    // These will create infinite loops at CSG.Tree()
  return fromPolygons(polygons).reTesselated()
}
1199253
public static int[] ind2subC(int[] shape, long index, long numIndices) {
        long denom = numIndices;
        int[] ret = new int[shape.length];
        for (int i = 0; i < shape.length; i++) {
            denom /= shape[i];
            if (index / denom >= Integer.MAX_VALUE)
                throw new IllegalArgumentException("Dimension can not be >= Integer.MAX_VALUE");
            ret[i] = (int) (index / denom);
            index %= denom;

        }
        return ret;
    }
967671
private Polygon getCellPolygon() {
        final Coordinate[] summits = new Coordinate[5];
        double x1 = minX + cellI * deltaX;
        double y1 = minY + cellJ * deltaY;
        double x2 = minX + (cellI + 1) * deltaX;
        double y2 = minY + (cellJ + 1) * deltaY;
        summits[0] = new Coordinate(x1, y1);
        summits[1] = new Coordinate(x2, y1);
        summits[2] = new Coordinate(x2, y2);
        summits[3] = new Coordinate(x1, y2);
        summits[4] = new Coordinate(x1, y1);
        final LinearRing g = GF.createLinearRing(summits);
        final Polygon gg = GF.createPolygon(g, null);
        cellI++;
        return gg;
    }
641204
def multigrid(bounds, points_count):
    
    if len(bounds)==1:
        return np.linspace(bounds[0][0], bounds[0][1], points_count).reshape(points_count, 1)
    x_grid_rows = np.meshgrid(*[np.linspace(b[0], b[1], points_count) for b in bounds])
    x_grid_columns = np.vstack([x.flatten(order='F') for x in x_grid_rows]).T
    return x_grid_columns
799403
def connectedVertices(self, index, returnIds=False):
        
        mesh = self.polydata()

        cellIdList = vtk.vtkIdList()
        mesh.GetPointCells(index, cellIdList)

        idxs = []
        for i in range(cellIdList.GetNumberOfIds()):
            pointIdList = vtk.vtkIdList()
            mesh.GetCellPoints(cellIdList.GetId(i), pointIdList)
            for j in range(pointIdList.GetNumberOfIds()):
                idj = pointIdList.GetId(j)
                if idj == index:
                    continue
                if idj in idxs:
                    continue
                idxs.append(idj)

        if returnIds:
            return idxs
        else:
            trgp = []
            for i in idxs:
                p = [0, 0, 0]
                mesh.GetPoints().GetPoint(i, p)
                trgp.append(p)
            return np.array(trgp)
251804
function buildEdgeVertex(index, x, y, tileBody, isGhost) {
	let vertex = {
		x: x,
		y: y,
		index: index,
		body: tileBody,
		isInternal: false,
		isGhost: isGhost,
		contact: null
	};
	
	vertex.contact = {
		vertex: vertex,
		normalImpulse: 0,
		tangentImpulse: 0
	};
	
	return vertex;
}
221843
def get_dual_variables(self):
        
        
        if self.problem:
            return self.problem.recover_dual_variables(self.lam*self.obj_sca,
                                                       self.nu*self.obj_sca,
                                                       self.mu*self.obj_sca,
                                                       self.pi*self.obj_sca)
        else:
            return None,None,None,None
918840
function(model, opt) {

        opt = opt || {};

        var inbound = opt.inbound;
        var outbound = opt.outbound;
        if ((inbound === undefined) && (outbound === undefined)) {
            inbound = outbound = true;
        }

        // the final array of connected link models
        var links = [];
        // a hash table of connected edges of the form: [edgeId] -> true
        // used for quick lookups to check if we already added a link
        var edges = {};

        if (outbound) {
            joint.util.forIn(this.getOutboundEdges(model.id), function(exists, edge) {
                if (!edges[edge]) {
                    links.push(this.getCell(edge));
                    edges[edge] = true;
                }
            }.bind(this));
        }
        if (inbound) {
            joint.util.forIn(this.getInboundEdges(model.id), function(exists, edge) {
                // skip links that were already added
                // (those must be self-loop links)
                // (because they are inbound and outbound edges of the same two elements)
                if (!edges[edge]) {
                    links.push(this.getCell(edge));
                    edges[edge] = true;
                }
            }.bind(this));
        }

        // if `deep` option is `true`, check also all the links that are connected to any of the descendant cells
        if (opt.deep) {

            var embeddedCells = model.getEmbeddedCells({ deep: true });

            // in the first round, we collect all the embedded elements
            var embeddedElements = {};
            embeddedCells.forEach(function(cell) {
                if (cell.isElement()) {
                    embeddedElements[cell.id] = true;
                }
            });

            embeddedCells.forEach(function(cell) {
                if (cell.isLink()) return;
                if (outbound) {
                    joint.util.forIn(this.getOutboundEdges(cell.id), function(exists, edge) {
                        if (!edges[edge]) {
                            var edgeCell = this.getCell(edge);
                            var sourceId = edgeCell.source().id;
                            var targetId = edgeCell.target().id;

                            // if `includeEnclosed` option is falsy, skip enclosed links
                            if (!opt.includeEnclosed
                                && (sourceId && embeddedElements[sourceId])
                                && (targetId && embeddedElements[targetId])) {
                                return;
                            }

                            links.push(this.getCell(edge));
                            edges[edge] = true;
                        }
                    }.bind(this));
                }
                if (inbound) {
                    joint.util.forIn(this.getInboundEdges(cell.id), function(exists, edge) {
                        if (!edges[edge]) {
                            var edgeCell = this.getCell(edge);
                            var sourceId = edgeCell.source().id;
                            var targetId = edgeCell.target().id;

                            // if `includeEnclosed` option is falsy, skip enclosed links
                            if (!opt.includeEnclosed
                                && (sourceId && embeddedElements[sourceId])
                                && (targetId && embeddedElements[targetId])) {
                                return;
                            }

                            links.push(this.getCell(edge));
                            edges[edge] = true;
                        }
                    }.bind(this));
                }
            }, this);
        }

        return links;
    }
1272809
private Vector2i getCellOf(Point3D coordinate) {
        int xCell = (int) ((coordinate.x - indexRegion.minX()) / x_size);
        int yCell = (int) ((coordinate.y - indexRegion.minY()) / y_size);
        return new Vector2i(xCell, yCell);
    }
179430
def decode_coords(ds, gridfile=None):
        
        extra_coords = set(ds.coords)
        for var in six.itervalues(ds.variables):
            if 'mesh' in var.attrs:
                mesh = var.attrs['mesh']
                if mesh not in extra_coords:
                    extra_coords.add(mesh)
                    try:
                        mesh_var = ds.variables[mesh]
                    except KeyError:
                        warn('Could not find mesh variable %s' % mesh)
                        continue
                    if 'node_coordinates' in mesh_var.attrs:
                        extra_coords.update(
                            mesh_var.attrs['node_coordinates'].split())
                    if 'face_node_connectivity' in mesh_var.attrs:
                        extra_coords.add(
                            mesh_var.attrs['face_node_connectivity'])
        if gridfile is not None and not isinstance(gridfile, xr.Dataset):
            gridfile = open_dataset(gridfile)
            ds.update({k: v for k, v in six.iteritems(gridfile.variables)
                       if k in extra_coords})
        if xr_version < (0, 11):
            ds.set_coords(extra_coords.intersection(ds.variables),
                          inplace=True)
        else:
            ds._coord_names.update(extra_coords.intersection(ds.variables))
        return ds
395673
def unpack_column_flat(self, value, components=None, offset=False, computed_type=None):
        
        if components:
            if isinstance(components, str):
                components = [components]
        else:
            components = self._dict.keys()

            # TODO: add this

        # we need to split the flat array by the lengths of each mesh
        N_lower = 0
        N_upper = 0
        offsetN = 0.0
        value_dict = {}
        for comp in components:
            if isinstance(self[comp], Meshes):
                # then we need to recursively extract to the underlying meshes
                # pass
                meshes = self[comp]._dict
            else:
                meshes = {comp: self[comp]}

            for c, mesh in meshes.items():
                if computed_type=='vertices' or (computed_type is None and mesh._compute_at_vertices):
                    N = mesh.Nvertices
                else:
                    N = mesh.Ntriangles
                N_upper += N
                value_dict[c] = value[N_lower:N_upper] - offsetN
                if offset:
                    offsetN += N
                N_lower += N

        return value_dict
896018
function subdivideLoop(cells, positions) {
  var edges       = nskel(cells, 1)
  var e_incidence = top.incidence(edges, cells)
  var dual        = top.dual(cells, positions.length)
  var npositions  = []
  var ncells      = []
  var e_indices   = new Array(edges.length)
  var v_indices   = new Array(positions.length)
  
  for(var i=0; i<e_indices.length; ++i) {
    e_indices[i] = -1
  }
  for(var i=0; i<v_indices.length; ++i) {
    v_indices[i] = -1
  }
  
  for(var f=0; f<cells.length; ++f) {
    var face = cells[f]
    
    for(var d=0; d<3; ++d) {
      var v = face[d]
      var u = face[(d+1)%3]
      e[0] = v
      e[1] = u
      var e_ptr = top.findCell(edges, e)
      var e_idx = e_indices[e_ptr]
      
      if(e_idx >= 0) {
        e_verts[d] = e_idx
      } else {
        //Compute edge-vertex
        var wing = e_incidence[e_ptr]
        var v0 = positions[u]
        var v1 = positions[v]
        var vertex = new Array(3)
        
        if(wing.length === 2) {
          var v2 = positions[opposite(u, v, cells[wing[0]])]
          var v3 = positions[opposite(u, v, cells[wing[1]])]
          
          for(var i=0; i<3; ++i) {
            vertex[i] = (3.0 * (v0[i] + v1[i]) + v2[i] + v3[i]) / 8.0
          }
        } else {
          for(var i=0; i<3; ++i) {
            vertex[i] = 0.5 * (v0[i] + v1[i])
          }
        }
        
        //Store vertex and continue
        e_indices[e_ptr] = e_verts[d] = npositions.length
        npositions.push(vertex)
      }
      
      if(v_indices[v] >= 0) {
        v_verts[d] = v_indices[v]
      } else {
        //Compute vertex-vertex weight
        
        //First, extract vertex neighborhood (slow and stupid here)
        var star = dual[v]
        var nbhd = [v]
        for(var i=0; i<star.length; ++i) {
          var tri = cells[star[i]]
          for(var j=0; j<3; ++j) {
            if(nbhd.indexOf(tri[j]) !== -1) {
              nbhd.push(tri[j])
            }
          }
        }
        
        //Next, compute weights
        var beta = (star.length === 3 ? 3.0/16.0 : 3.0/(8.0*star.length) )
        var center_weight = 1.0 - star.length * beta
        
        //Finally sum up weights
        var pos  = positions[v]
        var vertex = [0.1,0.1,0.1]
        for(var i=0; i<3; ++i) {
          vertex[i] = center_weight * pos[i]
        }
        for(var i=1; i<nbhd.length; ++i) {
          var p = positions[nbhd[i]]
          for(var j=0; j<3; ++j) {
            vertex[j] += beta * p[j]
          }
        }

        //Store result and continue        
        v_verts[d] = v_indices[v] = npositions.length
        npositions.push(vertex)
      }
    }
    
    //Add subdivided faces
    ncells.push([v_verts[0], e_verts[0], e_verts[2]],
                [e_verts[0], e_verts[1], e_verts[2]],
                [e_verts[0], v_verts[1], e_verts[1]],
                [e_verts[1], v_verts[2], e_verts[2]])
  }
  
  return { positions: npositions, cells: ncells }
}
769871
def toVerticesAndPolygons(self):
        
        offset = 1.234567890
        verts = []
        polys = []
        vertexIndexMap = {}
        count = 0
        for poly in self.polygons:
            verts = poly.vertices
            cell = []
            for v in poly.vertices:
                p = v.pos
                # use string key to remove degeneracy associated
                # very close points. The format %.10e ensures that
                # points differing in the 11 digits and higher are 
                # treated as the same. For instance 1.2e-10 and 
                # 1.3e-10 are essentially the same.
                vKey = '%.10e,%.10e,%.10e' % (p[0] + offset, 
                                              p[1] + offset, 
                                              p[2] + offset)
                if not vKey in vertexIndexMap:
                    vertexIndexMap[vKey] = len(vertexIndexMap)
                index = vertexIndexMap[vKey]
                cell.append(index)
                count += 1
            polys.append(cell)
        # sort by index
        sortedVertexIndex = sorted(vertexIndexMap.items(),
                                   key=operator.itemgetter(1))
        verts = []
        for v, i in sortedVertexIndex:
            p = []
            for c in v.split(','):
                p.append(float(c) - offset)
            verts.append(tuple(p))
        return verts, polys, count
1036156
def referenced_vertices(self):
        
        # no entities no reference
        if len(self.entities) == 0:
            return np.array([], dtype=np.int64)
        referenced = np.concatenate([e.points for e in self.entities])
        referenced = np.unique(referenced.astype(np.int64))

        return referenced
589237
private void init(DelaunayTriangulation delaunay, int xCellCount, int yCellCount, BoundingBox region) {
        indexDelaunay = delaunay;
        indexRegion = region;
        x_size = region.getWidth() / yCellCount;
        y_size = region.getHeight() / xCellCount;

        // The grid will hold a trinagle for each cell, so a point (x,y) will lie
        // in the cell representing the grid partition of region to a
        //  xCellCount on yCellCount grid
        grid = new Triangle[xCellCount][yCellCount];

        Triangle colStartTriangle = indexDelaunay.find(middleOfCell(0, 0));
        updateCellValues(0, 0, xCellCount - 1, yCellCount - 1, colStartTriangle);
    }