Welcome to the documentation!

Global intro

CristalX

Identification of individual grains in microscopic images

CristalX is a Python package that helps in the analysis of polycrystalline microstructures. Its name originates from the French word ‘cristal’, corresponding to the English word ‘crystal’.

Features

• Image segmentation to identify the grains in a microstructure

• Analysis tools for the segmented image

• Explicit geometrical representation of the grains

• Interacting with meshes created on the microstructure

• Mapping fields between a mesh and the grid of DIC measurements

• Simulation tools for the inverse problem arising from a combined numerical-experimental method (in progress …)

• Visualization and data exchange

Getting help

1. Read the documentation.

2. Check the existing issues. They may already provide an answer to you question.

3. Open a new issue.

Contributing

Read the docs/source/contributing.md file.

Citing CristalX

We have an article freely available on SoftwareX, showing the background and the design of CristalX.

When using CristalX in scientific publications, please cite the following paper:

• Csati, Z.; Witz, J.-F.; Magnier, V.; Bartali, A. E.; Limodin, N. & Najjar, D. CristalX: Facilitating simulations for experimentally obtained grain-based microstructures. SoftwareX, 2021, 14, 100669

BibTeX entry:

@Article{Csati2021,
  author    = {Zoltan Csati and Jean-Fran{\c{c}}ois Witz and Vincent Magnier and Ahmed El Bartali and Nathalie Limodin and Denis Najjar},
  journal   = {{SoftwareX}},
  title     = {{CristalX}: {F}acilitating simulations for experimentally obtained grain-based microstructures},
  year      = {2021},
  month     = jun,
  pages     = {100669},
  volume    = {14},
  doi       = {10.1016/j.softx.2021.100669},
}

Installation

CristalX is easy to install; you can get started in minutes.

Dependencies

CristalX is a written in pure Python, although it relies on packages that use other languages (mostly C and C++). However, it is not a problem from the user’s perspective as no manual compilation is needed. Those who are interested can see the dependencies in the environment.yml file in the root directory.

Whether you install CristalX by conda or try it online, all but one dependency is installed. In a common workflow, that single missing dependency will not affect you. For details, see …

Try it online

If you just want to get a taste of CristalX, you can try it online without installing anything.

Click here to jump to the root directory. Existing Jupyter notebooks are in the notebooks/ directory. You can modify them and create new ones. If you directly want to open the main tutorial, which describes a real-world application, click here.

Install it locally

If you want more control (debugging, inspecting variables, etc.) over CristalX or if you wish to contribute to the project, it is recommended to install it on your machine.

Obtain the source

If you are a user of CristalX, the best is to get the latest release:

• download it from GitHub

• or clone it with Git:

  git clone https://github.com/CsatiZoltan/CristalX.git
  cd CristalX
  git checkout v<version_number>
  

  where <version_number> is the version you want to use. E.g. if you want to use version 1.0.1, you need to type git checkout v1.0.1. See the available tags, corresponding to the published releases, for the possibilities.

  Note that in this case, you will be in a “detached HEAD” state, meaning that the HEAD does not point to a branch but to the specific tag. Any commit you make in this state will not be associated with a branch. Therefore, if you want to develop or contribute to CristalX, check out the master branch (see the next paragraph).

If you want to develop CristalX or simply want to have access to the latest features, you need to fetch the latest state:

• download it from GitHub

• or clone it with Git:

  git clone https://github.com/CsatiZoltan/CristalX.git
  

Install with conda

All you need to have is the conda package manager. Open a terminal (or a command prompt if you are under Windows) in the root directory and type

conda env create -f environment.yml

CristalX has been installed to a separate environment, so you can safely work inside it. Activate that environment:

conda activate CristalX

Once you have finished working with CristalX, either close the terminal or type

conda deactivate

to return to your default conda environment.

If you want to uninstall CristalX, make sure that the CristalX conda environment is not active and then

conda env remove -n CristalX

When uninstalled, conda env list will not show it.

Install with pip

If you do not have conda or you prefer pip, you can also install CristalX by typing

pip install -r requirements.txt

assuming that you have Python installed. It is highly recommended that you first create a new virtual environment to make sure you do not break your Python installation. However, an important component of CristalX (functions that rely on PythonOCC) will not be installed if you choose pip. The reason for this is that PythonOCC is not (yet) available on PyPI.

Documentation

The documentation for CristalX is available both locally and online. It is generated by Sphinx and consists of the code documentation and the various guides (including the one you are currently reading). If you want to extend the documentation, please read the corresponding section in the Developer guide.

Local documentation

It resides in the docs/ directory. Open a terminal there and type make html. The output will be written to docs/build/. Open docs/build/html/index.html to land on the home page of the documentation.

Online documentation

For every push event, the documentation is rebuilt by Read the Docs. The latest HTML version is available on https://cristalx.readthedocs.io/en/latest. Read the Docs is configured to build a PDF and EPUB output as well. All three formats can be downloaded for offline use if you click on the bottom of the left sidebar. For convenience, here are the direct links: HTML, PDF, EPUB.

Using the HTML documentation

For browsing on a computer, the HTML documentation is the most convenient and feature-rich.

The side-bar on the left provides quick navigation.

Code example boxes have two toggle buttons. Clicking on the prompt button >>> hides the prompt and the output. If you then click on the copy icon next to it, you can copy the contents of the box, making it easy to try the code chunks in your Python environment.

A direct link can be obtained for most sections, including the code documentation, by moving the cursor to the right of the text and then clicking on the chain icon.

When browsing the API documentation and want to look at the implementation, you can easily access the source code by clicking on the [source] tag.

How to use the codes

If you want to get to know the classes and functions, the best is to check the API documentation. Most components are amply documented and contain one or more examples. If you want to see them in action, check out our workflow for a real-world application. If you are curious about the underlying algorithms, consult with the Algorithms section.

Tutorials

The tutorials are available in two forms: Python scripts and Jupyter notebooks.

Python scripts

The scripts are intended to be run in batch mode. They need to be run in a specific order as the subsequent scripts rely on the output of the previous ones. These scripts are useful if you are already a bit familiar with CristalX or if you wish to adapt the scripts to your needs. Indeed, if you intend to extend the package, the safest way is to copy the relevant scripts and modify them. The scripts are named as run_moduleName, where moduleName is the name of the module in the grains/ directory that the script mainly relies on. In some way, run_moduleName acts as a demonstration of what the moduleName module is used for. The scripts are meant to be run as modules, i.e. navigate to the root of the project and type

python -m scripts.moduleName

Jupyter notebooks

Almost the same workflow is available as a single Jupyter notebook, located at the nootebooks/ directory. this allows you to interactively discover CristalX through an example. You can easily rerun parts of the code to see the effects of the parameters, and the rich output is embedded into the same document. If you are a novice Python user or just want to have a taste of CristalX, notebooks are the recommended way to get started.

Undo changes

When you experiment with CristalX, you will probably change parameter values. As the scripts communicate by reading and writing data, the original data that come with CristalX will be overwritten. There are multiple ways to undo the changes.

In Binder

As written in the Installation instructions, you can try CristalX online without the need to install anything on your computer. Then a separate virtual environment is created. Whatever changes you make there, they will not influence your local installation (if you have) or the data on the GitHub repository.

In a local installation

If you downloaded CristalX from GitHub, you can simply replace the new files with the original ones. In case you cloned CristalX with Git, you can easily discard the changes you make.

A detailed workflow

CristalX contains several modules that facilitate experimental and numerical works on polycrystalline microstructures. The usage of these modules is demonstrated on a complex example.

Problem formulation

Grain-based microstructures occur in nature but also develop during industrial processes. We are interested in how the size and the distribution of grains, and the material they are made up of influence the resistance of train wheels and axles to fatigue loading. The problem requires that

• the individual grains can be distinguished on a microscopic image

• a good quality mesh is generated on the microstructure

• the numerical solution obtained on this mesh can be compared to full-field measurements

The first two steps are discussed here.

Initial setup

First, we switch to the root directory of the project so that the package imports are consistent. Note that the relative path always works because the current directory in case of a Jupyter notebook is always the directory in which the notebook is.

import os
os.chdir('../')
main_dir = os.path.abspath(os.curdir)

Let us set the directory that will hold our data.

data_dir = os.path.abspath('data')
try:
    os.mkdir(data_dir)
except FileExistsError:
    pass

We now copy the necessary data files that come with CristalX into this freshly created directory.

from os.path import join
from shutil import copy
selected_files = ['1.png', '1_cropped.png', 'splinegons.png', '1_mesh_extended.npz']
for file in selected_files:
    copy(join('scripts/data', file), data_dir)

In interactive mode, debug information are displayed. They are not useful for presentation purposes as they pollute the rendered output. Similarly for runtime warnings.

import logging
logging.basicConfig(level=logging.ERROR)
import warnings
warnings.simplefilter('ignore')

For the figures to be visible in this notebook, we define a utility function to scale them.

import matplotlib.pyplot as plt
def scale_figure(dpi):
    plt.gcf().set_dpi(dpi)

Identify grains in a microstructure

The picture below was taken with a digital camera and shows the central part of a tensile specimen. In this section, we will identify the individual grains in the image.

from skimage.io import imshow
imshow(join(data_dir, '1.png'))
scale_figure(150)

So that the surrounding region does not interfere with the segmentation, we will work on a cropped region.

from skimage.io import imread
import numpy as np
image = join(data_dir, '1_cropped.png')

As you can see, the RGB image is represented as a 3D numpy array.

image_as_matrix = imread(image)
type(image_as_matrix), image_as_matrix.shape

(numpy.ndarray, (638, 589, 3))

However, you do not explicitly have to deal with this matrix as the Segmentation class will handle it.

The Segmentation class offers a set of methods to perform the consecutive steps of the segmentation workflow. We pass to the constructor the path to the image to be segmented. By default, the resulting images are shown at the end of each step. In batch mode, you probably want to set the interactive_mode optional parameter to False. But in this interactive document, we want to see the outputs.

from grains.segmentation import Segmentation
GS = Segmentation(image)

Image successfully loaded.

After the image has been loaded and stored in the GS object, we want to remove the noise. Since the interfaces among the grains are important for us, we use median filtering for smoothing as it preserves the contours. The smoothing is governed by the window size. The larger this number, the smoother the processed image is, and the higher the computational cost becomes.

filter_window_size = 5
filtered = GS.filter_image(filter_window_size)

Median filtering finished.

To improve the segmentation result later, we perform superpixel segmentation as an initial step. During this process, the image is subdivided into superpixels, groups of pixels that belong together based on colour, spatial distance or other properties. On this image, a superpixel algorithm called Quick Shift proved to be the best among the algorithms available in scikit-image.

segmented = GS.initial_segmentation(filtered)

Quick shift segmentation finished. Number of segments: 421

The result is indeed a segmented (also called labelled) image. Distinct positive integers, called labels, are associated to groups of pixels in the labelled image and each pixel belongs to one and only one such group.

segmented

array([[  8,   8,   8, ...,   2,   2,   2],
       [  8,   8,   8, ...,   2,   2,   2],
       [  8,   8,   8, ...,   2,   2,   2],
       ...,
       [400, 400, 400, ..., 414, 414, 414],
       [400, 400, 400, ..., 414, 414, 414],
       [400, 400, 400, ..., 414, 414, 414]])

The segmented image has the same size as the original image, but it is not an RGB image any more.

segmented.shape, image_as_matrix.shape

((638, 589), (638, 589, 3))

We can visualize a labelled image by associating different colours to the different labels.

from skimage.color import label2rgb
from numpy.random import random
nlabel = len(np.unique(segmented))
imshow(label2rgb(segmented, colors=random((nlabel, 3))))
scale_figure(150)

As the superpixel segmentation in the previous step resulted in an oversegmented image, the region adjacency graph is constructed and used to merge some of the neighbouring superpixels based on their similarity with respect to mean colour. Regions connected by edges with smaller weights than a prescribed threshold are combined. As you can see, the number of supersegments decreases.

cluster_merging_threshold = 7
reduced = GS.merge_clusters(segmented, threshold=cluster_merging_threshold)

Tiny clusters merged. Number of segments: 177

We need a few more steps before succesfully applying another segmentation technique to obtain the final segmented image. First, we detect the grain boundaries. They are shown superposed on the original image.

boundary = GS.find_grain_boundaries(reduced)

Grain boundaries found.

The boundary image is a binary image in which the True values indicate the boundaries among the labelled regions.

boundary

array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])

The we use thinning on the grain boundary to obtain a single-pixel wide skeleton.

skeleton = GS.create_skeleton(boundary)

Skeleton constructed.

If the automatic segmentation carried out so far is not good enough, the user can manually edit the skeleton as a graph in ImagePy. The combination of the automatic segmentation with human supervision is a powerful way to achieve good results in a relatively short amount of time. Here, we stay with the automatic method.

To recover the segments again, we use watershed segmentation on the skeleton. Had it been directly applied on the original image, the result of the watershed segmentation would have been an oversegmented image due to the noises that act as local minima. The success of the watershed segmentation depends on how well the catchment basins are identified, which are the locations where the flooding starts. The so-called marker-based watershed segmentation methods rely on markers (computed automatically or given by the user), i.e. the location of the catchment basins, as inputs. In the watershed_segmentation method of the Segmentation class, various transformations are used to obtain a desirable outcome.

watershed = GS.watershed_segmentation(skeleton)

Watershed segmentation finished. Number of segments: 247

To the left and to the right of the central zone, artificial “grains” appear. This is not a problem, they can either be merged, or left as it is and associate the same material to them.

Extending the image domain

The segmented image corresponds to the recrystallized part of the specimen. When we perform the traction test numerically we need a larger domain, the boundary on which the boundary conditions are applied. The size of this extended domain does not have to match with the physical size of the tensile specimen, but it has to be sufficiently large so that the far-field boundary conditions do no not influence the deformation state in the central region.

As a first step, we extend the domain by padding the corresponding numpy array. Since we want to handle the added domains (grains) on the left and on the right as separate grains, we associate different labels to them. There are about 250 grains in the central region, so choosing labels 500 and 501 ensures that these regions have unique labels.

The extended domain is constructed in such a way that it matches the region of the DIC measurements. This will be useful when comparing the results of the simulation with the experimental data.

from grains.simulation import change_domain
extended_image = change_domain(watershed, 0.4395, 0, 0, 0, 500)
extended_image = change_domain(extended_image, 0, 0.5568, 0, 0, 501)
imshow(extended_image)
scale_figure(150)

From image to geometry

In this section, we show how to represent the grains as an explicit geometry. See our paper for the motivation.

For grain representation, we can choose between polygon and splinegon approximations. We will base the meshing on the splinegon representation, but the polygons are also shown. The high-level functions polygonize and splinegonize wrap all the important algorithms. You can see from the code below that they share similar function signature. The only difference is that the spline parameters can be provided to the splinegonize function.

from grains.cad import polygonize, search_neighbor, splinegonize
polygons = polygonize(extended_image, search_neighbor(2, np.inf), connectivity=1)
splinegons, _ = splinegonize(extended_image, search_neighbor(2, np.inf), connectivity=1, degree_min=3, degree_max=3, continuity='C0', tol=1)

Let us plot the grains. For polygons:

from grains.cad import plot_polygons
plot_polygons(list(polygons.values()));
scale_figure(150)

Plotting the spline surfaces will bring up a window.

from grains.cad import plot_splinegons
from OCC.Display.SimpleGui import init_display
plot_splinegons(list(splinegons.values()), color=(0, 0, 1))

Layer manager created
Layer dimensions: 1024, 768

To allow reading this document statically, here is a screenshot of the window that popped up.

fig, ax = plt.subplots(dpi=500)
imshow(join(data_dir, 'splinegons.png'), ax=ax)
ax.set_axis_off()

Once we have the spline surfaces, we write them to a STEP file. This allows us to edit the geometry in a CAD program and to generate a mesh for the grains.

from grains.cad import regions2step
regions2step(list(splinegons.values()), join(data_dir, 'microstructure.stp'))

Repairing the geometry

For the sample microstructure, one grain was not identified, i.e. it contains a hole. The hole can be filled in by creating a new grain. We have a detailed guide on how to do this.

Mesh generation

The same guide linked in the previous section continues with the mesh generation. It happens that we constructed the conforming mesh in Salome, and then exported it to a .med file. Using the med module of CristalX, the triangular mesh cells for each grain and the boundary nodes have been extracted. We also have a tutorial on how to process a .med file.

Modifying the mesh

First, we obtain the mesh data that was saved using the med module.

mesh_file = join(data_dir, '1_mesh_extended.npz')
with np.load(mesh_file, allow_pickle=True) as mesh:
    nodes = mesh['nodes']
    elements = mesh['elements']
    element_groups = mesh['element_groups']
    node_groups = mesh['node_groups']
# Retrieve groups, which were stored in dictionaries (https://stackoverflow.com/a/40220343/4892892)
element_groups = element_groups.item()
node_groups = node_groups.item()

The array elements contains the label of the nodes in the mesh.

elements

array([[  603,   604,  4311],
       [  604,    19,   605],
       [  606,    18,  4312],
       ...,
       [13377, 13369, 13374],
       [13372, 13377, 13373],
       [13373, 13370, 13372]])

The array nodes holds the coordinates of each node of the mesh.

nodes

array([[2.50000000e-01, 2.59000000e+02],
       [8.00000000e+01, 2.59250000e+02],
       [2.31000000e+02, 2.59250000e+02],
       ...,
       [6.28062814e+02, 8.02235833e+02],
       [6.22855688e+02, 7.95968504e+02],
       [6.23756131e+02, 7.89612749e+02]])

The elements are available for each grain. For instance, grain 87 is labelled as Face_87 and it contains the following elements.

element_groups['Face_87']

array([23231, 23232, 23233, 23234, 23235, 23236, 23237, 23238, 23239,
       23240, 23241, 23242, 23243, 23244, 23245])

We mentioned earlier that some grains are artefacts of the watershed segmentation. In other words, they are outside the recrystallized central region. Therefore, those grains along with the two extensions (that were given the labels 500 and 501) were merged in Salome and a single mesh was constructed on the merged region. That region contains the most elements.

len(element_groups['homogeneous'])

3718

To be able to prescribe boundary conditions later, the boundary nodes were also exported from Salome.

node_groups

{'bottom': array([   8,   14,  448,  452,  467,  468,  470,  471,  472,  474,  475,
         477,  478,  479,  480,  481,  482,  495,  496,  497,  498,  499,
         557,  558,  559,  560,  561,  562,  563,  564,  565, 3673, 3674,
        3675, 3676, 3706, 3707, 3708, 3839, 3840, 3841, 3842, 3851, 3852,
        3861, 3862, 3863, 3864, 3865, 3866, 3867, 3874, 3880, 3881, 3882,
        3883, 3884, 3885, 3886, 3887, 3888, 3895, 3896, 3897, 3914, 3915,
        3922, 3923, 3924, 3925, 3926, 3930, 3931, 3932, 3933, 3943, 3944,
        3945, 3951, 3952, 3953, 3954, 3955, 3956, 3957, 3958, 3959, 3960,
        3961, 3962, 3963, 3964], dtype=int32),
 'left': array([  9,  10, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493,
        494], dtype=int32),
 'right': array([ 15,  16, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555,
        556], dtype=int32),
 'top': array([   0,   11,   17,   20,   21,   25,   28,   33,   34,   36,   44,
          45,   46,   52,   61,   63,  541,  542,  543,  544,  594,  595,
         596,  597,  598,  599,  600,  601,  602,  609,  653,  654,  655,
         656,  673,  674,  675,  676,  677,  678,  679,  702,  703,  704,
         705,  706,  707,  724,  725,  749,  750,  751,  759,  760,  761,
         762,  763,  764,  812,  813,  814,  815,  816,  817,  827,  828,
         829,  830,  831,  832,  848,  849,  850,  851,  896,  897,  898,
         920,  921, 1637, 1638, 1639, 1640], dtype=int32)}

Handling all these data individually is tedious. Hence, we created a mesh class that bundles them and defines operations on them. The TriMesh class is an abstraction for a triangular mesh. First, we initialize it with the elements and the nodes that were loaded from the 1_mesh_extended.npz file.

from grains.geometry import TriMesh
mesh = TriMesh(nodes, elements)

The TriMesh class has several options to plot the mesh, for now we consider the default plotting.

mesh.plot()
scale_figure(150)

Then the element and node sets are associated to the mesh object.

# Associate the element and node groups to it
for group_name, elements in element_groups.items():
    mesh.create_cell_set(group_name, elements)
for group_name, nodes in node_groups.items():
    mesh.create_vertex_set(group_name, nodes)

Mesh-based discretization methods often presume that the node ordering within an element is counter-clockwise. To make sure this is the case, we explicitly enforce it.

mesh.change_vertex_numbering('ccw', inplace=True);

The mesh was created on the geometry, which was constructed based on the image, i.e. in pixel units. For the simulation to be comparable with the experimental data (full-field measurement), the real size of the tensile speciment must be given. In the first image of this notebook, you can observe a ruler. With the help of it, we set the correspondance between the physical unit and the pixel unit: there are 29.55 pixels in 1 mm.

mesh.scale(1/29.55, inplace=True);

For convenience, we rotate the specimen so that its axis is parallel to the horizontal tensile loading.

from math import pi
mesh.rotate(-pi/2, inplace=True);

Let us check if this is what we wanted.

mesh.plot(vertex_legends=True)
scale_figure(150)

We are satisfied with the resulting transformations, so let us save them.

np.savez_compressed(join(data_dir, '1_mesh_extended_scaled.npz'),
                    nodes=mesh.vertices, elements=mesh.cells,
                    element_groups=mesh.cell_sets, node_groups=mesh.vertex_sets)

Conclusions

We successfully obtained a good quality mesh starting from an image.

In the near future, we will solve an inverse problem to identify characteristic parameters in the constitutive model. This optimization task requires the comparison of the numerical solution obtained on the mesh with the measured field values acquired through digital image correlation (DIC). Some of this work is already available in the dic and simulation modules.

This document is meant to serve as a global overview of what CristalX can be used for. Feel free to play with the parameters to investigate their effects. The functions and classes of CristalX provide many more functionalities. Browse the documentation to familiarize yourself with the details.

---

The initially created directory (data/) is now deleted.

from shutil import rmtree
try:
    rmtree(data_dir)
except FileNotFoundError:
    pass

Geometry and mesh processing in Salome

The cad module, as shown in the Algorithms section, allows to approximate a segmented image with planar spline surfaces (splinegons) and those surfaces can be written to a STEP file. Now, we demonstrate on the sample microstructure how to repair the geometry and generate a conforming mesh on the splinegons. All the manipulations are done in Salome 9.4.0.

Geometry

1. Import the geometry To import a STEP file (obtained by writing the splinegons into a STEP file), select the Geometry module and then File -> Import -> STEP.

2. Repair the geometry

   For the sample microstructure, one grain has not been identified, i.e. it contains a hole. To fill the hole, we perform the following steps (always select the object that was created the last time):

   1. Repair -> Close Contour and select the (in our sample) two boundary curves of the unidentified grain.

   2. Repair -> Suppress Holes and select the two boundary curves you selected in the previous step. Salome informs you that a face will be created in place of the hole.

   3. Repair -> Limit tolerance and set the tolerance to 1e-3.

   4. Repair -> Sewing and set 1e-2 for the tolerance.

3. Extract the grains We will need the mesh on each grain, therefore, we extract the grain faces by “exploding” the microstructure using New Entity -> Explode. In the dialog box, select the sewn geometry (end result of the geometry repairing workflow above) as Main Object and Face as Sub-shapes Type. Salome properly identifies the 250 sub-shapes, i.e. the grains.

4. Fix the “artificial” grains

   The sample microstructure describes the central part of a tensile specimen. However, to adhere to the Saint-Venant principle, the loading is exerted further from this central zone. Hence, we created long enough regions in the direction of the prescribed load by attaching two rectangular faces to the central region. These faces are not precise rectangles because the boundary of the microstructure does not consist of straight segments. Nevertheless, from now on, we will use the term rectangle to describe the extended region. To achieve perfect matching between the rectangle and the boundary of the microstructure, one side of the rectangle must contain the same lines as the boundary of the microstructure. The other three sides will be straight line segments.

   1. We want to extract the boundary lines of the microstructure. To do that, we explode the boundary grains into edges with New Entity -> Explode. The Main Object in the dialog box is a selected boundary grain, and the Sub-shape Type is Edge. Do this with all the boundary grains.

   2. Create the three sides of the rectangle

      1. Create the four vertices of the rectangle with New Entry -> Basic -> Point. Two of these points are the extremities of the microstructure, the other two are given such that

         • the rectangle’s three sides become parallel to the coordinate axes (so that boundary conditions are easily prescribed)

         • the length of the longer side of the rectangle is long enough compared to the microstructure (we applied a ratio of 5)

      2. Connect these vertices with three lines, which forms the three sides of the rectangle: New Entity -> Build -> Edge.

   3. Create the rectangle by connecting the three line segments constructed in the previous step with the boundary lines of the fourth side, obtained by exploding the boundary grains: New Entity -> Build -> Wire.

   4. Delete the original “artificial” grain as it will be replaced by the new one.

   5. Create the surface enclosed by the rectangle: New Entity -> Build -> Face, and select the wire created before. Perform steps 1-5 for the other rectangle (on the other side of the microstructure) as well.

   6. In order to create a conforming mesh on the whole domain (i.e. on the microstructure and on the two rectangles), we need to create a compound surface: New Entity -> Build -> Compound and select the grains plus the two rectangles.

   7. It is not enough to create the mesh, we also need to know which elements belong to which grains. To allow this, we explode the compound surface: select the compound object created in the previous step as the Main Object in New Entity -> Explode. The Sub-shape Type is Face, as before.

The geometry is now impeccable, let us start meshing.

Mesh

Select the Mesh module.

1. Create the mesh

   1. Choose Mesh -> Create Mesh from the menu. The geometry on which the mesh will be created is the compound surface.

   2. The mesh type in this study is Triangular and the algorithm is NETGEN 1D-2D. Choose NETGEN 2D Parameters as a hypothesis and experiment with the settings to suit your needs. Accept the changes.

   3. Choose Mesh -> Compute to generate the mesh with the chosen settings.

   If you want to alter the mesh, change the hypothesis and compute the new mesh.

2. Obtain the sub-meshes for the grains

   The segmented image contains regions that are not part of the recrystallized region, but they belong to the homogeneous region. We want to handle them in the same way as the two artificial grains. We do not need to merge the surfaces, just put them in the same element group. Therefore, we create

   • a common element group for the elements lying in the homogeneous region

   • one element group for each grain in the heterogeneous (recrystallized) region

   using Mesh -> Create Groups from Geometry and selecting all the 250 faces of the geometry. To form a single group for the heterogeneous region, merge the corresponding element groups in Mesh -> Union Groups. Once done, delete the groups that were used in forming the union.

3. We will explicitly need to prescribe boundary conditions on the left and right sides of the rectangular domain. For this, the nodes on those sides are collected with Mesh -> Create Group and set the Elements Type to Node. Select the nodes on the left edge and give the node set a name. Repeat it for the nodes on the right.

4. Export the mesh Click on File -> Export -> MED file. The sub-meshes (element groups, node groups) are also exported.

The mesh is ready for further processing. We have another guide that discusses how to handle the MED file from within Python.

Processing a .med file

After exporting the mesh from Salome to a MED file, we may want to perform certain operations on it. The MEDCoupling tool of Salome provides C++ and Python APIs for this purpose. However, that requires the user to

• have Salome installed as those APIs are available from the Salome kernel

• get to know the API

Moreover, it can happen that some mesh processing functionalities they may want to use does not exist. Since meshes consisting of cells of the same type (e.g. triangles) can be represented as homogeneous and contiguous arrays, converting the mesh from MED to numpy arrays seems a reasonable choice. This is what our med module does: it provides a thin wrapper around MEDCoupling to extract the mesh and the defined groups (cell and vertex groups) from the MED file and convert them to numpy arrays. This way, the user who deals with numerical modelling can implement their mesh processing algorithms based on numpy arrays, which is fast and straightforward. Furthermore, the person who performs the CAD operations and has Salome installed, can use our med module to export the mesh to numpy arrays so that the numerical analyst can directly work on it without having to have Salome installed and without any knowledge on the MEDCoupling API.

If you want to know more about the implementation details, read the documentation for the med module.

Using the module

To use our med module, access to the MEDLoader module from Salome is required. In the following, we assume that Salome has been installed and added to the path, so the command salome is available. If you

• want to use the Python REPL:

  $ salome shell -- python
  Python 3.6.5 (default, Dec 16 2019, 16:42:15) 
  [GCC 7.3.0] on linux
  Type "help", "copyright", "credits" or "license" for more information.
  

  Now, you have access to the MEDCoupling Python API. E.g.

  >>> from MEDLoader import MEDFileData
  

• want to use the PyCharm IDE

  $ salome shell "path_to_PyCharm/bin/pycharm.sh"
  

  This will start PyCharm. In the IDE (PyCharm in this example), set the interpreter to that of Salome’s Python. For me, it is located at BINARIES-UB18.04/Python/bin/python3, where UB18.04 refers to the fact that I downloaded a pre-compiled Salome binaries for Ubuntu 18.04.

To learn more about the salome command, read the manual.

Algorithms

The goal of this document is to provide more low-level algorithmic details than what is given in our paper. In fact, the material presented here is complimentary to the paper. The paper concentrates on the intuition behind the methods and shows the high-level structure, also citing the relevant literature.

From image to geometry

It is recommended to read Section 2.3 of the paper first and then this section.

The algorithm to construct a geometrical representation out of a segmented image is complex. To reduce the complexity, we created composable parts. We start with some notions below. Then the main steps of the algorithm is discussed in the upcoming subsections.

Grains can have different representations. When represented as an image (input in our algorithm), it is a set of pixels labelled with the same positive integer. A grain can be given different geometrical representations. It can be the assembly of primitive shapes, usually triangles. One can also describe a grain by its boundary. For example, piecewise linear boundary segments lead to a polygon, connecting spline curves lead to a splinegon. A grain can also be represented as a cycle of a graph.

We make distinction between topological and geometrical data. Topological data includes connectivity, neighbourhood and membership information. Coordinates of points are examples for geometrical data. Separating these two terms allows us to build abstractions in the code, making it understandable and extensible.

Skeleton of the image

Based on the segmented/labelled image, skan builds the skeleton network, as demonstrated on a sample image below. The different pixel colors represent different labels.

As a preprocessing step, we create an additional labelled region that surrounds the original labelled image. The reason is that this way, skan will create branches along the boundary of the original image, which will close the boundary grains. It also simplifies the algorithms we will use later because the boundary interfaces can be handled in the same way as the internal ones: every interface separates exactly two grains.

Which branches form a grain?

The question arises: given the branches, how to reconstruct the grains? First, let us try a graph theoretic approach. In the example below, we want to identify the three grains denoted by roman numerals. The network of branches can be represented as a graph, in which the vertices are the end points and the edges correspond to the branches.

Adequately selected cycles in this graph would give the grains. However, as discussed in the paper, determining the grains based on the graph exclusively would be very challenging because of the following characteristics.

• For general microstructures, the graph contains multiple edges (see the example above), which rules out many graph processing methods.

• Finding all the elementary cycles in the graph is too costly for graphs coming from realistic microstructures. Even if we found all the cycles, we would need a criterion to choose which ones correspond to grains. E.g. the cycle \( \bar{6}-\bar{2}-\bar{3}-\bar{5} \) does not encompass a single grain but the union of two grains.

• Another technique to find the correct cycles would be the minimum cycle basis. However, this basis is not unique so there is no guarantee that we find the adequate cycles.

The solution is to use both topological and geometrical information. The skeleton is superimposed on the labelled image (from which it was constructed by skan) and the labels around a given skeleton node or end point is detected. Staying with the example in the first subsection, the scheme is shown below.

In the following, we describe an algorithm to find out which two grains are incident to a branch (remember that always two grains neighbor a branch).

Algorithm

The neighbor search around every node of a branch is performed and the two most common labels are chosen. Since in the image representation the grain is a set of pixels having the same label, the selected two labels give the neighboring two grains to a branch.

Neighbor definitions

As described in the paper, certain scenarios necessitate to consider various neighborhood definitions. For those definitions, see the documentation of the grains.utils.neighborhood function.

Going back to the first configuration in the beginning of this section, the algorithm gives the following branch-grain connectivities.

\[\begin{split}\bar{1} &: [\mathrm{I}] \\ \bar{2} &: [\mathrm{II}] \\ \bar{3} &: [\mathrm{I}, \mathrm{II}] \\ \bar{4} &: [\mathrm{III}, \mathrm{II}] \\ \bar{5} &: [\mathrm{I}, \mathrm{III}] \\ \bar{6} &: [\mathrm{II}, \mathrm{I}] \\\end{split}\]

Inverting this relationship gives the grain-branch connectivities:

\[\begin{split}\mathrm{I} &: [\bar{1}, \bar{3}, \bar{5}, \bar{6}] \\ \mathrm{II} &: [\bar{2}, \bar{3}, \bar{6}, \bar{4}] \\ \mathrm{III} &: [\bar{5}, \bar{4}]\end{split}\]

The grain-branch connectivities are an intermediate representation (topological-geometrical). It is independent of how we geometrically represent a grain later.

Grains as oriented planar surfaces

The previous part of the reconstruction algorithm determined which branches bound a grain. In order to obtain a surface representation of a grain, the boundary must be oriented and hence the branches must be connected in the appropriate order. The following figure demonstrates for grain \( \mathrm{II} \) the working of a brute-force algorithm.

The branches are interlaced based on their common junctions. The arrows show which branches follow in order. Note that the default orientation of branches \( \bar{3} \) and \( \bar{4} \) needs to be swapped.

Finally, we arrive at a fully geometrical description because each grain is now given by a series of points (nodes and end points) along its boundary.

Geometrical representations of grains

Now that we have a list of points, we can build two geometrical representations of a grain. In the polygon representation, the list of points are the consecutive vertices of the polygon. Their coordinates are stacked, the first vertex being repeated to “close” the polygon. In the spline representation, the list of points on each branch act as a knots of a B-spline. Once the bounding splines have been constructed, the planar spline surface (splinegon) is spanned by those bounding splines.

Program design

CristalX is not a black-box library such as BLAS, neither is a GUI-based application intended for end-users. It is rather an easy-to-use and extensible set of Python codes that provide the basic functionalities that scientists can extend based on their needs. The following ideas were kept in mind while writing and maintaining CristalX.

• Driven by actual needs

Only implement features that are currently used. Adding extra features requires more testing, possibly more dependencies and therefore code bloat, and increases the cognitive load of the user. Instead, the emphasis is on creating a stable minimum core library that can be easily extended according to users’ demands. Consequently, application code is separated from the core modules.

• Build on well-established packages

We rely on the scientific Python stack: NumPy for array manipulations, SciPy for interpolation and some other computations, Matplotlib for visualization and scikit-image for image processing. This ensures interoperability with other scientific codes and that our software is hopefully bug-free.

• Minimize the dependencies

Rapid prototyping is essential in scientific code development and Python is an excellent choice to satisfy this requirement. At the same time, relying on fast libraries ensures that the computations are reasonably fast. The libraries mentioned in the previous point are easy to install, often already pre-installed in certain Python distributions.

• High-quality documentation

Future contributors will benefit from the rich documentation. Python doctests are extensively used, serving both as test cases and as examples of usage. The docstrings conform to the numpydoc style guide.

• We strive for decoupling the modules

Although part of the grains package, if the modules are independent, they can be reused in other projects too just by copy-pasting the required functions.

• Do not overuse classes

In the prototyping phase, prefer using free functions to methods. As an idea evolves, you will naturally find data and algorithms that belong together, and can refactor free functions into member functions of a class.

• Do not use deep hierarchies

Initially, stay away from excessive nesting to avoid fragmenting the code base. If the project grows big, you can still refactor the code by introducing deeper hierarchies. Deep nesting causes unnecessary cognitive load and it also makes the code more verbose at the caller’s site. Compare

import package1.package2.module

with

import package.module

• Gradually refactor code

As more and more features are added to the project, we will often find that similar tasks emerge in different contexts. It is a good time to think about how they can be generalized and to reconsider your model. This way, you will come up with utility functions best put into utils.py.

• Start writing code only after careful thinking

It is no point in writing code before you completely understand your problem domain you want to model. It is more efficient to build abstractions in your head or on paper, then to split it into modular chunks, and only after that start coding.

• Write the documentation before the code

If you document the function parameters and the return values in advance, as well as construct a doctest, you are enforced to think about the problem deeply and to create a good interface. Moreover, it guarantees that the documentation is not missing (what you would anyway have to write at some point, so why not at the beginning?).

• Give doctest-compatible examples

You hit two birds with the same stone: provide an example for the user and get some confidence that your code works as intended (at least for the particular example). As mentioned in the previous point, write them before the actual code implementation. Doctests do not replace careful testing.

• Keep the documentation as part of your code

The problem with wiki pages is that they are version controlled in a different Git repository. It makes it longer to change a hosting service (e.g. moving from GitHub to GitLab), you need to maintain two repositories, cannot change the documentation and the code in the same commit, and you have to rely on the rendering capabilities of the hosting service (e.g. GitHub cannot render math). It is therefore better to keep the documentation as part of your project in a dedicated directory (docs/ in our case), use a documentation generator (Sphinx in our case) and host it online (on Read the Docs in our case).

Coding conventions

Documentation

Writing documentation is necessary when you contribute to this project, either by writing code or by changing/extending the external documents. In this section, we give tips how you can write them. This is an opinionated topic and it lays down how it is currently done. Feel free to suggest new ideas.

Whenever you write the documentation, first test it locally before pushing changes: Read the Docs spends about 1000 seconds to build the documentation.

It took me a long time to experiment with the Sphinx settings that provide the output what you can see in the rendered documentation. Some notes concerning these efforts can be found in the Notes section. Another way to learn about Sphinx documentation is by reading the source of the existing documentation. On each page of the HTML documentation, the header contains a View page source hyperlink.

Code documentation

As shown in the figure, codes are written either as Python files or as Jupyter notebooks. With the proper extensions (see the docs/source/conf.py file), both of them are automatically included in the documentation by Sphinx. In what follows, we concentrate on documenting Python files.

The docstings are written in RST, following the numpydoc style, using the Napoleon Sphinx extension.

Provide doctests to demonstrate the use of the function you write and to provide minimal testing.

The line length is 100 characters, as defined in the /.editorconfig file.

External documentation

The external documentation is written in reStructuredText (RST) and in Markdown. Sphinx, by default, uses RST, extending it with more capabilities. However, recommonmark can parse Markdown files and automatically convert them to RST at documentation build time. This is a great help as more people are used to the simple Markdown than to the more complex (and more capable) RST. Actually, the current document you read was also written in Markdown. Of course, it is perfectly fine if you write the documentation exclusively in RST.

Currently, the structure of the documentation is written in RST, and most of the external documentation in Markdown.

Development

Intro …

General workflow

The development workflow can be followed in the following figure.

1. You write the Python source codes (.py) and the Jupyter notebooks (.ipynb) on your local machine. If you want, you can create the documentation locally using Sphinx.

2. When a logical unit has been finished, you commit the changes with git push. This will upload the new version of the modified files to a remote repository (currently GitHub). If you edit files in the remote repository and you want those changes to be present in your local copy, you can use git pull. For more details on Git, read the manual.

3. Several commit hooks are attached to the remote repository. When they detect a change, certain actions are activated. One the one hand, the online version of the documentation, hosted on Read the Docs will be updated. On the other hand, static analyzers will reanalyze the new version of the code. Some of them may create a pull request based on their recommendations.

Profiling

I use Pyinstrument for profiling the code. It comes with the CristalX installation.

Otherwise, you can install it with pip install pyinstrument or by conda after the conda-forge channel has been activated.

conda install -c conda-forge pyinstrument

As Pyinstrument has no dependencies, you can safely install it to your current environment. If you do not want to take a risk, create a new environment. With conda, you can do e.g.

conda create -n pyinstrument python=3.7 scipy matplotlib
conda activate pyinstrument
conda install -c conda-forge pyinstrument

The profiling module provides a wrapper around Pyinstrument. Put the code you want to profile in the profile context manager, e.g.

>>> import random
>>> from grains.profiling import profile
>>> with profile('html') as p:
...    for _ in range(1000000):
...        rand_num = random.uniform(1, 2.2)

For more details see the API reference.

Contributing

CristalX is an open-source project that welcomes contributions of any kind.

What can you help in?

The major fields in which you can help are the following (in increasing complexity).

Use cases

If you use CristalX in your project, let us know. Parts of your code, if they are general enough to be incorporated, could be included in CristalX. Make contact with us by opening an issue.

Documentation

If the documentation of some functions

• is missing

• is incomprehensible

• does not contain examples

• is not rendered properly

or if the code examples break, open a pull request. Similarly, open a pull request if the guides

• contain typos

• are incomprehensible

Code

• You can report bugs by opening an issue.

• If you want to help but do not know where to start, consider the currently open issues, especially the ones with the help wanted label.

• You can implement new features. First of all, contact us if you plan to work on a non-trivial feature. This will save work for you. Use the fork & PR workflow.

How to contribute

First, create a GitHub account. There are two ways to contribute.

Open an issue

To open a new issue, click on the green button on the Issues page.

Fork & pull request

At the moment, there are no code formatting guidelines, the best is to follow the formatting of the existing code.

1. Fork the GitHub repository by clicking on the Fork button at the top right corner.

   

2. Install CristalX locally by cloning your fork (git clone).

3. Create a new branch.

   Do not work on the default (master) branch but create a new feature branch, e.g.

   git checkout -b new_feature
   

4. Make your changes.

   Add or modify files and regularly commit your changes to your local clone with meaningful commit messages (git commit). Do not forget to test your code: if you do not provide unit tests, at least write doctests.

5. Push your changes to your remote fork on GitHub (git push).

6. Visit your forked repository on GitHub and click on the Pull request button. See the GitHub documentation for details.

   

7. Keep your fork up to date

   While you work on your forked repository, changes may be committed to the original (called upstream) repository. To make sure you keep your fork up to date with the upstream repository, follow the instructions in the GitHub documentation.

This was a basic overview, for more details check out the following documents:

• https://gist.github.com/Chaser324/ce0505fbed06b947d962

• https://github.com/susam/gitpr

Versioning

The versioning of CristalX starts with 1.0.0 and it somewhat follows the rules of Semantic Versioning 2.0.0, with the syntax MAJOR.MINOR.PATCH, where MAJOR introduces significant changes, MINOR comes with smaller changes, and PATCH provides a fixture or a tiny improvement either in the code or in the documentation.

Why not SemVer?

Semantic Versioning (SemVer) is a widely used versioning scheme, applicable for public APIs. Its purpose is to be rigorous on how to indicate when a bug fix, new features or incompatible changes in the public API are introduced. There is an excellent discussion about it here and a detailed guide here. Many criticize it for not being indicative about the rate of important changes. E.g. 1.8.5 –> 1.9.0 may include dozens of relevant improvements, while 1.9.5 –> 2.0.0 may merely be a simple clean-up that changes the public API. However, SemVer was never meant to be used for software version numbers.

CristalX is not a library, but a collection of tools that operate at a high level. Most of its functions are exposed to the user, except the ones marked with a leading underscore or two leading underscores. However, the public API would be a subset of these functions: e.g. the functions of the utils module are exposed but they are mostly intended to be used by other modules, not by the user. As CristalX is not a library but rather a tool for rapid prototyping, semantic versioning would not make much sense.

Why not CalVer?

Calendar Versioning (CalVer) is an alternative to SemVer. CristalX may not come with regular changes in the future or it will get updates at irregular intervals. We do not want to give the impression that relevant changes are introduced linearly in time. Moreover, the changelog and the time stamp in the git commits clearly show when new releases are published.

Our versioning

As there is no unconditionally best versioning system, we came up with our own, which seems to fit well for research code like CristalX. The starting point is SemVer with some differences. Our versioning is intended for humans. The MAJOR version in increased only for substantially new features. This is, of course, subjective but we want to avoid large MAJOR version numbers as e.g. in Firefox. The guideline to follow is that the novelty of a feature makes the MAJOR version increase, not the number of additions. Here are some examples. Assume that we are at 1.2.3. We fixed a set of related bugs in the code: 1.2.3 –> 1.2.4. Then we implemented several functionalities to speed up the code: 1.2.4 –> 1.3.0. A new module was created and several others were modified that allowed us to carry out groundbreaking research: 1.3.0 –> 2.0.0. An existing algorithm was improved to handle the corner cases: 2.0.0 –> 2.1.0.

Similarly to SemVer, when one of the digits increases, the ones right to it are set to zero (e.g. 1.0.1 –> 1.0.2, 1.0.3 –> 1.1.0, 1.8.6 –> 2.0.0).

We try to keep backward compatibility. Insignificant changes are postponed, and are included as part of a MINOR release. If you often feel the need to introduce changes to the function signatures, rather add new functions and give deprecation notices than remove or modify existing ones. This does not lead to code bloat in the long run because deprecated syntax is removed from time to time. To mark a function, a method or a class as deprecated, import the deprecated function from the deprecation package and follow its syntax.

The master branch of the Git repository contains the latest developments since the last published release. These unreleased modifications are allowed to contain incompatibilities (changes in the function signatures, etc.) compared to the latest release. This flexibility is essential for rapid prototyping in research codes like CristalX. However, these possible incompatibilities must be fixed for the next release, see the previous paragraph.

Connection with the Git workflow

The version numbers are reflected in the tag names. Your normal Git workflow stays the same: commit modifications and push them to the remote repository. When you want to mark a commit yet to be included in a certain version, type

git tag -m "Concise message" v<version_number>

where <version_number> has the form MAJOR.MINOR.PATCH. Note that it is preceded by the v letter, conventionally used for tags. As an example:

git tag -m "Initial release." v1.0.0

Keep the tag message short: the detailed changes since the previous version are collected in the changelog.

The tag is pushed by

git push --tags

Whenever you publish a tag, and hence update the changelog, also create a release for that tag on GitHub. Copy the changes the new version brings from the changelog to the description of the release. In this regard, we follow the way of JabRef. The zipped size of CristalX is quite small, so size constraint will not be a problem on GitHub.

Segmentation

This module contains the Segmentation class, responsible for the image segmentation of grain-based materials (rocks, metals, etc.)

Classes

Segmentation

Segmentation of grain-based microstructures

class grains.segmentation.Segmentation(image_location, save_location=None, interactive_mode=True)

Bases: object

Segmentation of grain-based microstructures

original_image

Matrix representing the initial, unprocessed image.

Type

ndarray

save_location

Directory where the processed images are saved

Type

str

create_skeleton(boundary_image)

Use thinning on the grain boundary image to obtain a single-pixel wide skeleton.

Parameters

boundary_image (bool ndarray) – A binary image containing the objects to be skeletonized.

Returns

skeleton (bool ndarray) – Thinned image.

filter_image(window_size, image_matrix=None)

Median filtering on an image. The median filter is useful in our case as it preserves the important borders (i.e. the grain boundaries).

Parameters

• window_size (int) – Size of the sampling window.

• image_matrix (3D ndarray with size 3 in the third dimension, optional) – Input image to be filtered. If not given, the original image is used.

Returns

filtered_image (3D ndarray with size 3 in the third dimension) – Filtered image, output of the median filter algorithm.

find_grain_boundaries(segmented_image)

Find the grain boundaries.

Parameters

segmented_image (ndarray) – Label image, output of a segmentation.

Returns

boundary (bool ndarray) – A bool ndarray, where True represents a boundary pixel.

initial_segmentation(*args)

Perform the quick shift superpixel segmentation on an image. The quick shift algorithm is invoked with its default parameters.

Parameters

*args (3D numpy array with size 3 in the third dimension) – Input image to be segmented. If not given, the original image is used.

Returns

segment_mask (ndarray) – Label image, output of the quick shift algorithm.

merge_clusters(segmented_image, threshold=5)

Merge tiny superpixel clusters. Superpixel segmentations result in oversegmented images. Based on graph theoretic tools, similar clusters are merged.

Parameters

• segmented_image (ndarray) – Label image, output of a segmentation.

• threshold (float, optional) – Regions connected by edges with smaller weights are combined.

Returns

merged_superpixels (ndarray) – The new labelled array.

save_array(filename, array)

Save an image as a numpy array. The array is saved in the standard numpy format, into the directory determined by the save_location attribute.

Parameters

• filename (str) – The array is saved under this name, with extension .npy

• array (ndarray) – An image represented as a numpy array.

save_image(filename, array, is_label_image=False)

Save an image as a numpy array. The array is saved in the standard numpy format, into the directory determined by the save_location attribute.

Parameters

• filename (str) – The array is saved under this name, with extension .npy

• array (ndarray) – An image represented as a numpy array.

• is_label_image (bool) – True if the array represents a labeled image.

watershed_segmentation(skeleton)

Watershed segmentation of a granular microstructure. Uses the watershed transform to label non-overlapping grains in a cellular microstructure given by the grain boundaries.

Parameters

skeleton (bool ndarray) – A binary image, the skeletonized grain boundaries.

Returns

segmented (ndarray) – Label image, output of the watershed segmentation.

Gala

A trimmed version of the Gala project (https://github.com/janelia-flyem/gala) with some additions (new function, added documentation for the existing ones). Gala is licensed by the Janelia Farm License: http://janelia-flyem.github.io/janelia_farm_license.html

Copyright 2012 HHMI. All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

Neither the name of HHMI nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Functions

imextendedmin(image, h[, connectivity])	Extended-minima transform.
hminima(a, thresh)	Suppress all minima that are shallower than thresh.
imhmin(a, thresh)	Suppress all minima that are shallower than thresh.
morphological_reconstruction(marker, mask[, …])	Perform morphological reconstruction of the marker into the mask.
regional_minima(a[, connectivity])	Find the regional minima in an ndarray.
complement(a)

grains.gala_light.complement(a)

grains.gala_light.hminima(a, thresh)

Suppress all minima that are shallower than thresh.

Parameters

• a (array) – The input array on which to perform hminima.

• thresh (float) – Any local minima shallower than this will be flattened.

Returns

out (array) – A copy of the input array with shallow minima suppressed.

grains.gala_light.imextendedmin(image, h, connectivity=1)

Extended-minima transform. The extended minima transform is the regional minima of the h-minima transform. The implementation follows the MATLAB function under the same name.

Parameters

• image (ndarray) – The input array on which to perform imextendedmin.

• h (float) – Any local minima shallower than this will be flattened.

• connectivity (int, optional) – Determines which elements are considered as neighbors of the central element. Elements up to a squared distance of connectivity from the center are considered neighbors. If connectivity=1, no diagonal elements are neighbors.

Returns

bool ndarray – True at places of the extended minima.

grains.gala_light.imhmin(a, thresh)

Suppress all minima that are shallower than thresh.

Parameters

• a (array) – The input array on which to perform hminima.

• thresh (float) – Any local minima shallower than this will be flattened.

Returns

out (array) – A copy of the input array with shallow minima suppressed.

grains.gala_light.morphological_reconstruction(marker, mask, connectivity=1)

Perform morphological reconstruction of the marker into the mask.

See the Matlab image processing toolbox documentation for details: http://www.mathworks.com/help/toolbox/images/f18-16264.html

grains.gala_light.regional_minima(a, connectivity=1)

Find the regional minima in an ndarray. As written in the MATLAB documentation of the imregionalmin function: “Regional minima are connected components of pixels with a constant intensity value, surrounded by pixels with a higher value.”

Analysis

This module contains the Analysis class, responsible for the analysis of segmented grain-based microstructures.

All the examples assume that the modules numpy and matplotlib.pyplot were imported as np and plt, respectively.

Classes

Analysis

Analysis of grain assemblies.

Functions

feret_diameter(prop)	Determines the maximum Feret diameter.
plot_prop(prop[, pixel_per_unit, show_axis])	Plots relevant region properties into a single figure.
plot_grain_characteristic(characteristic, …)	Plots the distribution of a given grain characteristic.
show_label_image(label_image[, alpha])	Displays a labeled image.
label_image_skeleton(label_image)	Skeleton of a labeled image.
thicken_skeleton(skeleton, thickness)	Thickens a skeleton by morphological dilation.
label_image_apply_mask(label_image, mask, value)	Changes parts of a labeled image to a given value.

class grains.analysis.Analysis(label_image, interactive_mode=False)

Bases: object

Analysis of grain assemblies.

original_image

Matrix representing the initial, unprocessed image.

Type

ndarray

save_location

Directory where the processed images are saved

Type

str

compute_properties()

Determines relevant properties of the grains. The area of each grain is determined in the units previously given in the set_scale method.

Parameters

• window_size (int) – Size of the sampling window.

• image_matrix (3D ndarray with size 3 in the third dimension, optional) – Input image to be filtered. If not given, the original image is used.

Returns

filtered_image (3D ndarray with size 3 in the third dimension) – Filtered image, output of the median filter algorithm.

set_scale(pixel_per_unit=1)

Defines a scale for performing computations in that unit. Image measures (area, diameter, etc.) are performed on a matrix corresponding to a label image. Therefore, the result of all the computations are obtained in pixel units. It is often of interest to access the results in physical units (mm, cm, inch, etc.). Manually converting pixels, pixel squares, etc. to pyhsical units, physical unit sqaures, etc. are tedious and error prone. Once the conversion between a pixel and a physical unit is given, all the subsequent calculations are performed in the desired physical unit.

Parameters

pixel_per_unit (float or int or scalar ndarray, optional) – Number of pixels contained in a certain unit. The default is 1, in which case all measurements are performed in pixel units.

Returns

None.

show_grains(grain_property=None)

Display the grains, optionally with a property superposed.

Parameters

grain_property ({None, ‘area’, ‘centroid’, ‘coordinate’,) – ‘equivalent_diameter’, ‘feret_diameter’, ‘label’}

optional

If not None, the selected property is shown on the grain as text.

Returns

None.

show_properties(gui=False)

Displays previously computed properties of the grains

Parameters

gui (bool, optional) – If true, the grain properties are shown in a GUI. If false, they are printed to stdout. The default is False. The GUI requires the dfgui modul, which can be obtained from https://github.com/bluenote10/PandasDataFrameGUI

Returns

None.

grains.analysis.feret_diameter(prop)

Determines the maximum Feret diameter.

Parameters

prop (RegionProperties) – Describes a labeled region.

Returns

max_feret_diameter (float) – Maximum Feret diameter of the region.

See also

skimage.measure.regionprops()

Measure properties of labeled image regions

Examples

>>> import numpy as np
>>> from skimage.measure import regionprops
>>> image = np.ones((2,2), dtype=np.int8)
>>> prop = regionprops(image)[0]
>>> feret_diameter(prop)
2.23606797749979

grains.analysis.label_image_apply_mask(label_image, mask, value)

Changes parts of a labeled image to a given value.

Convenience function, equivalent to label_image[mask] = value but the original array label_image is not overwritten.

Parameters

• label_image (ndarray) – Labeled input image, represented as a 2D numpy array of positive integers.

• mask (ndarray) – Boolean array of the same size as label_image, marking the pixels that will be replaced by value.

• value (int) – The masked pixels are replaced by this value.

Returns

ndarray – Copy of the input image, its selected pixels being replaced by the given value.

grains.analysis.label_image_skeleton(label_image)

Skeleton of a labeled image.

The skeleton of a labeled image is a single-pixel wide network that separates the labeled regions.

Parameters

label_image (ndarray) – Labeled input image, represented as a 2D numpy array of positive integers.

Returns

ndarray – A 2D bool numpy array having the same size as label_image, where True represents the skeleton pixels.

See also

thicken_skeleton()

grains.analysis.plot_grain_characteristic(characteristic, centers, interpolation='linear', grid_size=(100, 100), **kwargs)

Plots the distribution of a given grain characteristic.

One way to gain insight into a grain assembly is to plot the distribution of a certain grain property in the domain the grains occupy. In this function, for each grain, and arbitrary (scalar) quantity is associated to the center of the grain. In case of n grains, n data points span the interpolant and the given characteristic is interpolated on a grid of the AABB of the grain centers.

Parameters

• characteristic (ndarray) – Characteristic property, the distribution of which is sought. A 1D numpy array.

• centers (ndarray) – 2D numpy array with 2 columns, each row corresponding to a grain, and the two columns giving the Cartesian coordinates of the grain center.

• interpolation ({‘nearest’, ‘linear’, ‘cubic’}, optional) – Type of the interpolation for creating the distribution. The default is ‘linear’.

• grid_size (tuple of int, optional) – 2-tuple, the size of the grid on which the data is interpolated. The default is (100, 100).

Other Parameters

• center_marker (str, optional) – Marker indicating the center of the grains. The default is ‘P’. For a list of supported markers, see the documentation. If you do not want the centers to be shown, choose ‘none’.

• show_axis (bool, optional) – If True, the axes are displayed. The default is False.

Returns

None

See also

scipy.interpolate.griddata()

Notes

This function knows nothing about how the center of a grain is determined and what characteristic features a grain has. It only performs interpolation and visualization, hence decoupling the plotting from the actual representation of grains and their characteristics. For instance, a grain can be represented as a spline surface, as a polygon, as an assembly of primitives (often triangles), as pixels, just to mention some typical scenarios. Calculating the center of a grain depends on the grain representation at hand. Similarly, one can imagine various grain characteristics, such as area, diameter, Young modulus.

Examples

Assume that the grain centers are sampled from a uniformly random distribution on the unit square.

>>> n_data = 100
>>> points = np.random.random((n_data, 2))

The quantity we want to plot has a parabolic distribution with respect to the position of the grain centers.

>>> func = lambda x, y: 1 - (x-0.5)**2 - (y-0.5)**2
>>> plot_grain_characteristic(func(points[:, 0], points[:, 1]), points, center_marker='*')
>>> plt.show()

grains.analysis.plot_prop(prop, pixel_per_unit=1, show_axis=True)

Plots relevant region properties into a single figure. Four subfigures are created, giving the region’s

• image, its area and its center

• filled image, its area

• bounding box, its area

• convex image, its area

Parameters

• prop (RegionProperties) – Describes a labeled region.

• pixel_per_unit (float or int, optional) – Number of pixels contained in a certain unit. The default is 1, in which case all measurements are performed in pixel units.

Returns

fig (matplotlib.figure.Figure) – The figure object is returned in case further manipulations are necessary.

grains.analysis.show_label_image(label_image, alpha=1)

Displays a labeled image.

A random color is associated with each labeled region. If boundary pixels are present in the image, they are plotted in black.

Parameters

• label_image (ndarray) – Labeled input image, represented as a 2D numpy array of non-negative integers. The label 0 is assumed to denote a boundary pixel.

• alpha (float, optional) – Opacity of colorized labels. Must be within [0, 1].

Returns

None

grains.analysis.thicken_skeleton(skeleton, thickness)

Thickens a skeleton by morphological dilation.

Parameters

• skeleton (ndarray) – Skeleton of a binary image, represented as a bool 2D numpy array.

• thickness (int) – Thickness of the resulting boundaries.

Returns

ndarray – A 2D bool numpy array, where True represents the thickened skeleton.

See also

label_image_skeleton()

grains.analysis.truecolor2label(color_image)

Truecolor image into labeled image.

It is often the case that you need to deal with a labeled image that was saved as a truecolor image (e.g. RGB). A labeled region is then the set of pixels with the same colors.

Parameters

color_image (ndarray) – 3D array, the first two dimensions corresponding to the image pixels, the third one for the channels (e.g. RGB, HSV, CMYK, etc.).

Returns

ndarray – Labeled image, represented as a 2D numpy array of non-negative integers.

Meshing

Classes

SkeletonGeometry
QuadSkeletonGeometry
TriSkeletonGeometry
FixedDict	https://stackoverflow.com/a/14816446/4892892
OOF2

class grains.meshing.FixedDict(dictionary)

Bases: object

https://stackoverflow.com/a/14816446/4892892

class grains.meshing.OOF2

Bases: object

create_material(name)

Parameters

name (TYPE) – DESCRIPTION.

Raises

Exception – DESCRIPTION.

Returns

None.

create_microstructure(name=None)

Creates a microstructure from an image.

Parameters

name (str, optional) – Path to the image on which the microstucture is created, file extension included. If not given, the microstructure is given the same name as the input image.

Raises

Exception – DESCRIPTION.

Returns

None.

create_skeleton(nelem_x, nelem_y, geometry, name=None)

load_pixelgroups(microstructure_file)

Parameters

microstructure_file (str) – DESCRIPTION.

Returns

None.

materials2groups(materials, groups=None)

Parameters

• materials (list of str) – DESCRIPTION.

• groups (list of int, optional) – DESCRIPTION. The default is None.

Returns

None.

pixel2group()

read_image(label_image)

save_microstructure(name=None)

save_pixelgroups(name=None)

Parameters

name (str) – DESCRIPTION.

Returns

None.

script = []

show()

Returns

None.

write_script(name=None)

class grains.meshing.QuadSkeletonGeometry(leftright_periodicity=False, topbottom_periodicity=False)

Bases: grains.meshing.SkeletonGeometry

class grains.meshing.SkeletonGeometry(leftright_periodicity, topbottom_periodicity)

Bases: abc.ABC

class grains.meshing.TriSkeletonGeometry(leftright_periodicity=False, topbottom_periodicity=False, arrangement='conservative')

Bases: grains.meshing.SkeletonGeometry

grains.meshing.nt

alias of grains.meshing.modules

CAD

MED

Extracting and processing meshes from .med files. The functions were tested on the MEDCoupling API, version 9.4.0.

Todo

Support renumbering (https://docs.salome-platform.org/latest/dev/MEDCoupling/user/html/data_optimization.html).

Getting help:

• This module relies on the Python interface of MEDCoupling. Click here for the latest documentation.

• User’s manual for the Python interface

• To know more about the MED file format, which is a specialization of HDF5, see the documentation. For a discussion on the relation between the MED format and the APIS, see this page and that one.

• The definitions, such as group, used in this module are from the development guide.

• A (mostly English) tutorial for the Python interface to MEDCoupling is also useful. Particularly interesting are the mesh manipulation examples

• Main page of the documentation

Functions

read_mesh	Reads a mesh file in .med format.
get_nodes	Obtains the nodes and the node groups of a mesh.
get_elements	Obtains the elements for each group of a mesh.

grains.med.get_elements(mesh, numbering='global')

Obtains the elements for each group of a mesh.

Elements of the same dimension as the mesh are collected (e.g. faces for a 2D mesh).

Todo

put those elements that do not belong to any group into an automatically created group

Todo

support ordering elements in alphabetical order

Todo

implement the ‘global’ strategy

Parameters

• mesh (MEDFileUMesh) – Unstructured mesh.

• numbering ({‘global’}, optional) –

  Determines how to allocate element numbers in the mesh.

  ‘global’: numbers the elements without taking into account which group they belong to. Use this strategy if you are not sure whether an element belongs to more than one group. ‘group’: numbers the elements group-wise. This is much faster than the ‘global’ strategy, but use this option if you are sure that the groups of the mesh do not contain common elements.

  The default is ‘global’.

Returns

• elements (ndarray) – Element-node connectivities in a 2D numpy array, in which each row corresponds to an element and the columns are the nodes of the elements. It is assumed that all the elements have the same number of nodes.

• element_groups (dict) – The keys in the dictionary are the element group names, while the values are list of integers, giving the elements that belong to the particular group.

Warning

Currently, elements that do not fit into any groups are discarded.

See also

get_nodes(), change_node_numbering()

Notes

The element-node connectivities are read from the mesh. If you want to change the ordering of the nodes, use the change_node_numbering() function.

Both this and the get_nodes() function relies on getGroupsOnSpecifiedLev to obtain the groups based on a parameter, called meshDimRelToMaxExt. This parameter designates the relative dimension of the mesh entities whose IDs are required. If it is 1, it denotes the nodes. If 0, entities of the same dimension as the mesh are meant (e.g. group of volumes for a 3D mesh, or group of faces for a 2D mesh). When -1, entities of spatial dimension immediately below that of the mesh are collected (e.g. group of faces for a 3D mesh, or group of edges for a 2D mesh). For -2, entities of two dimensions below that of the mesh are fetched (e.g. group of edges for a 3D mesh).

grains.med.get_nodes(mesh)

Obtains the nodes and the node groups of a mesh.

Parameters

mesh (MEDFileUMesh) – Unstructured mesh.

Returns

• nodes (ndarray) – 2D numpy array with 2 columns, each row corresponding to a node, and the two columns giving the Cartesian coordinates of the nodes.

• node_groups (dict) – The keys in the dictionary are the node group names, while the values are list of integers, giving the nodes that belong to the particular group.

See also

get_elements(), getGroupArr

grains.med.read_mesh(filename)

Reads a mesh file in .med format. Only one mesh, the first one, is supported. However, that mesh can contain groups.

Parameters

filename (str) – Path to the mesh file.

Returns

MEDFileUMesh – Represents an unstructured mesh. For details, see the manual on https://docs.salome-platform.org/latest/dev/MEDCoupling/developer/classMEDCoupling_1_1MEDFileUMesh.html

Salome

The documentation generated from this file is available on https://cristalx.readthedocs.io/en/latest/api.html#module-grains.salome.

The aim of this module is to manage geometry and mesh operations on two-dimensional microstructures that tessellate a domain. This has the important consequence that the domain is assumed to consist of non-overlapping shapes that cover it. After generating a conforming mesh, each element belongs to a single group and no element exists outside the group. More precisely, elements that do not belong to any group are not taken into account by the Mesh class. They can still be accessed by the functions of Salome, but the use of such lower-level abstractions is against the philosophy of encapsulation this module provides.

The non-goals of this module include everything not related to the tessellation nature of 2D microstructures. The emphasis is on readability and not on speed.

This file can be used either as a module or as a script.

Using as a module

Developers, implementing new features, will use salome.py as a Python module. Although this module is part of the grains package in the CristalX project, it does not rely on the other modules of grains. This ensures that it can be used standalone and run as a script in the Salome environment. It only uses language constructs available in Python 3.6, and external packages shipped with Salome 9.4.0 onwards. To enable debugging, code completion and other useful development methods, consult with the documentation on … Most classes contain a protected variable that holds the underlying Salome object. Unless you debug, it is not necessary to directly deal with Salome objects programmatically.

Using as a script

When salome.py is run as a script, the contents in the if __name__ == "__main__": block is executed. Edit it to suit your needs. Similarly to the case when used as a module, the script can only be run from Salome’s own Python interpreter: either from the shell or from the GUI. To run it from the shell (including the GUI’s built-in Python command prompt), type

exec(open("<path_to_CristalX>/grains/salome.py", "rb").read())

You can also execute the script from the GUI by clicking on File ‣ Load Script…

Classes

Geometry	Represents the geometrical entities of a two-dimensional tessellated domain.
Face	Closed part of a plane.
Edge	A shape corresponding to a curve, and bounded by a vertex at each extremity.
Interface	An edge between two faces.
Mesh	Performs mesh manipulations on a tessellated geometry.
FaceMesh	Mesh on a face, part of the whole mesh.
InterfaceMesh	Mesh on an interface, part of the whole mesh.
CohesiveZone	Constructs zero-thickness elements along the interfaces.
GUI	Using GUI-related functionalities in Salome.

class grains.salome.CohesiveZone(mesh)

Bases: object

Constructs zero-thickness elements along the interfaces.

Parameters

mesh (Mesh) – Mesh into which the cohesive elements will be inserted.

_affected_elements(interface_mesh)

Face elements whose nodes must be renumbered when duplicating an interface mesh.

Parameters

interface_mesh (InterfaceMesh) – The original interface mesh that will be duplicated.

Returns

elements (set) – Elements of the mesh that require node renumbering.

See also

_enrich_interfaces()

_correct_junction_nodes()

Post-processing to handle inconsistent interface nodes at the junctions.

The interface-wise creation of new edge elements in _affected_elements() may result in edge element nodes that do not connect the opposite face element nodes on the two sides of the interface. This function checks the edge element nodes at the junctions and renumbers them so that they hold the same label as the face element nodes they connect to.

Returns

None.

See also

decouple_faces(), _affected_elements(), smeshBuilder.Mesh.FindCoincidentNodesOnPart(), smeshBuilder.Mesh.ChangeElemNodes()

_enrich_interfaces()

Inserts new interface elements and nodes into the mesh.

Although Salome has built-in functionality for duplicating nodes and creating elements, even accessible from the GUI with Modification -> Transformation -> Duplicate Nodes and/or Elements, it does not work with multiple intersecting interfaces or for closed interfaces. The reason is that the first step of the two-step procedure Salome performs fails in such situations. Therefore, this method uses a modified algorithm for the first step, and then calls the second step. These steps are the following:

1. Find the elements (called affected elements) in the mesh whose node numbers need to be changed due to the topological changes in the mesh caused by the introduction of new nodes.

2. The affected elements are fed to an existing function in Salome, which returns the 1D elements it creates from the duplicated nodes. The new interface mesh is stored in the CohesiveZone object.

Returns

None.

See also

decouple_faces(), _affected_elements(), smeshBuilder.Mesh.DoubleNodeElemGroups(), smeshBuilder.Mesh.MakeGroupByIds()

_generate_cohesive_element(bottom_element, top_element)

Creates a zero-thickness quadrilateral element.

Parameters

• bottom_element (int) – Edge element that will form the bottom edge of the cohesive element.

• top_element (int) – Edge element that will form the top edge of the cohesive element. It is assumed that the top element geometrically overlaps with the bottom element.

Returns

list of int – The four nodes of the cohesive element, numbered counter-clockwise. The node ordering adheres to the node numbering in Abaqus.

See also

Mesh.incident_elements(), Mesh.element_edge_normal(), smeshBuilder.Mesh.GetElemNodes(), smeshBuilder.Mesh.GetNodeXYZ()

create_cohesive_elements()

Creates zero-thickness quadrilateral elements along the interfaces.

It is necessary that the mesh has already been decoupled along the interfaces by the decouple_faces() method. That method introduced duplicated nodes and edge elements along the interfaces. The purpose of this method is to tie each interface (edge) element to its corresponding duplicate in order to form a four-noded zero-thickness element, referred to as cohesive element. The bottom edge of the new cohesive element corresponds to the original edge element, while its top edge is formed by the duplicated interface edge element.

Returns

cohesive_elements (list) – List of nodes that form the cohesive elements. The node numbering follows the node ordering of Salome, which is the same as the node ordering in Abaqus.

See also

decouple_faces(), _generate_cohesive_element(), smeshBuilder.Mesh.AddFace()

decouple_faces()

Decouples the face meshes along the interfaces.

The algorithm consists of two main steps. First, new interface meshes are created that overlap with the existing ones and contain independent nodes and interface elements. In the same step, the incident face mesh nodes are updated to reflect the topological changes. However, in this method, extra nodes are introduced at the junctions, leading to a kinematic inconsistency. Therefore, the extraneous interface mesh nodes are renumbered in the second step of the algorithm.

Returns

None.

See also

create_cohesive_elements()

class grains.salome.Edge(edge, name)

Bases: object

A shape corresponding to a curve, and bounded by a vertex at each extremity.

An Edge knows about the Face it is part of, and the faces neighboring it.

Parameters

• edge (GEOM_Object of shape type ‘EDGE’) – The main Salome object wrapped by this class.

• name (str) – Name of the edge.

See also

GEOM_Object, Shape type

length()

Length of the edge.

Returns

float – Length of the edge.

See also

geomBuilder.BasicProperties

class grains.salome.Face(face, name)

Bases: object

Closed part of a plane.

A Face knows about the Edge objects that bound it.

Parameters

• face (GEOM_Object of shape type ‘FACE’) – The main Salome object wrapped by this class.

• name (str) – Name of the face.

See also

GEOM_Object, Shape type

class grains.salome.FaceMesh(face_mesh, name, on_face)

Bases: object

Mesh on a face, part of the whole mesh.

Parameters

• face_mesh (smeshBuilder.Mesh.GroupOnGeom) – The main Salome object wrapped by this class.

• name (str) – Name of the face mesh.

• on_face (Face) – Geometrical face on which this mesh exists.

elements()

Retrieves the elements of the face mesh.

Returns

list of int – Elements belonging to the face mesh.

See also

SMESH.SMESH_IDSource.GetIDs()

nodes()

Retrieves the nodes of the face mesh.

Returns

list of int – Nodes belonging to the face mesh.

See also

SMESH.SMESH_GroupBase.GetNodeIDs()

class grains.salome.GUI

Bases: object

Using GUI-related functionalities in Salome.

Notes

A part of Salome’s GUI is exposed to Python. To get an idea of what is available, see https://docs.salome-platform.org/latest/gui/GUI/text_user_interface.html

exception SalomeNoDesktop

Bases: Exception

Raised when Salome is run without desktop, but a desktop functionality is invoked.

classmethod _get_component(obj)

Determines the component of an object.

This function maps a class of this module to the Salome module the class uses. For instance, class Face is mapped to ‘GEOM’.

Parameters

obj – Any object for which the component name is looked for.

Returns

str or None – The name of the component the object belongs to. If an object of an unsupported class is given, None is returned. For the list of supported classes, see the component_map member of GUI class.

See also

show()

static assert_salome_desktop()

Checks if Salome’s GUI is available, and raises an exception if it is not.

This function acts as a helper function when relying on Salome’s GUI.

Raises

SalomeNoDesktop – If Salome’s GUI is not available.

component_map = {<class 'grains.salome.Geometry'>: 'GEOM', <class 'grains.salome.Face'>: 'GEOM', <class 'grains.salome.Edge'>: 'GEOM', <class 'grains.salome.Interface'>: 'GEOM', <class 'grains.salome.Mesh'>: 'SMESH', <class 'grains.salome.FaceMesh'>: 'SMESH', <class 'grains.salome.InterfaceMesh'>: 'SMESH'}

static has_desktop()

Indicates if the Salome GUI is running.

Returns

bool – True if Salome’s GUI is available, False otherwise.

classmethod show(obj, show_only=False)

Shows objects in Salome’s GUI.

Todo

Support list of objects.

Parameters

• obj (iterable) – The object(s) to be shown in Salome. Objects of the following classes are supported: Geometry, Face, Edge, Interface, Mesh, FaceMesh, InterfaceMesh.

• show_only (bool, optional) – If True, the other objects are hidden. The default value is False.

Returns

None

Raises

• SalomeNoDesktop – If Salome’s GUI is not available.

• TypeError – If obj is not an object that can be displayed.

• ValueError – If the given object does not exist in the Salome study.

Examples

For example, you can display an interface mesh and a face mesh by calling

GUI.show([interface_mesh, face_mesh])

where interface_mesh and face_mesh are InterfaceMesh and FaceMesh objects respectively. This way of using the show method provides great flexibility as different types of objects can be handled at the same time.

static update_object_browser()

Refreshes Salome’s object browser.

Only makes sense if executed with the GUI enabled.

Returns

None

See also

has_desktop()

classmethod view(view='top')

Sets the viewpoint.

Parameters

view ({‘front’, ‘back’, ‘top’, ‘bottom’, ‘left’, ‘right’}, optional) – Position from which the scene is viewed. The default is ‘top’.

Returns

None

class grains.salome.Geometry(name='microstructure')

Bases: object

Represents the geometrical entities of a two-dimensional tessellated domain.

A Geometry object knows about the faces that tessellate the domain, and about the edges and interfaces that separate the faces.

Parameters

name (str, optional) – Name of the geometry.

See also

Face, Edge, Interface

_find_overlapping_edges()

Finds edges that are on top of each other.

Overlapping edges have the same length. Although its converse is not true in general, we will assume so. This part of the algorithm (i.e. deciding the overlapping edges) can later be refined.

Returns

overlapping_edges (list) – Each member of the list contains a list of (supposedly) two Edge objects.

See also

create_interfaces(), Edge.length()

static _has_smaller_ID(edges)

Selects the edge with a smaller ID.

When generating mesh with the NETGEN plugin of Salome, if several edges overlap, only the edge with the smallest ID holds a mesh. The purpose of this function is to find the edge with the smaller ID out of two overlapping edges.

Parameters

edges (list of Edge) – List of two Edge objects.

Returns

Edge – Either the first or the second element of the input list, depending on which of them has a smaller ID.

See also

create_interfaces()

Notes

For a detailed discussion on this highly important issue, see the corresponding forum thread.

create_interfaces()

Constructs unique interfaces that separate the faces.

Based on the edges (obtained by exploding the mesh), interfaces are created. Interfaces are unique separators of two neighboring faces. In other words,

• if the edge is a boundary edge, no interface is created,

• two neighboring faces have two overlapping edges, of which one is defined to be an interface.

It is assumed that the edges and faces of the geometry have already been obtained.

Returns

None.

See also

extract_faces(), extract_edges()

extract_edges()

Decomposes each face into edges.

This method must be called after the extract_faces() method, otherwise it has no effect.

Returns

None.

See also

extract_faces()

extract_faces()

Decomposes the geometry into faces.

Returns

None.

See also

extract_edges()

load(step_file)

Loads the geometry from a STEP file.

Parameters

step_file (str) – The STEP file containing the geometry.

Returns

None.

class grains.salome.Interface(edge, name, neighboring_faces)

Bases: object

An edge between two faces.

Similar to an Edge, but two neighboring Face objects share a common Interface. An Interface knows about the two faces that it separates.

Parameters

• edge (GEOM_Object of shape type ‘EDGE’) – The main Salome object wrapped by this class.

• name (str) – Name of the interface.

• neighboring_faces (list of Face) – The two neighboring faces.

See also

Edge, Face

length()

Length of the interface.

Returns

float – Length of the interface.

See also

geomBuilder.BasicProperties

class grains.salome.InterfaceMesh(interface_mesh, name, on_interface)

Bases: object

Mesh on an interface, part of the whole mesh.

Parameters

• interface_mesh (smeshBuilder.Mesh.GroupOnGeom) – The main Salome object wrapped by this class.

• name (str) – Name of the interface mesh.

• on_interface (Interface) – Interface on which this mesh exists.

elements()

Retrieves the elements of the interface mesh.

Returns

list of int – Elements belonging to the interface mesh.

See also

SMESH.SMESH_IDSource.GetIDs()

elements_by_nodes(nodes)

Connecting elements to given nodes.

Parameters

nodes (list of int) – Nodes for which we want to find the connecting elements.

Returns

list of int – Elements that are incident to the given nodes.

See also

smeshBuilder.Mesh.GetElementsByNodes()

endpoint_nodes()

Nodes at the extremities of the interface mesh.

Returns

ep_nodes (list of int) – Nodes of the end points of the interface on which the interface mesh is defined. If the interface is open, it has two end points. When closed, the two end points coincide and instead of the two coinciding nodes, a single node is returned.

nodes()

Retrieves the nodes of the interface mesh.

Returns

list of int – Nodes belonging to the interface mesh.

See also

SMESH.SMESH_GroupBase.GetNodeIDs()

class grains.salome.Mesh(geometry, name='Mesh')

Bases: object

Performs mesh manipulations on a tessellated geometry.

Parameters

• geometry (Geometry) – Geometry object on which the mesh exists.

• name (str, optional) – Name of the mesh.

See also

Geometry, FaceMesh, InterfaceMesh

class ElementType

Bases: enum.Enum

Subset of the element types recognized by Salome.

This enumeration is for convenience. Only those elements of Salome are considered that are relevant for the Mesh class.

See also

SMESH.ElementType

ALL

EDGE

FACE

NODE

element_edge_normal(element, edge)

Outward-pointing unit normal to an element edge.

The edge is assumed to be planar.

Parameters

• element (int) – ID of the element.

• edge (list of int) – Edge of the element for which the normal is to be found. The edge is given by its two nodes.

Returns

normal (tuple of float) – Outward-pointing unit normal.

See also

point_in_element()

generate()

Generates a mesh on the geometry.

Todo

Do not hardcode values and explain the need for consistent orientation.

generate_element_nodes(elements)

Nodes of selected elements, returned one at a time.

Parameters

elements (iterable) – Element IDs.

Yields

list of int – The first entry of the list is the element ID, the remaining entries are the node IDs of the element.

incident_elements(edge, element_type=None)

Searches for elements incident to an edge.

Parameters

• edge (list of int) – An edge of an element, given by its two nodes.

• element_type (Mesh.ElementType, optional) – Perform the search for the given element type only.

Returns

list of int – Element IDs that are incident to the given edge.

See also

smeshBuilder.Mesh.GetNodeInverseElements()

incident_face_mesh(interface_mesh)

Face meshes incident to an interface mesh.

Todo

Use this method in _affected_elements as well.

Parameters

interface_mesh (InterfaceMesh) – Interface mesh for which the connecting face meshes are sought.

Returns

face_mesh (list of FaceMesh) – Face meshes incident to an interface mesh.

Raises

Exception – If no face mesh is incident to the interface mesh.

obtain_face_meshes()

Retrieves the elements of the mesh on each face.

Returns

None.

obtain_interface_meshes()

Obtains the 1D interfacial mesh for each interface.

Returns

None.

one_ring(node, definition='connecting')

Elements around a node.

Parameters

• node (int) – Node ID for which the one-ring is searched.

• definition ({‘connecting’, ‘surrounding’}, optional) – What you mean by neighboring elements. See the notes below.

Returns

list – List of integers (element IDs).

Notes

One should make a distinction between elements connecting to a node and elements surrounding a node. For a mesh with no overlapping nodes, the two definitions give the same elements. However, if multiple nodes are located at the same geometrical point, it can happen that the incident elements are not connected to the same node.

point_in_element(element, point)

Checks whether a point is in an element.

This method is implemented for 2D meshes only.

Parameters

• element (int) – Element of the mesh.

• point (tuple of float) – Point coordinates (x,y).

Returns

bool – True if the given element contains the given point.

Raises

• Exception – If the mesh is not two-dimensional.

• ValueError – If the element does not exist in the mesh.

Notes

This method calls an efficient matplotlib function to determine whether a point is in a polygon. For alternative implementations, see this discussion.

See also

matplotlib.path.Path.contains_point()

Geometry

This module implements computational geometry algorithms, needed for other modules.

All the examples assume that the modules numpy and matplotlib.pyplot were imported as np and plt, respectively.

Classes

Mesh	Data structure for a general mesh.
TriMesh	Unstructured triangular mesh.
Polygon	Represents a polygon.

Functions

is_collinear(points[, tol])	Decides whether a set of points is collinear.
squared_distance(x, y)	Squared Euclidean distance between two points.
distance_matrix(points)	A symmetric square matrix, containing the pairwise squared Euclidean distances among points.
_polygon_area(x, y)	Computes the signed area of a non-self-intersecting, possibly concave, polygon.

class grains.geometry.Mesh(vertices, cells)

Bases: abc.ABC

Data structure for a general mesh.

This class by no means wants to provide intricate functionalities and does not strive to be efficient at all. Its purpose is to give some useful features the project relies on. The two main entities in the mesh are the vertices and the cells. They are expected to be passed by the user, so reading from various mesh files is not implemented. This keeps the class simple, requires few package dependencies and keeps the class focused as there are powerful tools to convert among mesh formats (see e.g. meshio).

Parameters

• vertices (ndarray) – 2D numpy array with 2 columns, each row corresponding to a vertex, and the two columns giving the Cartesian coordinates of the vertex.

• cells (ndarray) – Cell-vertex connectivities in a 2D numpy array, in which each row corresponds to a cell and the columns are the vertices of the cells. It is assumed that all the cells have the same number of vertices.

See also

change_vertex_numbering

Notes

Although not necessary, it is highly recommended that the local vertex numbering in the cells are the same, either clockwise or counter-clockwise. Some methods, such as get_boundary() even requires it. If you are not sure whether the cells you provide have a consistent numbering, it is better to renumber them by calling the change_vertex_numbering() method.

static _ismatrix(array)

Decides whether the input is a matrix.

Parameters

array (ndarray) – Numpy array to be checked.

Returns

bool – True if the input is a 2D array. Otherwise, False.

See also

_isvector()

static _isvector(array)

Decides whether the input is a vector.

Parameters

array (ndarray) – Numpy array to be checked.

Returns

bool – True if the input is a 1D array or if it is a column or row vector. Otherwise, False.

See also

_ismatrix()

associate_field(vertex_values, name='field')

Associates a scalar, vector or tensor field to the nodes.

Only one field can be present at a time. If you want to use a new field, call this method again with the new field values, which will replace the previous ones.

Parameters

• vertex_values (ndarray) – Field values at the nodes.

• name (str, optional) – Name of the field. If not given, it will be ‘field’.

Returns

None

create_cell_set(name, cells)

Forms a group from a set of cells.

Parameters

• name (str) – Name of the cell set.

• cells (list) – List of cells to be added to the set.

Returns

None

create_vertex_set(name, vertices)

Forms a group from a set of vertices.

Parameters

• name (str) – Name of the vertex set.

• vertices (list) – List of vertices to be added to the set.

get_boundary()

Extracts the boundary of the mesh.

It is expected that all the cells have the same orientation, i.e. the cell vertices are consistently numbered (either clockwise or counter-clockwise). See the constructor for details.

Returns

• boundary_vertices (ndarray) – Ordered 1D ndarray of vertices, the boundary vertices of the mesh.

• boundary_edges (dict) – The keys of the returned dictionary are 2-tuples, representing the two vertices of the boundary edges, while the values are the list of cells containing a particular boundary edge. The dictionary is ordered: the consecutive keys represent the consecutive boundary edges. Although a boundary edge is part of a single cell, that cell is given in a list so as to maintain the same format as the one used in the get_edges() method.

Notes

The reason why consistent cell vertex numbering is demanded is because in that case the boundary edges are oriented in such a way that the second vertex of a boundary edge is the first vertex of the boundary edge it connects to.

Examples

Let us consider the same example mesh as the one described in the get_edges() method.

>>> mesh = TriMesh(np.array([[0, 0], [1, 0], [2, 0], [0, 2], [0, 1], [1, 1]]),
...                np.array([[0, 1, 5], [4, 5, 3], [5, 4, 0], [2, 5, 1]]))

We extract the boundary of that mesh using

>>> bnd_vertices, bnd_edges = mesh.get_boundary()
>>> bnd_vertices
array([0, 1, 2, 5, 3, 4])
>>> bnd_edges
{(0, 1): [0], (1, 2): [3], (2, 5): [3], (5, 3): [1], (3, 4): [1], (4, 0): [2]}

get_edges()

Constructs edge-cell connectivities of the mesh.

The cells of the mesh do not necessarily have to have a consistent vertex numbering.

Returns

edges (dict) – The keys of the returned dictionary are 2-tuples, representing the two vertices of the edges, while the values are the list of cells containing a particular edge.

Notes

We traverse through the cells of the mesh, and within each cell the edges. The edges are stored as new entries in a dictionary if they are not already stored. Checking if a key exists in a dictionary is performed in O(1). The number of edges in a cell is independent of the mesh density. Therefore, the time complexity of the algorithm is O(N), where N is the number of cells in the mesh.

See also

get_boundary()

Examples

We show an example for a triangular mesh (as the Mesh class is abstract).

>>> mesh = TriMesh(np.array([[0, 0], [1, 0], [2, 0], [0, 2], [0, 1], [1, 1]]),
...                np.array([[0, 1, 5], [4, 5, 3], [5, 4, 0], [2, 5, 1]]))
>>> edges = mesh.get_edges()
>>> edges  
{(0, 1): [0], (1, 5): [0, 3], (5, 0): [0, 2], (4, 5): [1, 2], (5, 3): [1],
(3, 4): [1], (4, 0): [2], (2, 5): [3], (1, 2): [3]}
>>> mesh.plot(cell_labels=True, vertex_labels=True)
>>> plt.show()

class grains.geometry.Polygon(vertices)

Bases: object

Represents a polygon.

This class works as expected as long as the given polygon is simple, i.e. it is not self-intersecting and does not contain holes.

A simple class that has numpy and matplotlib (for the visualization) as the only dependencies. It does not want to provide extensive functionalities (for those, check out Shapely). The polygon is represented by its vertices, given in a consecutive order.

Parameters

vertices (ndarray) – 2D numpy array with 2 columns, each row corresponding to a vertex, and the two columns giving the Cartesian coordinates of the vertex.

Raises

• Exception – If all the vertices of the polygon lie along the same line. If the polygon is not given in R^2.

• ValueError – If the polygon does not have at least 3 vertices.

Examples

Try to give a “polygon”, in which all vertices are collinear

>>> poly = Polygon(np.array([[0, 0], [1, 1], [2, 2]]))  
Traceback (most recent call last):
...
Exception: All vertices are collinear. Not a valid polygon.

Now we give a valid polygon:

>>> pentagon = Polygon(np.array([[2, 1], [0, 0], [0.5, 3], [-1, 4], [3, 5]]))

Use Python’s print function to display basic information about a polygon:

>>> print(pentagon)
A non-convex polygon with 5 vertices, oriented clockwise.

area()

Signed area of the polygon.

The signed area is computed by the shoelace formula 1

\[A = \frac{1}{2}\sum\limits_{i=1}^N (x_i y_{i+1} - x_{i+1} y_i)\]

Returns

float – Signed area.

References

1

http://paulbourke.net/geometry/polygonmesh/centroid.pdf

Examples

>>> poly = Polygon(np.array([[0, 0], [1, 0], [1, 1], [-1, 1]]))
>>> poly.area()
1.5
>>> poly = Polygon(np.array([[-1, 1], [1, 1], [1, 0], [0, 0]]))
>>> poly.area()
-1.5

centroid()

Centroid of the polygon.

The centroid is computed according to the following formula 2

\[ \begin{align}\begin{aligned}C_x = \frac{1}{6A}\sum\limits_{i=1}^N (x_i + x_{i+1})(x_i y_{i+1} - x_{i+1} y_i)\\C_y = \frac{1}{6A}\sum\limits_{i=1}^N (y_i + y_{i+1})(x_i y_{i+1} - x_{i+1} y_i)\end{aligned}\end{align} \]

where \(A\) is the signed area determined by the area() method and \(x_i, y_i\) are the vertex coordinates with \(x_{N+1}=x_1\) and \(y_{N+1}=y_1\).

Returns

tuple – 2-tuple, the coordinates of the centroid.

References

2

http://paulbourke.net/geometry/polygonmesh/centroid.pdf

Examples

>>> poly = Polygon(np.array([[0, 0], [0, 1], [1, 1], [1, 0]]))
>>> poly.centroid()
(0.5, 0.5)
>>> poly = Polygon(np.array([[2, 1], [0, 0], [0.5, 3], [-1, 4], [3, 5]]))
>>> poly.centroid()  
(1.254..., 2.807...)

diameter(definition='set')

Diameter of the polygon.

Multiple definitions are supported for the diameter:

• Diameter of a set. The polygon is considered as a set \(A\) of points comprised

of the polygon vertices. Let \((X,d)\) be a metric space. The diameter of the set is defined as

(1)\[\mathrm{diam}(A) = \sup\{ d(x,y)\ |\ x, y \in A\}.\]

Here, the Euclidean metric is used.

• Equivalent diameter. Diameter of the circle of the same area as that of the polygon.

Parameters

definition ({‘set’, ‘equivalent’}, optional) – The default is ‘set’.

Returns

float – Diameter of the polygon, based on the chosen definition.

Notes

When definition is ‘set’, computing the diameter by (1) is equivalent to determining the distance of the furthest points in the convex hull of \(A\). Therefore, the diameter will always be the distance between two points on the convex hull of \(A\). Then for each vertex of the hull finding which other hull vertex is farthest away from it, the rotating caliper algorithm can be used. Our brute-force method is simpler as it needs neither the convex hull nor the rotating caliper algorithm: all the pairwise distances among the polygon vertices are computed and the largest one is chosen. Pair of points a maximum distance apart. Since the polygons in our applications do not have that many vertices, this simplistic approach is a viable alternative.

Examples

>>> poly = Polygon(np.array([[2, 5], [0, 1], [4, 3], [4, 5]]))
>>> poly.diameter('set')  
5.6568542...
>>> poly.diameter('equivalent')  
3.1915382...
>>> poly = Polygon(np.array([[2, 1], [3, -4], [-1, -1], [-4, -2], [-3, 0]]))
>>> poly.diameter('set')  
7.2801098...
>>> poly.diameter('equivalent')  
4.2967398...

is_convex()

Decides whether the polygon is convex.

Returns

bool – True if the polygon is convex. False otherwise.

Notes

The algorithm works by checking if all pairs of consecutive edges in the polygon are either all clockwise or all counter-clockwise oriented. This method is valid only for simple polygons. The implementation follows this code, extended for the case when two consecutive edges are collinear. If the polygon was not simple, a more complicated algorithm would be needed, see e.g. here.

Examples

A triangle is always convex:

>>> poly = Polygon(np.array([[1, 1], [0, 1], [0, 0]]))
>>> poly.is_convex()
True

Let us define a concave deltoid:

>>> poly = Polygon(np.array([[-1, -1], [0, 1], [1, -1], [0, 5]]))
>>> poly.is_convex()
False

Give a polygon that has two collinear edges:

>>> poly = Polygon(np.array([[0.5, 0], [1, 0], [1, 1], [0, 1], [0, 0]]))
>>> poly.is_convex()
True

orientation()

Orientation of the polygon.

Returns

str – ‘cw’ if the polygon has clockwise orientation, ‘ccw’ if counter-clockwise.

plot(*args, **kwargs)

Plots the polygon.

Parameters

ax (matplotlib.axes.Axes, optional) – The Axes instance the polygon resides in. The default is None, in which case a new Axes within a new figure is created.

Other Parameters

• show_axes (bool, optional) – If True, the coordinate system is shown. The default is True.

• vertex_labels (bool, optional) – If True, vertex labels are shown. The default is False.

• args, kwargs (optional) – Additional arguments and keyword arguments to be specified. Those arguments are the ones supported by matplotlib.axes.Axes.plot().

Returns

None

Examples

Consider the pentagon used in the example of the constructor. Plot it in black with red diamond symbols representing its vertices. Moreover, display the vertex numbers and do not show the coordinate system.

>>> pentagon = Polygon(np.array([[2, 1], [0, 0], [0.5, 3], [-1, 4], [3, 5]]))
>>> pentagon.plot('k-d', vertex_labels=True, markerfacecolor='r', show_axes=False)
>>> plt.show()

plot_options = {'ax': None, 'show_axes': True, 'vertex_labels': False}

class grains.geometry.TriMesh(vertices, cells)

Bases: grains.geometry.Mesh

Unstructured triangular mesh.

Vertices and cells are both stored as numpy arrays. This makes the simple mesh manipulations easy and provides interoperability with the whole scientific Python stack.

Parameters

• vertices (ndarray) – 2D numpy array with 2 columns, each row corresponding to a vertex, and the two columns giving the Cartesian coordinates of the vertex.

• cells (ndarray) – Cell-vertex connectivities in a 2D numpy array, in which each row corresponds to a cell and the columns are the vertices of the cells. It is assumed that all the cells have the same number of vertices.

cell_area(cell)

Computes the area of a cell.

Parameters

cell (int) – Cell label.

Returns

area (float) – Area of the cell.

See also

cell_set_area()

cell_set_area(cell_set)

Computes the area of a cell set.

Parameters

cell_set (str) – Name of the cell set.

Returns

area (float) – Area of the cell set.

See also

cell_area()

cell_set_to_mesh(cell_set)

Creates a mesh from a cell set.

The cell orientation is preserved. I.e. if the cells had a consistent orientation (clockwise or counter-clockwise), the cells of the new mesh inherit this property.

Parameters

cell_set (str) – Name of the cell set being used to construct the new mesh. The cell set must be present in the cell_sets member variable of the current mesh object.

Returns

TriMesh – A new TriMesh object, based on the selected cell set of the original mesh.

Notes

The implementation is based on https://stackoverflow.com/a/13572640/4892892.

Examples

>>> mesh = TriMesh(np.array([[0, 0], [1, 0], [0, 1], [1, 1]]),
...                np.array([[0, 1, 2], [1, 3, 2]]))
>>> mesh.create_cell_set('set', [1])
>>> new_mesh = mesh.cell_set_to_mesh('set')
>>> new_mesh.cells  # note that the vertices have been relabelled
array([[0, 2, 1]])
>>> new_mesh.vertices
array([[1, 0],
       [0, 1],
       [1, 1]])
>>> new_mesh.plot(cell_labels=True, vertex_labels=True)
>>> plt.show()

change_vertex_numbering(orientation, inplace=False)

Changes cell vertex numbering.

Parameters

• orientation ({‘ccw’, ‘cw’}) – Vertex numbering within a cell, either ‘ccw’ (counter-clockwise, default) or ‘cw’ (clock-wise).

• inplace (bool, optional) – If True, the vertex ordering is updated in the mesh object. The default is False.

Returns

reordered_cells (ndarray) – Same format as the cells member variable, with the requested vertex ordering.

Notes

Supposed to be used with planar P1 or Q1 finite s.

Examples

>>> mesh = TriMesh(np.array([[1, 1], [3, 5], [7,3]]), np.array([0, 1, 2]))
>>> mesh.change_vertex_numbering('ccw')
array([[2, 1, 0]])

plot(*args, **kwargs)

Plots the mesh.

Parameters

ax (matplotlib.axes.Axes, optional) – The Axes instance the mesh resides in. The default is None, in which case a new Axes within a new figure is created.

Other Parameters

• cell_sets, vertex_sets (bool, optional) – If True, the cell/vertex sets (if exist) are highlighted in random colors. The default is True.

• cell_legends, vertex_legends (bool, optional) – If True, cell/vertex set legends are shown. The default is False. For many sets, it is recommended to leave these options as False, otherwise the plotting becomes very slow.

• cell_labels, vertex_labels (bool, optional) – If True, cell/vertex labels are shown. The default is False. Recommended to be left False in case of many cells/vertices. Cell labels are positioned in the centroids of the cells.

• args, kwargs (optional) – Additional arguments and keyword arguments to be specified. Those arguments are the ones supported by matplotlib.axes.Axes.plot().

Returns

None

Notes

If you do not want to plot the cells, only the vertices, pass the '.' option, e.g.:

mesh.plot('k.')

to plot the vertices in black. Here, mesh is a TriMesh object.

Examples

A sample mesh is constructed by creating uniformly randomly distributed points on the rectangular domain [-1, 1] x [1, 2]. These points will constitute the vertices of the mesh, while its cells are the Delaunay triangles on the vertices.

>>> from grains.geometry import TriMesh
>>> msh = TriMesh(*TriMesh.sample_mesh(1))

The cells are drawn in greeen, in 3 points of line width, and the vertices of the mesh are shown in blue.

>>> msh.plot('go-', linewidth=3, markerfacecolor='b', vertex_labels=True)
>>> plt.show()

(Source code, png, hires.png, pdf)

Notes

The plotting is done by calling triplot(), which internally makes a deep copy of the triangles. This increases the memory usage in case of many elements.

plot_field(component, *args, show_mesh=True, **kwargs)

Plots a field on the mesh.

The aim of this method is to support basic post-processing for finite element visualization. Only the basic contour plot type is available. For vector or tensor fields, the components to be plotted must be chosen. For faster and more comprehensive plotting capabilities, turn to well-established scientific visualization software, such as ParaView or Mayavi. Another limitation of the plot_field() method is that field values are assumed to be associated to the vertices of the mesh, which restricts us to \(P1\) Lagrange elements.

Parameters

• component (int) – Positive integer, the selected component of the field to be plotted. Components are indexed from 0.

• show_mesh (bool, optional) – If True, the underlying mesh is shown. The default is True.

• ax (matplotlib.axes.Axes, optional) – The Axes instance the plot resides in. The default is None, in which case a new Axes within a new figure is created.

Other Parameters

See them described in the :meth:`plot` method.

Returns

None

See also

plot()

Examples

The following example considers the same type of mesh as in the example shown for plot().

>>> msh = TriMesh(*TriMesh.sample_mesh(1))

We pretend that the field is an analytical function, evaluated at the vertices.

>>> field = lambda x, y: 1 - (x + y**2) * np.sign(x)
>>> field = field(msh.vertices[:, 0], msh.vertices[:, 1])

We associate this field to the mesh and plot it with and without the mesh

>>> msh.associate_field(field, 'analytical field')
>>> _, (ax1, ax2) = plt.subplots(1, 2)
>>> msh.plot_field(0, 'bo-', ax=ax1, linewidth=1, markerfacecolor='k')
>>> msh.plot_field(0, ax=ax2, show_mesh=False)
>>> plt.show()

plot_options = {'ax': None, 'cell_labels': False, 'cell_legends': False, 'cell_sets': True, 'vertex_labels': False, 'vertex_legends': False, 'vertex_sets': True}

rotate(angle, point=(0, 0), inplace=False)

Rotates a 2D mesh about a given point by a given angle.

Parameters

• angle (float) – Angle of rotation, in radians.

• point (list or tuple, optional) – Coordinates of the point about which the mesh is rotated. If not given, it is the origin.

• inplace (bool, optional) – If True, the vertex positions are updated in the mesh object. The default is False.

Returns

rotated_vertices (ndarray) – Same format as the vertices member variable, with the requested rotation.

Notes

Rotating a point \(P\) given by its coordinates in the global coordinate system as \(P(x,y)\) around a point \(A(x,y)\) by an angle \(\alpha\) is done as follows.

1. The coordinates of \(P\) in the local coordinate system, the origin of which is \(A\), is expressed as

\[P(x',y') = P(x,y) - A(x,y).\]

2. The rotation is performed in the local coordinate system as \(P'(x',y') = RP(x',y')\), where \(R\) is the rotation matrix:

\[\begin{split}R = \begin{pmatrix} \cos\alpha & -\sin\alpha \\ \sin\alpha & \hphantom{-}\cos\alpha \end{pmatrix}.\end{split}\]

3. The rotated point \(P'\) is expressed in the original (global) coordinate system:

\[P'(x,y) = P'(x',y') + A(x,y).\]

static sample_mesh(sample, param=100)

Provides sample meshes.

Parameters

• sample (int) – Integer, giving the sample mesh to be considered. Possibilities:

• param – Parameters to the sample meshes. Possibilities:

Returns

• nodes (ndarray) – 2D numpy array with 2 columns, each row corresponding to a vertex, and the two columns giving the Cartesian coordinates of the vertices.

• cells (ndarray) – Cell-vertex connectivity in a 2D numpy array, in which each row corresponds to a cell and the columns are the vertices of the cells. It is assumed that all the cells have the same number of vertices.

scale(factor, inplace=False)

Scales the geometry by modifying the coordinates of the vertices.

Parameters

• factor (float) – Each vertex coordinate is multiplied by this non-negative number.

• inplace (bool, optional) – If True, the vertex positions are updated in the mesh object. The default is False.

Returns

None.

write_inp(filename)

Writes the mesh to an Abaqus .inp file.

Parameters

filename (str) – Path to the mesh file.

Returns

None

grains.geometry._polygon_area(x, y)

Computes the signed area of a non-self-intersecting, possibly concave, polygon.

Directly taken from http://rosettacode.org/wiki/Shoelace_formula_for_polygonal_area#Python

Parameters

x, y (list) – Coordinates of the consecutive vertices of the polygon.

Returns

float – Area of the polygon.

Warning

If numpy vectors are passed as inputs, the resulting area is incorrect! WHY?

Notes

The code is not optimized for speed and for numerical stability. Intended to be used to compute the area of finite element cells, in which case the numerical stability is not an issue (unless the cell is degenerate). As this function is called possibly as many times as the number of cells in the mesh, no input checking is performed.

Examples

>>> _polygon_area([0, 1, 1], [0, 0, 1])
0.5

grains.geometry.distance_matrix(points)

A symmetric square matrix, containing the pairwise squared Euclidean distances among points.

Parameters

points (ndarray) – 2D numpy array with 2 columns, each row corresponding to a point, and the two columns giving the Cartesian coordinates of the points.

Returns

dm (ndarray) – Distance matrix.

See also

squared_distance()

Examples

>>> points = np.array([[1, 1], [3, 0], [-1, -1]])
>>> distance_matrix(points)
array([[ 0.,  5.,  8.],
       [ 5.,  0., 17.],
       [ 8., 17.,  0.]])

grains.geometry.is_collinear(points, tol=None)

Decides whether a set of points is collinear.

Works in any dimensions.

Parameters

• points (ndarray) – 2D numpy array with N columns, each row corresponding to a point, and the N columns giving the Cartesian coordinates of the point.

• tol (float, optional) – Tolerance value passed to numpy’s matrix_rank function. This tolerance gives the threshold below which SVD values are considered zero.

Returns

bool – True for collinear points.

See also

numpy.linalg.matrix_rank()

Notes

The algorithm for three points is from Tim Davis.

Examples

Two points are always collinear

>>> is_collinear(np.array([[1, 0], [1, 5]]))
True

Three points in 3D which are supposed to be collinear (returns false due to numerical error)

>>> is_collinear(np.array([[0, 0, 0], [1, 1, 1], [5, 5, 5]]), tol=0)
False

The previous example with looser tolerance

>>> is_collinear(np.array([[0, 0, 0], [1, 1, 1], [5, 5, 5]]), tol=1e-14)
True

grains.geometry.squared_distance(x, y)

Squared Euclidean distance between two points.

For points \(x(x_1, ..., x_n)\) and \(y(y_1, ... y_n)\) the following metric is computed

\[\sum\limits_{i=1}^n (x_i - y_i)^2\]

Parameters

x, y (ndarray) – 1D numpy array, containing the coordinates of the two points.

Returns

float – Squared Euclidean distance.

See also

distance_matrix()

Examples

>>> squared_distance(np.array([0, 0, 0]), np.array([1, 1, 1]))
3.0

Abaqus

Warning

This module will substantially be rewritten. See this issue.

This module allows to create and manipulate Abaqus input files through the Abaqus keywords, thereby providing automation. Note that it is not intended to be a complete API to Abaqus. If you want fine control over the whole Abaqus ecosystem, consult with the Abaqus Scripting Reference Guide (ASRG). However, ASRG needs Abaqus to be installed, moreover, you must use the Python interpreter embedded into Abaqus. That version of Python is very old even in the latest versions of Abaqus. Furthermore, if you need to use custom Python packages for your work, chances are high that they will not work with the embedded interpreter, and may even crash the installation. To use this module, no Abaqus installation is needed. In fact, only functions from the Python standard library are used.

The documentation of Abaqus version 2017 is hosted on the following website: https://abaqus-docs.mit.edu/2017/English/SIMACAEEXCRefMap/simaexc-c-docproc.htm. Throughout the documentation of this module (grains.abaqus), we will make references to that website. If the links cease to exist, please let me know by opening an issue. Alternatively, once you have registered, you can browse the documentation on the official website.

Classes

Geometry	Geometrical operations on the mesh.
Material	Adds, removes, modifies materials.
Procedure	Handling analysis steps during the simulation.

Functions

extract(keyword)	Obtains Abaqus keyword and its parameters.
validate_file(file, caller)	Input or output file validation.

class grains.abaqus.Geometry

Bases: object

Geometrical operations on the mesh.

_Geometry__format()

Formats the material data in the Abaqus .inp format.

The internal representation of the material data in converted to a string understood by Abaqus.

abaqus_formatlist

List of strings, each element of the list corresponding to a line (with

line ending) in the Abaqus .inp file. In case of no

material, an empty list is returned.

The output is a list so that further concatenation operations are easy. If you want a string, merge the elements of the list:

output = ‘’.join(output)

This is what the show method does.

read(inp_file)

Reads material data from an Abaqus .inp file.

Parameters

inp_file (str) – Abaqus input (.inp) file to be created.

Returns

None.

Notes

• This method is designed to read material data. Although the logic could be used to process other properties (parts, assemblies, etc.) in an input file, they are not yet implemented in this class.

• This method assumes that the input file is valid. If it is, the material data can be extacted. If not, the behavior is undefined: the program can crash or return garbage. This is by design: the single responsibility principle dictates that the validity of the input file must be provided by other methods. If the input file was generated from within Abaqus CAE, it is guaranteed to be valid. The write method of this class also ensures that the resulting input file is valid. This design choice also makes the program logic simpler. For valid syntax in the input file, check the Input Syntax Rules section in the Abaqus user’s guide.

• To read material data from an input file, one has to identify the structure of .inp files in Abaqus. Abaqus is driven by keywords and corresponding data. For a list of accepted keywords, consult the Abaqus Keywords Reference Guide. There are three types of input lines in Abaqus:

  • keyword line: begins with a star, followed by the name of the keyword. Parameters, if any, are separated by commas and are given as parameter-value pairs. Keywords and parameters are not case sensitive. Example:

    *ELASTIC, TYPE=ISOTROPIC, DEPENDENCIES=1

    Some keywords can only be defined once another keyword has already been defined. E.g. the keyword ELASTIC must come after MATERIAL in a valid .inp file.

  • data line: immediately follows a keyword line. All data items must be separated by commas. Example:

    -12.345, 0.01, 5.2E-2, -1.2345E1

  • comment line: starts with ** and is ignored by Abaqus. Example:

    ** This is a comment line

• Internally, the materials are stored in a dictionary. It holds the material data read from the file. The keys in this dictionary are the names of the materials, and the values are dictionaries themselves. Each such dictionary stores a behavior for the given material. E.g. an elastoplastic material is governed by an elastic and a plastic behavior. The parameters for each behavior are stored in a list.

scale(factor)

Scales the geometry by modifying the coordinates of the nodes.

Parameters

factor (float) – Each nodal coordinate is multiplied by this non-negative number.

Returns

None.

Notes

The modification happens in-place.

write(output_file=None)

Writes material data to an Abaqus .inp file.

Parameters

output_file (str, optional) – Output file name to write the modifications into. If not given, the original file name is appended with ‘_mod’.

Returns

None.

Notes

• If the output file name is the same as the input file, the original .inp file will be overwritten. This is strongly not recommended.

• The whole content of the original input file is read to memory. It might be a problem for very large .inp files. In that case, a possible implementation could be the following:

  1. Remove old material data

  2. Append new material data to the proper position in the file

  Appending is automatically done at the end of the file. Moving the material description to the end of the file is not possible in general because defining materials cannot be done from any module, i.e. the *MATERIAL keyword cannot follow an arbitrary keyword. In this case, Abaqus throws an AbaqusException with the following message:

  It can be suboption for the following keyword(s)/level(s): model

class grains.abaqus.Material(from_Abaqus=False)

Bases: object

Adds, removes, modifies materials. Requirements: be able to - create an empty .inp file, containing only the materials - add materials to an existing .inp file

TODO: only document public attributes! .. attribute:: materials

materials, their behaviors and their parameters are stored here. Intended for internal representation. To view them, use the show method, to write them to an input file, use the write method.

type

dict

read(inp_file)

Reads material data from an Abaqus .inp file.

write(output_file=None)

Writes material data to an Abaqus .inp file.

remove(inp_file, output_file=None)

Removes material definitions from an Abaqus .inp file.

create(inp_file)

Creates empty Abaqus .inp file.

show()

Shows material data as it appears in the Abaqus .inp file.

Notes

The aim of this class is to programmatically add, remove, modify materials in a format understandable by Abaqus. This class does not target editing materials through a GUI (that can be done in Abaqus CAE).

_Material__format()

Formats the material data in the Abaqus .inp format.

The internal representation of the material data in converted to a string understood by Abaqus.

abaqus_formatlist

List of strings, each element of the list corresponding to a line (with

line ending) in the Abaqus .inp file. In case of no

material, an empty list is returned.

The output is a list so that further concatenation operations are easy. If you want a string, merge the elements of the list:

output = ‘’.join(output)

This is what the show method does.

static _Material__isnumeric(x)

Decides if the input is a scalar number.

Parameters

x (any type) – Input to be tested.

Returns

bool – True if the given object is a scalar number.

add_linearelastic(material, E, nu)

Adds linear elastic behavior to a given material.

Parameters

• material (str) – Name of the material the behavior belongs to.

• E (int, float) – Young’s modulus.

• nu (int, float) – Poisson’s ratio.

Returns

None.

add_material(material)

Defines a new material by its name.

Parameters

material (str) – Name of the material to be added. A material can have multiple behaviors (e.g. elastoplastic).

Returns

None.

add_plastic(material, sigma_y, epsilon_p)

Adds metal plasticity behavior to a given material.

Parameters

• material (str) – Name of the material the behavior belongs to.

• sigma_y (int, float) – Yield stress.

• epsilon_p (int, float) – Plastic strain.

Returns

None.

static add_sections(inp_file, output_file=None)

Adds section definitions to an Abaqus .inp file.

Defines properties for elements by associating materials to them. The element set containing the elements for which the material behavior is being defined is assumed to have the same name as that of the material. E.g. if materials with names mat-1 and mat-2 exist, element sets with names mat-1 and mat-2 must also exist. If such element sets do not exist, Abaqus will throw a warning and the section assignment will not be successful.

Parameters

• inp_file (str) – Abaqus .inp file from which the materials should be removed.

• output_file (str, optional) – Output file name to write the modifications into. If not given, the original file name is appended with ‘_mod’.

Returns

None.

Notes

If fine control is required for associating custom material names to custom element set names, that can be done from the Abaqus GUI. The purpose of this method is automation for large number of element sets, each associated to a (possibly distinct) material. In that case, custom element set names are not reasonable any more, and having the same names for the element sets and for the materials is completely meaningful.

create(inp_file)

Creates empty Abaqus .inp file.

Parameters

inp_file (str) – Abaqus input file to be created. If an extension is not given, the default .inp is used.

Returns

None.

read(inp_file)

Reads material data from an Abaqus .inp file.

Parameters

inp_file (str) – Abaqus input (.inp) file to be created.

Returns

None.

Notes

• This method is designed to read material data. Although the logic could be used to process other properties (parts, assemblies, etc.) in an input file, they are not yet implemented in this class.

• This method assumes that the input file is valid. If it is, the material data can be extacted. If not, the behavior is undefined: the program can crash or return garbage. This is by design: the single responsibility principle dictates that the validity of the input file must be provided by other methods. If the input file was generated from within Abaqus CAE, it is guaranteed to be valid. The write method of this class also ensures that the resulting input file is valid. This design choice also makes the program logic simpler. For valid syntax in the input file, check the Input Syntax Rules section in the Abaqus user’s guide.

• To read material data from an input file, one has to identify the structure of .inp files in Abaqus. Abaqus is driven by keywords and corresponding data. For a list of accepted keywords, consult the Abaqus Keywords Reference Guide. There are three types of input lines in Abaqus:

  • keyword line: begins with a star, followed by the name of the keyword. Parameters, if any, are separated by commas and are given as parameter-value pairs. Keywords and parameters are not case sensitive. Example:

    *ELASTIC, TYPE=ISOTROPIC, DEPENDENCIES=1

    Some keywords can only be defined once another keyword has already been defined. E.g. the keyword ELASTIC must come after MATERIAL in a valid .inp file.

  • data line: immediately follows a keyword line. All data items must be separated by commas. Example:

    -12.345, 0.01, 5.2E-2, -1.2345E1

  • comment line: starts with ** and is ignored by Abaqus. Example:

    ** This is a comment line

• Internally, the materials are stored in a dictionary. It holds the material data read from the file. The keys in this dictionary are the names of the materials, and the values are dictionaries themselves. Each such dictionary stores a behavior for the given material. E.g. an elastoplastic material is governed by an elastic and a plastic behavior. The parameters for each behavior are stored in a list.

static remove(inp_file, output_file=None)

Removes material definitions from an Abaqus .inp file.

Parameters

• inp_file (str) – Abaqus .inp file from which the materials should be removed.

• output_file (str, optional) – Output file name to write the modifications into. If not given, the original file name is appended with ‘_mod’.

Returns

None.

show()

Shows material data as it appears in the Abaqus .inp file.

Returns

None.

write(output_file=None)

Writes material data to an Abaqus .inp file.

Parameters

output_file (str, optional) – Output file name to write the modifications into. If not given, the original file name is appended with ‘_mod’.

Returns

None.

Notes

• If the output file name is the same as the input file, the original .inp file will be overwritten. This is strongly not recommended.

• The whole content of the original input file is read to memory. It might be a problem for very large .inp files. In that case, a possible implementation could be the following:

  1. Remove old material data

  2. Append new material data to the proper position in the file

  Appending is automatically done at the end of the file. Moving the material description to the end of the file is not possible in general because defining materials cannot be done from any module, i.e. the *MATERIAL keyword cannot follow an arbitrary keyword. In this case, Abaqus throws an AbaqusException with the following message:

  It can be suboption for the following keyword(s)/level(s): model

class grains.abaqus.Procedure(from_Abaqus=False)

Bases: object

Handling analysis steps during the simulation.

TODO: only document public attributes!

steps

analysis steps, each step collecting the necessary information. Intended for internal representation. To view them, use the show() method, to write them to an input file, use the write() method.

Type

dict

_Procedure__format()

Formats the data in the Abaqus .inp format.

The internal representation of the steps data in converted to a string understood by Abaqus.

abaqus_formatlist

List of strings, each element of the list corresponding to a line (with

line ending)

in the Abaqus .inp file. In case of no data, an empty list is returned.

The output is a list so that further concatenation operations are easy. If you want a string, merge the elements of the list: output = ''.join(output) This is what the show() method does.

add_analysis(step, time_period=1.0, initial_increment=None, min_increment=0, max_increment=None)

Adds an analysis type to a given step.

Note

Currently only static stress/displacement analysis is supported, indicated by the STATIC Abaqus keyword. No optional parameters are supported.

Todo

Check user inputs. Define exceptions (ValueError class) with the error messages Abaqus CAE throws.

Parameters

• step (str) – Name of the step the analysis type belongs to.

• initial_increment (float, optional) – Initial time increment. This value will be modified as required if the automatic time stepping scheme is used. If this entry is zero or is not specified, a default value that is equal to the total time period of the step is assumed.

• time_period (float, optional) – Time period of the step. The default is 1.0.

• min_increment (float, optional) – Only used for automatic time incrementation. If ABAQUS/Standard finds it needs a smaller time increment than this value, the analysis is terminated. If this entry is zero, a default value of the smaller of the suggested initial time increment or 1e–5 times the total time period is assumed.

• max_increment (float, optional) – Maximum time increment allowed. Only used for automatic time incrementation. If this value is not specified, no upper limit is imposed, i.e. max_increment = time_period.

Returns

None

Examples

>>> proc = Procedure()
>>> proc.create_step('Step-1', max_increments=1000)
>>> proc.add_analysis('Step-1', initial_increment=0.001, min_increment=1e-8)

add_boundary_condition(name, step, nodes, first_dof, last_dof=None, magnitude=0.0)

Adds boundary condition to a given step.

Boundary conditions can be prescribed on node sets, using the Abaqus keyword BOUNDARY. Optional parameters are not supported.

Multiple boundary conditions can be given on the same node set. It is the responsibility of the user to ensure that the constraints are compatible.

Note

Currently, boundary conditions can only be defined in a user-defined step, and not in the initial step.

Parameters

• name (str) – Name of the boundary condition to be added. Boundary condition names must be unique.

• step (str) – Name of the step the boundary condition belongs to.

• nodes (str or int) – If a string, the label of the node set, if a positive integer, the label of a single node on which the boundary condition is prescribed.

• first_dof (int) – First degree of freedom constrained.

• last_dof (int, optional) – Last degree of freedom constrained. Not necessary to be given if only one degree of freedom is being constrained.

• magnitude (float) – Magnitude of the variable. If this magnitude is a rotation, it must be given in radians. The default value is 1.0.

Returns

None

Examples

Let us consider a two-dimensional structure, fixed on a part of its boundary (called ‘left’), and displaced on another part (called ‘right’). First, we create an analysis step (named ‘Step-1’).

>>> proc = Procedure()
>>> proc.create_step('Step-1', max_increments=1000)

Now, we prescribe the zero displacements.

>>> proc.add_boundary_condition('fixed', 'Step-1', 'left', first_dof=1, last_dof=2)

Since we did not specify the magnitude of the displacement components, they are zero by default. The other boundary condition prescribes different values for the horizontal and vertical components. Hence, they are given separately.

>>> proc.add_boundary_condition('pulled_right', 'Step-1', 'right', first_dof=1, magnitude=2)
>>> proc.add_boundary_condition('fixed_right', 'Step-1', 'right', first_dof=2, magnitude=0)

Note that we can also prescribe boundary conditions at a particular node, e.g.

>>> proc.add_boundary_condition('at_node', 'Step-1', 1, first_dof=2, magnitude=-1.2)

where we set the vertical displacement component to -1.2 at node 1.

create_step(name, nlgeom=False, max_increments=100)

Defines a new step.

A step is the fundamental part of the simulation workflow in Abaqus. Among the available options, only the NLGEOM and INC options are supported currently.

Todo

Check in the write() method whether every step contains an analysis, as this is required by Abaqus.

Parameters

• name (str) – Name of the analysis step to be added. Step names must be unique.

• nlgeom (bool, optional) – Omit this parameter or set to False to perform a geometrically linear analysis during the current step. Set it to True to indicate that geometric nonlinearity should be accounted for during the step. Once the nlgeom option has been switched on, it will be active during all subsequent steps in the analysis. The default is False.

• max_increments (int) – The analysis will stop if the maximum number of increments is exceeded before the complete solution for the step has been obtained. The default is 100.

Returns

None

read(inp_file)

Reads procedure data from an Abaqus .inp file.

Parameters

inp_file (str) – Abaqus input (.inp) file to be created.

Returns

None.

Notes

• This method is designed to read material data. Although the logic could be used to process other properties (parts, assemblies, etc.) in an input file, they are not yet implemented in this class.

• This method assumes that the input file is valid. If it is, the material data can be extacted. If not, the behavior is undefined: the program can crash or return garbage. This is by design: the single responsibility principle dictates that the validity of the input file must be provided by other methods. If the input file was generated from within Abaqus CAE, it is guaranteed to be valid. The write method of this class also ensures that the resulting input file is valid. This design choice also makes the program logic simpler. For valid syntax in the input file, check the Input Syntax Rules section in the Abaqus user’s guide.

• To read material data from an input file, one has to identify the structure of .inp files in Abaqus. Abaqus is driven by keywords and corresponding data. For a list of accepted keywords, consult the Abaqus Keywords Reference Guide. There are three types of input lines in Abaqus:

  • keyword line: begins with a star, followed by the name of the keyword. Parameters, if any, are separated by commas and are given as parameter-value pairs. Keywords and parameters are not case sensitive. Example:

    *ELASTIC, TYPE=ISOTROPIC, DEPENDENCIES=1

    Some keywords can only be defined once another keyword has already been defined. E.g. the keyword ELASTIC must come after MATERIAL in a valid .inp file.

  • data line: immediately follows a keyword line. All data items must be separated by commas. Example:

    -12.345, 0.01, 5.2E-2, -1.2345E1

  • comment line: starts with ** and is ignored by Abaqus. Example:

    ** This is a comment line

• Internally, the materials are stored in a dictionary. It holds the material data read from the file. The keys in this dictionary are the names of the materials, and the values are dictionaries themselves. Each such dictionary stores a behavior for the given material. E.g. an elastoplastic material is governed by an elastic and a plastic behavior. The parameters for each behavior are stored in a list.

show()

Shows the text for the step module, as it appears in the Abaqus .inp file.

Returns

None.

write(output_file=None)

Writes step definition data to an Abaqus .inp file.

Parameters

output_file (str, optional) – Output file name to write the modifications into. If not given, the original file name is appended with ‘_mod’.

Returns

None.

Notes

• If the output file name is the same as the input file, the original .inp file will be overwritten. This is strongly not recommended.

• The whole content of the original input file is read to memory. It might be a problem for very large .inp files. In that case, a possible implementation could be the following:

  1. Remove old material data

  2. Append new material data to the proper position in the file

  Appending is automatically done at the end of the file. Moving the material description to the end of the file is not possible in general because defining materials cannot be done from any module, i.e. the *MATERIAL keyword cannot follow an arbitrary keyword. In this case, Abaqus throws an AbaqusException with the following message:

  It can be suboption for the following keyword(s)/level(s): model

grains.abaqus.extract(keyword)

Obtains Abaqus keyword and its parameters.

Parameters

keyword (str) –

Some examples:

‘*Elastic, type=ORTHOTROPIC’ ‘*Damage Initiation, criterion=HASHIN’ ‘*Nset, nset=Set-1, generate’

Returns

separated (list) – DESCRIPTION.

grains.abaqus.validate_file(file, caller)

Input or output file validation.

Parameters

• file (str) – Existing Abaqus .inp file or a new .inp file to be created.

• caller ({‘read’, ‘write’, ‘create’}) – Method name that called this function.

Returns

None.

DIC

This module is used to process field values obtained by digital image correlation (DIC). Plotting, comparison with numerical solutions, etc. are implemented.

Classes

DIC	Performs operations on digital image correlation data.

Functions

plot_strain(strain[, minval, maxval, …])

Plots a scalar strain field.

class grains.dic.DIC(u, v)

Bases: object

Performs operations on digital image correlation data.

Parameters

u, v (ndarray) – The two components of the displacement field, discretized on a rectangular grid.

Raises

ValueError – If the displacement components are not discretized on the same grid or if the grid size is not at least 3-by-3. This latter requirement is not a problem in practice (which camera is not capable of shooting photos in 9 pixels resolution?), while it allows the evaluation of the numerical derivatives with second order accuracy on the boundaries too.

Notes

Throughout the class methods, the following terms are used. When not precised, image refers to the DIC data plotted as an image. We can perceive the image as a pixel grid, in which the vertices are the centres of the pixels. An image coordinate system \((X,Y)\) can be defined on this grid such that the vertices have integer coordinates, the origin \(A\) is at the top left hand corner, and the \(Y\)-axis points towards the right, while the \(X\)-axis points to the bottom. When the experimental result is compared with a numerical solution (assumed to be available at the nodes of a mesh), one needs to map values defined on the DIC grid onto the nodes of the mesh, and vice versa. The mesh exists on the physical domain, for which a physical coordinate system \((x,y)\) is attached to. Hence, there is a need to express \((X,Y)\) in terms of \((x,y)\). Let \(A(a_x,a_y)\) be the origin of the \((X,Y)\) coordinate system given in terms of the \(x,y\) coordinates, and \(s[\mathrm{pixel/mm}]\) is the scale for the DIC image. The coordinate transformation is then given by

\[\begin{split}\begin{gather} x = a_x + \frac{Y}{s} \\ y = a_y - \frac{X}{s} \end{gather}\end{split}\]

The DIC grid in the physical coordinate system \((x,y)\) is called physical grid.

It should be noted that the special alignment of \((x,y)\) with respect to \((X,Y)\) is not a major constraint. In practical applications (when a Cartesian coordinate system is used), the orientation of \((x,y)\) in this way is a convention.

static equivalent_strain(strain_tensor, strain_measure)

Computes a scalar quantity from the strain tensor.

Parameters

• strain_tensor (ndarray) – Components of the strain tensor \(\boldsymbol{\varepsilon}\) as an \(m \times n \times 3\) array. The first two dimensions correspond to the grid points they are determined at, the third dimension gives the components of the tensor in the following order: \(\varepsilon_{11}, \varepsilon_{12}, \varepsilon_{22}\).

• strain_measure ({‘von Mises’}) – One of the following strain measures.

  • ‘von Mises’

    \[\varepsilon_M := \sqrt{\frac{2}{3} \boldsymbol{\varepsilon}^{\mathrm{dev}} : \boldsymbol{\varepsilon}^{\mathrm{dev}}} = \sqrt{\frac{2}{3} \left( \boldsymbol{\varepsilon}:\boldsymbol{\varepsilon} - \frac{1}{3}(\mathrm{tr} \boldsymbol{\varepsilon})^2 \right)}\]

    where \(\boldsymbol{\varepsilon}^{\mathrm{dev}}\) is the deviatoric part of the strain tensor. This can further be simplified (see the notes below) to

    \[\varepsilon_M = \frac{2}{3}\sqrt{\varepsilon_{11}^2 + \varepsilon_{22}^2 + 3\varepsilon_{12}^2 - \varepsilon_{11}\varepsilon_{22}}\]

Returns

ndarray – Equivalent strain at the grid points, given as an \(m \times n\) numpy array.

See also

strain()

Notes

Under plane stress conditions, \(\varepsilon_{33}\) is not zero, but it is computed as

\[\varepsilon_{33} = -\frac{\nu}{E}(\sigma_{11} + \sigma_{22})\]

Since the stresses are not known from the DIC, one must settle with plane strain conditions where \(\varepsilon_{33} = 0\). This is what we follow in the DIC class. We remark that in stereo-DIC multiple cameras are used, which allows the measurement of \(\varepsilon_{33}\).

plot_displacement(component, ax=None)

Plots a component of the displacement field.

Parameters

• component ({1, 2}) – Component to be plotted, where 1 corresponds to the first, 2 corresponds to the second component of the displacement field.

• ax (matplotlib.axes.Axes, optional) – The Axes instance the grid resides in. The default is None, in which case a new Axes within a new figure is created.

Returns

None

See also

plot_strain()

plot_physicalgrid(ax=None)

Plots the DIC grid in the physical coordinate system.

This method is only available once the relation between the image coordinate system and the physical coordinate system has been set up by the set_transformation() method.

Parameters

ax (matplotlib.axes.Axes, optional) – The Axes instance the grid resides in. The default is None, in which case a new Axes within a new figure is created.

Returns

ax (matplotlib.axes.Axes)

See also

plot_pixelgrid()

Notes

See in the plot_pixelgrid() method.

plot_pixelgrid(ax=None)

Plots the DIC grid in the image coordinate system.

Parameters

ax (matplotlib.axes.Axes, optional) – The Axes instance the grid resides in. The default is None, in which case a new Axes within a new figure is created.

Returns

ax (matplotlib.axes.Axes)

See also

plot_physicalgrid()

Notes

For a DIC image with \(m \times n\) number of pixels, the number of grid lines is \((m + 1) + (n + 1)\). Plotting all these lines would not only slow down the program, it would also make the grid practically indistinguishable from a filled rectangle. Therefore, for high resolution images, only grid lines around the boundary of the image are plotted. The target resolution above which this strategy is used can be given in the class constructor.

plot_strain(component, minval=None, maxval=None, colorbar=True)

Plots a component of the infinitesimal strain tensor.

The partial derivatives of the displacement field are computed with numerical differentiation.

Parameters

• component (tuple, {(1,1), (1,2), (2,1), (2,2)}) – Component to be plotted, where

  • (1,1) denotes \(\varepsilon_{11}\)

  • (1,2) and (2,1) denote \(\varepsilon_{12} = \varepsilon_{21}\)

  • (2,2) denotes \(\varepsilon_{22}\)

  for the infinitesimal strain tensor

  \[\begin{split}\varepsilon = \begin{pmatrix} \varepsilon_{11} & \varepsilon_{12} \\ \varepsilon_{21} & \varepsilon_{22} \end{pmatrix}\end{split}\]

• minval, maxval (float, optional) – Set the color scaling for the image by fixing the values that map to the colormap color limits. If minval is not provided, the limit is determined from the minimum value of the data. Similarly for maxval.

• colorbar (bool, optional) – If True, a horizontal colorbar is placed below the strain map. The default is True.

Returns

None

See also

plot_displacement(), numpy.gradient()

Deprecated since version 1.1.0: This will be removed in 1.3.0. This method will be modified to perform the plotting only. The computation of the various strain measures will be done in the strain() method. The function that replaces this function will keep the same name and is currently implemented outside this class.

plot_superimposedmesh(mesh, *args, **kwargs)

Plots the DIC grid with a mesh superimposed on it.

Since the finite element mesh represents a physical domain, the DIC grid is plotted in the physical coordinate system, which is determined by the set_transformation() method.

Note

Maybe no need to couple this module with the geometry module just for the sake of this function. It is probably easier to simply pass the nodes and the elements, and call triplot.

Parameters

• mesh (Mesh) – A subclass of grains.geometry.Mesh to be plotted along with the physical DIC grid.

• args (optional) – The same as for matplotlib.axes.Axes.plot(). The settings given here influence the mesh plotting.

• kwargs (optional) – The same as for matplotlib.axes.Axes.plot(). The settings given here influence the mesh plotting.

Returns

ax (matplotlib.axes.Axes) – The Axes object on which the plot is drawn.

See also

grains.geometry.TriMesh.plot()

Examples

Let us create a grid and a random mesh (the same as in the Examples section of the project_onto_mesh()).

>>> from grains.dic import DIC
>>> from grains.geometry import TriMesh
>>> x_grid, y_grid = np.mgrid[-1:1:100j, 1:2:50j]
>>> exact_solution = lambda x, y: 1 - (x + y**2) * np.sign(x)
>>> grid = DIC(np.random.rand(*np.shape(x_grid.T)), np.random.rand(*np.shape(x_grid.T)))
>>> grid.set_transformation((-1, 2), 50)
>>> n_nodes = 100  # modify this to see how good the interpolated solution is
>>> msh = TriMesh(*TriMesh.sample_mesh(1, n_nodes))
>>> grid.plot_superimposedmesh(msh, linewidth=3, markerfacecolor='b')
>>> plt.show()

(Source code, png, hires.png, pdf)

project_onto_grid(nodes, nodal_values, method='linear')

Project a numerically computed field onto the DIC grid.

Todo

Use the example code of this method to create the plotting method of this class. When done, rewrite the example with the methods.

Parameters

• nodes (ndarray) – 2D numpy array with 2 columns, each row corresponding to a node of the mesh, and the two columns giving the Cartesian coordinates of the nodes.

• nodal_values (ndarray) – 1D numpy array (vector), the numerically computed field at the nodes of the mesh.

• method ({‘nearest’, ‘linear’, ‘cubic’}, optional) – Type of the interpolation. The default is linear.

Returns

ndarray – 2D numpy array, the projected field values at the physical DIC grid points.

See also

project_onto_mesh(), scipy.interpolate.griddata()

Examples

The DIC grid and the nodal finite element values (interpolated in-between, which is an errenous way to visualize discontinuous functions) can be seen in the upper left and lower left figures, respectively.

>>> from grains.dic import DIC
>>> from grains.geometry import TriMesh
>>> x_grid, y_grid = np.mgrid[-1:1:100j, 1:2:50j]
>>> exact_solution = lambda x, y: 1 - (x + y**2) * np.sign(x)
>>> grid = DIC(exact_solution(x_grid, y_grid).T, exact_solution(x_grid,y_grid).T)
>>> grid.set_transformation((-1, 2), 50)

The finite element solution is obtained at the nodes (upper right figure) of the mesh. Here, we assumed that the nodes are sampled from a uniformly random distribution on \([-1,1] \times [1, 2]\).

>>> n_nodes = 100  # modify this to see how good the interpolated solution is

For the sake of this example, we also assume that the nodal values are “exact”, i.e. they admit the function values at the nodes. In reality, of course, this will not be the case, but this allows us to investigate the effect moving from the FE mesh to the DIC grid.

>>> mesh = TriMesh(*TriMesh.sample_mesh(1, n_nodes))
>>> fe_values = exact_solution(mesh.vertices[:, 0], mesh.vertices[:, 1])
>>> mesh.associate_field(fe_values)
>>> interpolated = grid.project_onto_grid(mesh.vertices, fe_values, method='linear')

The FE field available at the nodes are interpolated at the DIC grid points, as shown in the lower right figure. Note that no extrapolation is performed, so values are not computed at the grid points lying outside the convex hull of the finite element nodes.

>>> ax = []
>>> ax.append(plt.subplot(221))
>>> plt.plot(x_grid, y_grid, 'k.', ms=1)
>>> plt.title('DIC grid')
>>> ax.append(plt.subplot(222))
>>> plt.plot(mesh.vertices[:, 0], mesh.vertices[:, 1], 'k.', ms=1)
>>> plt.title('FE nodes')
>>> ax.append(plt.subplot(223))
>>> mesh.plot_field(0, ax=ax[-1])
>>> plt.title('FE field')
>>> ax.append(plt.subplot(224))
>>> plt.imshow(interpolated, extent=(-1, 1, 1, 2), origin='lower', vmin=-4, vmax = 5)
>>> plt.title('FE field interpolated at the DIC grid')
>>> plt.tight_layout()
>>> for a in ax:
...     a.set_aspect('equal', 'box')
>>> plt.show()

(Source code, png, hires.png, pdf)

You are invited to modify this example to simulate a finer mesh by increasing n_nodes. Also try different interpolation techniques.

project_onto_mesh(nodes, method='linear')

Project the experimental displacement field onto a mesh.

Parameters

nodes (ndarray) – 2D numpy array with 2 columns, each row corresponding to a node of the mesh, and the two columns giving the Cartesian coordinates of the nodes.

Returns

• u_nodes, v_nodes (ndarray) – 1D numpy arrays (vectors), the components of the projected displacement field at the given nodes.

• method ({‘nearest’, ‘linear’, ‘cubic’}, optional) – Type of the interpolation. The default is linear.

See also

project_onto_grid()

Examples

Suppose we have measured field data, obtained by DIC. The experimental field values are obtained as an image. The pixel centers of this image form a grid, each grid point having an associated scalar value, the value of the measured field. In this example, we pretend that the measured field is given as a function. The DIC grid and the measured field values can be seen in the upper left and lower left figures, respectively.

>>> from grains.dic import DIC
>>> from grains.geometry import TriMesh
>>> x_grid, y_grid = np.mgrid[-1:1:100j, 1:2:50j]
>>> exact_solution = lambda x, y: 1 - (x + y**2) * np.sign(x)
>>> grid = DIC(exact_solution(x_grid, y_grid).T, exact_solution(x_grid,y_grid).T)

The finite element solution is obtained at the nodes (upper right figure) of the mesh. Here, we assumed that the nodes are sampled from a uniformly random distribution on \([-1,1] \times [1, 2]\).

>>> n_nodes = 100  # modify this to see how good the interpolated solution is
>>> mesh = TriMesh(*TriMesh.sample_mesh(1, n_nodes))

For the sake of this example, we also assume that the nodal values are “exact”, i.e. they are the function values at the nodes. In reality, of course, this will not be the case, but this allows us to investigate the effect moving from the FE mesh to the DIC grid.

>>> grid.set_transformation((-1, 2), 50)
>>> fe_values = exact_solution(mesh.vertices[:, 0], mesh.vertices[:, 1])
>>> interpolated = grid.project_onto_mesh(mesh.vertices, 'nearest')

The FE solution available at the nodes are interpolated at the DIC grid points, as shown in the lower right figure. Interpolation with continuous functions cannot resolve well the discontinuity present in the “exact solution”, i.e. in the measurement data. A discontinuous manufactured solution was intentionally prepared to illustrate this.

>>> ax = []
>>> ax.append(plt.subplot(221))
>>> grid.plot_physicalgrid(ax[0])
>>> plt.title('DIC grid')
>>> ax.append(plt.subplot(222))
>>> mesh.plot('k.', ax=ax[1], markersize=1)
>>> plt.title('FE nodes')
>>> ax.append(plt.subplot(223))
>>> plt.imshow(exact_solution(x_grid, y_grid).T, extent=(-1, 1, 1, 2), origin='lower',
...            vmin=-4, vmax = 5)
>>> plt.title('Exact solution on the DIC grid')
>>> ax.append(plt.subplot(224))
>>> interpolated[0][np.isnan(interpolated[0])] = 0
>>> mesh.associate_field(interpolated[0])
>>> mesh.plot_field(0, show_mesh=False, ax=ax[3])
>>> plt.title('FE solution interpolated at the DIC grid')
>>> plt.show()

(Source code, png, hires.png, pdf)

You are invited to modify this example to simulate a finer mesh by increasing n_nodes. Also try different interpolation techniques.

set_transformation(origin, pixels_per_physicalunit)

Sets the transformation rule between the pixel and the physical coordinate systems.

To determine the position and the size of the DIC image in the physical space, a transformation rule needs to be given, as described in the Notes section of the DIC class.

Parameters

• origin (tuple) – 2-tuple of float, the coordinates of the origin of the DIC grid (upper left corner) in the physical coordinate system.

• pixels_per_physicalunit (float)

Returns

None

strain(strain_measure)

Computes the strain tensor from the displacement vector.

The symmetric strain tensor \(\boldsymbol{\varepsilon}\) with components

\[\begin{split}\boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_{11} & \varepsilon_{12} \\ \varepsilon_{21} & \varepsilon_{22} \end{pmatrix}\end{split}\]

is computed for the given strain measure. The partial derivatives of the displacement field \((u,v)\), available on an \(m \times n\) grid, are computed with numerical differentiation with second order accuracy.

Parameters

strain_measure ({‘infinitesimal’, ‘Green-Lagrange’}) – One of the following strain measures.

• ‘infinitesimal’

  \[\begin{split}\varepsilon_{11} = \frac{\partial u}{\partial x} \\ \varepsilon_{12} = \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\ \varepsilon_{22} = \frac{\partial v}{\partial y} \\\end{split}\]

• ‘Green-Lagrange’

  \[\begin{split}\varepsilon_{11} = \frac{1}{2} \left( 2\frac{\partial u}{\partial x} + \left(\frac{\partial u}{\partial x}\right)^2 + \left(\frac{\partial v}{ \partial x}\right)^2 \right) \\ \varepsilon_{12} = \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} + \frac{\partial u}{\partial x} \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \frac{\partial v}{\partial y} \right) \\ \varepsilon_{22} = \frac{1}{2} \left( 2\frac{\partial v}{\partial y} + \left(\frac{\partial u}{\partial y}\right)^2 + \left(\frac{\partial v}{ \partial y}\right)^2 \right) \\\end{split}\]

Returns

strain_tensor (ndarray) – Components of the strain tensor as an \(m \times n \times 3\) array. The first two dimensions correspond to the grid points they are determined at, the third dimension gives the components of the tensor in the following order: \(\varepsilon_{11}, \varepsilon_{12}, \varepsilon_{22}\).

See also

equivalent_strain(), plot_strain(), numpy.gradient()

Notes

The derivatives of the displacement field are computed in the the physical coordinate system.

Examples

Consider a rectangular body with a linear displacement function in horizontal direction and no displacement vertically. The deformation of the body is such that only the horizontal strain is nonzero.

>>> u = np.array([[0, 1e-3, 2e-3, 3e-3], [0, 1e-3, 2e-3, 3e-3], [0, 1e-3, 2e-3, 3e-3]])
>>> v = np.zeros((3,4))
>>> d = DIC(u, v)
>>> E = d.strain('infinitesimal')
>>> E[:, :, 0]
array([[0.001, 0.001, 0.001, 0.001],
       [0.001, 0.001, 0.001, 0.001],
       [0.001, 0.001, 0.001, 0.001]])
>>> np.allclose(E[:, :, 2], 0)
True
>>> np.allclose(E[:, :, 1], 0)
True

The Green-Lagrange strain tensor is slightly different. Comparing with the infinitesimal strain tensor shows how good the assumption of small deformations is.

>>> E_gl = d.strain('Green-Lagrange')
>>> E_gl[:, :, 0] - E[:, :, 0]
array([[5.e-07, 5.e-07, 5.e-07, 5.e-07],
       [5.e-07, 5.e-07, 5.e-07, 5.e-07],
       [5.e-07, 5.e-07, 5.e-07, 5.e-07]])

The other two strain components remain zero.

>>> np.allclose(E_gl[:, :, 1:2], 0)
True

grains.dic.plot_strain(strain, minval=None, maxval=None, colorbar=True, label='')

Plots a scalar strain field.

Parameters

• strain (ndarray) – Scalar strain field sampled on an m-by-n grid, given as a 2D numpy array.

• minval, maxval (float, optional) – Set the color scaling for the image by fixing the values that map to the colormap color limits. If minval is not provided, the limit is determined from the minimum value of the data. Similarly for maxval.

• colorbar (bool, optional) – If True, a horizontal colorbar is placed below the strain map. The default is True.

• label (str, optional) – Label describing the strain field. LaTeX code is also accepted, e.g. r’$varepsilon_{yy}$’. If not given, no text is displayed.

Returns

None

See also

DIC.strain()

Deprecated since version 1.1.0: This will be removed in 1.3.0. This function will be a static method of the DIC class. See also the deprecation warning of DIC.plot_strain().

Simulation

Functions

data_Pierre	Yield stresses and average grain diameters from Pierre’s thesis.
hallpetch_constants	Determines the two Hall-Petch constants.
hallpetch	Computes the yield stress from the Hall-Petch relation.
hallpetch_plot	Plots the Hall-Petch formula for given grain sizes and yield stresses.
change_domain	Extends or crops the domain an image fills.
nature_of_deformation	Characterizes the intergranular/intragranular deformations.

grains.simulation.change_domain(image, left, right, bottom, top, padding_value=nan)

Extends or crops the domain an image fills.

The original image is extended, cropped or left unchanged in the left, right, bottom and top directions by padding the corresponding 2D or 3D array with a given value or removing existing elements. Non-integer image length is avoided by rounding up. This choice prevents ending up with an empty image during cropping or no added pixel during extension.

Parameters

• image (ndarray) – 2D or 3D numpy array representing the original image.

• left, right (float) – When positive, the domain is extended, when negative, the domain is cropped by the given value relative to the image width.

• bottom, top (float) – When positive, the domain is extended, when negative, the domain is cropped by the given value relative to the image height.

• padding_value (float, optional) – Value for the added pixels. The default is numpy.nan.

Returns

changed_image (ndarray) – 2D or 3D numpy array representing the modified domain.

Examples

Crop an image at the top, extended it at the bottom and on the left, and leave it unchanged on the right. Note the rounding for non-integer pixels.

>>> import numpy as np
>>> image = np.array([[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]])
>>> modified = change_domain(image, 0.5, 0, 1/3, -1/3, 0)
>>> image  # no in-place modification, the original image is not overwritten
array([[0.1, 0.2, 0.3],
       [0.4, 0.5, 0.6],
       [0.7, 0.8, 0.9]])
>>> modified
array([[0. , 0. , 0.4, 0.5, 0.6],
       [0. , 0. , 0.7, 0.8, 0.9],
       [0. , 0. , 0. , 0. , 0. ]])

grains.simulation.data_Pierre()

Yield stresses and average grain diameters from Pierre’s thesis.

Returns

• sigma_y (list of floats) – Yield stresses.

• d (list of floats) – Diameters of the grains.

grains.simulation.hallpetch(sigma_0, k, d)

Computes the yield stress from the Hall-Petch relation.

Parameters

• sigma_0 (float) – Starting stress for dislocation movement (material constant).

• k (float) – Strengthening coefficient (material constant).

• d (float or list of floats) – Diameter of the grain.

Returns

sigma_y (float or list of floats) – Yield stress.

grains.simulation.hallpetch_constants(sigma_y, d)

Determines the two Hall-Petch constants. Given available measurements for the grains sizes and the yield stresses, the two constants in the Hall-Petch formula are computed.

Parameters

• sigma_y (list of floats) – Yield stresses.

• d (list of floats) – Diameters of the grains.

Returns

• sigma_0 (float) – Starting stress for dislocation movement (material constant).

• k (float) – Strengthening coefficient (material constant).

Notes

If two data points are given in the inputs (corresponding to two measurements), the output parameters have unique values:

k = (sigma_y[0] - sigma_y[1]) / (1/sqrt(d[0]) - 1/sqrt(d[1])) sigma_0 = sigma_y[0] - k/sqrt(d[0])

If there are more than two measurements, the resulting linear system is overdetermined. In both cases, the outputs are determined using least squares fitting.

grains.simulation.hallpetch_plot(sigma_y, d, units=('MPa', 'mm'))

Plots the Hall-Petch formula for given grain sizes and yield stresses.

Parameters

• sigma_y (ndarray) – Yield stresses.

• d (ndarray) – Grain diameters.

• units (2-tuple of str, optional) – Units for the yield stress and the grain diameters. The default is (“MPa”,”mm”).

Returns

fig (matplotlib.figure.Figure) – The figure object is returned in case further manipulations are necessary.

grains.simulation.nature_of_deformation(microstructure, strain_field, interface_width=3, visualize=True)

Characterizes the intergranular/intragranular deformations.

To decide whether the strain localization is intergranular (happens along grain boundaries, also called interfaces) or intragranular in a given microstructure, the strain field is projected on the microstructure. Here, by strain field we mean a scalar field, often called equivalent strain that is derived from a strain tensor.

It is irrelevant for this function whether the strain field is obtained from a numerical simulation or from a (post-processed) full-field measurement. All what matters is that the strain field be available on a grid of the same size as the microstructure.

The strain field is assumed to be localized on an interface if its neighborhood, with band width defined by the user, contains higher strain values than what is outside the band (i.e. the grain interiors). A too large band width identifies small grains to have boundary only, without any interior. This means that even if the strain field in reality localizes inside such small grains, the localization is classified as intergranular. However, even for extreme deformations, one should not expect that the strain localizes on an interface with a single-point width. Moreover, using a too small band width is susceptible to the exact position of the interfaces, which are extracted from the grain microstructure. A judicial balance needs to be achieved in practice.

Parameters

• microstructure (ndarray) – Labelled image corresponding to the segmented grains in a grain microstructure.

• strain_field (ndarray) – Discrete scalar field of the same size as the microstructure.

• interface_width (int, optional) – Thickness of the band around the interfaces.

• visualize (bool, optional) – If True, three plots are created. Two of them show the deformation field within the bands and outside the bands. They are linked together, so when you pan or zoom on one, the other plot will follow. The third plot contains two histograms on top of each other, giving the frequency of the strain values within the bands and outside the bands.

Returns

• boundary_strain (ndarray) – Copy of strain_field, but values outside the band are set to NaN.

• bulk_strain (ndarray) – Copy of strain_field, but values within the band are set to NaN.

• bands (ndarray) – Boolean array of the same size as strain_field, with True values corresponding to the band.

See also

grains.dic.DIC.strain()

Computes a strain tensor from the displacement field.

grains.dic.DIC.equivalent_strain()

Extracts a scalar quantity from a strain tensor.

matplotlib.pyplot.hist()

Plots a histogram.

Notes

1. From the modelling viewpoint, it is important to know whether the strain localizes to the grain boundaries or it is dominant within the grains as well. In the former case, simplifications in the models save computational time in the simulations.

2. In dynamics, the evolution of the strain field is relevant. E.g. an initially intergranular deformation can turn into diffuse localization that occurs within the grains as well. In that case, a strain field must be obtained at each time step, and this function can be called for each such instance.

Examples

The following figure was created by this function with visualize set to True and band_width chosen to be 3.

Utilities

This module provides general utility functions used by the grains package. The specific helper functions reside in the proper module. For example, a function that works on a general list goes here, but a computational geometry algorithm goes to the geometry module. The functions in the utils module can be interesting for other projects too, partly because of their general scope, and partly because of the few dependencies.

Functions

duplicates	Set of duplicate elements in a sequence.
toggle	Return True for False values and False for True values in a list.
index_list	Index a list by another list.
flatten_list	Merge a list of lists to a single list.
argsorted	Return the indices that would sort a list or a tuple.
map_inplace	Apply a function to each member of an iterable in-place.
non_unique	Finds indices of non-unique elements in a 1D or 2D ndarray.
parse_kwargs	Compares keyword arguments with defaults.
compress	Creates a zip archive from a single file.
decompress	Decompresses the contents of a zip archive into the current directory.
decompress_inmemory	Decompresses the contents of a zip archive into a dictionary.
neighborhood	Neighboring points to a grid point.

grains.utils.argsorted(sequence, reverse=False)

Return the indices that would sort a list or a tuple.

Implementation is taken from https://stackoverflow.com/a/6979121/4892892.

Parameters

• sequence (list, tuple) – Input sequence in which the sorted indices to be found.

• reverse (bool) – If set to True, then the elements are sorted as if each comparison was reversed.

Returns

list – List of indices that would sort the input list/tuple.

See also

sorted(), numpy.argsort()

Examples

>>> argsorted([2, 1.1, 1.1])
[1, 2, 0]
>>> argsorted([2, 1.1, 1.1], reverse=True)
[0, 1, 2]
>>> argsorted(())
[]

grains.utils.compress(filename, level=9)

Creates a zip archive from a single file.

Parameters

• filename (str) – Name of the file to be compressed.

• level (int, optional) – Level of compression. Integers 0 through 9 are accepted. The default is 9.

Returns

None.

See also

zipfile.ZipFile

grains.utils.decompress(filename, path=None)

Decompresses the contents of a zip archive into the current directory.

Parameters

• filename (str) – Name of the zip archive.

• path (str, optional) – Directory to extract to. The default is the directory the function is called from.

See also

zipfile.ZipFile

grains.utils.decompress_inmemory(filename)

Decompresses the contents of a zip archive into a dictionary.

Parameters

filename (str) – Name of the zip archive.

Returns

data (dict) – The keys of the dictionary are the compressed file names (without extension), the corresponding values are their contents.

See also

zipfile.ZipFile

grains.utils.duplicates(sequence)

Set of duplicate elements in a sequence.

Parameters

sequence (sequence types (list, tuple, string, etc.)) – Sequence possibly containing repeating elements.

Returns

set – Set of unique values.

Notes

Copied from https://stackoverflow.com/a/9836685/4892892

Examples

Note that the order of the elements in the resulting set does not matter.

>>> a = [1, 2, 3, 2, 1, 5, 6, 5, 5, 5]  # list
>>> duplicates(a)
{1, 2, 5}
>>> a = (1, 1, 0, -1, -1, 0)  # tuple
>>> duplicates(a)
{0, 1, -1}
>>> a = 'abbcdkc'  # string
>>> duplicates(a)
{'c', 'b'}

grains.utils.flatten_list(nested_list)

Merge a list of lists to a single list.

Parameters

nested_list (list) – List containing other lists.

Returns

list – Flattened list.

Notes

• Only a single level (i.e. list of lists) is handled, see the second example.

• Several methods, such as list comprehension, monoid and loops, are proposed in https://stackoverflow.com/questions/952914/how-to-make-a-flat-list-out-of-list-of-lists. Here, the list comprehension approach is used.

Examples

>>> nested_list = [['some'], ['items']]
>>> flatten_list(nested_list)
['some', 'items']
>>> multiply_nested_list = [[['item'], 'within', 'item']]
>>> flatten_list(multiply_nested_list)
[['item'], 'within', 'item']

grains.utils.index_list(lst, indices)

Index a list by another list.

Parameters

• lst (list) – List to be indexed.

• indices (list) – Indices of the original list that will form the new list.

Returns

list – Members of lst, selected by indices.

Examples

>>> index_list(['c', ['nested', 'list'], 13], [1, 2])
[['nested', 'list'], 13]

grains.utils.map_inplace(function, __iterable)

Apply a function to each member of an iterable in-place.

Parameters

• function (function object) – Function to be applied to the entries of the iterable.

• __iterable (iterable) – Iterable.

Notes

Comprehensions or functional tools work on iterators, thereby not modifying the original container (https://stackoverflow.com/a/4148523/4892892). For in-place modification, the conventional for loop approach is used (https://stackoverflow.com/a/4148525/4892892).

Examples

>>> lst = ['2', 2]; func = lambda x: x*2
>>> map_inplace(func, lst); lst
['22', 4]
>>> lifespan = {'cat': 15, 'dog': 12}; die_early = lambda x: x/2
>>> map_inplace(die_early, lifespan); lifespan
{'cat': 7.5, 'dog': 6.0}

grains.utils.neighborhood(center, radius, norm, method='ball', bounds=None)

Neighboring points to a grid point.

Given a point in a subspace of \(\mathbb{Z}^n\), the neighboring points are determined based on the specified rules.

Todo

Currently, the neighbors are deterministically but not systematically ordered. Apart from testing purposes, this does not seem to be a big issue. Still, a logical ordering is desirable. E.g. ordering in increasing coordinate values, first in the first dimension and lastly in the last dimension.

Parameters

• center (tuple of int) – Point \(x_0\) around which the neighborhood is searched. It must be an n-tuple of integers, where \(n\) is the spatial dimension.

• radius (int, positive) – Radius of the ball or sphere in which the neighbors are searched. If the radius is 1, the immediate neighbors are returned.

• norm ({1, inf}) – Type of the vector norm \(\| x - x_0 \|\), where \(x\) is a point whose distance is searched from the center \(x_0\).

  Type	Name	Definition
  1	1-norm	\(\| x \|_1 = \sum_{i=1}^n |x_i|\)
  numpy.inf	maximum norm	\(\| x \|_{\infty} = \max\{ |x_1|, \ldots, |x_n| \}\)

  where inf means numpy’s np.inf object.

• method ({‘ball’, ‘sphere’}, optional) – Specifies the criterion of how to decide whether a point \(x\) is in the neighborhood of \(x_0\). The default is ‘ball’.

  For ‘ball’:

  \[\| x - x_0 \| \leq r\]

  For ‘sphere’:

  \[\| x - x_0 \| = r\]

  where \(r\) is the radius passed as the radius parameter and the type of the norm is taken based on the norm input parameter.

• bounds (list of tuple, optional) – Restricts the neighbors within a box. The dimensions of the n-dimensional box are given as a list of 2-tuples: [(x_1_min, x_1_max), … , (x_n_min, x_n_max)]. The default value is an unbounded box in all dimensions. Use np.inf to indicate unbounded values.

Returns

neighbors (tuple of ndarray) – Tuple of length n, each entry being a 1D numpy array: the integer indices of the points that are in the neighborhood of center.

Notes

1. The von Neumann neighborhood with range \(r\) is a special case when radius=r, norm=1 and method='ball'.

2. The Moore neighborhood with range \(r\) is a special case when radius=r, norm=np.inf and method='ball'.

Examples

Find the Moore neighborhood with range 2 the point (1) on the half-line [0, inf).

>>> neighborhood((1,), 2, np.inf, bounds=[(0, np.inf)])
(array([0, 1, 2, 3]),)

Find the von Neumann neighborhood with range 2 around the point (2, 2), restricted on the domain [0, 4] x [0, 3].

>>> neighborhood((2, 2), 2, 1, bounds=[(0, 4), (0, 3)])
(array([2, 1, 2, 3, 0, 1, 2, 3, 4, 1, 2, 3]), array([0, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3]))

Find the Moore neighborhood with range 1 around the point (0, -4) such that the neighbors lie on the half-plane [0, 2] x (-inf, -4].

>>> neighborhood((0, -4), 1, np.inf, bounds=[(0, 2), (-np.inf, -4)])
(array([0, 1, 0, 1]), array([-5, -5, -4, -4]))

Find the sphere of radius 2, measured in the 1-norm, around the point (-1, 0, 3), within the half-space {(x,y,z) in Z^3 | y>=0}.

>>> neighborhood((-1, 0, 3), 2, 1, 'sphere', [(-np.inf, np.inf), (0, np.inf), (-np.inf, np.inf)])  
(array([-3, -2, -2, -1, -1,  0,  0,  1, -2, -1, -1,  0, -1]),
 array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2]),
 array([3, 2, 4, 1, 5, 2, 4, 3, 3, 2, 4, 3, 3]))

grains.utils.non_unique(array, axis=None)

Finds indices of non-unique elements in a 1D or 2D ndarray.

Parameters

• array (ndarray) – Array in which the non-unique elements are searched.

• axis ({None, 0, 1}, optional) – The axis to operate on. If None, array will be flattened. If an integer, the subarrays indexed by the given axis will be flattened and treated as the elements of a 1-D array with the dimension of the given axis. Object arrays or structured arrays that contain objects are not supported if the axis kwarg is used. The default is None.

Returns

• nonunique_values (list) – Unique (individual, row or column) entries.

• nonunique_indices (list) – Each element of the list corresponds to non-unique elements, whose indices are given in a 1D numpy array.

Examples

In a 1D array, the repeated values and their indices are found by

>>> val, idx = non_unique(np.array([1, -1, 0, -1, 2, 5, 0, -1]))
>>> val
[-1, 0]
>>> idx
[array([1, 3, 7]), array([2, 6])]

In the matrix below, we can see that rows 0 and 2 are identical, as well as rows 1 and 4.

>>> val, idx = non_unique(np.array([[1, 3], [2, 4], [1, 3], [-1, 0], [2, 4]]), axis=0)
>>> val
[array([1, 3]), array([2, 4])]
>>> idx
[array([0, 2]), array([1, 4])]

By transposing the matrix above, the same holds for the columns.

>>> val, idx = non_unique(np.array([[1, 2, 1, -1, 2], [3, 4, 3, 0, 4]]), axis=1)
>>> val
[array([1, 3]), array([2, 4])]
>>> idx
[array([0, 2]), array([1, 4])]

If the dimensions along which to find the duplicates are not given, the input is flattened and the indexing happens in C-order (row-wise).

>>> val, idx = non_unique(np.array([[1, 2, 1, -1, 2], [3, 4, 3, 0, 4]]))
>>> val
[1, 2, 3, 4]
>>> idx
[array([0, 2]), array([1, 4]), array([5, 7]), array([6, 9])]

grains.utils.parse_kwargs(kwargs, defaults)

Compares keyword arguments with defaults.

Allows processing keyword arguments supplied to a function by the user by comparing them with an admissible set of options defined by the developer. There are three cases:

1. The keyword given by the user is present in the set the developer provides. Then the value belonging to the keyword overrides the default value.

2. The keyword given by the user is not present in the set the developer provides. In this case, the unrecognized keyword, along with its value, is saved separately.

3. The keyword existing in the set the developer provides is not given by the user. Then the default value is used.

Parameters

• kwargs (dict) – Keyword arguments (parameter-value pairs) passed to a function.

• defaults (dict) – Default parameter-value pairs.

Returns

• parsed (dict) – Dictionary with the same keys as defaults, the parsed parameter-value pairs.

• unknown (dict) – Dictionary containing the parameter-value pairs not present in defaults.

Notes

The default values, given in the input dictionary defaults, are never overwritten.

Examples

>>> default_options = {'opt1': 1, 'opt2': 'string', 'opt3': [-1, 0]}
>>> user_options = {'opt3': [2, 3, -1], 'opt2': 'string',  'opt4': -2.14}
>>> parsed_options, unknown_options = parse_kwargs(user_options, default_options)
>>> parsed_options
{'opt1': 1, 'opt2': 'string', 'opt3': [2, 3, -1]}
>>> unknown_options
{'opt4': -2.14}

grains.utils.toggle(lst)

Return True for False values and False for True values in a list.

Parameters

lst (list) – An arbitrary list, possibly containing other lists.

Returns

list – Element-wise logical not operator applied on the input list.

Notes

Solution taken from https://stackoverflow.com/a/51122372/4892892.

Examples

>>> toggle([True, False])
[False, True]
>>> toggle(['h', 0, 2.3, -2, 5, []])
[False, True, False, False, False, True]

Profiling

This module implements profiling facilities so that user code can be tested for speed.

Functions

profile([output_format, output_filename, …])

Profiles a block of code with Pyinstrument.

grains.profiling.profile(output_format='html', output_filename=None, html_open=False)

Profiles a block of code with Pyinstrument.

Intended to be used in a with statement.

Parameters

• output_format ({‘html’, ‘text’}, optional) – Shows the result either as text or in an HTML file. If output_format = 'html', the file is saved according to the output_filename parameter. The default is ‘html’.

• output_filename (str, optional) – Only taken into account if output_format = 'html'. If not given (default), the html output is saved to the same directory the caller resides. The name of the html file is the same as that of the caller.

• html_open (bool, optional) – Only taken into account if output_format = 'html'. If True, the generated HTML file is opened in the default browser. The default is False.

Notes

In the implementation, we move two levels up in the stack frame, one for exiting the context manager and one for exiting this generator. This assumes that profile() was called as a context manager. As a good practice, provide the output_filename input argument.

Examples

Measure the time needed to generate 1 million uniformly distributed random numbers.

>>> import random
>>> with profile('html') as p:
...    for _ in range(1000000):
...        rand_num = random.uniform(1, 2.2)

Index

Index

Changelog

All notable changes to this project will be documented in this file. The format is based on Keep a Changelog, and this project does not adhere to Semantic Versioning. We refer to GitHub issues by their numbers. If there is no issue associated to the change, the commit hashes implementing the change are linked.

Unreleased

1.1.0 - 2021-03-04

Added

• Grains embedded into other grains are now identified in the splinegon creation algorithm. 4cef7ec, 6c4ad0a

• Computation of the infinitesimal and Green-Lagrange strain tensors for DIC displacement field. 7eaa0b5, 9e5466e

• Computation of the equivalent von Mises strain. 2b9b0be

• Characterization of the strain localization (intergranular/intragranular) in the microstructure. a74f954, 13576f7

• The README file shows how to cite our paper. #30

• Added a new segmented microstructure (6b539ba) and the corresponding DIC measurements (3c845e7).

Deprecated

• The content of dic.DIC.plot_strain will be replaced with that of dic.plot_strain in version 1.3.0. The latter function will be removed. b8633c2

1.0.1 - 2020-11-19

Added

• Document on describing the versioning scheme we will follow. e74474c

• Added a changelog to the project. b48821b, f6ca9c1

Removed

• The documentation no longer shows the recent git commits. 1fd3a50

1.0.0 - 2020-11-16

This is the initial release of CristalX.

License

GNU LESSER GENERAL PUBLIC LICENSE

Version 3, 29 June 2007

Copyright (C) 2007 Free Software Foundation, Inc. <https://fsf.org/> Everyone is permitted to copy and distribute verbatim copies of this license document, but changing it is not allowed.

This version of the GNU Lesser General Public License incorporates

the terms and conditions of version 3 of the GNU General Public License, supplemented by the additional permissions listed below.

0. Additional Definitions.

As used herein, “this License” refers to version 3 of the GNU Lesser

General Public License, and the “GNU GPL” refers to version 3 of the GNU General Public License.

“The Library” refers to a covered work governed by this License,

other than an Application or a Combined Work as defined below.

An “Application” is any work that makes use of an interface provided

by the Library, but which is not otherwise based on the Library. Defining a subclass of a class defined by the Library is deemed a mode of using an interface provided by the Library.

A “Combined Work” is a work produced by combining or linking an

Application with the Library. The particular version of the Library with which the Combined Work was made is also called the “Linked Version”.

The “Minimal Corresponding Source” for a Combined Work means the

Corresponding Source for the Combined Work, excluding any source code for portions of the Combined Work that, considered in isolation, are based on the Application, and not on the Linked Version.

The “Corresponding Application Code” for a Combined Work means the

object code and/or source code for the Application, including any data and utility programs needed for reproducing the Combined Work from the Application, but excluding the System Libraries of the Combined Work.

1. Exception to Section 3 of the GNU GPL.

You may convey a covered work under sections 3 and 4 of this License

without being bound by section 3 of the GNU GPL.

2. Conveying Modified Versions.

If you modify a copy of the Library, and, in your modifications, a

facility refers to a function or data to be supplied by an Application that uses the facility (other than as an argument passed when the facility is invoked), then you may convey a copy of the modified version:

a) under this License, provided that you make a good faith effort to ensure that, in the event an Application does not supply the function or data, the facility still operates, and performs whatever part of its purpose remains meaningful, or

b) under the GNU GPL, with none of the additional permissions of this License applicable to that copy.

3. Object Code Incorporating Material from Library Header Files.

The object code form of an Application may incorporate material from

a header file that is part of the Library. You may convey such object code under terms of your choice, provided that, if the incorporated material is not limited to numerical parameters, data structure layouts and accessors, or small macros, inline functions and templates (ten or fewer lines in length), you do both of the following:

a) Give prominent notice with each copy of the object code that the Library is used in it and that the Library and its use are covered by this License.

b) Accompany the object code with a copy of the GNU GPL and this license document.

4. Combined Works.

You may convey a Combined Work under terms of your choice that,

taken together, effectively do not restrict modification of the portions of the Library contained in the Combined Work and reverse engineering for debugging such modifications, if you also do each of the following:

a) Give prominent notice with each copy of the Combined Work that the Library is used in it and that the Library and its use are covered by this License.

b) Accompany the Combined Work with a copy of the GNU GPL and this license document.

c) For a Combined Work that displays copyright notices during execution, include the copyright notice for the Library among these notices, as well as a reference directing the user to the copies of the GNU GPL and this license document.

4. Do one of the following:

   0) Convey the Minimal Corresponding Source under the terms of this License, and the Corresponding Application Code in a form suitable for, and under terms that permit, the user to recombine or relink the Application with a modified version of the Linked Version to produce a modified Combined Work, in the manner specified by section 6 of the GNU GPL for conveying Corresponding Source.

   1) Use a suitable shared library mechanism for linking with the Library. A suitable mechanism is one that (a) uses at run time a copy of the Library already present on the user’s computer system, and (b) will operate properly with a modified version of the Library that is interface-compatible with the Linked Version.

e) Provide Installation Information, but only if you would otherwise be required to provide such information under section 6 of the GNU GPL, and only to the extent that such information is necessary to install and execute a modified version of the Combined Work produced by recombining or relinking the Application with a modified version of the Linked Version. (If you use option 4d0, the Installation Information must accompany the Minimal Corresponding Source and Corresponding Application Code. If you use option 4d1, you must provide the Installation Information in the manner specified by section 6 of the GNU GPL for conveying Corresponding Source.)

5. Combined Libraries.

You may place library facilities that are a work based on the

Library side by side in a single library together with other library facilities that are not Applications and are not covered by this License, and convey such a combined library under terms of your choice, if you do both of the following:

a) Accompany the combined library with a copy of the same work based on the Library, uncombined with any other library facilities, conveyed under the terms of this License.

b) Give prominent notice with the combined library that part of it is a work based on the Library, and explaining where to find the accompanying uncombined form of the same work.

6. Revised Versions of the GNU Lesser General Public License.

The Free Software Foundation may publish revised and/or new versions

of the GNU Lesser General Public License from time to time. Such new versions will be similar in spirit to the present version, but may differ in detail to address new problems or concerns.

Each version is given a distinguishing version number. If the

Library as you received it specifies that a certain numbered version of the GNU Lesser General Public License “or any later version” applies to it, you have the option of following the terms and conditions either of that published version or of any later version published by the Free Software Foundation. If the Library as you received it does not specify a version number of the GNU Lesser General Public License, you may choose any version of the GNU Lesser General Public License ever published by the Free Software Foundation.

If the Library as you received it specifies that a proxy can decide

whether future versions of the GNU Lesser General Public License shall apply, that proxy’s public statement of acceptance of any version is permanent authorization for you to choose that version for the Library.