diff --git a/meshmode/interop/FInAT/lagrange_element.py b/meshmode/interop/FInAT/lagrange_element.py index b3d4c50f977e9e970a2ac39a563b55889fad4310..4b3d977c04ba8a0db2517f78290e7b7147d275b0 100644 --- a/meshmode/interop/FInAT/lagrange_element.py +++ b/meshmode/interop/FInAT/lagrange_element.py @@ -22,7 +22,6 @@ THE SOFTWARE. import numpy as np import numpy.linalg as la -import six from meshmode.interop import ExternalImportHandler from meshmode.interop.fiat import FIATSimplexCellImporter @@ -44,6 +43,8 @@ class FinatLagrangeElementImporter(ExternalImportHandler): :param finat_element: A FInAT element of type :class:`finat.fiat_elements.Lagrange` or :class:`finat.fiat_elements.DiscontinuousLagrange` + (or :class:`finat.spectral.GaussLegendre` or + :class:`finat.spectral.GaussLobattoLegendre` in 1D) which uses affine mapping of the basis functions (i.e. ``finat_element.mapping`` must be ``"affine"``) @@ -51,6 +52,7 @@ class FinatLagrangeElementImporter(ExternalImportHandler): :raises TypeError: If :param:`finat_element` is not of type :class:`finat.fiat_elements.Lagrange` or :class:`finat.fiat_elements.DiscontinuousLagrange` + (or the 1D classes listed above) :raises ValueError: If :param:`finat_element` does not use affine mappings of the basis functions @@ -91,39 +93,33 @@ class FinatLagrangeElementImporter(ExternalImportHandler): if they have not already been computed. """ if self._unit_nodes is None or self._unit_vertex_indices is None: + # FIXME : This should work, but uses some private info # {{{ Compute unit nodes - node_nr_to_coords = {} - unit_vertex_indices = [] - - # FIXME : This works, but is very ad-hoc. It is probably better - # to get permission from the FInAT people to reach into - # the fiat element and get the nodes explicitly - # Get unit nodes - for dim, element_nrs in six.iteritems( - self.data.entity_support_dofs()): - for element_nr, node_list in six.iteritems(element_nrs): - # Get the nodes on the element (in meshmode reference coords) - pts_on_element = self.cell_importer.make_points( - dim, element_nr, self.data.degree) - # Record any new nodes - i = 0 - for node_nr in node_list: - if node_nr not in node_nr_to_coords and \ - i < len(pts_on_element): - node_nr_to_coords[node_nr] = pts_on_element[i] - i += 1 - # If is a vertex, store the index - if dim == 0: - unit_vertex_indices.append(node_nr) - - # store vertex indices - self._unit_vertex_indices = np.array(sorted(unit_vertex_indices)) - - # Convert unit_nodes to array, then change to (dim, nunit_nodes) - # from (nunit_nodes, dim) - unit_nodes = np.array([node_nr_to_coords[i] for i in - range(len(node_nr_to_coords))]) - self._unit_nodes = unit_nodes.T.copy() + + # each point evaluator is a function p(f) evaluates f at a node, + # so we need to evaluate each point evaluator at the identity to + # recover the nodes + point_evaluators = self.data._element.dual.nodes + unit_nodes = [p(lambda x: x) for p in point_evaluators] + unit_nodes = np.array(unit_nodes).T + self._unit_nodes = \ + self.cell_importer.affinely_map_firedrake_to_meshmode(unit_nodes) + + # Is this safe?, I think so bc on a reference element + close = 1e-8 + # Get vertices as (dim, nunit_vertices) + unit_vertices = np.array(self.data.cell.vertices).T + unit_vertices = \ + self.cell_importer.affinely_map_firedrake_to_meshmode(unit_vertices) + self._unit_vertex_indices = [] + for n_ndx in range(self._unit_nodes.shape[1]): + for v_ndx in range(unit_vertices.shape[1]): + diff = self._unit_nodes[:, n_ndx] - unit_vertices[:, v_ndx] + if np.max(np.abs(diff)) < close: + self._unit_vertex_indices.append(n_ndx) + break + + self._unit_vertex_indices = np.array(self._unit_vertex_indices) # }}} @@ -135,9 +131,10 @@ class FinatLagrangeElementImporter(ExternalImportHandler): def unit_vertex_indices(self): """ - :return: An array of shape *(dim+1,)* of indices + :return: A numpy integer array of indices so that *self.unit_nodes()[self.unit_vertex_indices()]* - are the vertices of the reference element. + are the nodes of the reference element which coincide + with its vertices (this is possibly empty). """ self._compute_unit_vertex_indices_and_nodes() return self._unit_vertex_indices @@ -203,12 +200,12 @@ class FinatLagrangeElementImporter(ExternalImportHandler): def make_resampling_matrix(self, element_grp): """ - :param element_grp: A - :class:`meshmode.discretization.InterpolatoryElementGroupBase` whose - basis functions span the same space as the FInAT element. - :return: A matrix which resamples a function sampled at - the firedrake unit nodes to a function sampled at - *element_grp.unit_nodes()* (by matrix multiplication) + :param element_grp: A + :class:`meshmode.discretization.InterpolatoryElementGroupBase` whose + basis functions span the same space as the FInAT element. + :return: A matrix which resamples a function sampled at + the firedrake unit nodes to a function sampled at + *element_grp.unit_nodes()* (by matrix multiplication) """ from meshmode.discretization import InterpolatoryElementGroupBase assert isinstance(element_grp, InterpolatoryElementGroupBase), \ diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py index cf0c1622b24669b249527022470a8b524d0fd91f..c6c0ca2ca4e4241a0baa99860a083214645b4de5 100644 --- a/meshmode/interop/fiat/simplex_cell.py +++ b/meshmode/interop/fiat/simplex_cell.py @@ -105,30 +105,18 @@ class FIATSimplexCellImporter(ExternalImportHandler): self._mat, self._shift = get_affine_mapping(reference_vertices, self._unit_vertices) - def make_points(self, dim, entity_id, order): + def affinely_map_firedrake_to_meshmode(self, points): """ - Constructs a lattice of points on the *entity_id*th facet - of dimension *dim*. - - Args are exactly as in - :meth:`fiat.FIAT.reference_element.Cell.make_points` - (see `FIAT docs `_), - but the unit nodes are (affinely) mapped to :mod:`modepy` + Map points on the firedrake reference simplex to + :mod:`modepy` `unit coordinates `_. - :arg dim: Dimension of the facet we are constructing points on. - :arg entity_id: identifier to determine which facet of - dimension *dim* to construct the points on. - :arg order: Number of points to include in each direction. - - :return: an *np.array* of shape *(dim, npoints)* holding the - coordinates of each of the ver + :arg points: *n* points on the reference simplex + as a numpy array of shape *(dim, n)* + :return: A numpy array of shape *(dim, n)* wherein the + firedrake refernece simplex has been affinely mapped + to the modepy reference simplex """ - points = self.data.make_points(dim, entity_id, order) - if not points: - return points - points = np.array(points) - # Points is (npoints, dim) so have to transpose - return (np.matmul(self._mat, points.T) + self._shift[:, np.newaxis]).T + return np.matmul(self._mat, points) + self._shift[:, np.newaxis] # }}} diff --git a/meshmode/interop/firedrake/function_space_shared_data.py b/meshmode/interop/firedrake/function_space_shared_data.py index 52c0e6d2f547dba44eae1e54f4130dd59613f507..9bd03f27a088041249736315c5c3a7b45ec41bc4 100644 --- a/meshmode/interop/firedrake/function_space_shared_data.py +++ b/meshmode/interop/firedrake/function_space_shared_data.py @@ -146,7 +146,6 @@ def reordering_array(mesh_importer, key, fspace_data): (order.shape[0]//nunit_nodes, nunit_nodes) + order.shape[1:]) flip_mat = finat_element_importer.flip_matrix() - print("ORIENT=", mesh_importer._orient) reorder_nodes(mesh_importer.orientations(), new_order, flip_mat, unflip=firedrake_to_meshmode) new_order = new_order.flatten() diff --git a/meshmode/interop/firedrake/mesh_geometry.py b/meshmode/interop/firedrake/mesh_geometry.py index 8ce88f266fdd54738107e1bb75f94ad9bb22fe78..da965cbe4148b7f264622b5a83b99d91d930bd11 100644 --- a/meshmode/interop/firedrake/mesh_geometry.py +++ b/meshmode/interop/firedrake/mesh_geometry.py @@ -180,6 +180,14 @@ class FiredrakeMeshGeometryImporter(ExternalImportHandler): fspace_importer = self.coordinates_importer.function_space_importer() finat_element_importer = fspace_importer.finat_element_importer + from finat.fiat_elements import Lagrange + from finat.spectral import GaussLobattoLegendre + if not isinstance(finat_element_importer.data, + (Lagrange, GaussLobattoLegendre)): + raise TypeError("Coordinates must live in a continuous space " + " (Lagrange or GaussLobattoLegendre), not %s" + % type(finat_element_importer.data)) + # Convert cell node list of mesh to vertex list unit_vertex_indices = finat_element_importer.unit_vertex_indices() cfspace = self.data.coordinates.function_space()