From 70657ca9e56c3fea77d3efd3b336797e8cd91ad2 Mon Sep 17 00:00:00 2001
From: Ben Sepanski <ben_sepanski@baylor.edu>
Date: Thu, 28 May 2020 11:41:52 -0500
Subject: [PATCH] Added an importer for FIAT's simplex cell
---
meshmode/interop/fiat/simplex_cell.py | 123 ++++++++++++++++++++++++++
1 file changed, 123 insertions(+)
create mode 100644 meshmode/interop/fiat/simplex_cell.py
diff --git a/meshmode/interop/fiat/simplex_cell.py b/meshmode/interop/fiat/simplex_cell.py
new file mode 100644
index 00000000..51c61a04
--- /dev/null
+++ b/meshmode/interop/fiat/simplex_cell.py
@@ -0,0 +1,123 @@
+from warnings import warn # noqa
+import numpy as np
+import numpy.linalg as la
+
+from meshmode.mesh.interop import ExternalImporter
+
+
+__doc__ = """
+Interoperators
+--------------
+.. autoclass:: SimplexCellAnalog
+ :members:
+"""
+
+# {{{ Compute an affine mapping from given input/outputs
+
+def get_affine_mapping(reference_vects, vects):
+ r"""
+ Returns (mat, shift),
+ a matrix *mat* and vector *shift* which maps the
+ *i*th vector in *reference_vects* to the *i*th vector in *vects by
+
+ ..math::
+
+ A ri + b -> vi, \qquad A = mat, b = shift
+
+ :arg reference_vects: An np.array of *n* vectors of dimension *ref_dim*
+ :arg vects: An np.array of *n* vectors of dimension *dim*, with
+ *ref_dim* <= *dim*.
+
+ NOTE : Should be shape *(ref_dim, nvectors)*, *(dim, nvectors)* respectively.
+
+ *mat* will have shape *(dim, ref_dim)*, *shift* will have shape *(dim,)*
+ """
+ # Make sure both have same number of vectors
+ ref_dim, num_vects = reference_vects.shape
+ assert num_vects == vects.shape[1]
+
+ # Make sure d1 <= d2 (see docstring)
+ dim = vects.shape[0]
+ assert ref_dim <= dim
+
+ # If there is only one vector, set M = I, b = vect - reference
+ if num_vects == 1:
+ mat = np.eye(dim, ref_dim)
+ shift = vects[:, 0] - np.matmul(mat, reference_vects[:, 0])
+ else:
+ ref_span_vects = reference_vects[:, 1:] - reference_vects[:, 0, np.newaxis]
+ span_vects = vects[:, 1:] - vects[:, 0, np.newaxis]
+ mat = la.solve(ref_span_vects, span_vects)
+ shift = -np.matmul(mat, reference_vects[:, 0]) + vects[:, 0]
+
+ return mat, shift
+
+# }}}
+
+
+# {{{ Interoperator for FIAT's reference_element.Simplex
+
+class FIATSimplexCellImporter(ExternalImporter):
+ """
+ Importer for a :mod:`FIAT` simplex cell.
+ There is no data created from the :attr:`from_data`. Instead,
+ the data is simply used to obtain FIAT's reference
+ nodes according to :mod:`modepy`'s reference coordinates
+ using :meth:`make_modepy_points`.
+
+ .. attribute:: from_data
+
+ An instance of :class:`fiat.FIAT.reference_element.Simplex`.
+
+ .. attribute:: to_data
+
+ :raises ValueError: If accessed.
+ """
+ def __init__(self, cell):
+ """
+ :arg cell: a :class:`fiat.FIAT.reference_element.Simplex`.
+ """
+ # Ensure this cell is actually a simplex
+ from FIAT.reference_element import Simplex
+ assert isinstance(cell, Simplex)
+
+ super(SimplexCellInteroperator, self).__init__(cell)
+
+ # Stored as (dim, nunit_vertices)
+ from modepy.tools import unit_vertices
+ self._unit_vertices = unit_vertices(cell.get_dimension()).T
+
+ # Maps FIAT reference vertices to :mod:`meshmode`
+ # unit vertices by x -> Ax + b, where A is :attr:`_mat`
+ # and b is :attr:`_shift`
+ reference_vertices = np.array(cell.vertices).T
+ self._mat, self._shift = get_affine_mapping(reference_vertices,
+ self._unit_vertices)
+
+ def make_modepy_points(self, dim, entity_id, order):
+ """
+ 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 <fiat.FIAT.reference_element.Cell.make_points>`_),
+ but the unit nodes are (affinely) mapped to :mod:`modepy`
+ `unit coordinates <https://documen.tician.de/modepy/nodes.html>`_.
+
+ :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
+ """
+ points = self.from_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
+
+# }}}
--
GitLab