Add new capabilities

- general: switched "active" and "passive" for orientation convention (see
  incompatibility note below), changed convention for HCP directions.
- module -T: added -format svg, added -n ... -morpho lamellae and (the
  equivalent) -morpho lamellar(n=...)", fixed regularization in the case of
  curved domain edges, improved -domain and -loadtess to load CAD files (tess,
  obj, ovm), added tessellation key length(d_x,d_y,d_z).
- module -M: improved meshing of 2D tessellations with curved domain edges,
  made minor bug fixes and improvements.
- module -S: made minor bug fixes.

Author: rquey

VERSIONS

@@ -1,8 +1,25 @@
-New in 4.9.1-6 (16 Sep 2024):
-- module -T: added -n ... -morpho lamellae and (the equivalent)
-  -morpho lamellar(n=...)", made minor fix to regularization of
-  2D cut domain, added tessellation key length(d_x,dd_y,d_z).
-- made minor fixes and improvements.
+New in 4.10.0 (01 Oct 2024):
+- general: switched "active" and "passive" for orientation convention (see
+  incompatibility note below), changed convention for HCP directions.
+- module -T: added -format svg, added -n ... -morpho lamellae and (the
+  equivalent) -morpho lamellar(n=...)", fixed regularization in the case of
+  curved domain edges, improved -domain and -loadtess to load CAD files (tess,
+  obj, ovm), added tessellation key length(d_x,d_y,d_z).
+- module -M: improved meshing of 2D tessellations with curved domain edges,
+  made minor bug fixes and improvements.
+- module -S: made minor bug fixes.
+
+* Incompatible changes in -T:
+  - switched "active" and "passive" for orientation convention.  The default
+    convention is now "passive".  The orientation data themselves remain the
+    same.  Bumped version numbers to 3.5 for tess, 2.2 for tesr, 2.3 for msh
+    and 1.1 for sim.  When an older tess, tesr, msh or sim is parsed, the
+    defined convention is changed from active to passive and vice versa.  This
+    ensures no change in Neper's own behavior.  Postprocessors that are
+    sensitive to the change of convention should be updated accordingly.
+  - changed convention for HCP directions (see the documentation).
+  - changed -domain "cell(<tess_file>,<cell_id>)" into -domain
+    "file(<tess_file_>),cell(<cell_id>)".
 
 New in 4.9.0 (15 May 2024):
 - module -T: added -per 1 for -morpho cube, added -ori ks (Kurdjumov-Sachs)

doc-dev/index.rst

@@ -13,6 +13,7 @@ Neper's Developer Documentation
    smoothing.rst
    1dlaguerre.rst
    ellipse.rst
+   odfsampling.rst
    gpl.rst
    fdl.rst
 

doc-dev/odfsampling.rst

@@ -0,0 +1,117 @@
+.. _odfsampling:
+
+Random sampling of ODFs
+=======================
+
+Given an ODF defined on a mesh, either on elements (piecewise constant) or on nodes (piecewise linear), how can we sample it randomly?
+
+Case of a uniform ODF
+---------------------
+
+If the ODF is uniform (no texture), then we can choose the orientations without taking it into account, from random numbers (:math:`n_1`, :math:`n_2`, :math:`n_3`, :math:`n_4 \in [0,\,1]`):
+
+- Euler-Bunge angles: :math:`\varphi_1=2\,\pi\,n_1`, :math:`\varphi_2= \hbox{acos} (2\,n_2-1)`, :math:`\varphi_3=2\,\pi\,n_3`
+- unit quaternion: if :math:`{{n_1}^2+{n_2}^2+{n_3}^2+{n_4}^2} \leq 1` then accept orientation :math:`q_i=n_i`
+- homochoric vector:  if :math:`\sqrt{{n_1}^2+{n_2}^2+{n_3}^2} \leq 1` then accept orientation :math:`x_i=n_i`
+
+Case of a non-uniform ODF
+-------------------------
+
+.. note::
+
+  Reminder on the expressions of the volume element, :math:`dg`:
+
+  1. Rodrigues space:
+
+    - :math:`dg = \left(\frac{\rho}{1+\rho^2}\right)^2 \, d\rho \sin{\chi} \, d\chi \, d\zeta` (polar coordinates)
+    - :math:`dg = \left(\frac{1}{1+\rho^2}\right)^2 \, dr_1 \, dr_2 \, dr_3` (rectangular coordinates)
+
+    - By definition, :math:`\int_g dg = \pi^2` (with :math:`\rho\in\left[0,\,+\infty\right[`, :math:`\chi\in\left[0,\,\pi\right]` and :math:`\zeta\in\left[0,\,2\,\pi\right]` or :math:`\left\{r_1,\, r_2,\, r_3\right\} \in \mathbb{R}^3`).  If integration is done on a fundamental region, then the integral is equal to :math:`\pi^2/n_c`, where :math:`n_c` is the multiplicity (24 for cubic, etc.).
+
+  2. Euler space (Bunge convention):
+
+    - :math:`dg = \sin{\phi} \, d\varphi_1 \, d\phi \, d\varphi_2`
+
+    - By definition, :math:`\int_{\varphi_1=0}^{2\,\pi}  \int_{\phi=0}^{\pi/2} \int_{\varphi_2=0}^{\pi/2} dg = \pi^2`
+
+    - :math:`\sqrt{\hbox{det}(g)} = \frac{1}{8\,\pi^2} \sin{\phi}`, where :math:`g` is the metric tensor
+
+  3. Homochoric space:
+
+ The integral of the ODF, :math:`f(g)`, over full orientation space is therefore equal to :math:`\pi^2`.  The integral over the cubic fundamental region is equal to :math:`\pi^2/24`.
+
+If the ODF is non-uniform, the frequency of orientations must adhere to the intensity of the ODF (by definition).  One may distinguish random sampling and uniform sampling, however.  Uniform sampling (just as we implemented it for the case of a uniform ODF) attempts to fit the ODF at best given the number of orientations available.  Random sampling (just as we apply it to the case of a uniform ODF) follows the ODF from a statistical point of view, but does not attempt to fit it at best; it retains a random character.  This is the case developed here.
+
+Let the ODF be defined as
+
+.. math::
+
+  \begin{equation}
+    f(\boldsymbol{x}) = \frac{1}{V} \, \frac{dV(\boldsymbol{x})}{d\boldsymbol{x}},
+  \end{equation}
+
+where :math:`\boldsymbol{x}` is the orientation.
+
+Rejection sampling (inefficient)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Let :math:`f_\text{max}` be the maximum of :math:`f`.  The procedure is then the following:
+
+  1. Pick a random orientation, :math:`\boldsymbol{x}`, following a uniform ODF;
+  2. Generate a random number :math:`t\in[0,\,1]`;
+  3. if :math:`t < f(\boldsymbol{x}) / f_\text{max}`, then accept orientation.
+
+Terminate when all orientations have been generated.
+
+The efficiency of the method scales with :math:`1/f_\text{max}`; that is, the stronger the texture is, the more orientations need to be generated and tested to attain the desired number of orientations.  One may also note that the evaluation of :math:`f` at a particular orientation may be costly (for large meshes), which makes the process relatively inefficient.
+
+Inverse transform sampling (100% efficient)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. note:: This is the method implemented in Neper.
+
+The ODF is defined over (multidimensional) orientation space; however, the discretization of orientation space (its mesh elements) provides a unidimensional ODF function that depends on the orientations (the numbers of the mesh elements), :math:`f_u (x)`, where :math:`x` is an integer than varies from 1 to :math:`N` (the total number of elemental orientations).  (Note that the respective positions if the elements have no influence.) The "frequency of occurrence" of an orientation represented by element :math:`x`, which can be seen as a *probability*, :math:`p_x`, is given by:
+
+.. math::
+
+   \begin{equation}
+     p_x = \frac{n_c}{\pi^2} \, v_x \, f_u(x)
+   \end{equation}
+
+where :math:`v_x` is the volume of orientation space that orientation :math:`x` represents (taking into account the orientation multiplicity), and :math:`f_u(x)` its ODF value, and :math:`n_c` is the multiplicity of the space (24 for cubic).
+
+:math:`v_x` is computed by integration of :math:`dg` over the associated element.  In practice, it is taken as the volume of the element multiplied by :math:`(1 / (1 + \rho^2))^2`, in Rodrigues space (:math:`\rho = \sqrt{{r_1}^2 + {r_2}^2 + {r_3}^2}` assumed constant). So, :math:`p_x` sums to 1.
+
+Formally, the inverse transform sampling is based on the fact that, for any random variable :math:`X \in \mathbb{R}`, the random variable :math:`F_X^{-1}(U)` has the same probability density function as :math:`X`, where :math:`F_X^{-1}` is the generalized inverse of the cumulative distribution function :math:`F_X` of :math:`X` and :math:`U` is uniform on :math:`[0,\,1]`.
+
+In the present case, discrete (random) variables are discussed. The cumulative distribution function can be expressed as
+
+.. math::
+
+   \begin{equation}
+     P(x) = \sum_{j\leq x} p_j
+   \end{equation}
+
+which is a step function.  As a step function does not have an inverse, the concept of generalized inverse is of particular importance.
+
+The generalized inverse of :math:`P` is:
+
+.. math::
+
+   \begin{equation}
+     P^-(y) = \hbox{inf}\left\{x \in \mathbb{R} : P(x) \geq y\right\}, y \in \mathbb{R}
+   \end{equation}
+
+Note that entire ranges of values of :math:`y` now correspond to a single value of :math:`x`.
+
+The procedure then is to draw a random number between 0 and 1, :math:`y`, and compute the value of :math:`P^-(y)`.  This value is the element from which we generate an orientation.
+
+To complete the process, for any selected element, :math:`y`, the orientation is chosen randomly within the element.  To do so, we pick a random position inside the element, following a uniform space distribution, using Rocchini and Cignoni's method (https://doi.org/10.1080/10867651.2000.10487528), which provides the orientation.
+
+This method is direct and does not involve any testing or rejection.  It is *almost* correct in the sense that steps 1 and 2 are correct, but step 3 assumes a uniform density of the space inside an element.  This is reasonable assumption for properly conditioned spaces (such as Rodrigues space) and properly refined meshes, and is the same assumption as the one made to discretize the ODF itself.
+
+Useful references:
+  - https://link.springer.com/article/10.1007/s00186-013-0436-7
+  - Morawiec's book for the rejection method
+  - L13-Grain_Bndries_RFspace-15Mar16.pdf for the volume element
+  - find a book in statistics for rejection and inverse sampling methods, which are standard

doc/conf.py

@@ -11,8 +11,8 @@
 import sphinx_rtd_theme
 
 project = u'Neper'
-version = u'4.9.1-6'
-release = u'4.9.1-6'
+version = u'4.10.0'
+release = u'4.10.0'
 author = u'Romain Quey'
 copyright = u'Romain Quey'
 language = 'en'

doc/exprskeys.rst

@@ -640,10 +640,36 @@ Key                        Descriptor                               Number of pa
 :data:`quaternion`         quaternion                               4
 ========================== ======================================== =============================
 
-The convention can be added to the descriptor, either :data:`active` or :data:`passive`, as in :data:`rodrigues:active`.  When no convention is provided, :data:`active` is assumed.
+The convention can be added to the descriptor, either :data:`active` or :data:`passive`, as in :data:`rodrigues:active`.  When no convention is provided, :data:`passive` is assumed.
 
 Some options can take parameter values as argument, in which case the orientation must be expressed as :data:`<descriptor>(<parameter1>,<parameters2>,...)`. An example is :data:`rodrigues(0.1,0.2,0.3)`.
 
+.. _orientation_convention:
+
+Orientation Convention
+~~~~~~~~~~~~~~~~~~~~~~
+
+The crystal coordinate systems are attached to the crystal lattice as illustrated below, in the case of cubic and hexagonal symmetries:
+
+.. image:: imgs/directionconvention.png
+
+The so-called "passive" orientation convention is used by default, which is based on the rotation of the sample coordinate system into the crystal coordinate system.  Under this convention, the values of all :ref:`rotation_and_orientation_descriptors` are provided below for a simple (but representative) configuration that corresponds to a rotation of 30° about :math:`X_s`:
+
+.. image:: imgs/orientationconvention.png
+
+======================================== =============================
+**Descriptor**                           **Value** (angles in degrees)
+Rodrigues vector                         :math:`(0.267949192,\,0,\,0)`
+Euler angles (Bunge convention)          :math:`(0,\,30,\,0)`
+Euler angles (Kocks convention)          :math:`(270,\,30,\,90)`
+Euler angles (Roe convention)            :math:`(270,\,30,\,90)`
+rotation matrix                          :math:`\left(\begin{array}{ccc}1 & 0 & 0 \\ 0 & 0.866025404 & 0.5 \\ 0 & -0.5 & 0.866025404\\\end{array}\right)`
+rotation axis / angle pair               :math:`(1,\,0,\,0) / 30`
+quaternion                               :math:`(0.965925826,\,0.258819045,\,0,\,0)`
+======================================== =============================
+
+The values of the orientation descriptors under the "active" convention are obtained by taking the opposite rotation.
+
 .. _ideal_orientations:
 
 Ideal Orientations
@@ -667,7 +693,7 @@ Keys are available for ideal orientations (lowercased is accepted):
 :data:`Copper2`          :math:`(\overline{1}\,1\,2)[1\,\overline{1}\,1]`
 ======================== ====================================
 
-When loading orientations from an external file, use :data:`file(<file_name>[,des=<descriptor>])` where the orientation descriptor is among those listed above and is :data:`rodrigues:active` by default.
+When loading orientations from an external file, use :data:`file(<file_name>[,des=<descriptor>])` where the orientation descriptor is among those listed above and is :data:`rodrigues:passive` by default.
 
 .. _orientation_fibers:
 

doc/fileformat.rst

@@ -78,8 +78,10 @@ and
     *edge
      <total_number_of_dom_edges>
      <dom_edge_id> <number_of_dom_vertices> [<dom_ver_1> <dom_ver_2>]
+                   <dom_edge_type>
+                   <number_of_params> <dom_edge_param1> <dom_edge_param2> ...
                    <dom_edge_label>
-                   <number_of_dom_tess_edges> <edge_1> <edge_2> ...
+                   <number_of_dom_tess_edges> <dom_tess_edge_1> <dom_tess_edge_2> ...
      ...
     *face
      <total_number_of_dom_faces>
@@ -88,8 +90,7 @@ and
                    <dom_face_type>
                    <number_of_params> <dom_face_param1> <dom_face_param2> ...
                    <dom_face_label>
-                   <number_of_dom_tess_faces>
-                   <dom_tess_face_1> <dom_tess_face_2> ...
+                   <number_of_dom_tess_faces> <dom_tess_face_1> <dom_tess_face_2> ...
      ...
    **periodic
     *general
@@ -127,7 +128,7 @@ where (with identifiers being integer numbers),
 
 - :data:`**format` denotes the beginning of the format field.
 
-- :data:`<format>` is the file format, currently `2.0` (character string).
+- :data:`<format>` is the file format, currently `3.5` (character string).
 
 - :data:`**general` denotes the beginning of the general information field.
 
@@ -257,6 +258,12 @@ where (with identifiers being integer numbers),
 
 - :data:`<dom_ver_#>` are identifiers of the domain vertices of a domain edge or face.
 
+- :data:`<dom_edge_type>` is the type of an edge.
+
+- :data:`<number_of_params>` is the number of parameters.
+
+- :data:`<dom_edge_param#>` are the parameters of a domain edge.
+
 - :data:`<dom_edge_label>` is the label of a domain edge, formatted as :data:`x0y0`, :data:`x0y1`, :data:`x1y0`, ... (for a cubic domain).
 
 - :data:`<number_of_dom_tess_edges>` is the number of tessellation edges of a domain edge.
@@ -273,8 +280,6 @@ where (with identifiers being integer numbers),
 
 - :data:`<dom_face_type>` is the type of a face, among `plane`, `cylinder` or `sphere`.
 
-- :data:`<number_of_params>` is the number of parameters of a domain face.
-
 - :data:`<dom_face_param#>` are the parameters of a domain face.  For a planar face, they are the parameters of the equation of the face, listed in the order :data:`<face_eq_d>`, :data:`<face_eq_a>`, :data:`<face_eq_b>` and :data:`<face_eq_c>`.  For a cylindrical face, they are the coordinates of the base point, the axis and the radius.  For a spherical face, they are the coordinates of the center and the radius.
 
 - :data:`<dom_face_label>` is the label of a domain face, formatted as :data:`x0`, :data:`x1`, :data:`y0`, :data:`y1`, :data:`z0` or :data:`z1` (for a cubic domain).  For a cylindrical domain, it is formatted as :data:`z0`, :data:`z1`, :data:`f1`, :data:`f2`, ... Otherwise, it is one of :data:`f1`, :data:`f2`, ...  Edge and vertex labels are derived accordingly.
@@ -401,7 +406,7 @@ where
 
 - :data:`**format` denotes the beginning of the format field.
 
-- :data:`<format>` is the file format, currently `2.0` (character string).
+- :data:`<format>` is the file format, currently `2.2` (character string).
 
 - :data:`**general` denotes the beginning of the general information field.
 
@@ -572,7 +577,7 @@ where
 
 - :data:`$MeshVersion` denotes the beginning of a mesh version field.
 
-- :data:`<mesh_version>` is the mesh file version (currently :data:`2.2.3`).
+- :data:`<mesh_version>` is the mesh file version (currently :data:`2.3`).
 
 - :data:`$EndMeshVersion` denotes the end of a mesh version field.
 
@@ -769,7 +774,7 @@ where
 
 Results can have integer values, real values, vectorial values or tensorial values. In the result files, values for the different entities (nodes, elements, etc.) are written on successive lines, with components written on successive columns (space delimited). The components of a vector, :data:`v`, are written as :data:`v1` :data:`v2` :data:`v3`. The components of a symmetrical tensor, :data:`t`, are written using Voigt notation, as :data:`t11` :data:`t22` :data:`t33` :data:`t23` :data:`t31` :data:`t12`. The components of a skwe-symmetrical tensor, :data:`t`, are written using :data:`t12` :data:`t13` :data:`t23`. The components of a non-symmetrical tensor, :data:`t`, are written as :data:`t11` :data:`t12` :data:`t13` :data:`t21` :data:`t22` :data:`t23` :data:`t31` :data:`t32` :data:`t33`.
 
-The directory also contains a hidden file, :file:`.sim`, containing information on the simulation and the content of the simulation directory.  This file is only for internal use and is formatted as follows:
+The directory also contains a hidden file, :file:`.sim`, containing information on the simulation and the content of the simulation directory.  This file is only for internal use and is formatted as follows (the sections depend on the actual content of the simulation directory and most of them are optional):
 
 .. code-block:: plain
 
@@ -785,6 +790,8 @@ The directory also contains a hidden file, :file:`.sim`, containing information
      <msh_file>
     *ori
      <ori_file>
+    *opt
+     <opt_file>
     *bcs
      <bcs_file>
     *phase

doc/imgs/cubicconvention.asy

@@ -0,0 +1,15 @@
+settings.outformat="pdf";
+
+unitsize(3cm);
+
+draw(shift((0.15,0.15))*scale(1.8,1.4)*polygon(4),invisible);
+
+draw(scale(0.8)*polygon(4));
+draw(circle((0,0),0.025), linewidth(0.8));
+
+// First graph
+draw((0,0)--(1,0), Arrow);
+draw((0,0)--(0,1), Arrow);
+label("$X_c$", (1,0), N);
+label("$Y_c$", (0,1), E);
+label("$Z_c$", (0,0), SW);

doc/imgs/directionconvention.png

@@ -0,0 +1,15 @@
+settings.outformat="pdf";
+
+unitsize(3cm);
+
+draw(shift((0.15,0.15))*scale(1.8,1.4)*polygon(4),invisible);
+
+draw(scale(0.8)*polygon(6));
+draw(circle((0,0),0.025), linewidth(0.8));
+
+// First graph
+draw((0,0)--(1,0), Arrow);
+draw((0,0)--(0,1), Arrow);
+label("$X_c$", (1,0), N);
+label("$Y_c$", (0,1), E);
+label("$Z_c$", (0,0), SW);