Periodic microstructures and bug fixes

.github/workflows/python_package.yaml

@@ -8,6 +8,7 @@ jobs:
 
     runs-on: ${{ matrix.os }}
     strategy:
+      fail-fast: false
       matrix:
         python-version: [3.9, '3.10', '3.11', '3.12']
         os: ['ubuntu-latest', 'macos-latest', 'windows-latest']
@@ -30,7 +31,7 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         pip install setuptools wheel
-        pip install flake8 pytest==6.2.5 pytest-cov coveralls
+        pip install flake8 pytest pytest-cov coveralls
     - name: Install package requirements
       run: pip install -r requirements.txt
     - name: Install package

.gitignore

@@ -6,6 +6,10 @@
 test_reports/
 .vscode/
 docs/build*/
+# generated by the documentation build (the examples are run by sphinx-gallery)
+docs/source/auto_examples/
+docs/source/sphinx_gallery/*.csv
+src/microstructpy/examples/*/
 
 ###############################################################################
 #                                                                             #

.readthedocs.yaml

@@ -10,7 +10,7 @@ formats: all
 build:
   os: ubuntu-22.04
   tools:
-    python: "3.8"
+    python: "3.10"
   apt_packages:
     - freeglut3-dev
 
@@ -24,5 +24,6 @@ sphinx:
 python:
     install:
     - requirements: docs/requirements.txt
+    - requirements: requirements.txt
     - method: pip
       path: .

CHANGELOG.rst

@@ -6,6 +6,181 @@ 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`_.
 
+Unreleased
+----------
+Added
+'''''
+- Periodic microstructures, in 2D and 3D: the ``<periodic>`` field of the
+  domain (or the ``periodic`` argument of ``cli.run``, ``SeedList.position``
+  and ``PolyMesh.from_seeds``) selects the periodic axes. Seeds crossing a
+  periodic face are placed without overlapping the opposite side, the
+  Laguerre tessellation is periodic across those faces (cells crossing a
+  face are cut and their pieces tile the domain), and the triangular,
+  tetrahedral and raster meshes have matching nodes on opposite faces
+  with the quality and size settings (``min_angle``, ``max_volume``, the
+  ``max_volume`` of each phase, ``max_edge_length``) acting as on
+  non-periodic meshes: in 2D the cells next to the periodic faces are
+  copied outside the faces while meshing, so that Triangle refines both
+  faces the same way; in 3D the facets are triangulated with the points
+  TetGen adds on them in a first pass, identically on opposite faces, and
+  the mesh is built again with the facets fixed. The pairs of periodic
+  points/nodes and facets
+  are stored in the meshes and their text files; the Abaqus output has a
+  node set per periodic face in matching order; the verification unwraps
+  grains that are split by the faces. The examples ``periodic_2D.xml``,
+  ``periodic_3D.xml`` and ``periodic_tiling.py`` demonstrate periodic
+  microstructures, ``pbx_2D.xml`` and ``pbx_3D.xml`` a periodic
+  particulate composite (crystalline inclusions in a binder), and
+  ``pbx_interface_2D.xml`` and ``pbx_interface_3D.xml`` meshes refined at
+  the grain boundaries. Cells of the same amorphous phase
+  that touch across a periodic face are merged into one region, like cells
+  that share a facet, and the merged region is labelled with the smallest
+  seed number among its cells by every mesher and writer. The element
+  attributes and the facets of periodic meshes are computed from the
+  geometry of the polymesh, since TetGen can leave sub-faces of a facet
+  unmarked when it may not modify the boundary and its region attributes
+  then leak between cells. gmsh is not supported for periodic meshes.
+- ``periodic_margin`` (a setting, and an argument of ``SeedList.position``
+  and ``cli.run``): the minimum distance between the surface of a seed and
+  a periodic face. A seed that ends within the margin of a face, or crosses
+  it by less, is placed elsewhere, since it would leave a thin piece of its
+  grain on the opposite face and elements much smaller than the target size
+  of the mesh there. ``auto`` uses half the target edge length of the mesh,
+  or an eighth of the size of the smallest seed if that is smaller (the
+  smallest seed needs about four elements across it, and cannot satisfy a
+  margin larger than itself).
+- ``PolyMesh.from_seeds(edge_opt=True, periodic_margin=...)``: the edge
+  optimization of periodic meshes treats the thickness of each piece of a
+  cell at a periodic face as a feature like an edge, and moves the seeds of
+  the pieces thinner than the margin (and of their neighbors) normal to the
+  face until the piece reaches the margin or the cell no longer crosses the
+  face. The CLI passes its ``periodic_margin`` setting on. A trial of the
+  optimization is now kept when the shortest feature that it changes gets
+  longer (for the shortest edge of the mesh, the criterion is unchanged),
+  a target that does not improve in ``n_iter`` trials is left alone and the
+  next one is taken, and the seeds moved across a periodic face are wrapped
+  back into the domain. The CLI writes and plots the seeds again after the
+  optimization, so that the seed files match the polygonal mesh. The
+  periodic examples use ``periodic_margin`` ``auto`` and ``edge_opt``.
+  With ``min_angle`` (the minimum angle of the mesh to be built, passed by
+  the CLI from ``mesh_min_angle``), the corners of the cells at the
+  periodic faces narrower than that angle are features too, since the
+  mesher cannot reach the minimum angle there and fills them with shells of
+  very small elements: the seeds on both sides of the facet are moved along
+  it to open the corner.
+- When the nodes on the periodic faces of a 2D mesh do not match after the
+  first pass, the next pass starts from all the points of the mesh (those on
+  a periodic face and on its image merged and put on both faces), so that
+  Triangle only refines it around the merged points, instead of meshing the
+  cells again with the points on the faces only, which split the narrow
+  corners of the cells again at every pass, down to very small elements.
+  The copies of the cells at the corners of the domain, used while meshing,
+  were open on one side and partly discarded by Triangle.
+
+Fixed
+'''''
+- Seed generation is reproducible: the RNG seed chain no longer depends on
+  the (hash-randomized) iteration order of the phase keywords, and
+  ``SeedList.from_info`` and ``cli.run`` no longer modify the ``rng_seeds``
+  and ``filetypes`` arguments (or their mutable defaults).
+- ``<dist_type> cdf </dist_type>`` inputs are no longer distorted when the
+  x-values in the CSV file are not evenly spaced (``density=False`` is now
+  passed to ``scipy.stats.rv_histogram``); ``pdf`` is accepted as an alias
+  of ``histogram``.
+- 3D ``mesh_max_volume`` and per-phase ``max_volume`` are now honored by
+  TetGen; in 2D a per-phase ``max_volume`` larger than the global value is
+  no longer capped, and an infinite ``mesh_max_volume`` is no longer passed
+  to Triangle as the (mis-parsed) switch ``ainf``.
+- ``Ellipsoid.approximate`` mapped the axes incorrectly for the ordering
+  c >= a >= b, so those grains were tessellated with the wrong orientation;
+  the b >= c >= a ordering is now sorted explicitly as well.
+- Seeds read back from ``seeds.txt`` can be repositioned; ellipsoid seeds
+  with a rotation sequence are written in a form that can be read back.
+- Cells that intersect a circular or elliptical domain without having a
+  vertex inside it are no longer dropped, cells cut twice by the boundary
+  are clipped correctly, and the stored areas of clipped cells are correct
+  (``PolyMesh.volumes``, ``verification.volume_fractions``).
+- ``_segment_cross`` no longer hangs for large coordinate values;
+  ``sample_pos_within`` raises instead of looping forever when the position
+  distribution does not cover the domain.
+- ``cli.plot_tri`` no longer hangs in 3D when a void grain touches the
+  boundary of the domain.
+- Relative ``<filename>`` and ``<directory>`` paths inside repeated tags
+  (e.g. several ``<material>`` blocks) are resolved relative to the input
+  file; a top-level ``<include>`` no longer discards materials; values such
+  as ``true_cdf.csv`` no longer cause infinite recursion; ``inf`` is parsed
+  as a float.
+- Verification: ``angle_rad`` inputs are no longer replaced by a uniform
+  distribution, ``<orientation> random </orientation>`` and vector-valued
+  parameters (``side_lengths``, ``axes``) no longer crash, unknown phase
+  fields are ignored, the caller's phases are not modified.
+- ``RasterMesh``: elements are counter-clockwise / right-handed (valid for
+  Abaqus CPS4/C3D8), facets and their attributes are correct, ``vtk`` and
+  ``abaqus`` output work (including with voids), 3D plotting works.
+- ``TriMesh.write``: valid ``.ele``/``.edge``/``.face`` files, Abaqus
+  exterior surface unions reference only defined surfaces, full-precision
+  points in text files, no dangling headers for meshes without attributes.
+- ``PolyMesh.from_seeds(edge_opt=True)`` leaves the seed list in the
+  accepted state (positions and breakdowns consistent) and is quiet unless
+  ``verbose``; ``PolyMesh.write(format='poly')`` writes the file;
+  ``PolyMesh.__eq__`` is silent and no longer cubic.
+- ``Ellipse(axes=...)``, ``Ellipse(matrix=...)``, ``Ellipsoid(c=..,
+  ratio_bc=..)``, ``Square.area_expectation(side_lengths=...)``, the
+  ``*_expectation`` methods with numpy scalars, ``Sphere.plot`` and 3D
+  ``PolyMesh.plot``/``SeedList.plot_breakdown`` on a fresh figure,
+  ``Rectangle.within`` for rotated rectangles, ``reflect`` for ellipses
+  and ellipsoids, single-material ``color_by`` settings, numpy arrays as
+  per-item plot keywords.
+
+Changed
+'''''''
+- pyvoro is installed from the ``pyvoro-rimoli`` package instead of
+  ``pyvoro-mmalahe``. Both provide the same ``pyvoro`` module, but
+  pyvoro-mmalahe bundles Voro++ 0.4.6, whose radical (Laguerre)
+  tessellation can return a cell uncut: the cells then overlap, and
+  Triangle and TetGen can crash on the resulting polygonal mesh.
+  pyvoro-rimoli bundles the current Voro++, where this is fixed, and has
+  wheels for Linux, macOS and Windows. Uninstall pyvoro-mmalahe before
+  upgrading (``pip uninstall pyvoro-mmalahe``), since the two packages
+  install the same files.
+- The continuous integration installs the current pytest and no longer
+  installs tox from the requirements: the pinned tox 3.14 forced an old
+  pluggy that the current pytest-cov cannot load, so no test could run.
+  The jobs of the test matrix no longer cancel each other on a failure.
+- Read the Docs builds the documentation with Python 3.10 instead of 3.8,
+  which pyvoro-rimoli and the current versions of other dependencies do not
+  support, and installs ``requirements.txt`` like the documentation check
+  of the continuous integration.
+- The facets of the triangular and tetrahedral meshes created from a
+  polygonal mesh are sorted (nodes in ascending order within a facet,
+  facets in lexicographic order), whatever the mesher. Triangle and TetGen
+  list the edges/faces of a mesh in an order, and with an orientation, that
+  vary from one run to the next, so the mesh files of otherwise identical
+  runs differed in the order of their facets.
+- ``max_edge_length`` (``mesh_max_edge_length``) acts in 3D on the triangles
+  of the grain boundaries: when it is set, the facets of the polyhedral mesh
+  are triangulated to that edge length (with Triangle, minimum angle 20
+  degrees) before TetGen meshes the cells, for periodic and non-periodic
+  meshes alike, so that the elements can be smaller at the interfaces than
+  inside the grains (see the ``pbx_interface_3D.xml`` example). The
+  geometric tests of a mesh against its polymesh use the non-planarity of
+  the facets (from the snapping of the points to the periodic faces) as
+  their tolerance, and each element is assigned to the cell in which its
+  centroid is deepest.
+- The overlap tolerance fit ``rtol='fit'`` uses the coefficients published
+  in Hart and Rimoli, CMAME 370 (2020) 113242, Eqs. (14) and (15). For very
+  wide size distributions this allows less overlap than before (2D
+  asymptote 0.18 instead of 0.36), so some seeds of high-cv inputs may be
+  rejected during placement.
+- ``Ellipsoid.limits`` is exact for rotated ellipsoids (it was sampled).
+- A ``Seed`` created with a ``position`` (or a geometry with a center) has
+  its breakdown at that position; the geometry center is no longer reset
+  to the origin.
+- ``Ellipse``, ``Ellipsoid`` and ``NBox`` geometries compare equal when
+  their parameters are equal.
+- Unused sampling helpers were removed from ``seeding.seedlist``.
+
+
 `1.5.9`_ - 2023-10-05
 --------------------------
 Added

README.rst

@@ -44,6 +44,9 @@ To install MicroStructPy, download it from PyPI using::
 
 If there is an error with the install, try ``pip install pybind11`` first,
 then install MicroStructPy.
+When upgrading from an earlier version, run ``pip uninstall pyvoro-mmalahe``
+first: MicroStructPy now uses ``pyvoro-rimoli``, which installs the same
+``pyvoro`` module.
 
 
 MicroStructPy can also be installed from source::

docs/source/cli/domain.rst

@@ -1,3 +1,5 @@
+.. _cli_domain:
+
 ====================================
 ``<domain>`` - Microstructure Domain
 ====================================
@@ -298,3 +300,65 @@ Below are some example square domain definitions.
             <center> 5, 0 </center>
         </domain>
     </input>
+
+Periodicity
+^^^^^^^^^^^
+
+Rectangular domains (rectangle, square, box, cube) can be made periodic with
+the ``<periodic>`` field, which lists the periodic axes. Seeds that cross a
+periodic face are placed so that they do not overlap seeds on the opposite
+side, the tessellation is periodic across those faces (a grain that crosses a
+face continues on the opposite side), and the triangular mesh has matching
+nodes on opposite faces. The pairs of periodic nodes are written with the
+meshes, and the Abaqus output contains a node set per periodic face, in
+matching order.
+
+.. code-block:: XML
+
+    <?xml version="1.0" encoding="UTF-8"?>
+    <!-- Example periodic domains -->
+    <input>
+        <domain>
+            <shape> square </shape>
+            <side_length> 10 </side_length>
+            <!-- periodic in both directions -->
+            <periodic> True </periodic>
+        </domain>
+
+        <domain>
+            <shape> rectangle </shape>
+            <side_lengths> 10, 5 </side_lengths>
+            <!-- periodic in x only -->
+            <periodic> x </periodic>
+        </domain>
+
+        <domain>
+            <shape> cube </shape>
+            <side_length> 10 </side_length>
+            <!-- periodic in x and z, free in y -->
+            <periodic> xz </periodic>
+        </domain>
+    </input>
+
+A seed that barely crosses a periodic face, or ends just inside it, leaves a
+thin piece of its grain on the opposite face and very small elements there;
+the ``periodic_margin`` setting rejects such positions and, with
+``edge_opt``, moves the seeds of the cells that still leave such pieces (see
+:ref:`cli_settings`).
+
+The mesher must be Triangle/TetGen for periodic microstructures (with gmsh the
+nodes on opposite faces are not guaranteed to match), and the mesh size of a
+raster mesh must divide the domain length along the periodic axes.
+The mesh quality and size settings (``mesh_min_angle``, ``mesh_max_volume``,
+the ``max_volume`` of each phase and ``mesh_max_edge_length``) apply to
+periodic meshes as to non-periodic ones. In 2D, the cells next to the
+periodic faces are copied outside the faces while meshing, so that Triangle
+refines both faces of a pair the same way; if some nodes on the faces have
+no image, the mesh is built again from all of its points, with the points
+of a periodic face and of its image merged and put on both faces, so that
+Triangle only refines it around those points. In 3D, the
+mesh is built once as usual, then the facets are triangulated with the points
+TetGen added on them (the facets on opposite periodic faces with the points
+of both) and the mesh is built again with the facets fixed; the elements next
+to the periodic faces are slightly more numerous and of slightly lower
+quality than in a non-periodic mesh.