diff --git a/meshmode/interop/firedrake/connection.py b/meshmode/interop/firedrake/connection.py index 9720751a6d87e1e606647e6d73ef720ef49aa58c..e0de6fc9fcf3b6c64a9919b89ba979dffb85c8cd 100644 --- a/meshmode/interop/firedrake/connection.py +++ b/meshmode/interop/firedrake/connection.py @@ -53,7 +53,7 @@ def _reorder_nodes(orient, nodes, flip_matrix, unflip=False): :arg orient: An array of shape *(nelements)* of orientations, >0 for positive, <0 for negative - :arg nodes: a *(nelements, nunit_nodes)* or shaped array of nodes + :arg nodes: a *(nelements, nunit_nodes)* shaped array of nodes :arg flip_matrix: The matrix used to flip each negatively-oriented element :arg unflip: If *True*, use transpose of *flip_matrix* to @@ -647,29 +647,42 @@ class ToFiredrakeConnection(FiredrakeConnection): fd_mesh, fd_cell_order, perm2cells = \ export_mesh_to_firedrake(discr.mesh, group_nr, comm) fspace = FunctionSpace(fd_mesh, 'DG', el_group.order) + # get firedrake unit nodes and map onto meshmode reference element + dim = fspace.mesh().topological_dimension() + fd_ref_cell_to_mm = get_affine_reference_simplex_mapping(dim, True) + fd_unit_nodes = get_finat_element_unit_nodes(fspace.finat_element) + fd_unit_nodes = fd_ref_cell_to_mm(fd_unit_nodes) + # To get the meshmode to firedrake node assocation, we need to handle # local vertex reordering and cell reordering. # - # Get a copy of cell_node_list with local vertex reordering undone - reordered_cell_node_list = np.copy(fspace.cell_node_list) + # **_cell_node holds the node nrs in shape *(ncells, nunit_nodes)* + fd_cell_node = fspace.cell_node_list + mm_cell_node = el_group.view(np.arange(discr.nnodes)) + reordering_arr = np.arange(el_group.nnodes, dtype=fd_cell_node.dtype) for perm, cells in six.iteritems(perm2cells): - perm_inv = tuple(np.argsort(perm)) + # reordering_arr[i] should be the fd node corresponding to meshmode + # node i + # + # The jth meshmode cell corresponds to the fd_cell_order[j]th + # firedrake cell. If *nodeperm* is the permutation of local nodes + # applied to the *j*th meshmode cell, the firedrake node + # fd_cell_node[fd_cell_order[j]][k] corresponds to the + # mm_cell_node[j, nodeperm[k]]th meshmode node. + # + # Note that the permutation on the unit nodes may not be the + # same as the permutation on the barycentric coordinates (*perm*). + # Importantly, the permutation is derived from getting a flip + # matrix from the Firedrake unit nodes, not necessarily the meshmode + # unit nodes + # flip_mat = get_simplex_element_flip_matrix(el_group.order, - el_group.unit_nodes, - perm_inv) - reordered_cell_node_list[cells] = \ - np.einsum("ij,jk->ik", - fspace.cell_node_list[cells], - np.rint(flip_mat)) - - # making reordering_arr using the discr then assign using a view of it - reordering_arr = np.arange(discr.nnodes, - dtype=fspace.cell_node_list.dtype) - by_cell_view = el_group.view(reordering_arr) - by_cell_view = reordered_cell_node_list[fd_cell_order] # noqa : F841 - # we only want to keep the part relevant to el_group, the rest of the - # array was just there so we could use *el_group.view* - reordering_arr = reordering_arr[el_group.node_nr_base:el_group.nnodes] + fd_unit_nodes, + np.argsort(perm)) + flip_mat = np.rint(flip_mat).astype(np.int32) + fd_permuted_cell_node = np.matmul(fd_cell_node[fd_cell_order[cells]], + flip_mat.T) + reordering_arr[mm_cell_node[cells]] = fd_permuted_cell_node super(ToFiredrakeConnection, self).__init__(discr, fspace, diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py index 928ba1ba8fc5d629c862ada2fa5da2ea4ac14581..25304e8f71f05d08da2f2a4a3ffbdb1cd2f52a90 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -75,14 +75,14 @@ def _get_firedrake_nodal_info(fdrake_mesh_topology, cells_to_use=None): # If you don't understand dmplex, look at the PETSc reference # here: https://cse.buffalo.edu/~knepley/classes/caam519/CSBook.pdf # used to get topology info - # FIXME... not sure how to get around the private access + # FIXME : not sure how to get around the private access top_dm = top._topology_dm # Get range of dmplex ids for cells, facets, and vertices f_start, f_end = top_dm.getHeightStratum(1) v_start, v_end = top_dm.getDepthStratum(0) - # FIXME... not sure how to get around the private accesses + # FIXME : not sure how to get around the private accesses # Maps dmplex vert id -> firedrake vert index vert_id_dmp_to_fd = top._vertex_numbering.getOffset @@ -678,8 +678,9 @@ def export_mesh_to_firedrake(mesh, group_nr=None, comm=None): *mesh.groups[group_nr].vertex_indices*) corresponds to the *fdrake_cell_ordering[i]*th :mod:`firedrake` cell * *perm2cell* is a dictionary, mapping tuples to - lists of meshmode element indices. Each meshmode element index - appears in exactly one of these lists. The corresponding + 1-D numpy arrays of meshmode element indices. + Each meshmode element index + appears in exactly one of these arrays. The corresponding tuple describes how firedrake reordered the local vertex indices on that cell. In particular, if *c* is in the list *perm2cell[p]* for a tuple *p*, then @@ -714,7 +715,7 @@ def export_mesh_to_firedrake(mesh, group_nr=None, comm=None): if comm is None: from pyop2.mpi import COMM_WORLD comm = COMM_WORLD - # FIXME : this is a private member..., and there are some below + # FIXME : not sure how to get around the private accesses import firedrake.mesh as fd_mesh plex = fd_mesh._from_cell_list(group.dim, cells, coords, comm) # Nb : One might be tempted to pass reorder=False and thereby save some @@ -758,6 +759,8 @@ def export_mesh_to_firedrake(mesh, group_nr=None, comm=None): perm = tuple(np.argsort(mm_order)[np.argsort(np.argsort(fdrake_order))]) perm2cells.setdefault(perm, []) perm2cells[perm].append(mm_cell_id) + perm2cells = {perm: np.array(cells) + for perm, cells in six.iteritems(perm2cells)} # Now make a coordinates function from firedrake import VectorFunctionSpace, Function @@ -791,8 +794,8 @@ def export_mesh_to_firedrake(mesh, group_nr=None, comm=None): flip_mat = get_simplex_element_flip_matrix(group.order, group.unit_nodes, perm) - flipped_group_nodes[:, cells, :] = \ - np.einsum("ijk,ke->ije", group.nodes[:, cells, :], flip_mat) + flipped_group_nodes[:, cells, :] = np.matmul(group.nodes[:, cells, :], + flip_mat.T) # reorder to firedrake cell ordering reordered_cell_node_list = coords_fspace.cell_node_list[cell_index_mm2fd]