Motivation

By allowing the degradation along the interfaces, one can better capture the localization, which is often intergranular in polycrystalline materials. The cohesive zone model is simple to implement for non-moving interfaces, once we have access to the 1D elements on the interfaces.

Useful links

Cohesive zone modelling in Abaqus

Here come the links I bookmarked...

Workflow

The following tasks need to be done:

• Fork Phon and make it work with Python 3.
  • Issue with non-mutable OrderedDict: KristofferC/Phon#9 (comment)
  • Use list(dictionary) instead of dictionary.keys() to loop through and modify dictionary.
  • Support 2D nodes (maybe it's only needed during the .inp file writing)
  • Support CPS3 element types for Abaqus.
• Prepare the med module to be able to extract edges and edge groups from a .med file.
• Create a function/class in the simulation module (not in the geometry module to keep it uncoupled) that will manage the whole cohesive zone insertion.
  • Create a wrapper class for Phon in simulation (if it proves to be a general solution, I can still refactor it later). Going through and understand everything in Phon would take too much time to worth it for now. Rather create a function/class that wraps it. That wrapper should contain a converter to create Mesh object, then to call the create_cohesive_elements function, and finally convert it back to my format. I should not use Phon's capability to write to an input file because it deletes important information, such as material definitions. Some links to start with:
    • https://phon.readthedocs.io/en/latest/examples/insert_interface.html
    • https://phon.readthedocs.io/en/latest/phon.mesh_objects.html#module-phon.mesh_objects.mesh
• Document the workflow by
  • creating an example in a Jupyter notebook
  • write about the extraction of the 1D interface mesh in https://cristalx.readthedocs.io/en/latest/salome.html

Only the very high level modules (e.g. simulation) should import other modules. For other improvements, see Better Code Hub.

Currently, the DIC class knows about the displacement field at a given time step only. On the other hand, one is often interested in the evolution of the displacement field. It is clear that loading the displacement measurement data for each time step is not viable in general on personal computers¹. Moreover, the member functions of the class would become more verbose because the user would need to give at which time step a given operation should be applied. Since the class cannot know in advance how the series of displacement fields are stored, it is the responsibility of the user to pass it to DIC. The DIC class should have a method that loads the new field.

¹E.g. in my current dataset, the DIC data is given as an HDF5 file with the size of 38 GB.

Create an Abaqus class with the following desired capabilities:

• do not rely on the Salome API so that it can be factored out later to a standalone module (see also #28)
• additional functionalities required for running an Abaqus simulation (materials, section definitions, etc).
  This might seem unimportant at first sight (why would one do it, when they can be set from within Abaqus?), however, setting a different material property for each grain in a microstructure is tedious.

Also, create a script that shows the use of these.

Those in which it makes sense to store the performed operations and the possible warnings for later retrieval.

• FixedDict in meshing.py
• (maybe) plot_prop in analysis.py
• save_image in segmentation.py
• common methods in the Material and Geometry classes in abaqus.py
• simulation.py is not intended to be run, factor out the runnable components to a script

• Sourcery

  [![Sourcery](https://img.shields.io/badge/Sourcery-enabled-brightgreen)](https://sourcery.ai)`
  

  

While asking for something else, I got an answer to visualize a mesh with its partitions exploded. It is easy to implement and helps in visually distinguishing small partitions, which would remain indistinguishable due to the neighboring larger partitions.

Similarly as done for matplotlib: https://github.com/CsatiZoltan/Polycrystalline-microstructures/blob/89cdfed050a268a5904ab126dc42efbf258e1fea/grains/analysis.py#L154

The center of the grain can be found by computing the center of mass, assuming that the grain is homogeneous (refer to core_shape_properties.py):

from OCC.GProp import GProp_GProps
from OCC.BRepGProp import brepgprop_SurfaceProperties
...
# Compute inertia properties
props = GProp_GProps()
brepgprop_SurfaceProperties(region, props)  # region is a *TopoDS_Face* object
area = props.Mass()  # assuming that the region is homogeneous
cog = props.CentreOfMass()
cog_x, cog_y, cog_z = cog.Coord()

More and more functions work on labeled images in the codebase. Currently, these are

• show_label_image
• label_image_skeleton
• thicken_skeleton
• label_image_apply_mask
• some methods of the Analysis class could also come here

in the analysis module and

• Segmentation.save_image
• Segmentation.save_array

in the segmentation module.

It is time to collect them in a class. We will call this class LabeledImage¹ and it should have, at first thought, the following methods:

• __init__(self, label_image)
• find_skeleton
• thicken_skeleton
• apply_mask(self, mask, value)
• change_label(self, old_labels, new_labels)
• show(self, color='random', labels=False) ²
• save_png(self, filename)
• save_numpy(self, filename)
• _validate(label_image)
• __str__(self)

and with the following data members:

• n_label

Once the class has been created, deprecate the old functions and also mention this class creation principal as an example under the "Do not overuse classes" bullet point of the program design.

¹ MATLAB calls it "label image". However, the word "label" is also a verb and we want to avoid the misunderstanding that we perform labeling as done in computer vision. We also prefer "labeled" to "labelled" as American English is used in the rest of CristalX.
² Other values for the color parameter:

• 'seeded_random', see #45
• 'optimal', see #9

E.g. change_node_numbering() in the See also section of grains.med.get_elements.

For the list of changes, see here.

When plotting the same grain configuration twice, the displayed colors are different because colors are randomly allocated to each grain. This is not only annoying when making comparisons, but it can also deceive people who do not know about the random color allocation.
The solution is to use a given seed for the random number generator. Then the same colors are used for subsequent plotting.

from numpy.random import default_rng
rng = default_rng(seed)  # uses the PCG64 generator
vals = rng.random((m, n))

See also #9.

Implement it in the analysis module as show_label_image.

PyCharm complains about it.

Although this project works with Binder, the notebook is launched quite slowly. Try to speed it up: https://discourse.jupyter.org/t/how-to-reduce-mybinder-org-repository-startup-time/4956

Thorough checking for input types clutters the code and is not Pythonic. The user is expected to read the docstring of the parameters and pass the proper type.
Later, the type checking could be enforced by using static analyzers, such as mypy.

https://github.com/KristofferC/Phon/blob/1b09eb2af26151e6d5cb36f74c4de8c83a103bf4/docs/conf.py#L113-L115

See https://www.youtube.com/watch?v=XI58hlGA7Is. Although I came up with most of them by myself, it's worth to follow.

Depending on the input image, different strategies may work better. Implement the strategies seen in this talk (source code here) and extend them with the one currently being the most promising in my case.

Useful utilities for positioning the text can be found in core_geometry_utils.py of the PythonOCC-examples project for version 0.18 and in ShapeFactory.py in the new versions.
The text itself can be displayed using the DisplayMessage method, see core_topology_glue.py.

E.g. on Read the Docs.

When referring to other functions in the docstrings, use (see here)

:py:func:`function_name`

so that the generated Sphinx documentation turn them into links.

Some issues in the sibling repository are not connected to the GUI but to the core. Resolve those issues.

• CsatiZoltan/GrainSegmentation#2
• CsatiZoltan/GrainSegmentation#3
• CsatiZoltan/GrainSegmentation#6
• CsatiZoltan/GrainSegmentation#9
• CsatiZoltan/GrainSegmentation#11

In cad.plot_polygon, indicate the orientation by placing arrows on the sides: https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.pyplot.arrow.html.

von Neumann, Moore, and their hollow versions. See also #27

A recurring task in this project is to interpolate scattered data on a grid. The function griddata returns nan for the locations being outside the convex hull of the input data, for all but the nearest neighbor interpolant. A brute-force method to obtain interpolated values at those external points too is to perform two interpolations: one with the requested interpolant (linear, cubic, etc.) and one with the nearest neighbor. Then use the nearest neighbour interpolated values for the data points where the other interpolants return nan. A possibly faster way could be to perform the nearest neighbor interpolation only at the external points. Two approaches in mind:

• https://stackoverflow.com/questions/35807321/scipy-interpolation-with-masked-data
• masked arrays in Numpy.

Follow the numpydoc style to document modules, classes, etc. The newly written codes usually adhere to this guideline, but the old ones do not.

Currently, random colors are generated when saving a label image.
https://github.com/CsatiZoltan/Polycrystalline-microstructures/blob/20e3eb63df753251911fbe1c2b8ae410cb0c8c7d/grains/segmentation.py#L265
An ongoing issue deals with it: scikit-image/scikit-image#4507
Ideally, colors should be chosen so that nearby labelled regions are easily distinguishable. See the above issue for details.
This functionality is already implemented in ImageJ: https://imagej.net/Glasbey

If we publish this work, show a citing request on the main page of the documentation, such as this.
See also here.

Version numbers and tags are used in more and more places. It is easy to forget about increasing them all. At least create a script that makes the necessary changes, but the best would be a pre-commit hook.

Places where the version numbers change:

• conf.py
• grains.__init__.py
• CHANGELOG.md
• git tags
• GitHub releases

This project was initialized with the MIT license. As I plan to use GPLv3 licensed code, the whole project will have to use either GPLv3 or a compatible license. More on the GPLv3 license: https://choosealicense.com/licenses/gpl-3.0/.

I am not sure about the implications of re-licensing this project. Will it make others more difficult to use my code?

https://github.com/CsatiZoltan/GrainSegmentation

Processing an Abaqus input file (.inp) is more complicated than it seemed. The current implementation used nested dictionaries to extract data from a .inp file. For deeply nested dictionaries, it hampers the readability and it is difficult to insert a new level somewhere in the middle. Another problem is that many Abaqus commands are similar but still need slightly different inputs. The many edge cases render the current approach inflexible for extension.

The .inp files have a set of syntax rules and keywords. This makes it a domain specific language. Hence, an AST can be used to analyze it. Recognizing this gave the idea to represent the Abaqus module hierarchy as a tree. E.g. if the .inp file contains

* Material, name=mat-1
** possible comment line
* Elastic
** possible comment line
data line
* Plastic
data line 1
** possible comment
data line 2
* Material, name=mat-2
...

it could be represented as

Root
|--- ...
|--- Material (object containing the name: 'mat-1')
      |--- Elastic (object containing the elastic properties)
      |--- Plastic (object containing the plastic properties)
|--- Material (object containing the name: 'mat-2')
...

Being a tree data structure, including new commands or changing the hierarchy is easy. E.g. if we want the Material command to be part of the Property Abaqus module, we insert a Property object under the Root and set the Material objects to be children of the Property (singleton) object.

Once we have the tree, we can fill it in with data coming from the input file. We define parent-children relationships so that the added nodes (commands) are inserted to the correct position of the tree. For each command, create a class

• which serves as a container holding the data read from the .inp file
• checks for admissibility
• the object of which is a node of the tree

Manipulations (modifying element connectivities, adding materials, etc.) are easy to carry out in this intermediate representation, no need to deal with the text representation. Writing the (possibly modified) state into file is straightforward by traversing the tree.

anytree, written in pure Python, is a promising package, in which custom objects can form the nodes of the tree.

A MED file (.med) is an HDF5 file in disguise. Processing a MED file to extract mesh information requires the MEDCoupling module. I couldn't compile a standalone MEDCoupling interface. What I currently do is

1. Download Salome.
2. Set up my IDE to use Salome's built-in Python interpreter.
3. Use that interpreter when working with MED mesh processing.

Comments for the workflow above:

1. This is not a problem because I anyway use Salome to repair the geometry and to export it to .med.
2. It has its quirks, but I created a short script for it.
3. This point is inconvenient. The reason is that when I use the Python interpreter of Salome
   • I cannot install packages that need compilation because Python links with the wrong dynamically linked libraries
   • even if I could install packages to Salome's Python, those packages can crash the original environment (happened to me)
   • I am restricted to the preinstalled versions of numpy and matplotlib

So the goal is to use load and process the .med file using h5py, which has very few dependencies. As a side effect, the user of my code does not need the several GB sized Salome so that they can handle the .med file. The structure of a .med file can easily be visualized with ViTables.

Commit b45154e includes some temporary measures, used for debugging.

When a large grain contains one or more grains, the large grain has to be modified to exclude the small grains that act as holes. This makes the large grain a non-simply connected domain.
This modification can be carried out as

• Boolean operation (set difference)
• boundary representation

Currently, I plan to support the Boolean solution for the splinegon representation only, as that's what I actively work with. However, a general implementation would be nice, which does not need to have information about the actual geometrical representation. I.e. it would work for the polygon representation too. That would however require several changes in the code:

• The cad module would import the geometry module for the sake polygon-based operations.
• The region_branches variable will no longer be a simple list, but a list of lists in order to incorporate the holes
• The geometry.Polygon class must be modified so that holes can be represented, which needs
  • updating the methods for calculating the area, the equivalent diameter, etc.
  • the data structure will not remain a simple nx2 NumPy array, but a list of NumPy arrays, that describe the outer boundary and the internal boundaries as well

The strategy implemented in 6c4ad0a is very simple and general, but it is very slow. The profiling results are below. They were executed on another microstructure, also contained in the zipped file.
measurements.zip

Create a function in the geometry module that displays the geometry on top of the label image it was constructed from. This is useful for presenting the main purpose of the module and for visual debugging too. May get ideas from skan.draw.overlay_skeleton_2d_class.

Different images often require different meshing techniques. Once experimented in the OOF2 GUI, one would like to save the corresponding OOF2 commands. Create a method add_command in class OOF2.

See https://docs.pyvista.org/index.html and its source code.

Physics

From the modelling viewpoint, it is of importance to know whether the strain localizes to the grain boundaries (also called interfaces) or it is dominant within the grains as well. Indeed, the cohesive zone/band model can largely speed up the computations in the first case. To decide whether the localization is intergranular or intragranular in a given microstructure, we can make use of experimental data (if available).
From now on, we will assume that thanks to full-field measurement (e.g. digital image correlation - DIC), we have access to the strain field at every point¹ in the microstructure.

Algorithm

Things to do:

1. Define a band on the interfaces. The thickness of the band is a parameter of the model. ce6c499

2. Be able to compute

   • a strain measure (tensor) from the displacement field 7eaa0b5
   • an equivalent strain from the strain tensor 2b9b0be

   on the whole domain (where the displacement field exists)

3. Take the grain microstructure and the (equivalent) strain field and

   • localize the strain field on the band a74f954
   • possibly return a histogram of the magnitude of the equivalent strain in the band a74f954
   • determine what portion of the large strain values lie in the band, compared to the rest of the domain (i.e. the internal regions of the grains) a74f954

4. The ratio of inter/intragranular deformation should be

   • executed for different time steps, as the deformation evolves 13576f7; and
   • the data above should be visualized somehow 13576f7

5. Once the methodology works properly, create a Jupyter notebook from it.

Expected results

For many microstructures, the strain localization is dominantly intergranular. If that is the case,

1. Device a model to allow simplifications (e.g. no need to account plastic deformation in the grain interiors)
2. Allow extracting the generated mesh on the interfaces (they will be a set of 1D elements)

Footnotes

¹ I.e. in every pixel. Use subpixel interpolation if you need a value at an arbitrary position.

Commit dc50d41 introduced a deprecation life-cycle to CristalX.
Managing deprecated functionalities, and the fact in which release it was marked as deprecated and in which future release it is planned to be removed, becomes cumbersome as CristalX grows.
Therefore, create a function in the grain package (in __init__.py) to fetch these data. I found 3 main possibilities:

1. Analyze the source code with the inspect built-in Python module
   "Method 2" in https://stackoverflow.com/a/5910893/4892892
2. Create an additional decorator to decorate the decorator you want to track
   "Method 3" in https://stackoverflow.com/a/5910893/4892892. In my case, I would need to decorate the deprecated decorator from the deprecation package.
3. Analyze the AST
   Provided by https://stackoverflow.com/a/9580006/4892892, this is also a source parsing method, as the first item.

Links to the documentation do not work. Get the new links.

Create image formats that the PDF backend can display.
See https://stackoverflow.com/a/45970146/4892892