Add raster mesh - hotfix (#46)

addresses #44

Author: kip-hart

CHANGELOG.rst

@@ -6,6 +6,13 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog`_,
 and this project adheres to `Semantic Versioning`_.
 
+`1.5.1`_ - 2021-09-09
+--------------------------
+Fixed
+'''''''
+- Plaid issue in 2D raster mesh array output.
+- Initialization and plotting of 3D raster meshes.
+
 `1.5.0`_ - 2021-09-08
 --------------------------
 Added
@@ -243,7 +250,8 @@ Added
 
 .. LINKS
 
-.. _`Unreleased`: https://github.com/kip-hart/MicroStructPy/compare/v1.5.0...HEAD
+.. _`Unreleased`: https://github.com/kip-hart/MicroStructPy/compare/v1.5.1...HEAD
+.. _`1.5.1`: https://github.com/kip-hart/MicroStructPy/compare/v1.5.0...v1.5.1
 .. _`1.5.0`: https://github.com/kip-hart/MicroStructPy/compare/v1.4.10...v1.5.0
 .. _`1.4.10`: https://github.com/kip-hart/MicroStructPy/compare/v1.4.9...v1.4.10
 .. _`1.4.9`: https://github.com/kip-hart/MicroStructPy/compare/v1.4.8...v1.4.9

docs/source/api/meshing/index.rst

@@ -3,11 +3,11 @@ microstructpy.meshing
 
 .. automodule:: microstructpy.meshing
 
-The meshing module contains two mesh classes,
-the :class:`.PolyMesh` and the :class:`.TriMesh`.
+The meshing module contains three mesh classes,
+the :class:`.PolyMesh`, the :class:`.RasterMesh` and the :class:`.TriMesh`.
 The polygonal mesh contains a 2D or 3D tessellation of the microstructure
-domain, while the triangular mesh is more suitable for direct numerical
-simulation (finite element analysis).
+domain, while the raster aand triangular meshes are more suitable for
+direct numerical simulation (finite element analysis).
 
 .. only:: html
 
@@ -16,4 +16,5 @@ simulation (finite element analysis).
 .. toctree::
 
     polymesh
+    rastermesh
     trimesh

docs/source/api/meshing/rastermesh.rst

@@ -0,0 +1,11 @@
+.. _api_meshing_rastermesh:
+
+.. currentmodule:: microstructpy.meshing
+
+microstructpy.meshing.RasterMesh
+================================
+
+.. autoclass:: RasterMesh
+    :members:
+    :undoc-members:
+    :inherited-members:

src/microstructpy/__init__.py

@@ -4,4 +4,4 @@
 import microstructpy.seeding
 import microstructpy.verification
 
-__version__ = '1.5.0'
+__version__ = '1.5.1'

src/microstructpy/meshing/trimesh.py

@@ -578,31 +578,7 @@ def plot(self, index_by='element', material=[], loc=0, **kwargs):
             xlim = [float('inf'), -float('inf')]
             ylim = [float('inf'), -float('inf')]
         if n_dim == 2:
-            simps = np.array(self.elements)
-            pts = np.array(self.points)
-            xy = pts[simps, :]
-
-            plt_kwargs = {}
-            for key, value in kwargs.items():
-                if type(value) in (list, np.array):
-                    plt_value = []
-                    for e_num, e_att in enumerate(self.element_attributes):
-                        if index_by == 'element':
-                            ind = e_num
-                        elif index_by == 'attribute':
-                            ind = int(e_att)
-                        else:
-                            e_str = 'Cannot index by {}.'.format(index_by)
-                            raise ValueError(e_str)
-                        v = value[ind]
-                        plt_value.append(v)
-                else:
-                    plt_value = value
-                plt_kwargs[key] = plt_value
-
-            pc = collections.PolyCollection(xy, **plt_kwargs)
-            ax.add_collection(pc)
-            ax.autoscale_view()
+            _plot_2d(ax, self, index_by, **kwargs)
         else:
             if n_obj > 0:
                 zlim = ax.get_zlim()
@@ -804,7 +780,7 @@ def from_polymesh(cls, polymesh, mesh_size, phases=None):
         elem_atts = np.full(cens.shape[0], -1)
         for r_num, region in enumerate(polymesh.regions):
             # A. Create a bounding box
-            r_kps = np.unique([polymesh.facets[f] for f in region])
+            r_kps = np.unique([k for f in region for k in polymesh.facets[f]])
             r_pts = p_pts[r_kps]
             r_mins = r_pts.min(axis=0)
             r_maxs = r_pts.max(axis=0)
@@ -1260,8 +1236,12 @@ def as_array(self, element_attributes=True):
 
         Args:
             element_attributes (bool): *(optional)* Flag to return element
-            attributes in the array. Set to True return attributes and set to
-            False to return element indices. Defaults to True.
+                attributes in the array. Set to True return attributes and
+                set to False to return element indices. Defaults to True.
+
+        Returns:
+            numpy.ndarray: Array of values of element atttributes, or indices.
+
         """
         # 1. Convert 1st node of each element into array indices
         pts = np.array(self.points)
@@ -1270,7 +1250,7 @@ def as_array(self, element_attributes=True):
 
         corner_pts = pts[np.array(self.elements)[:, 0]]
         rel_pos = corner_pts - mins
-        elem_tups = (rel_pos / sz).astype(int)
+        elem_tups = np.round(rel_pos / sz).astype(int)
 
         # 2. Create array full of -1 values
         inds_maxs = elem_tups.max(axis=0)
@@ -1290,7 +1270,119 @@ def as_array(self, element_attributes=True):
     # ----------------------------------------------------------------------- #
     # Plot Function                                                           #
     # ----------------------------------------------------------------------- #
-    # Inherited from TriMesh
+    def plot(self, index_by='element', material=[], loc=0, **kwargs):
+        """Plot the mesh.
+
+        This method plots the mesh using matplotlib.
+        In 2D, this creates a :class:`matplotlib.collections.PolyCollection`
+        and adds it to the current axes.
+        In 3D, it creates a
+        :meth:`mpl_toolkits.mplot3d.axes3d.Axes3D.voxels` and
+        adds it to the current axes.
+        The keyword arguments are passed though to matplotlib.
+
+        Args:
+            index_by (str): *(optional)* {'element' | 'attribute'}
+                Flag for indexing into the other arrays passed into the
+                function. For example,
+                ``plot(index_by='attribute', color=['blue', 'red'])`` will plot
+                the elements with ``element_attribute`` equal to 0 in blue, and
+                elements with ``element_attribute`` equal to 1 in red.
+                Defaults to 'element'.
+            material (list): *(optional)* Names of material phases. One entry
+                per material phase (the ``index_by`` argument is ignored).
+                If this argument is set, a legend is added to the plot with
+                one entry per material. Note that the ``element_attributes``
+                must be the material numbers for the legend to be
+                formatted properly.
+            loc (int or str): *(optional)* The location of the legend,
+                if 'material' is specified. This argument is passed directly
+                through to :func:`matplotlib.pyplot.legend`. Defaults to 0,
+                which is 'best' in matplotlib.
+            **kwargs: Keyword arguments that are passed through to matplotlib.
+
+        """
+        n_dim = len(self.points[0])
+        if n_dim == 2:
+            ax = plt.gca()
+        else:
+            ax = plt.gcf().gca(projection=Axes3D.name)
+        n_obj = _misc.ax_objects(ax)
+        if n_obj > 0:
+            xlim = ax.get_xlim()
+            ylim = ax.get_ylim()
+        else:
+            xlim = [float('inf'), -float('inf')]
+            ylim = [float('inf'), -float('inf')]
+        if n_dim == 2:
+            _plot_2d(ax, self, index_by, **kwargs)
+        else:
+            if n_obj > 0:
+                zlim = ax.get_zlim()
+            else:
+                zlim = [float('inf'), -float('inf')]
+
+            inds = self.as_array(element_attributes=index_by=='attribute')
+            plt_kwargs = {}
+            for key, value in kwargs.items():
+                if type(value) in (list, np.array):
+                    plt_value = np.empty(inds.shape, dtype=object)
+                    for i, val_i in enumerate(value):
+                        plt_value[inds == i] = val_i
+                else:
+                    plt_value = value
+                plt_kwargs[key] = plt_value
+
+            # Scale axes
+            pts = np.array(self.points)
+            mins = pts.min(axis=0)
+            sz = self.mesh_size
+            pt_tups = np.round((pts - mins) / sz).astype(int)
+            maxs = pt_tups.max(axis=0)
+            grids = np.indices(maxs + 1, dtype=float)
+            for pt, pt_tup in zip(pts, pt_tups):
+                for i, x in enumerate(pt):
+                    grids[i][tuple(pt_tup)] = x
+            ax.voxels(*grids, inds >= 0, **plt_kwargs)
+
+        # Add legend
+        if material and index_by == 'attribute':
+            p_kwargs = [{'label': m} for m in material]
+            for key, value in kwargs.items():
+                if type(value) not in (list, np.array):
+                    for kws in p_kwargs:
+                        kws[key] = value
+
+                for i, m in enumerate(material):
+                    if type(value) in (list, np.array):
+                        p_kwargs[i][key] = value[i]
+                    else:
+                        p_kwargs[i][key] = value
+
+            # Replace plural keywords
+            for p_kw in p_kwargs:
+                for kw in _misc.mpl_plural_kwargs:
+                    if kw in p_kw:
+                        p_kw[kw[:-1]] = p_kw[kw]
+                        del p_kw[kw]
+            handles = [patches.Patch(**p_kw) for p_kw in p_kwargs]
+            ax.legend(handles=handles, loc=loc)
+
+        # Adjust Axes
+        mins = np.array(self.points).min(axis=0)
+        maxs = np.array(self.points).max(axis=0)
+        xlim = (min(xlim[0], mins[0]), max(xlim[1], maxs[0]))
+        ylim = (min(ylim[0], mins[1]), max(ylim[1], maxs[1]))
+        if n_dim == 2:
+            plt.axis('square')
+            plt.xlim(xlim)
+            plt.ylim(ylim)
+        elif n_dim == 3:
+            zlim = (min(zlim[0], mins[2]), max(zlim[1], maxs[2]))
+            ax.set_xlim(xlim)
+            ax.set_ylim(ylim)
+            ax.set_zlim(zlim)
+            _misc.axisEqual3D(ax)
 
 
 def facet_check(neighs, polymesh, phases):
@@ -1794,3 +1886,31 @@ def _facet_in_normal(pts, cen_pt):
     vn *= sgn  # flip so center is inward
     un = vn / np.linalg.norm(vn)
     return un, pts.mean(axis=0)
+
+
+def _plot_2d(ax, mesh, index_by, **kwargs):
+    simps = np.array(mesh.elements)
+    pts = np.array(mesh.points)
+    xy = pts[simps, :]
+
+    plt_kwargs = {}
+    for key, value in kwargs.items():
+        if type(value) in (list, np.array):
+            plt_value = []
+            for e_num, e_att in enumerate(mesh.element_attributes):
+                if index_by == 'element':
+                    ind = e_num
+                elif index_by == 'attribute':
+                    ind = int(e_att)
+                else:
+                    e_str = 'Cannot index by {}.'.format(index_by)
+                    raise ValueError(e_str)
+                v = value[ind]
+                plt_value.append(v)
+        else:
+            plt_value = value
+        plt_kwargs[key] = plt_value
+
+    pc = collections.PolyCollection(xy, **plt_kwargs)
+    ax.add_collection(pc)
+    ax.autoscale_view()