Polymesh exterior normals, output plot fixes, faster seed positioning (#23)

* stop tracking duplicate aluminum micro file

* facets on the exterior of the polymesh now have outward normals (#16)

* updated changelog

* changed plot limits, trying to get 3D axes equal

* updated trimesh plotting in CLI

* fixed typo in order of lines for empirical PDFs

* added option to maximize the minimum edge length in a polygonal mesh

* simplified seedlist indexing

* optimized seedlist positioning method

uses BFS in aabb tree for overlap detection and samples 100 positions at a time

* fixed the logo script

* fixed 2d axes behavior

* simplified foam example

* updated changelog

* added missing zlim

* updated sphinx version

* indented changelog

* removed unneeded version string functions

* bumped version

* simplified version getting

Author: kip-hart

CHANGELOG.rst

@@ -6,6 +6,25 @@ 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.3.0`_ - 2020-06-25
+--------------------------
+Added
+'''''
+- Option to reduce the presence of short edges in polygonal meshes.
+
+Changed
+'''''''
+- Optimized seed positioning algorithm by using breadth-first search
+  in the AABB tree.
+- Facets in polygonal meshes are now always defined with a positive
+  outward normal vector.
+
+Fixed
+'''''
+- Plotting of 3D meshes. 
+- Documentation for empirical PDFs.
+- Minor errors in examples.
+
 `1.2.2`_ - 2020-05-14
 --------------------------
 Fixed
@@ -88,7 +107,8 @@ Added
 
 .. LINKS
 
-.. _`Unreleased`: https://github.com/kip-hart/MicroStructPy/compare/v1.2.2...HEAD
+.. _`Unreleased`: https://github.com/kip-hart/MicroStructPy/compare/v1.3.0...HEAD
+.. _`1.3.0`: https://github.com/kip-hart/MicroStructPy/compare/v1.2.2...v1.3.0
 .. _`1.2.2`: https://github.com/kip-hart/MicroStructPy/compare/v1.2.1...v1.2.2
 .. _`1.2.1`: https://github.com/kip-hart/MicroStructPy/compare/v1.2.0...v1.2.1
 .. _`1.2.0`: https://github.com/kip-hart/MicroStructPy/compare/v1.1.2...v1.2.0

docs/requirements.txt

@@ -1,5 +1,5 @@
 matplotlib
 numpy==1.17.3
 pybind11==2.4.3
-sphinx==2.2.1
+sphinx==3.1.1
 sphinx-gallery==0.4.0

docs/source/cli/material.rst

@@ -150,12 +150,12 @@ the complete list of available distributions and their input parameters.
 
 In the case that the distribution type is "pdf" then the only other field
 should be ``filename``.
-For a PDF, the file should contain two lines: the first has the (n+1) bin
-locations and the second has the (n) bin heights.
+For a PDF, the file should contain two lines: the first has the (n) bin
+heights and the second has the (n+1) bin locations.
 A PDF file could contain, for example::
 
-    1, 2, 2.5
     0.5, 1
+    1, 2, 2.5
 
 For a CDF, the file should have two columns: the first being the size and the
 second being the CDF value.

docs/source/cli/settings.rst

@@ -40,6 +40,9 @@ quality, among other things. The default settings are:
 
             <rtol> fit </rtol>
 
+            <edge_opt> False </edge_opt>
+            <edge_opt_n_iter> 100 </edge_opt_n_iter>
+
             <mesh_max_volume> inf </mesh_max_volume>
             <mesh_min_angle> 0 </mesh_min_angle>
             <mesh_max_edge_length> inf </mesh_max_edge_length>
@@ -244,6 +247,40 @@ an rtol value that maximizes the fit between input and output distributions.
 Acceptable values of rtol are 0 to 1 inclusive, though rtol below 0.2 will
 likely result in long runtimes.
 
+edge_opt
+--------
+
+The edge_opt field provides the option to maximize the shortest edge in the
+polygonal/polyhedral mesh.
+The default is ``<edge_opt> False </edge_opt>``, which skips the optimization
+process.
+This optimization is performed by making small adjustments to the positions of
+seeds surrounding the shortest edge, assessing if the change created an
+improvement, then either a) attempting a different change for the same edge if
+there was not improvement or b) moving on to the new shortest edge.
+The optimization algorithm exits when ``edge_opt_n_iter`` iterations have been
+performed on the same edge.
+
+This flag is useful if the polygonal/polyhedral or triangular/tetrahedral are
+used in numerical simulations, such as finite element analysis.
+A high ratio of longest edge to shortest edge leads to a high ratio in maximum
+to minimum eigenvalue in FEA stiffness matrices, which can create problems for
+the FEA solver.
+Setting ``edge_opt`` to ``True`` will reduce short edges in the polygonal mesh,
+which translates into reduced short edges in the triangular mesh.
+This optimization process, however, will increase the time to generate a
+polygonal mesh.
+To track the progress of the optimizer, set ``verbose`` to ``True``.
+
+edge_opt_n_iter
+---------------
+
+This field specifies how many times the optimizer should attempt to increase
+the length of the shortest edge in the polygonal mesh.
+The default is ``<edge_opt_n_iter> 100 </edge_opt_n_iter>``, which limits the
+optimizer to 100 attempts per edge.
+This field is ignored if ``edge_opt`` is set to ``False``.
+
 mesh_max_volume
 ---------------
 

requirements.txt

@@ -1,4 +1,4 @@
-aabbtree==2.3.1
+aabbtree==2.5.0
 matplotlib==3.0.2
 pybind11==2.4.3
 MeshPy==2018.2.1

setup.py

@@ -13,7 +13,6 @@
 from setuptools import setup
 
 desc = 'Microstructure modeling, mesh generation, analysis, and visualization.'
-vs_fname = ('src/microstructpy', '_vs.py')
 
 
 def read(*fname):
@@ -26,18 +25,12 @@ def find_version(*fname):
         if line.startswith('__version__') and '=' in line:
             ver_str = line.split('=')[-1].strip().strip('\"').strip('\'')
             break
-
-    for line in read(*vs_fname).split('\n'):
-        if '+' in line:
-            tag = line.split('+')[-1].strip().strip('\"').strip('\'')
-            ver_str += '+' + tag
-            break
     return ver_str
 
 
 setup(
     name='microstructpy',
-    version=find_version('src/microstructpy', '__init__.py'),
+    version=find_version('src', 'microstructpy', '__init__.py'),
     license='MIT License',
     description=desc,
     long_description=read('README.rst'),
@@ -82,8 +75,8 @@ def find_version(*fname):
         'multi-sphere'
     ],
     install_requires=[
-        'aabbtree',
-        'matplotlib>=2.2.0,<3.1.0',
+        'aabbtree>=2.5.0',
+        'matplotlib>=3.0.0',
         'pybind11',  # must come before meshpy for successful install
         'lsq-ellipse',
         'meshpy>=2018.2.1',

src/microstructpy/__init__.py

@@ -3,7 +3,5 @@
 import microstructpy.meshing
 import microstructpy.seeding
 import microstructpy.verification
-from microstructpy._vs import _ver_str
 
-__version__ = '1.2.2'
-__version__ = _ver_str(__version__)
+__version__ = '1.3.0'

src/microstructpy/_misc.py

@@ -19,6 +19,12 @@
 
 mpl_plural_kwargs = {'edgecolors', 'facecolors', 'linewidths', 'antialiaseds',
                      'offsets'}
+plt_3d_adj = {
+    'left': 0.4,
+    'right': 1,
+    'bottom': 0,
+    'top': 0.8,
+}
 
 
 # --------------------------------------------------------------------------- #
@@ -163,3 +169,21 @@ def tangent_sphere(points, radii=None, simplices=None):
     # return results
     spheres = np.hstack((posc, rc.reshape(-1, 1)))
     return np.squeeze(spheres)
+
+
+def axisEqual3D(ax):
+    '''From stackoverflow: https://stackoverflow.com/a/19248731'''
+    extents = np.array([getattr(ax, 'get_{}lim'.format(d))() for d in 'xyz'])
+    sz = extents[:, 1] - extents[:, 0]
+    centers = np.mean(extents, axis=1)
+    maxsize = max(abs(sz))
+    r = maxsize/2
+    for ctr, dim in zip(centers, 'xyz'):
+        getattr(ax, 'set_{}lim'.format(dim))(ctr - r, ctr + r)
+
+
+def ax_objects(ax):
+    n = 0
+    for att in ['collections', 'images', 'lines', 'patches', 'texts']:
+        n += len(getattr(ax, att))
+    return n

src/microstructpy/_vs.py

@@ -174,7 +174,8 @@ def read_input(filename):
 
 
 def run(phases, domain, verbose=False, restart=True, directory='.',
-        filetypes={}, rng_seeds={}, plot_axes=True, rtol='fit',
+        filetypes={}, rng_seeds={}, plot_axes=True, rtol='fit', edge_opt=False,
+        edge_opt_n_iter=100,
         mesh_max_volume=float('inf'), mesh_min_angle=0,
         mesh_max_edge_length=float('inf'), verify=False, color_by='material',
         colormap='viridis', seeds_kwargs={}, poly_kwargs={}, tri_kwargs={}):
@@ -243,10 +244,18 @@ def run(phases, domain, verbose=False, restart=True, directory='.',
 
             The default value is ``'fit'``, which uses the mean and variance
             of the size distribution to estimate a value for rtol.
+        edge_opt (bool): *(optional)* This option will maximize the minimum
+            edge length in the PolyMesh. The seeds associated with the
+            shortest edge are displaced randomly to find improvement and
+            this process iterates until `n_iter` attempts have been made
+            for a given edge. Defaults to False.
+        edge_opt_n_iter (int): *(optional)* Maximum number of iterations per
+            edge during optimization. Ignored if `edge_opt` set to False.
+            Defaults to 100.
         mesh_max_volume (float): *(optional)* The maximum volume (area in 2D)
             of a mesh cell in the triangular mesh. Default is infinity,
             which turns off the maximum volume quality setting.
-            Value should be stritly positive.
+            Value should be strictly positive.
         mesh_min_angle (float): *(optional)* The minimum interior angle,
             in degrees,  of a cell in the triangular mesh. For 3D meshes,
             this is the dihedral angle between faces of the tetrahedron.
@@ -384,7 +393,8 @@ def run(phases, domain, verbose=False, restart=True, directory='.',
         if verbose:
             print('Creating polygon mesh.')
 
-        pmesh = PolyMesh.from_seeds(seeds, domain)
+        pmesh = PolyMesh.from_seeds(seeds, domain, edge_opt, edge_opt_n_iter,
+                                    verbose)
 
     # Write polymesh
     poly_types = filetypes.get('poly', [])
@@ -614,7 +624,6 @@ def plot_seeds(seeds, phases, domain, plot_files=[], plot_axes=True,
             ax.get_yaxis().set_visible(False)
         else:
             ax._axis3don = False
-    fig.add_axes(ax)
 
     # Plot seeds
     edge_kwargs.setdefault('edgecolors', {2: 'k', 3: 'none'}[n_dim])
@@ -624,9 +633,26 @@ def plot_seeds(seeds, phases, domain, plot_files=[], plot_axes=True,
     else:
         seeds.plot(facecolors=seed_colors, **edge_kwargs)
 
+    # Crop to Domain
+    d_lims = domain.limits
+    if n_dim == 2:
+        plt.axis('square')
+        plt.xlim(d_lims[0])
+        plt.ylim(d_lims[1])
+    elif n_dim == 3:
+        plt.gca().set_xlim(d_lims[0])
+        plt.gca().set_ylim(d_lims[1])
+        plt.gca().set_zlim(d_lims[2])
+
+        _misc.axisEqual3D(plt.gca())
+
     # Save plot
     for fname in plot_files:
-        plt.savefig(fname, bbox_inches='tight', pad_inches=0)
+        if n_dim == 3:
+            fig.subplots_adjust(**_misc.plt_3d_adj)
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0.15)
+        else:
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0)
 
     plt.close('all')
 
@@ -747,7 +773,11 @@ def plot_poly(pmesh, phases, plot_files=['polymesh.png'], plot_axes=True,
 
     # save plot
     for fname in plot_files:
-        plt.savefig(fname, bbox_inches='tight', pad_inches=0)
+        if n_dim == 3:
+            fig.subplots_adjust(**_misc.plt_3d_adj)
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0.15)
+        else:
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0)
     plt.close('all')
 
 
@@ -848,7 +878,6 @@ def plot_tri(tmesh, phases, seeds, pmesh, plot_files=[], plot_axes=True,
             ax.get_yaxis().set_visible(False)
         else:
             ax._axis3don = False
-    fig.add_axes(ax)
 
     # Determine which facets are visible
     vis_regions = set()
@@ -888,18 +917,27 @@ def plot_tri(tmesh, phases, seeds, pmesh, plot_files=[], plot_axes=True,
     # plot triangle mesh
     edge_kwargs.setdefault('linewidths', {2: 0.5, 3: 0.1}[n_dim])
     edge_kwargs.setdefault('edgecolors', 'k')
-    if given_names and color_by in ('material', 'material number'):
+    if color_by in ('material', 'material number'):
         n = len(phases)
+        if n_dim == 2:
+            cs = seed_colors
+        else:
+            cs = facet_colors
+            cs.append('none')
+
         cs = [_phase_color_by(i, phases, color_by, colormap) for i in range(n)]
         cs.append('none')
 
         old_e_att = np.copy(tmesh.element_attributes)
         old_f_att = np.copy(tmesh.facet_attributes)
         tmesh.element_attributes = [seeds[i].phase for i in old_e_att]
-        tmesh.facet_attributes = facet_phases
+        tmesh.facet_attributes = [facet_phases[a] for a in old_f_att]
 
-        tmesh.plot(facecolors=cs, index_by='attribute', material=phase_names,
-                   **edge_kwargs)
+        if given_names:
+            tmesh.plot(facecolors=cs, index_by='attribute',
+                       material=phase_names, **edge_kwargs)
+        else:
+            tmesh.plot(facecolors=cs, index_by='attribute', **edge_kwargs)
 
         tmesh.element_attributes = old_e_att
         tmesh.facet_attributes = old_f_att
@@ -909,7 +947,11 @@ def plot_tri(tmesh, phases, seeds, pmesh, plot_files=[], plot_axes=True,
 
     # save plot
     for fname in plot_files:
-        plt.savefig(fname, bbox_inches='tight', pad_inches=0)
+        if n_dim == 3:
+            fig.subplots_adjust(**_misc.plt_3d_adj)
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0.15)
+        else:
+            plt.savefig(fname, bbox_inches='tight', pad_inches=0)
 
     plt.close('all')
 
@@ -943,8 +985,8 @@ def dict_convert(dictionary, filepath='.'):
     Second, if the value of ``dist_type`` is ``histogram``, then the remaining
     key should also be ``filename`` and its value should be the path to a CSV
     file.
-    For the histogram, the first row of this CDF should be the *n+1* bin
-    locations and the second row should be the *n* bin heights.
+    For the histogram, the first row of this CDF should be the *n* bin heights
+    and the second row should be the *n+1* bin locations.
 
     Additionally, if a key in the dictionary contains ``filename`` or
     ``directory`` and the value associated with that key is a relative path,