def _refine(node_coords, cells_nodes, edge_nodes, cells_edges):
num_nodes = len(node_coords)<EOL>num_new_nodes = len(edge_nodes)<EOL>node_coords.resize(num_nodes + num_new_nodes, <NUM_LIT:3>, refcheck=False)<EOL>new_node_gid = num_nodes<EOL>num_edges = len(edge_nodes)<EOL>num_cells = len(cells_nodes)<EOL>assert num_cells == len(cells_edges)<EOL>num_new_edges = <NUM_LIT:2> * num_edges + <NUM_LIT:3> * num_cells<EOL>new_edges_nodes = numpy.empty(num_new_edges, dtype=numpy.dtype((int, <NUM_LIT:2>)))<EOL>new_edge_gid = <NUM_LIT:0><EOL>num_new_cells = <NUM_LIT:4> * num_cells<EOL>new_cells_nodes = numpy.empty(num_new_cells, dtype=numpy.dtype((int, <NUM_LIT:3>)))<EOL>new_cells_edges = numpy.empty(num_new_cells, dtype=numpy.dtype((int, <NUM_LIT:3>)))<EOL>new_cell_gid = <NUM_LIT:0><EOL>is_edge_divided = numpy.zeros(num_edges, dtype=bool)<EOL>edge_midpoint_gids = numpy.empty(num_edges, dtype=int)<EOL>edge_newedges_gids = numpy.empty(num_edges, dtype=numpy.dtype((int, <NUM_LIT:2>)))<EOL>for cell_id, cell in enumerate(zip(cells_edges, cells_nodes)):<EOL><INDENT>cell_edges, cell_nodes = cell<EOL>local_edge_midpoint_gids = numpy.empty(<NUM_LIT:3>, dtype=int)<EOL>local_edge_newedges = numpy.empty(<NUM_LIT:3>, dtype=numpy.dtype((int, <NUM_LIT:2>)))<EOL>local_neighbor_midpoints = [[], [], []]<EOL>local_neighbor_newedges = [[], [], []]<EOL>for k, edge_gid in enumerate(cell_edges):<EOL><INDENT>edgenodes_gids = edge_nodes[edge_gid]<EOL>if is_edge_divided[edge_gid]:<EOL><INDENT>local_edge_midpoint_gids[k] = edge_midpoint_gids[edge_gid]<EOL>local_edge_newedges[k] = edge_newedges_gids[edge_gid]<EOL><DEDENT>else:<EOL><INDENT>node_coords[new_node_gid] = <NUM_LIT:0.5> * (<EOL>node_coords[edgenodes_gids[<NUM_LIT:0>]] + node_coords[edgenodes_gids[<NUM_LIT:1>]]<EOL>)<EOL>local_edge_midpoint_gids[k] = new_node_gid<EOL>new_node_gid += <NUM_LIT:1><EOL>edge_midpoint_gids[edge_gid] = local_edge_midpoint_gids[k]<EOL>new_edges_nodes[new_edge_gid] = numpy.array(<EOL>[edgenodes_gids[<NUM_LIT:0>], local_edge_midpoint_gids[k]]<EOL>)<EOL>new_edge_gid += <NUM_LIT:1><EOL>new_edges_nodes[new_edge_gid] = numpy.array(<EOL>[local_edge_midpoint_gids[k], edgenodes_gids[<NUM_LIT:1>]]<EOL>)<EOL>new_edge_gid += <NUM_LIT:1><EOL>local_edge_newedges[k] = [new_edge_gid - <NUM_LIT:2>, new_edge_gid - <NUM_LIT:1>]<EOL>edge_newedges_gids[edge_gid] = local_edge_newedges[k]<EOL>is_edge_divided[edge_gid] = True<EOL><DEDENT>edgenodes_lids = [<EOL>numpy.nonzero(cell_nodes == edgenodes_gids[<NUM_LIT:0>])[<NUM_LIT:0>][<NUM_LIT:0>],<EOL>numpy.nonzero(cell_nodes == edgenodes_gids[<NUM_LIT:1>])[<NUM_LIT:0>][<NUM_LIT:0>],<EOL>]<EOL>local_neighbor_midpoints[edgenodes_lids[<NUM_LIT:0>]].append(<EOL>local_edge_midpoint_gids[k]<EOL>)<EOL>local_neighbor_midpoints[edgenodes_lids[<NUM_LIT:1>]].append(<EOL>local_edge_midpoint_gids[k]<EOL>)<EOL>local_neighbor_newedges[edgenodes_lids[<NUM_LIT:0>]].append(local_edge_newedges[k][<NUM_LIT:0>])<EOL>local_neighbor_newedges[edgenodes_lids[<NUM_LIT:1>]].append(local_edge_newedges[k][<NUM_LIT:1>])<EOL><DEDENT>new_edge_opposite_of_local_node = numpy.empty(<NUM_LIT:3>, dtype=int)<EOL>for k in range(<NUM_LIT:3>):<EOL><INDENT>new_edges_nodes[new_edge_gid] = local_neighbor_midpoints[k]<EOL>new_edge_opposite_of_local_node[k] = new_edge_gid<EOL>new_edge_gid += <NUM_LIT:1><EOL><DEDENT>new_cells_nodes[new_cell_gid] = local_edge_midpoint_gids<EOL>new_cells_edges[new_cell_gid] = new_edge_opposite_of_local_node<EOL>new_cell_gid += <NUM_LIT:1><EOL>for k in range(<NUM_LIT:3>):<EOL><INDENT>new_cells_nodes[new_cell_gid] = numpy.array(<EOL>[<EOL>cells_nodes[cell_id][k],<EOL>local_neighbor_midpoints[k][<NUM_LIT:0>],<EOL>local_neighbor_midpoints[k][<NUM_LIT:1>],<EOL>]<EOL>)<EOL>new_cells_edges[new_cell_gid] = numpy.array(<EOL>[<EOL>new_edge_opposite_of_local_node[k],<EOL>local_neighbor_newedges[k][<NUM_LIT:0>],<EOL>local_neighbor_newedges[k][<NUM_LIT:1>],<EOL>]<EOL>)<EOL>new_cell_gid += <NUM_LIT:1><EOL><DEDENT><DEDENT>return node_coords, new_cells_nodes, new_edges_nodes, new_cells_edges<EOL>
Canonically refine a mesh by inserting nodes at all edge midpoints
    and make four triangular elements where there was one.
    This is a very crude refinement; don't use for actual applications.