# Mesh-based graph convolutional neural networks for modeling materials with
microstructure

Ari Frankel, Cosmin Safta, Coleman Alleman, Reese Jones111Corresponding
author<EMAIL_ADDRESS>
Sandia National Laboratories, Livermore, CA 94551

###### Abstract

Predicting the evolution of a representative sample of a material with
microstructure is a fundamental problem in homogenization. In this work we
propose a graph convolutional neural network that utilizes the discretized
representation of the initial microstructure directly, without segmentation or
clustering. Compared to feature-based and pixel-based convolutional neural
network models, the proposed method has a number of advantages: (a) it is deep
in that it does not require featurization but can benefit from it, (b) it has
a simple implementation with standard convolutional filters and layers, (c) it
works natively on unstructured and structured grid data without interpolation
(unlike pixel-based convolutional neural networks), and (d) it preserves
rotational invariance like other graph-based convolutional neural networks. We
demonstrate the performance of the proposed network and compare it to
traditional pixel-based convolution neural network models and feature-based
graph convolutional neural networks on multiple large datasets.

##  1 Introduction

Predicting the evolution of a system with a complex initial state represents a
wide class of physical problems of scientific and technological interest. For
instance, simulating the evolution of materials with complex microstructure is
necessary for predicting the behavior of highly engineered materials [1, 2, 3,
4, 5, 6]. With the advent of machine learning for physical applications and
the availability of considerable experimental and high-fidelity simulation
data, models and architectures for these and related applications have begun
to arise [7, 8, 9, 10]. These models can be used for a number of tasks such as
sub-grid accurate constitutive modeling [7], material design by structure
property exploration [11], and uncertainty quantification of materials with
high intrinsic variability [12]. For this work, we are interested in
predicting the evolution of the physical response of a sample given its
initial state and a history of loading. For this class of problems, we assume
the initial state can be represented as a field or collection of fields
captured in a multispectral/multichannel image and this image is data on a
structured grid or an unstructured mesh.

In the ever expanding field of machine learning (ML) [13, 14, 15], there are
many methods suitable to the task of supervised learning where the objective
is to represent an input-output map to high fidelity. Neural networks (NN)
[16, 17] are a particularly versatile sub-category of machine learning
techniques suitable for regression tasks. They can be designed to be smooth,
expressive models of physical behavior and have been shown to be effective in
reducing complex processes to low-dimensional latent spaces [18, 19]. In
particular, convolutional neural networks (CNNs) [20, 21, 22, 23] are an
efficient version of fully connected networks applied to image data. Their
main advantage is the reduction of the weight space needed to encode images to
a manageable size by exploiting spatial correlation, locality, and similarity
with a compact kernel filter. They have been spectacularly successful in a
variety of image (and time series) applications [24, 25, 26].

Graph convolutional neural networks (GCNNs) [27, 28] can be seen as
counterparts to CNNs which operate on topologically related, as opposed to
spatially or temporally related, data. Bruna et al. [27] recognized that many
of the properties that make CNNs effective on gridded data could be translated
to graph-based data, such as data on unstructured meshes. They developed both
local (e.g. pooling) and spectral/global (e.g. convolution in the Fourier
domain) operations on graph data. They made the connection with Fourier bases
through the graph Laplacian (of a binary adjacency matrix) to translate the
convolution operation to graphs. In another arena (graph wavelets/graph signal
processing), some of the graph convolution developments were preceded by
Hammond et al. [29] who developed the mathematics of spectral analysis and
filtering on the more general context of kernel-weighted graphs. In order to
surmount the global and expensive nature of applying filters in the spectral
domain, where an eigen-decomposition is required, Defferrard et al. [30]
(ChebNet) introduced (spatially) compact polynomial filters. In particular,
they approximated the action of a general spectral filter with a truncated
Chebychev expansion of the eigenvalue matrix of the graph Laplacian,
leveraging the fact that repeated application of these filters ($k$ times)
leads to $k$-hop diffusion of information. Kipf and Welling [31] took these
developments to their logical conclusion with the Graph Convolutional Network
(GCN). They truncated the Chebychev expansion to first order and relied on
deep/multi-layer networks to build expressive representations.222App. A
provides a brief synopsis of the GCN which we base our developments on. For a
more complete overview of GCNNs see reviews by Wu et al. [32]. Zhang et al.
[33], and Zhou et al. [34]. Note that related kernel-based NN methods, such as
the work of Trask et al. [35, 36] exists and have been shown to be effective
in physical applications. Also, Bronstein et al. [28] provides an insightful
perspective on applications of deep learning to non-Euclidean data, graphs,
and manifolds as well as the relevant mathematics.

One of the issues with using convolutional NNs for physical problems is the
need to preserve fundamental physical relationships such as frame
indifference, where scalars are unchanged (invariant) to rotational changes in
observer and tensorial quantities rotate in a corresponding way (equivariance)
with rotation of the observation coordinate frame. The advantages of
preserving these spatial symmetries in filter based NNs was recognized early
on [37]. By construction CNNs embed shift/translation invariance where every
local neighborhood is processed in a similar manner. Traditional convolutional
filters operate strictly on structured grids of pixelated images and generally
do not preserve rotational symmetries but provide a richer set of filters than
methods that do. Numerous efforts have been made to endow CNNs with rotational
invariance. Dielemann et al. [38] augmented the layers of a CNN to add
rotations of the image by $\pi/4$ and inversions. Cohen and Welling [39]
outlined the mathematics whereby convolutional network will be equivariant
with respect to any group, including rotation and reflections. Worrall et al.
[40] developed filters based on circular harmonics. Chidester et al. [41] used
a (discrete) Fourier transform to embed rotational invariance by filter
augmentation. Using concepts from abstract algebra, Kondor and Trivedi [42]
proved that convolutional structure is a necessary and sufficient condition
for equivariance to the action of a compact symmetry group. This finding
serves as a requirement for CNN to be equivariant. Recently, Finzi et al. [43]
provided a significant extension of Cohen and Welling treatment of small,
discrete symmetry groups to continuous (Lie) groups, e.g. rotations in 3D
(the special orthogonal group SO(3)). In contrast to pixel-based CNNs, the
existing graph based filters can operate on unstructured spatial data by use
of user-defined neighbor attributions but lose some spatial information in the
process. Furthermore, graph-based convolutional networks can have an inherent
rotational invariance (assuming the node data is invariant/equivariant) since
the representation has been lifted out of its spatial embedding.

Currently, there are many applications of pixel-based CNNs to physics and
mechanics problems that have data on a structured grid. Some of the early work
focused on classification and simple property prediction. Chowdhury et al.
[44] used a CNN to classify microstructures. Lubbers et al. [45] developed a
CNN based method for inferring low-dimensional microstructure representations
and used the model for microstructure generation. DeCost et al. [46] also
applied CNNs to microstructure representations and used the model to connect
processing to structure. Contemporaneous with these developments, Kondo et al.
[47] used a CNN to predict ionic conductivity based on microstructure. Many
publications have followed these initial applications, some have targeted
evolution prediction tasks. Frankel and collaborators devised a hybrid between
a CNN and recurrent NN (RNN) network to predict the elastic-plastic stress
response of polycrystalline samples [7] and have used a convolutional long
short-term memory architecture (convLSTM) [48], which integrates CNN and RNN
aspects into a single architecture, to predict the evolution of the stress
field [8]. The latter work [8] was based on convLSTM network of Shi et al.
[49] which was developed for atmospheric predictions and provides a framework
for representing the solutions of time dependent partial differential
equations. On the other hand, we are aware of only a few applications of GCNNs
to physics or mechanics to date. Vlassis et al. [9] employed a feature-based
graph neural network to model the elastic interaction of grains in
polycrystalline samples using an adaptation of the CNN-RNN network in Ref.
[7]. Chen et al. [50] employed a GCNN trained to sparse (and unstructured)
diffusion data to model human tissue.

In contrast to these approaches, the proposed graph based convolutional neural
network processes the structured or unstructured microstructural images
directly and in an invariant manner. The adjacency matrix, which
conveys/defines neighbors and spatial information, is based on the topology of
the image data itself. The formulation does not require an obvious
segmentation of the microstructure, which may be occluded by noise in real
data. Most importantly it does not require featurization of the multi-
channel/hyperspectral image data [51, 52]. Nevertheless, as we will show, it
can benefit from obvious features but no feature engineering is needed to
obtain good accuracy. The main contributions of the work are: a generalization
of graph/pixel CNNs for predicting homogenized response of samples with
complex microstructure, and a demonstrative comparison of performance relative
to existing methods.

In Sec. 2 we describe the physical problem which is a homogenization of
physical response suitable for sub-grid/multi-scale applications [53, 7]. In
particular we apply the methods to homogenization of the evolution of stress
of a sample volume where internal state determined by microstructure. In Sec.
3 we describe the proposed neural network architecture and relate it to
traditional CNNs and feature-engineered GCNNs. In Sec. 4 we focus on
comparing the performance of pixel-based convolutional neural network (CNN), a
feature-based “reduced” graph convolutional neural network (rGCNN), and the
proposed “direct” graph-based convolutional neural network (dGCNN) which
operates directly on structured or unstructured image data like a CNN without
the need for featurization or segmentation required by the rGCNN. To assess
the performance of these models, we employ two exemplars of the homogenization
problem: (a) the prediction of the stress evolution in a porous metal and (b)
the prediction of the stress evolution in a polycrystalline material. With
these physical problems we created four datasets using: (a) an ensemble of
three-dimensional (3D) realizations of the porous metal, (b) an ensemble of
two-dimensional (2D) realizations of polycrystals, (c) a 3D ensemble with low
variance in the grain sizes, and (d) a 3D ensemble with high variance in the
grain sizes, all of which are subjected to a tension deformation process. The
first dataset is used to demonstrate the efficacy of the proposed mesh-based
GCNNs on unstructured meshes; the remainder are used to compare CNNs to GCNNs
on structured meshes. We exploit the lower computational cost of the 2D
polycrystalline dataset in a deeper exploration of variants and hyper-
parameters than would be possible with the 3D data. In Sec. 5 we summarize
the findings of exploring data and parameter efficiencies, architecture
variations, and adjacency manipulation, and discuss ongoing work.

##  2 Physical problem

Predicting the physical response of a sample given a complex initial state is
representative of a general class of problems in homogenization [54, 55]. In
particular, we focus on the prediction of the evolution of the volume average
of stress ${\boldsymbol{\sigma}}$ in a representative volume $V$ :

$\bar{{\boldsymbol{\sigma}}}(t)=\frac{1}{V}\int{\boldsymbol{\sigma}}\left({\boldsymbol{\epsilon}}(t),{\boldsymbol{\phi}}(\mathbf{X})\right)\,\mathrm{d}\mathbf{X}\
,$ (1)

given a microstructural field ${\boldsymbol{\phi}}(\mathbf{X})$ observed at
time $t=0$ and a time-dependent, homogeneous loading determined by the imposed
strain ${\boldsymbol{\epsilon}}(t)$. This class of problems is the basis for
multiscale models [56], material structure optimization [57], and material
variability and uncertainty quantification [12]. The microstructural field
${\boldsymbol{\phi}}(\mathbf{X})$ characterizes location-dependent
inhomogeneity that influences the state of the material, where $\mathbf{X}$ is
the position vector in the reference configuration of the sample. Examples of
${\boldsymbol{\phi}}$ include phase in a multiphase composites [58], elastic
modulus in a material with inclusions [59] or pores [60], and the local defect
density in a defective material [61].

Rotational equivariance requires

$\mathbf{Q}\boxtimes\bar{{\boldsymbol{\sigma}}}(t,{\boldsymbol{\phi}})=\frac{1}{V}\int{\boldsymbol{\sigma}}\left(\mathbf{Q}\boxtimes{\boldsymbol{\epsilon}}(t),\mathbf{Q}\boxtimes{\boldsymbol{\phi}}(\mathbf{X})\right)\,\mathrm{d}\mathbf{X}\
,$ (2)

where $\mathbf{Q}$ is an orthogonal tensor (rotation) and $\boxtimes$ is the
Kronecker product, which is defined as
$\mathbf{Q}\boxtimes\mathbf{A}=\mathbf{Q}\mathbf{A}\mathbf{Q}^{T}=\sum_{i,j}A_{ij}\mathbf{Q}\mathbf{e}_{i}\otimes\mathbf{Q}\mathbf{e}_{j}$
for a second order tensor $\mathbf{A}$. In effect this is a requirement that
rotation of the inputs by $\mathbf{Q}$ must lead to a corresponding rotation
of the output by $\mathbf{Q}$. It has fundamental consequences for the form
that the function
${\boldsymbol{\sigma}}({\boldsymbol{\epsilon}},{\boldsymbol{\phi}})$ is
allowed to take [62, 53].

###  2.1 Exemplars

We employ two exemplars to test and demonstrate performance of the NN
architectures: (a) a porous metal and (b) a polycrystalline metal. The
mechanical response in both systems is complex. They undergo an elastic-to-
plastic transition with loading and heterogeneous deformation due to the
microstructure. For simplicity and data storage/memory considerations we focus
on the primary component, $\sigma(t)$, of the volume averaged stress
$\bar{{\boldsymbol{\sigma}}}$ for both exemplars. Each dataset was created
with standard, well-documented software packages [63, 64, 65].

####  2.1.1 Porous plasticity

In the porous metal exemplar, ${\boldsymbol{\phi}}(\mathbf{X})$ is the local
density field with ${\boldsymbol{\phi}}(\mathbf{X})=0$ in the pores and equal
to the density of the metal elsewhere, which we normalize to one. Ref. [66]
describes a similar material model. Here aluminum was chosen as a
representative material.

The metal response follows from a widely-employed $J_{2}$ elastic-plastic
model [67] where the stress $\mathbf{S}$ is given by a linear elastic rule:

$\mathbf{S}=\mathbb{C}:\mathbf{E}_{e}$ (3)

Here “:” is a double inner product that allows the 4th order elastic modulus
tensor $\mathbb{C}$ to map the elastic strain $\mathbf{E}_{e}$ to the stress
$\mathbf{S}$. Note that $\mathbf{E}_{e}$ is distinct from the applied strain
${\boldsymbol{\epsilon}}(t)$ driving the evolution of the sample. For an
isotropic material like aluminum the components of $\mathbb{C}$ reduce to

$[\mathbb{C}]_{ijkl}=\frac{E}{(1+\nu)}\left(\frac{\nu}{(1-2\nu)}\delta_{ij}\delta_{kl}+\frac{1}{2}(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})\right)$
(4)

which depends only on Young’s modulus $E=$ 59.2 GPa and Poisson’s ratio $\nu=$
0.33.

The plastic flow is derived from the von Mises yield condition

$\sigma_{\text{vm}}(\mathbf{S})-\check{\sigma}(\epsilon_{p})\leq 0$ (5)

which limits the elastic regime to a convex region in stress space and offsets
the elastic strain $\mathbf{E}_{e}$ from the total strain. Here
$\sigma_{\text{vm}}=\sqrt{\frac{3}{2}\mathbf{s}\cdot\mathbf{s}}$ is the von
Mises stress where
$\mathbf{s}=\mathbf{S}-{\operatorname{tr}}(\mathbf{S})\mathbf{I}$ is the
deviatoric part of $\mathbf{S}$, and $\epsilon_{p}$ is the equivalent plastic
which is a measure of the accumulated plastic strain computed from the plastic
velocity gradient $\mathbf{L}_{p}$. The yield limit $\check{\sigma}$ is given
by a Voce hardening law

$\check{\sigma}=Y-H\exp(-\alpha\epsilon_{p})$ (6)

with parameters: initial yield $Y=$ 200.0 MPa, hardening $H=$ 163.6 MPa, and
saturation exponent $\alpha=$ 73.3.

Realizations were created by a random placement scheme of spherical voids in
the sample cube with constraints on pore overlap with other pores and the
sample boundary [68]. This process created unit cells with mean porosity
0.09$\pm$0.03 following a beta distribution and at most 20 pores per cell. The
cubic samples were on the order of 1.53 mm3 with pore radius $\approx$ 150
$\mu$m (refer to Fig. 2, which will be discussed in Sec. 2.2). Pores in each
of the 1121 realizations we created were explicitly meshed and resulted in
unstructured discretizations with 14,640 to 101,360 elements.

Each realization was subjected to quasi-static uniaxial tension up to 20%
engineering strain with minimal Dirichlet boundary conditions (no lateral
constraints, uniform displacement on the ends). These simulations were
performed with Sierra [65]. From these simulations we extracted microstructure
${\boldsymbol{\phi}}(\mathbf{X})$, applied strain $\epsilon(t)$, volume-
averaged stress $\bar{{\boldsymbol{\sigma}}}(t)$ data to demonstrate the
efficacy of mesh-based GCNNs in the Results section.

####  2.1.2 Crystal plasticity

Our second exemplar of the homogenization problem, Eq. (1), used a crystal
plasticity (CP) constitutive model where ${\boldsymbol{\phi}}(\mathbf{X})$ is
a field of crystal orientations associated with grains. Although each grain
has a relatively simple response, the collective behavior is difficult to
predict without a detailed simulation since each grain influences its
neighbors [7]. For this exemplar, steel was chosen as a representative
material.

The response of each crystal follows an elastic-viscoplastic constitutive
relation based on well-known meso-scale models [69, 70, 71, 72, 73, 74, 75].
For the crystal elasticity, we employed the same linear stress model, Eq. (3),
as in the porous metal exemplar albeit with a different elastic modulus tensor
$\mathbb{C}$. In this case ferrous (face centered) cubic symmetry for
$\mathbb{C}$ was assumed, and the independent components of the elastic
modulus tensor $\mathbb{C}$ were those of steel:
$C_{11},C_{12},C_{44}=\\{204.6,137.7,126.2\\}$ GPa. The overall response
reflects the anisotropy of each grain; however, the response of polycrystals
with random orientations $\phi(\mathbf{X})$ was determined by the collective
response (which tends to isotropy with large sample sizes). In each crystal,
plastic flow

$\mathbf{L}_{p}=\sum_{\alpha}\dot{\gamma}_{\alpha}\mathbf{s}_{\alpha}\otimes\mathbf{n}_{\alpha}$
(7)

can occur on any of the 12 face-centered cubic slip planes, where
$\mathbf{L}_{p}$ is the plastic velocity gradient, $\dot{\gamma}_{\alpha}$ is
the slip rate, $\mathbf{s}_{\alpha}$ is the slip direction, and
$\mathbf{n}_{\alpha}$ is the slip plane normal. We employed a common power-law
form for the slip rate relation

$\dot{\gamma}_{\alpha}=\dot{\gamma}_{0}\left|\frac{\tau_{\alpha}}{g_{\alpha}}\right|^{m-1}\tau_{\alpha}\
,$ (8)

driven by the shear stress $\tau_{\alpha}$ resolved on slip system $\alpha$.
The reference slip rate was chosen to be $\dot{\gamma}_{0}$ = 1.0 s-1, the
rate sensitivity exponent was $m=20$, the initial slip is set to zero, and the
slip resistance $g_{\alpha}$ was given the initial value $g_{\alpha}$ = 122.0
MPa [53]. The slip resistance evolved according to [76, 77]

$\dot{g}_{\alpha}=(H-R_{d}g_{\alpha})\sum_{\alpha}|\dot{\gamma}_{\alpha}|$ (9)

where the hardening modulus was chosen to be $H=355.0$ MPa and the recovery
constant was $R_{d}=2.9$. Both Eq. (8) and Eq. (9) are integrated with
standard implicit numerical integrators as part of a global equilibrium
solution algorithm. See Ref. [53] for additional details.

For this exemplar, we created multiple sets of
$\\{{\boldsymbol{\phi}}(\mathbf{X}),\epsilon(t);\sigma(t)\\}$ data to train
and compare the NN models described in the next section: (a) a 2D dataset
consisting of 12,000 realizations [7], (b) a 3D dataset consisting of 10,000
realizations with low variance in the grain sizes, and (c) a corresponding
10,000 realization 3D dataset with high variance in the number of grains per
realization. The nominal sample length for each realization was 1 $\mu$m (a
realization is shown in Fig. 4a, which will be discussed in the next section)
. The variance in the grain sizes is directly related to the variety of grain
topologies in the particular ensemble, this aspect will be used in
explorations described in the Results section.

Although the GCNN method can be applied to unstructured grids and complex
geometries of different sizes, the CNN cannot without interpolation or some
other intermediate data processing. In this exemplar we chose structured
computational grids to facilitate direct comparison of pixel and graph based
methods. The 2D dataset was computed on a 32$\times$32 structured FE mesh and
output over 31 time steps (max strain 0.3%) ; and the two 3D CP datasets used
a 25$\times$25$\times$25 mesh and output over 51 steps (max strain 0.4%). Each
polycrystal realization was subjected to quasi-static uniaxial tension at a
constant engineering strain-rate of $\dot{\epsilon}=1$ s-1 with minimal
Dirichlet boundary conditions. The time evolution of each system was observed
over a limited number of time-steps that covered the physics of interest: the
transition from elastic to full plastic flow.

Realizations of the microstructure ${\boldsymbol{\phi}}(\mathbf{X})$ consisted
of a crystal orientation vector field that encodes the rotation of a crystal
in a reference orientation to that in the polycrystal [7]. The orientation
vector ${\boldsymbol{\phi}}$ is the unit eigenvector $\mathbf{p}$ of the
rotation tensor $\mathbf{R}$, which takes $\mathbb{C}$ from a canonical
orientation to that of a particular grain $\mathbf{R}\boxtimes\mathbb{C}$,
scaled by the rotation angle $\theta$ around that axis

${\boldsymbol{\phi}}=\theta\mathbf{p}\ \ \text{such that}\ \
\mathbf{R}\mathbf{p}=\mathbf{p}\ \text{and}\ \|\mathbf{p}\|=1,$ (10)

where $\theta$ can be obtained from the non-unitary eigenvalues of
$\mathbf{R}$. Refer to Ref. [7] for more details. The computational cell for
each realization is a cube, as can be seen in Fig. 4a, which is partitioned
into sub-regions, called grains, with distinct
${\boldsymbol{\phi}}(\mathbf{X})$. The sub-regions evoke a natural topology
for a grain-based graph [9]. All polycrystal realizations where created with
Dream3D [63] using spatial correlations to obtain a reasonable number of
grains and angle distribution functions that gave a uniform texture. The 2D
simulations were run with Albany [64] and 3D simulations were run with Sierra
[65]. See Refs. [53, 78, 79] for related efforts.

###  2.2 System characterization and response

In each exemplar, the non-linear plasticity model and heterogeneous
microstructure evoke a complex response to loading
${\boldsymbol{\epsilon}}(t)$. The local stress fields reflect internal
inhomogeneities at the pores or the grain boundaries, as characterized by
${\boldsymbol{\phi}}(\mathbf{X})$, and display large gradients at elastic-
plastic transitions. Spatial averaging to extract the system response
$\bar{{\boldsymbol{\sigma}}}(t)$ does some smoothing of the evolution but the
range of microstructures evokes a distribution of responses. The behavior is
generally similar across the ensemble of realizations but variations are
difficult to predict from simple statistics such as mean grain or pore size.

####  2.2.1 Porous plasticity

Fig. 1 shows the distribution of sample porosities over an ensemble of 1121
realizations. The distribution follows the target beta distribution with 0.09
mean porosity and a standard deviation of 0.03.

Due to variation in the number of voids, their size, and their placement
relative to each other and inside the geometry, the response varies as the
pores decrease the material stiffness and neighboring pores increase local
stress concentrations. The three representative realizations shown in Fig. 2
at the final strain of 20% display significant heterogeneities in their
deformation due to the dense packing of large pores. The middle realization,
in particular, clearly shows a plastic localization plane due to a collection
of pores leading to a weak section in the sample. The stress fields for each
of the samples display correspondingly large local variations and gradients.
The average stress histories, $\\{\sigma(t)\\}$, shown in Fig. 3 display a
variation of 22%. After rescaling by a mixture rule based on the solid
fraction, 10% variation remains. This indicates that the details of the pore
configurations control a significant portion of the plastic response.

Figure 1: Porosity distribution.  Figure 2: Pore realizations showing exterior
and interior surface mesh at 20 % strain colored by tensile stress (blue: $<$
0, red: 700 MPa). The original, undeformed configuration is outlined.  Figure
3: Porous elastic-plastic stress response. Color distinguishes the 64
realizations shown.

####  2.2.2 Crystal plasticity

The CP microstructure also evokes a complex and evolving stress field, as Fig.
4b–d illustrates.

Descriptive statistics of the 2D and 3D ensembles are shown in Fig. 5. The
distribution of grain densities of each realization (the reciprocal of the
number of grains in the particular realization) is relatively broad and long-
tailed in the high variance 3D dataset compared to that of the low variance 3D
dataset. The grain distribution of the 2D dataset has more compact support
relative to the high variance 3D ensemble but is significantly wider than the
low variance 3D ensemble, thus representing an intermediate distribution. The
distribution of the individual grain volumes over all realizations show
similar trends. The peak width of the high variance 3D dataset, however, is
relatively narrow compared to the high variance 3D dataset. The 2D data has a
pronounced tail and indistinct peak which indicates a wide variance in grain
sizes across the ensemble. Fig. 6 illustrates how the variance of the inputs
${\boldsymbol{\phi}}(\mathbf{X})$ is reflected in the variance of the output
$\sigma(t)$.

These datasets are particularly challenging to represent, relative to existing
work, since each realization had a distinct topology/grain assignment and
texture. Previous work [7, 9] employed a limited number of grain topologies
(the tiling of the domain by distinct subregions) with unique texture
assignments (the specific orientations assigned to the grains). Specifically,
here each dataset had on the order of 10,000 unique topologies, whereas in
Refs. [7, 9] there were on the order of 10–100 topologies. It is plausible, in
a small dataset, when the number of convolutional filters approaches the
number of microstructural topologies a simpler learning process results since
each filter can specialize to a particular topology. In this study this is not
the case since each of the 104 samples had both a unique grain structure and
grain orientation (texture). Clearly, with datasets of this size and variety,
the filters must learn generalized, predictive features.

Figure 4: Polycrystal orientation (a) colored by the first Euler angle, and
stress states: (b) elastic, (c) transition, (d) plastic, colored by
$\sigma_{11}$ with the same scale for all three panels (blue: 0 MPa, red: 250
MPa).  Figure 5: Comparative statistics for 2D and 3D high and low variance CP
ensembles. Grain densities per realization and volumes per grain have been
referenced to unit cell volumes in 2D and 3D, respectively.  Figure 6: CP
stress response for 2D (top), 3D low grain size variance (middle), and 3D high
variance (bottom) ensembles. Color distinguishes the 64 realizations shown.

##  3 Neural network architecture

The overall neural network (NN) architecture to model the problem of interest
Eq. (1) is analogous to the hybrid network from our previous work [7, 8] and
also used in Ref. [9]. It is illustrated in Fig. 7 and consists of two main
components: (a) a convolutional neural network (CNN or GCNN, orange) to
process the spatially complex initial state ${\boldsymbol{\phi}}(\mathbf{X})$
and (b) a recurrent neural network (RNN, blue) to evolve the quantity of
interest $\sigma(t)$ given the time-dependent loading $\epsilon(t)$ as an
input. As can be seen in Fig. 7, the output features of the CNN become inputs
to the RNN along with the time-dependent loading. This hybrid NN has the
potential to be quite complex with many hyper–parameters, e.g. a different
kernel width for each convolution layer and different number of nodes for each
dense layer.

In this work we simplify the NN architecture from that we previously employed
to facilitate exploration and comparison of CNN and GCNN approaches. This is
not particularly constraining since there are redundancies in the
approximation power of such a network. We determine the shape of the network
diagrammed in Fig. 7 with three hyper-parameters: (a) the number of filters,
$N_{\text{f}}$, applied in parallel to node data; (b) the number of
convolutional layers, $N_{\text{c}}$; and (c) the number of densely connected
layers, $N_{\text{d}}$, post-convolution used in the reduction of image to
(hidden) features. In the following we will use the abbreviation
$({N_{\text{f}}}\text{:}{N_{\text{c}}}\text{:}{N_{\text{d}}})$ to refer to
architectures defined by these three parameters, for instance
$({4}\text{:}{2}\text{:}{1})$ is a network with 4 filters, 2 convolutional
layers, and 1 dense layer. Note the global (average) pooling after the
convolutional layers reduces the output of the image processing CNN to
$N_{\text{f}}$ outputs and hence determines with width of the entire
convolutional component including the densely connected layers. Note, in our
previous work [7] we used an encoder for this task (c). Fig. 8 illustrates the
details of a typical convolutional component in the overall architecture
(yellow in Fig. 7). A batch normalization layer (not shown) is inserted
between each convolutional layer to aid in conditioning the output and
training [80]. To facilitate comparison of the CNN and GCNN approaches we fix
the kernel width of pixel convolutional layer to 3 to make an analog of a
graph where connections are nearest neighbors for the physical problem. The
particular RNN we employ is the well-known long-short-term memory unit (LSTM)
[81]. In preliminary studies we also tried another standard RNN, the gated
recurrent unit (GRU) [82], which achieved marginally worse accuracy on
average. The RNN used $\tanh$ activations and all other layers employ
rectifying linear units (ReLUs) for their nonlinear activation functions. We
also tried $\tanh$ activations for the entire network but that configuration
choice performed relatively poorly [83, 84]. Lastly we apply the standard
technique of using a linear mixing layer (no nonlinearity, only affine
transformation) just prior to output of interest.

The primary variation in the networks we explore is in the convolutional unit
(yellow in the schematic Fig. 7 and shown in detail in Fig. 8) and in
particular the construction of the convolution filter, which is our focus.

###  3.1 Proposed graph structure

At the level of digitized data in a pixel-based CNN, the microstructural input
is values of ${\boldsymbol{\phi}}$ at pixels (or a cell or an element) that
are addressed in a grid-wise fashion. In general the field
${\boldsymbol{\phi}}(\mathbf{X})$ can be represented as a pixelized image on a
structured grid or on an unstructured mesh. For a graph representation, such
as that used in Vlassis et al. [9], the reduction of a clearly segmentable
microstructure, such as that illustrated in Fig. 4a, leads to nodes
representing homogeneous regions and graph edges encoding adjacency between
regions. This reduction loses spatial information such as the shape of the
regions in the clustering/aggregation to nodes. Hence, that framework
typically requires enrichment of the node data $\mathsf{x}$ by featurization,
i.e. picking measures/statistics that quantify the information lost in the
reduction. In that approach, the number of nodes $N_{\text{nodes}}$ is
variable, and equal to the number of regions in the particular sample. In that
context, and in general, the node data $\mathsf{x}$ is number of nodes
$N_{\text{nodes}}$ by number of input features $N_{\text{f}}$ derived from
${\boldsymbol{\phi}}(\mathbf{X})$. Some features are obvious, such as the
value of ${\boldsymbol{\phi}}$ in the clustered region represented by the node
and the volume of the region, and others are not. It is easy to see that these
features can result from preconceived filters and clustering operations
applied to image ${\boldsymbol{\phi}}(\mathbf{X})$. We will refer to this
feature-based, reduced graph convolutional network as a rGCNN.

To avoid nebulous task of featurization, we propose a direct graph CNN (dGCNN)
where $\mathsf{x}$ is identically ${\boldsymbol{\phi}}(\mathbf{X})$ at the
cell/element centers of the computational grid but flattened (in an arbitrary
order) to fit in the graph convolution paradigm. The graph convolutions are
permutationally invariant so ordering does not affect output. In the proposed
network $N_{\text{nodes}}$ is the number of pixels or unstructured elements in
the image and the edges are derived directly from the mesh topology, i.e.
element/pixel neighbors are graph neighbors. This approach has qualitative
advantage of being a graph-based representation, which has intrinsic
invariance properties, while working directly on the structured or
unstructured data ${\boldsymbol{\phi}}(\mathbf{X})$ not a segmented or
clustered version of it. This representation does not preclude the use of
derived, informative features to boost accuracy by adding them to
$\mathsf{x}$, as we will demonstrate.

The proposed architecture was implemented with Spektral [86] and TensorFlow
[85].

###  3.2 Comparison of pixel and graph based convolution

To understand the difference between a pixel-based and a graph-based
convolution, let us examine a single convolutional filter. Note that the
proposed architecture shown in Fig. 7 and Fig. 8 has several filters. In
general, a convolutional filter has trainable weights $\mathsf{W}$ and bias
$\mathsf{b}$. The weight matrix $\mathsf{W}$ effects an affine transform by
matrix multiplication that mixes input features $\mathsf{x}$ into output
features $\mathsf{y}$ and $\mathsf{b}$ provides an offset to tune activation
of the subsequently applied nonlinearity $\mathsf{y}=f(\mathsf{x})$. In a
traditional grid/pixel convolution (CNN), the node features are addressable by
grid index, for instance in 2D:

$\mathsf{y}_{(i,j)}=\sum_{k^{\prime},l^{\prime}}\mathsf{W}_{(k^{\prime},l^{\prime})}\mathsf{x}_{(i+k^{\prime},j+l^{\prime})}+b\
,$ (11)

where $i,j$ are indices over the pixelated image/discretized field,
$k^{\prime},l^{\prime}$ indices over the convolution/filter which is
$N_{\text{kernel}}\times N_{\text{kernel}}$ in size, with
$N_{\text{kernel}}\ll N_{\text{nodes}}$ being the kernel width. To see the
similarities with graph convolution we recast this convolutional
multiplication in Eq. (11) by mapping both $(i,j)$ and
$(k^{\prime},l^{\prime})$ to single indices $I$ and $L$ by imposing a fixed,
arbitrary ordering across the kernel $\mathsf{W}$ (refer to Fig. 9). For the
flattened input $\mathsf{x}$, the corresponding output $\mathsf{y}$ is

$\mathsf{y}_{I}=\sum_{J}\left[\sum_{K^{\prime}}\mathsf{W}_{K^{\prime}}\mathsf{A}_{IJ}^{(K^{\prime})}\right]\mathsf{x}_{J}+\mathsf{b}_{K^{\prime}}\
,$ (12)

where multiplication by a matrix $\mathsf{A}^{(K^{\prime})}$ provides a
masking operation translating the global indices to local dependencies. Here
$\mathsf{A}^{(K^{\prime})}$ is a global $N_{\text{nodes}}\times
N_{\text{nodes}}$ adjacency matrix for each entry $K^{\prime}$ in the
convolutional kernel (i.e. each pixel under the filter kernel is treated
uniquely); and $\mathsf{A}_{IJ}^{(K^{\prime})}=1$ if $I$ and $J$ are neighbors
(by some definition e.g. shared face, shared node in the computational mesh
or distance) and 0 otherwise. The definition of neighbors determines the
direct interactions, in rough analogy to choosing the kernel width in a pixel-
based convolutional filter.

Likewise, for a graph convolutional network (GCN) layer [31], the convolution
operation takes the form:

$\mathsf{y}_{I}=\sum_{J}W\mathsf{A}_{IJ}\mathsf{x}_{J}+b\ ,$ (13)

where the adjacency $\mathsf{A}$ plays the same masking/connectivity role as
$\mathsf{A}^{(K^{\prime})}$ and an ordering of the input data $\mathsf{x}$ is
necessary to associate the indices $I$ and $J$ with particular cells. Again
$\mathsf{A}_{IJ}=1$ if $I$ and $J$ are neighbors by some definition and 0
otherwise. Based on the derivation of the GCN, refer to App. A, self-loops are
added such that $\mathsf{A}_{II}=1$. For each filter there are two trainable
parameters, $W$ and $b$, versus $k^{d}$ for a pixel-based CNN filter (where
$d$ is the spatial dimension and $k$ is the kernel width). These graph-based
filters have permutational invariance by construction since all neighbors have
same weight. They are also typically normalized by the number of neighbors,
for instance in a GCN, the binary adjacency matrix $\tilde{\mathsf{A}}$ is
normalized by the degree matrix
$\mathsf{D}_{IJ}=\sum_{I}\tilde{\mathsf{A}}_{IJ}\delta_{IJ}$

${\mathsf{A}}=\mathsf{D}^{-1/2}\tilde{\mathsf{A}}\mathsf{D}^{-1/2}$ (14)

to convey average neighborhood information and improve the conditioning of
$\mathsf{A}$.

Finally, in both the pixel and graph convolutional cases an activation
function $f$ is applied element-wise to the resulting node data $\mathsf{y}$

$\mathsf{x}^{\dagger}=f\left(\mathsf{y}(\mathsf{x})\right)$ (15)

to obtain the input $\mathsf{x}^{\dagger}$ to the next layer from that of the
present layer $\mathsf{x}$.

Clearly the adjacency $\mathsf{A}$ plays a crucial role in convolution by
defining what constitutes influence between nodes. In physical problems
typically influence is local and decays with distance. For some applications,
such as electrostatic interactions the interactions do not decay quickly with
distance and hence are long-ranged [87]. Most other applications have
relatively short-ranged influence, such as contact/interface interactions. As
mentioned, the source field for the microstructure is pixelated image grid or
computational mesh which have their own obvious topology and neighbors. In the
proposed direct graph we define adjacency by the pixel/cell neighbors of the
image grid/mesh not by the data ${\boldsymbol{\phi}}(\mathbf{X})$ on it. When
reduced to a graph in this fashion, where nodes are the pixels/cells, the only
topologically distinct nearest neighbors are: (a) those that share an edge
with the particular cell or (b) only a vertex. In CNNs compact kernels are
preferable since they limit the number of parameters that need to be trained.
For GCNNs, and in particular the GCN, neighbors are not distinct and hence
there is only a single weight. As in CNNs longer range influence can be
captured by multiple layers applied in sequence.

The Kipf and Welling [31] GCN is well-known to be an effective (non-spectral)
convolutional filter. By comparing it to standard CNNs we employed in our
previous work [78] we devised a few variants based on manipulating the node
adjacencies and associated weights. Fig. 9 shows the variety of rotationally
and permutationally invariant filters we explored. Note these filters assume
that the neighbors are at equivalent distances from the node at the center and
the cells are comparable sizes. We generalize this type of filter to multiple
adjacencies $\mathsf{A}^{(K^{\prime})}$ with their associated weights
$\mathsf{W}_{K^{\prime}}$ via:

$\mathsf{y}_{I}=\sum_{J}\left[\sum_{K^{\prime}}\mathsf{W}_{K^{\prime}}\mathsf{A}_{IJ}^{(K^{\prime})}\right]\mathsf{x}_{J}\
.$ (16)

The pixel-based filter in a standard CNN treats every pixel in the kernel
independently (refer to Fig. 9a) and retains a sense of how pixels located
relative to the central pixel. This richness has the side effect of not
satisfying the symmetries required by invariance. To elaborate with an
illustration, if the data
$\\{{\boldsymbol{\phi}},{\boldsymbol{\epsilon}};{\boldsymbol{\sigma}}\\}$ is
rotated by $\pi/2$ around a coordinate axis the data at pixels are uniquely
mapped by $\mathsf{Q}$ to new indices without interpolation
$\mathsf{x}^{*}=\mathsf{Q}\mathsf{x}$. In this simple case the linear
transformation

$\mathsf{Q}\mathsf{y}_{I}\equiv\mathsf{Q}\sum_{J}\left[\sum_{K^{\prime}}\mathsf{W}_{K^{\prime}}\mathsf{A}_{IJ}^{(K^{\prime})}\right]\mathsf{x}_{J}\neq\sum_{J}\left[\sum_{K^{\prime}}\mathsf{W}_{K^{\prime}}\mathsf{A}_{IJ}^{(K^{\prime})}\right]^{*}\mathsf{Q}\mathsf{x}_{J}$
(17)

in Eq. (12) applied to the inputs $\mathsf{x}_{J}=\phi(\mathbf{X}_{J})$ will
not produce a rotated output since image $\mathsf{x}=\phi(\mathbf{X})$ rotates
but the filter $\left[\mathsf{W}_{K^{\prime}}\mathsf{A}^{(K^{\prime})}\right]$
for the $K^{\prime}$-th adjacency does not. A graph treatment, unlike a pixel
convolution, effectively allows the transformation $\mathsf{Q}$ to commute
with the adjacency,
$\mathsf{Q}\mathsf{A}^{(K^{\prime})}\mathsf{x}=\left[\mathsf{A}^{(K^{\prime})}\right]^{*}\mathsf{Q}\mathsf{x}$,
since the node association of the graph adjacency is invariant to spatial
transformations. Refer to Refs. [42, 43] for a more formal treatment. The
equivariance problem remains for the action of the weights
$\mathsf{W}_{K^{\prime}}$.

Permutational invariance of the weights is one means of solving this issue.
Strictly speaking, the weights of all of the neighbors of the central pixel
are required to be identical for permutational invariance (and for invariance
with respect to arbitrary rotations). We implemented filter variants where
edge and vertex weights are used exclusively (Fig. 9f uses edge neighbors
whereas Fig. 9e employs only node neighbors) or given independent weights,
such as the “$\ast$” pattern in Fig. 9b. These filters are in contrast to the
CNN filter, shown in Fig. 9a, which has no inherent symmetries. Note there is
some evidence that with sufficient data CNNs learn rotational invariance [88]
but the benefits of a smaller parameter space and a compact representation
that satisfies this constraint exactly are clear. In the GCN, as designed by
Kipf and Welling [31], the weight of the center (shown in gray in Fig. 9f vs.
colored in Fig. 9c) was chosen to be the same as the neighbors (refer to App.
A). In this work we also tried variants where the self-weight is independent
of the neighbor weight (Fig. 9c vs. Fig. 9f). As in Eq. (12), each
independent, trainable weight $\mathsf{W}_{K^{\prime}}$ can be associated with
a different pre-determined, binary adjacency matrix
$\mathsf{A}^{(K^{\prime})}$, or, equivalently, the weighted adjacency
$\sum_{K^{\prime}}\mathsf{A}^{(K^{\prime})}_{IJ}\mathsf{W}_{K^{\prime}}$ can
be considered the trainable entity.

Figure 7: Hybrid neural network with CNN (yellow) and RNN (blue) components to
transform the inputs (red): spatially-dependent
${\boldsymbol{\phi}}(\mathbf{X})$ microstructure and time-dependent loading
history ${\boldsymbol{\epsilon}}(t)$, to the output
${\boldsymbol{\sigma}}({\boldsymbol{\phi}},t)$ (cyan).  Figure 8: A (4:2:2)
convolutional unit. $\mathsf{A}^{(i)}$ is the adjacency matrix for the $i$-th
neighbors (0 is self), $N_{\text{f}}$ is the number of filters, and $f$ is the
activation function for the layer. Note the last dense layer is linear (no
non-linear activation function).

(a) pixel

(b) $\ast$

(c) #

(d) O

(e) X

(f) +

Figure 9: Convolutional filters in 2D: (a) CNN where all neighbors have
independent weights (denoted by colors); and alternative GCN-like filters that
preserve invariant outputs: (b) $\ast$: edges and nodes have separate weights
(c) #: all edges have a single weight. (d) O: all neighbors have a single
weight. (e) X: all node neighbors have single weight. (f) +: all edges have a
single weight. Edge neighbors share an edge/face, whereas node neighbors only
share a vertex/node. Color: trainable, gray: determined by trainable, white:
not used.

##  4 Results

To establish the efficacy of the proposed GCNN-RNN architecture applied
directly to unstructured mesh data without featurization we demonstrate its
performance on the porous material dataset. Since hyper-parameter optimization
and other computationally intensive tasks are more feasible with the 2D
dataset, we use it explore a variety of hyper-parameters and architecture
choices with the 2D dataset. After determining what graph-based filters
perform well and have good accuracy per parameter and dataset size, we
demonstrate the proposed architecture on the two 3D CP datasets.

The data was conditioned to aid training. The input data
${\boldsymbol{\phi}}(\mathbf{X})$ and $\epsilon(t)$ were normalized by their
maximum values since both had lower bounds of zero. The output data
$\sigma({\boldsymbol{\phi}}_{a},t)$ was transformed to the difference
$\Delta\sigma({\boldsymbol{\phi}}_{a},t)$ between the data
$\sigma({\boldsymbol{\phi}}_{a},t)$ and its mean trend
$\langle\sigma\rangle=\frac{1}{N}\sum_{a}{\sigma}({\boldsymbol{\phi}}_{a},t)$
($N$ is the number of realizations) and normalized by the standard deviation
of $\sigma({\boldsymbol{\phi}}_{a},t)$ over time $t$. We chose to train the
model response $\hat{\sigma}({\boldsymbol{\phi}}_{a},t)$ to the difference
from the mean trend to emphasize the variation in response between
microstructures ${\boldsymbol{\phi}}_{a}$. We evaluated the performance of
each of the convolutional networks primarily with the root mean squared error
(RMSE)

$\varepsilon(t)=\frac{1}{\max\langle\sigma\rangle}\sqrt{{\sum_{a}(\hat{\sigma}({\boldsymbol{\phi}}_{a},t)-\sigma({\boldsymbol{\phi}}_{a},t))^{2}}}$
(18)

relative to the maximum $\sigma$ over the dataset, and the (Pearson)
correlation coefficient

$C(t)=\frac{\sum_{a}\Delta\hat{\sigma}({\boldsymbol{\phi}}_{a},t)\,\Delta\sigma({\boldsymbol{\phi}}_{a},t)}{\sqrt{\sum_{a}\Delta\hat{\sigma}({\boldsymbol{\phi}}_{a},t)^{2}\
\sum_{a}\Delta{\sigma}({\boldsymbol{\phi}}_{a},t)^{2}}}$ (19)

of the normalized data. In Eq. (18) and Eq. (19) the sum is over all
realizations $a$ in the test set, and the NN model response is denoted as
$\hat{\sigma}({\boldsymbol{\phi}}_{a},t)$. We used a 70/10/20
train/validation/test split for the smaller porous material dataset and a
80/10/10 split for the larger CP datasets.

Recall the abbreviated architecture naming convention
$({N_{\text{f}}}\text{:}{N_{\text{c}}}\text{:}{N_{\text{d}}})$, which will be
used throughout this section.

###  4.1 Demonstration on unstructured mesh data

We used the smaller 3D porous material dataset to demonstrate that GCNNs
applied directly to unstructured mesh data are effective at representing
homogenized material behavior. For this study we used a (32:4:2) convolutional
unit with “O” type filters that treat all nearest neighbor elements equally.
The field data ${\boldsymbol{\phi}}(\mathbf{X})$, in this case, is the binary
density field (0: void, 1: metal) on the native unstructured computational
mesh. Since element and systems sizes varied across the ensemble, we augmented
${\boldsymbol{\phi}}(\mathbf{X})$ with the volumes. The stress response
$\bar{\sigma}(t)$ is comprised of 400 tensile loading steps to reach 20%
strain. Since each realization has a different mesh and each adjacency matrix
is large (on the order of 105-106 rows and columns) albeit sparse, we trained
the network by evaluating and updating the network weights one sample at a
time, i.e. a batch size of 1.

The predictions for 8 randomly selected trajectories shown in Fig. 10 are
smooth and display minimal errors that are fairly uniform across the
evolution. The initial elastic regime is well captured, as is the ultimate
(peak) strength and subsequent plastic flow. Over the entire test set of 224
samples, the normalized mean RMSE was 0.00409 with 0.00052 standard deviation,
and the mean correlation was 0.995. Clearly the architecture shown in Fig. 7
with graph convolution layers applied to unstructured field data is an
effective model of this microstructural response.

Figure 10: Comparison of true (lines) and predictions (points with
corresponding color) for 8 realizations of the porous metal data.

###  4.2 Efficacy of the components of hybrid network

In this study we vary the three chosen architecture hyper-parameters:
$N_{\text{f}}$, $N_{\text{c}}$ and $N_{\text{d}}$ for three architectures: (a)
a CNN, (b) the proposed dGCNN, and (c) a rGCNN endowed with angle and volume
features. Both GCNN architectures employ the standard GCN filter, with a “+”
pattern and dependent center weight illustrated in Fig. 9f. All models are
trained to the 2D CP dataset.

Fig. 11 shows the correlations $C(t)$ and errors $\varepsilon(t)$ over time
for the three architectures for a range of filters
$N_{\text{f}}\in\\{1,\ldots,6\\}$ for 1 or 2 convolutional layers
($N_{\text{c}}=\\{1,2\\}$) and one dense layer ($N_{\text{d}}=1$). Note that
two dense layers, $N_{\text{d}}=2$, produced similar results. Generally the
correlation of all three hybrid CNN-RNNs is better earlier in the process,
while the error tends to peak early on near the end of the elastic regime
where the response variance is highest. Referring to Fig. 6a, we observe that
elastic-to-plastic transition occurs around 0.1% and this transition is
apparent in Fig. 11 where the correlation appears to transition between two
plateaus. All architectures improve with more filters although there is a
clearly a limit to the improvement, which suggests a small number of relevant
features. It is also apparent that more filters are needed to capture the
later plastic regime accurately than the initial elastic regime.

The simplest CNN models (fewest parameters) that have the best performance
are: $({4}\text{:}{1}\text{:}{1})$ with 269 parameters,
$({3}\text{:}{2}\text{:}{1})$ with 271 parameters (shown in Fig. 11), and
$({3}\text{:}{1}\text{:}{2})$ with 187 parameters,
$({2}\text{:}{2}\text{:}{2})$ with 151 (not shown). This demonstrates the
fungibility of the nodes in the network. The proposed dGCNN with
$N_{\text{c}}=1$ with either $N_{\text{d}}=\\{1,2\\}$ does well in the elastic
regime for $N_{\text{f}}>3$ but achieves no greater than $0.7$ correlation for
plastic for all cases with $N_{\text{f}}<9$; however, with a second
convolutional layer three filters are sufficient to achieve the accuracy of
the best CNN architectures. This indicates that second nearest neighbors are
required to represent the plastic flow well using the standard GCN filter in
the dGCNN architecture. This finding will be revisited in Sec. 4.3. Lastly,
all the variants of rGCNN with angle and volume node features gave similar
(poor) performance. The rGCNN clearly require more features for improvement.
This finding will be expanded on in Sec. 4.4.

Table 1 compares the number of trainable parameters for the three types of
convolutional networks. Clearly, the additional weights in the pixel-based
convolutional layers incur a cost in complexity and training. Also the
parameter complexity of the dGCNN is essentially equivalent to the rGCNN since
the same filters are being used on different graphs and data. The data is
certainly different, with the rGCNN storing more features on a more compact
adjacency than the dGCNN; with sparse storage this is not a significant
advantage for moderately sized meshes.

| CNN | dGCNN | rGCNN
---|---|---|---
$({1}\text{:}{1}\text{:}{1})$ | 41 | 26 | 27
$({2}\text{:}{1}\text{:}{1})$ | 99 | 75 | 71
$({4}\text{:}{1}\text{:}{1})$ | 269 | 209 | 213
$({6}\text{:}{1}\text{:}{1})$ | 511 | 421 | 427
$({1}\text{:}{2}\text{:}{1})$ | 55 | 32 | 33
$({2}\text{:}{2}\text{:}{1})$ | 145 | 69 | 85
$({4}\text{:}{2}\text{:}{1})$ | 433 | 245 | 249
$({6}\text{:}{2}\text{:}{1})$ | 865 | 487 | 493
Table 1: Parameter counts for architectures with GCN filters applied to 2D CP
data

(a) CNN

(b) dGCNN

(c) rGCNN

Figure 11: Architecture comparison (2D CP data): (a) CNN, (b) dGCNN: GCN
applied directly grid, (c) rGCNN: GCN applied to clustered/segmented data.
Left panels: 1 convolutional layer, 1 dense layer; right panels: 2
convolutional layers, 1 dense layer.

###  4.3 Comparison of convolutional filters

Motivated by the fact that a CNN can achieve good performance with only one
convolutional layer we tried richer variants of the standard GCN filter (where
the self-weight is set equal to the neighbor weight). Using the patterns shown
in Fig. 9, we explored their relative performance with one convolutional layer
($N_{\text{c}}=1$) and four filters ($N_{\text{f}}=4$). Fig. 12a shows that
patterns +, X, O, which only have one trainable weight and have interactions
with edge, vertex-only, and edge+vertex neighbors, respectively, have
comparable and less than satisfactory performance. Even the $\ast$ pattern,
where edge and node neighbors are given separate weights (2 independent
weights) has an inferior performance to the CNN (which has 9 independent
weights). Only the # pattern, where the center pixel is given an independent
weight from the vertex neighbors, has performance on par with the CNN. Given
this finding we endowed each of the basic patterns {+, X, O, $\ast$} with an
independent central weight. As shown in Fig. 12b, this was sufficient to
improve the performance of all but the X (vertex only neighbors) to be
comparable with the CNN. It is physically plausible that that edge neighbors
have a stronger influence than vertex neighbors and it appears that they are
required for the filters to learn predictive interactions.

(a)

(b)

Figure 12: Filter comparison (2D CP data): (a) patterns illustrated Fig. 9 and
(b) edge,vertex, and both filters augmented with an independent center weight.

###  4.4 Selected features

In Sec. 4.2 we observed that the feature-dependent rGCNN formulation had sub-
par performance especially in the plastic regime that was not improved with a
more complex GCNN component. Fig. 13 shows the performance of a
$({4}\text{:}{2}\text{:}{1})$ rGCNN network does improve with an expanded
feature set. Here, in addition to the orientation angle and the volume
associated with the clusters/grains represented by the graph nodes (2 total
features), we added the surface area of each grain (3 total features) and the
area of the grain that is on the surface of the cell (4 total features) to the
node features. The improvement with these additional features is marginal.
However, if we allow for an independent self-weight by changing the standard
GCN pattern to the # pattern, as in the previous section, the performance is
dramatically improved and the improvement with additional features is more
distinct. Although these networks have considerably smaller adjacency matrices
than dGCNNdue to the clustering of the pixels, the performance is sub-par,
particularly in the plastic regime. This serves as an illustration of the
difficulty of improvement by feature selection, as opposed to deep
learning.333Note only elastic response was modeled in Ref. [9].

Figure 13: Comparison of rGCNN with increasing number of features (2D CP
data). Features: angle, volume, area, surface area. Graph labels correspond to
2: {angle, volume}, 3: {angle, volume, area}, and 4: {angle, volume, area,
surface area}. Dashed lines: self-weight equal to neighbor weight; solid
lines: self-weight independent of neighbor weight.

###  4.5 Data efficiency

As a last trial with the 2D CP dataset, we investigated how efficient the best
networks are with smaller datasets. Here we compared (a)
$({4}\text{:}{1}\text{:}{1})$ CNN with 269 parameters, (b)
$({4}\text{:}{2}\text{:}{1})$ dGCNN with the + pattern (GCNN+) and 317
parameters, and (c) $({4}\text{:}{1}\text{:}{1})$ dGCNN with the # pattern
(GCNN#) and 249 parameters. The training set was reduced from 80% of the
12,000 realizations (9600) by a fraction that ranged from 0.01 to 1.0. The
test set was a fixed 20% (2400) of the full realizations and the results where
averaged over 9 trials. Fig. 14 shows that the majority of learning
(improvement in accuracy) occurs by the time the training size is
approximately equal to the number of parameters. After the step-down in error
(at 0.04 of the total training set for the CNN, at 0.02 for the GCNN+, and at
0.05 GCNN#) the improvement is relatively slow but steady. This data
demonstrates that these small networks can be effective at the prediction of
the homogenized response task, Eq. (1), with much smaller data sets.

Figure 14: 2D CP data: data efficiency comparison for
$({4}\text{:}{1}\text{:}{1})$ CNN with 269 parameters,
$({4}\text{:}{2}\text{:}{1})$ dGCNN+ 317 parameters, and
$({4}\text{:}{1}\text{:}{1})$ dGCNN# with 249 parameters. Note $1-C(t)$ is
being plotted in the upper panel.

###  4.6 Boosting with preconceived features

Now that we have discovered effective adjacencies and guidance on architecture
hyper-parameters, we turn to using the 3D CP datasets. Motivated by the fact
that some features of the image ${\boldsymbol{\phi}}(\mathbf{X})$ have obvious
bearing on the output due to physical reasoning, we boosted the dGCNN with
some of the features we employed with the solely feature-based rGCNN. The
proposed architecture can accommodate pre-selected features by simply
augmenting the image/cell microstructural field ${\boldsymbol{\phi}}$ with
additional channels. Fig. 15 shows the effect of adding a channel with the
volume fraction of the associated grain to each pixel. Clearly there is a
distinct and uniform benefit; however, it is somewhat marginal due to the fact
that the selected feature is likely, at least partially, redundant/correlated
with the output of the trainable filters. Additional means of augmenting with
more global data, such as the average grain size or equivalently the grain
density, via inputs concatenated to the pooling layer (gray in Fig. 7) output
going to the dense layers (green in Fig. 7) would also likely prove
beneficial.

Figure 15: The effect of boosting the microstructural input field with with
volume fraction.

###  4.7 Generalizability

Using CNN, GCNN+, and GCNN# with 32 filters, 1 or 2 convolutional layers, and
1 dense, layer Fig. 16 shows that the proposed architecture with the # filter
with an independent central weight can outperform a corresponding CNN and
GCNN+. The benefits of the more complex $({32}\text{:}{2}\text{:}{1})$
configuration over the $({32}\text{:}{1}\text{:}{1})$ appear to be most
significant for the two dGCNNs. It is also apparent that all the types of
convolutional neural networks perform better, at least in terms of
correlation, on the higher variance dataset than on the lower variance
dataset.

Now focusing on the GCNN#, Fig. 17 shows the distribution of RMSE errors (Eq.
(18)) is approximately Gaussian with some outlier errors above 5% for the high
variance dataset and 3% for the low variance dataset. A direct comparison of
the true and predicted values over a sequence of strains, shown in Fig. 18,
indicates that the GCNN# overpredicts values near the mean which may be due it
being harder to distinguish near-mean response microstructures from those that
produce extreme/outlier responses.

Following this conjecture, Fig. 19 illustrates that training on the low
variance ensemble and testing on the high variance ensemble (which also has a
different mean) does relatively poorly compared to the reverse. It appears
that the network generalizes well to different distributions of inputs if they
are in the span of the training set i.e. it does well at interpolation and
less well at extrapolation to potential out-of-distribution samples.

(a) low variance

(b) high variance

Figure 16: 3D data: performance of $({32}\text{:}{1}\text{:}{1})$ and
$({32}\text{:}{2}\text{:}{1})$ configurations of a CNN, dGCNN+, and dGCNN# for
the low and high variance datasets (2:GCNN+ denotes a
$({32}\text{:}{2}\text{:}{1})$ direct graph CNN using a “+” pattern).

(a) low variance

(b) high variance

Figure 17: 3D data: cumulative distribution of error for various strains
($({32}\text{:}{2}\text{:}{1})$ dGCNN#).

(a) low variance

(b) high variance

Figure 18: 3D data: distribution of true and predicted stress values for
various strains ($({32}\text{:}{2}\text{:}{1})$ dGCNN#).  Figure 19: 3D data:
generalizability results for train:test pairs ($({32}\text{:}{2}\text{:}{1})$
dGCNN#).

##  5 Discussion

The response of microstructures for use in multiscale simulations, structure-
property investigations, and uncertainty quantification can be accurately
modeled with graphs. The proposed formulation used the topology of the data
discretization directly instead of a segmentation or clustering of the image
data. This aspect should have particular advantages for image data where the
segmentation is not obvious, hard to compute, or is obscured by noise.
Furthermore, it has a simple implementation and avoids the need for feature
engineering, but can benefit from it. The architecture draws on both purely
graph-based networks and permutationally invariant convolutional filters. We
demonstrated that endowing the widely-used GCN filter [31] with an independent
self-weight (as suggested by the reduction of the ChebNet) can significantly
improve accuracy without adding additional layers and their parameters. The
independent self-weight allows for differencing the node data of the self and
its neighbors instead of only averaging. This can be seen as giving the filter
the ability to infer edge features between the central pixel and its
neighbors. For physical problems driven by gradients this change to the filter
is important. We also found that pixel edge neighbors are more crucial for a
predictive model than vertex-only neighbors. Lastly we were able to
demonstrate that small, efficient graph convolutional networks can be
effective at the task of predicting the homogenized evolution of complex
microstructure. This has significant applications in sub-grid constitutive
models in large scale simulations, structure-property property investigations,
and material uncertainty quantification.

An apparent downside of the proposed approach is the graph and its adjacency
grows with resolution of image (number of pixels/elements). This issue is
partially offset by sparse storage of the adjacency matrix, in general, and
largely ameliorated by data that is on the same discretization. In future work
we will investigate low-rank approximations to the adjacency matrix [89, 90,
91, 92, 93], dimensionality reduction techniques [94, 95], and the use of
graph auto-encoders [96, 97, 98, 99] to reduce the mesh-based graphs in-line.
We are also pursuing the larger topic of processing image with multi-
resolution filters [100], e.g. spanning the pixel to the cluster level.

## Acknowledgments

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Advanced Scientific Computing Research program. Sandia
National Laboratories is a multimission laboratory managed and operated by
National Technology and Engineering Solutions of Sandia, LLC., a wholly owned
subsidiary of Honeywell International, Inc., for the U.S. Department of
Energy’s National Nuclear Security Administration under contract DE-
NA-0003525. This paper describes objective technical results and analysis. Any
subjective views or opinions that might be expressed in the paper do not
necessarily represent the views of the U.S. Department of Energy or the United
States Government.

## References

  * [1] T Kraft, M Rettenmayr, and HE Exner.  An extended numerical procedure for predicting microstructure and microsegregation of multicomponent alloys.  Modelling and Simulation in Materials Science and Engineering, 4(2):161, 1996.
  * [2] Xiaolei Yin, Wei Chen, Albert To, Cahal McVeigh, and Wing Kam Liu.  Statistical volume element method for predicting microstructure–constitutive property relations.  Computer methods in applied mechanics and engineering, 197(43-44):3516–3529, 2008.
  * [3] Somnath Ghosh and Dennis Dimiduk.  Computational methods for microstructure-property relationships, volume 1.  Springer, 2011.
  * [4] Ole Stenzel, Omar Pecho, Lorenz Holzer, Matthias Neumann, and Volker Schmidt.  Predicting effective conductivities based on geometric microstructure characteristics.  AIChE Journal, 62(5):1834–1843, 2016.
  * [5] Yulan Li, Shenyang Hu, Xin Sun, and Marius Stan.  A review: applications of the phase field method in predicting microstructure and property evolution of irradiated nuclear materials.  npj Computational Materials, 3(1):1–17, 2017.
  * [6] Carl Herriott and Ashley D Spear.  Predicting microstructure-dependent mechanical properties in additively manufactured metals with machine-and deep-learning methods.  Computational Materials Science, 175:109599, 2020.
  * [7] Ari L Frankel, Reese E Jones, Coleman Alleman, and Jeremy A Templeton.  Predicting the mechanical response of oligocrystals with deep learning.  Computational Materials Science, 169:109099, 2019.
  * [8] Ari Frankel, Kousuke Tachida, and Reese Jones.  Prediction of the evolution of the stress field of polycrystals undergoing elastic-plastic deformation with a hybrid neural network model.  Machine Learning: Science and Technology, 2020.
  * [9] Nikolaos N Vlassis, Ran Ma, and WaiChing Sun.  Geometric deep learning for computational mechanics Part I: Anisotropic hyperelasticity.  Computer Methods in Applied Mechanics and Engineering, 371:113299, 2020.
  * [10] Anup Pandey and Reeju Pokharel.  Machine learning enabled surrogate crystal plasticity model for spatially resolved 3d orientation evolution under uniaxial tension.  arXiv preprint arXiv:2005.00951, 2020.
  * [11] Juhwan Noh, Jaehoon Kim, Helge S Stein, Benjamin Sanchez-Lengeling, John M Gregoire, Alan Aspuru-Guzik, and Yousung Jung.  Inverse design of solid-state materials via a continuous representation.  Matter, 1(5):1370–1384, 2019.
  * [12] M Khalil, GH Teichert, C Alleman, NM Heckman, RE Jones, K Garikipati, and BL Boyce.  Modeling strength and failure variability due to porosity in additively manufactured metals.  Computer Methods in Applied Mechanics and Engineering, 373:113471, 2021.
  * [13] Christopher M Bishop.  Pattern recognition and machine learning.  springer, 2006.
  * [14] Trevor Hastie, Robert Tibshirani, Jerome Friedman, and James Franklin.  The elements of statistical learning: data mining, inference and prediction.  The Mathematical Intelligencer, 27(2):83–85, 2005.
  * [15] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [16] John J Hopfield.  Neural networks and physical systems with emergent collective computational abilities.  Proceedings of the national academy of sciences, 79(8):2554–2558, 1982.
  * [17] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep Learning.  MIT Press, 2016.  http://www.deeplearningbook.org.
  * [18] Kookjin Lee and Kevin Carlberg.  Deep conservation: A latent-dynamics model for exact satisfaction of physical conservation laws.  arXiv preprint arXiv:1909.09754, 2019.
  * [19] Lawson Fulton, Vismay Modi, David Duvenaud, David IW Levin, and Alec Jacobson.  Latent-space dynamics for reduced deformable simulation.  In Computer graphics forum, volume 38, pages 379–391. Wiley Online Library, 2019.
  * [20] Yann LeCun, Bernhard Boser, John S Denker, Donnie Henderson, Richard E Howard, Wayne Hubbard, and Lawrence D Jackel.  Backpropagation applied to handwritten zip code recognition.  Neural computation, 1(4):541–551, 1989.
  * [21] Keiron O’Shea and Ryan Nash.  An introduction to convolutional neural networks.  arXiv preprint arXiv:1511.08458, 2015.
  * [22] Saad Albawi, Tareq Abed Mohammed, and Saad Al-Zawi.  Understanding of a convolutional neural network.  In 2017 International Conference on Engineering and Technology (ICET), pages 1–6. IEEE, 2017.
  * [23] Zhuwei Qin, Fuxun Yu, Chenchen Liu, and Xiang Chen.  How convolutional neural network see the world-a survey of convolutional neural network visualization methods.  arXiv preprint arXiv:1804.11191, 2018.
  * [24] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [25] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  Advances in neural information processing systems, 25:1097–1105, 2012.
  * [26] Weibo Liu, Zidong Wang, Xiaohui Liu, Nianyin Zeng, Yurong Liu, and Fuad E Alsaadi.  A survey of deep neural network architectures and their applications.  Neurocomputing, 234:11–26, 2017.
  * [27] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  arXiv preprint arXiv:1312.6203, 2013.
  * [28] Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst.  Geometric deep learning: going beyond euclidean data.  IEEE Signal Processing Magazine, 34(4):18–42, 2017.
  * [29] David K Hammond, Pierre Vandergheynst, and Rémi Gribonval.  Wavelets on graphs via spectral graph theory.  Applied and Computational Harmonic Analysis, 30(2):129–150, 2011\.
  * [30] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  arXiv preprint arXiv:1606.09375, 2016.
  * [31] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907, 2016.
  * [32] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  IEEE transactions on neural networks and learning systems, 2020\.
  * [33] Ziwei Zhang, Peng Cui, and Wenwu Zhu.  Deep learning on graphs: A survey.  IEEE Transactions on Knowledge and Data Engineering, 2020.
  * [34] Jie Zhou, Ganqu Cui, Shengding Hu, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph neural networks: A review of methods and applications.  AI Open, 1:57–81, 2020.
  * [35] Nathaniel Trask, Ravi G Patel, Ben J Gross, and Paul J Atzberger.  Gmls-nets: A framework for learning from unstructured data.  arXiv preprint arXiv:1909.05371, 2019.
  * [36] Nathaniel Trask, Andy Huang, and Xiaozhe Hu.  Enforcing exact physics in scientific machine learning: a data-driven exterior calculus on graphs.  arXiv preprint arXiv:2012.11799, 2020.
  * [37] Henry A Rowley, Shumeet Baluja, and Takeo Kanade.  Rotation invariant neural network-based face detection.  In Proceedings. 1998 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No. 98CB36231), pages 38–44. IEEE, 1998.
  * [38] Sander Dieleman, Jeffrey De Fauw, and Koray Kavukcuoglu.  Exploiting cyclic symmetry in convolutional neural networks.  In International conference on machine learning, pages 1889–1898. PMLR, 2016.
  * [39] Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In International conference on machine learning, pages 2990–2999. PMLR, 2016.
  * [40] Daniel E Worrall, Stephan J Garbin, Daniyar Turmukhambetov, and Gabriel J Brostow.  Harmonic networks: Deep translation and rotation equivariance.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5028–5037, 2017.
  * [41] Benjamin Chidester, Minh N Do, and Jian Ma.  Rotation equivariance and invariance in convolutional neural networks.  arXiv preprint arXiv:1805.12301, 2018.
  * [42] Risi Kondor and Shubhendu Trivedi.  On the generalization of equivariance and convolution in neural networks to the action of compact groups.  In International Conference on Machine Learning, pages 2747–2755. PMLR, 2018.
  * [43] Marc Finzi, Samuel Stanton, Pavel Izmailov, and Andrew Gordon Wilson.  Generalizing convolutional neural networks for equivariance to lie groups on arbitrary continuous data.  In International Conference on Machine Learning, pages 3165–3176. PMLR, 2020.
  * [44] Aritra Chowdhury, Elizabeth Kautz, Bülent Yener, and Daniel Lewis.  Image driven machine learning methods for microstructure recognition.  Computational Materials Science, 123:176–187, 2016.
  * [45] Nicholas Lubbers, Turab Lookman, and Kipton Barros.  Inferring low-dimensional microstructure representations using convolutional neural networks.  Physical Review E, 96(5):052111, 2017.
  * [46] Brian L DeCost, Toby Francis, and Elizabeth A Holm.  Exploring the microstructure manifold: image texture representations applied to ultrahigh carbon steel microstructures.  Acta Materialia, 133:30–40, 2017.
  * [47] Ruho Kondo, Shunsuke Yamakawa, Yumi Masuoka, Shin Tajima, and Ryoji Asahi.  Microstructure recognition using convolutional neural networks for prediction of ionic conductivity in ceramics.  Acta Materialia, 141:29–38, 2017.
  * [48] Shi Xingjian, Zhourong Chen, Hao Wang, Dit-Yan Yeung, Wai-Kin Wong, and Wang-Chun Woo.  Convolutional LSTM network: A machine learning approach for precipitation nowcasting.  In Advances in neural information processing systems, pages 802–810, 2015.
  * [49] Xingjian Shi, Zhourong Chen, Hao Wang, Dit Yan Yeung, Wai Kin Wong, and Wangchun Woo.  Convolutional lstm network: A machine learning approach for precipitation nowcasting.  Advances in Neural Information Processing Systems, 2015.
  * [50] Geng Chen, Yoonmi Hong, Yongqin Zhang, Jaeil Kim, Khoi Minh Huynh, Jiquan Ma, Weili Lin, Dinggang Shen, Pew-Thian Yap, UNC/UMN Baby Connectome Project Consortium, et al.  Estimating tissue microstructure with undersampled diffusion data via graph convolutional neural networks.  In International Conference on Medical Image Computing and Computer-Assisted Intervention, pages 280–290. Springer, 2020.
  * [51] Miguel A Bessa, R Bostanabad, Zeliang Liu, A Hu, Daniel W Apley, C Brinson, Wei Chen, and Wing Kam Liu.  A framework for data-driven analysis of materials under uncertainty: Countering the curse of dimensionality.  Computer Methods in Applied Mechanics and Engineering, 320:633–667, 2017.
  * [52] M Mozaffar, R Bostanabad, W Chen, K Ehmann, Jian Cao, and MA Bessa.  Deep learning predicts path-dependent plasticity.  Proceedings of the National Academy of Sciences, 116(52):26414–26420, 2019.
  * [53] Reese Jones, Jeremy A Templeton, Clay M Sanders, and Jakob T Ostien.  Machine learning models of plastic flow based on representation theory.  Computer Modeling in Engineering & Sciences, pages 309–342, 2018\.
  * [54] Toshio Mura.  Micromechanics of defects in solids.  Springer Science & Business Media, 2013.
  * [55] Siavouche Nemat-Nasser and Muneo Hori.  Micromechanics: overall properties of heterogeneous materials.  Elsevier, 2013.
  * [56] Patrizia Trovalusci, Danilo Capecchi, and Giuseppe Ruta.  Genesis of the multiscale approach for materials with microstructure.  Archive of Applied Mechanics, 79(11):981–997, 2009.
  * [57] Chau Le, Tyler E Bruns, and Daniel A Tortorelli.  Material microstructure optimization for linear elastodynamic energy wave management.  Journal of the Mechanics and Physics of Solids, 60(2):351–378, 2012\.
  * [58] Jérémie Bouquerel, Kim Verbeken, and BC De Cooman.  Microstructure-based model for the static mechanical behaviour of multiphase steels.  Acta Materialia, 54(6):1443–1456, 2006.
  * [59] Charles Roduit, Serguei Sekatski, Giovanni Dietler, Stefan Catsicas, Frank Lafont, and Sandor Kasas.  Stiffness tomography by atomic force microscopy.  Biophysical journal, 97(2):674–677, 2009.
  * [60] Nathan M Heckman, Thomas A Ivanoff, Ashley M Roach, Bradley H Jared, Daniel J Tung, Harlan J Brown-Shaklee, Todd Huber, David J Saiz, Josh R Koepke, Jeffrey M Rodelas, et al.  Automated high-throughput tensile testing reveals stochastic process parameter sensitivity.  Materials Science and Engineering: A, 772:138632, 2020.
  * [61] M Čekada, P Panjan, D Kek-Merl, M Panjan, and G Kapun.  Sem study of defects in pvd hard coatings.  Vacuum, 82(2):252–256, 2007.
  * [62] Miroslav Silhavy.  The mechanics and thermodynamics of continuous media.  Springer Science & Business Media, 2013.
  * [63] Michael A Groeber and Michael A Jackson.  DREAM.3D.  http://dream3d.bluequartz.net, 2019.
  * [64] Andrew G Salinger, Roscoe A Bartlett, Andrew M Bradley, Qiushi Chen, Irina P Demeshko, Xujiao Gao, Glen A Hansen, Alejandro Mota, Richard P Muller, Erik Nielsen, et al.  Albany multiphysics code.  https://github.com/gahansen/Albany, 2019.
  * [65] James R Stewart, H Carter Edwards, and et al.  Sierra mechanics.  https://www.sandia.gov/ASC/integrated_codes.html, 2020.
  * [66] Francesco Rizzi, Mohammad Khalil, Reese E Jones, Jeremy A Templeton, Jakob T Ostien, and Brad L Boyce.  Bayesian modeling of inconsistent plastic response due to material variability.  arXiv preprint arXiv:1809.01009, 2018.
  * [67] Jacob Lubliner.  Plasticity theory.  Courier Corporation, 2008.
  * [68] JA Brown, JD Carroll, B Huddleston, Z Casias, and KN Long.  A multiscale study of damage in elastomeric syntactic foams.  Journal of materials science, 53(14):10479–10498, 2018.
  * [69] Geoffrey Ingram Taylor.  The mechanism of plastic deformation of crystals. Part I. Theoretical.  Proceedings of the Royal Society of London. Series A, 145(855):362–387, 1934.
  * [70] E Kroner.  On the plastic deformation of polycrystals.  Acta Metallurgica, 9(2):155–161, 1961.
  * [71] JFW Bishop and Rodney Hill.  XLVI. A theory of the plastic distortion of a polycrystalline aggregate under combined stresses.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42(327):414–427, 1951.
  * [72] JFW Bishop and Rodney Hill.  CXXVIII. A theoretical derivation of the plastic properties of a polycrystalline face-centred metal.  The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 42(334):1298–1307, 1951.
  * [73] Jean Mandel.  Généralisation de la théorie de plasticité de WT Koiter.  International Journal of Solids and structures, 1(3):273–295, 1965\.
  * [74] Paul R Dawson.  Computational crystal plasticity.  International journal of solids and structures, 37(1-2):115–130, 2000.
  * [75] Franz Roters, Philip Eisenlohr, Luc Hantcherli, Denny Dharmawan Tjahjanto, Thomas R Bieler, and Dierk Raabe.  Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications.  Acta Materialia, 58(4):1152–1211, 2010.
  * [76] UF Kocks.  Laws for work-hardening and low-temperature creep.  Journal of engineering materials and technology, 98(1):76–85, 1976\.
  * [77] H Mecking, UF Kocks, and H Fischer.  Hardening, recovery, and creep in fcc mono-and polycrystals.  In Presented at the 4th Intern. Conf. on Strength of Metals and Alloys, Nancy, 30 Aug.-3 Sep. 1976, 1976.
  * [78] Ari L Frankel, Reese E Jones, Coleman Alleman, and Jeremy A Templeton.  Predicting the mechanical response of oligocrystals with deep learning.  Computational Materials Science, 169:109099, 2019.
  * [79] Ari Frankel, Kousuke Tachida, and Reese Jones.  Prediction of the evolution of the stress field of polycrystals undergoing elastic-plastic deformation with a hybrid neural network model.  Machine Learning: Science and Technology, 1(3):035005, 2020.
  * [80] Johan Bjorck, Carla Gomes, Bart Selman, and Kilian Q Weinberger.  Understanding batch normalization.  arXiv preprint arXiv:1806.02375, 2018.
  * [81] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [82] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio.  Learning phrase representations using rnn encoder-decoder for statistical machine translation.  arXiv preprint arXiv:1406.1078, 2014.
  * [83] Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pages 249–256. JMLR Workshop and Conference Proceedings, 2010.
  * [84] Xavier Glorot, Antoine Bordes, and Yoshua Bengio.  Deep sparse rectifier neural networks.  In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 315–323. JMLR Workshop and Conference Proceedings, 2011.
  * [85] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: An end-to-end open source machine learning platform.  https://www.tensorflow.org/, 2021.
  * [86] Daniele Grattarola.  Spektral: a python library for graph deep learning.  https://graphneural.network/, 2021.
  * [87] Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Multipole graph neural operator for parametric partial differential equations.  arXiv preprint arXiv:2006.09535, 2020.
  * [88] Facundo Quiroga, Franco Ronchetti, Laura Lanzarini, and Aurelio F Bariviera.  Revisiting data augmentation for rotational invariance in convolutional neural networks.  In International Conference on Modelling and Simulation in Management Sciences, pages 127–141. Springer, 2018.
  * [89] Berkant Savas and Inderjit S Dhillon.  Clustered low rank approximation of graphs in information science applications.  In Proceedings of the 2011 SIAM International Conference on Data Mining, pages 164–175. SIAM, 2011.
  * [90] Emile Richard, Pierre-André Savalle, and Nicolas Vayatis.  Estimation of simultaneously sparse and low rank matrices.  arXiv preprint arXiv:1206.6474, 2012.
  * [91] Vadim Lebedev, Yaroslav Ganin, Maksim Rakhuba, Ivan Oseledets, and Victor Lempitsky.  Speeding-up convolutional neural networks using fine-tuned cp-decomposition.  arXiv preprint arXiv:1412.6553, 2014.
  * [92] Cheng Tai, Tong Xiao, Yi Zhang, Xiaogang Wang, et al.  Convolutional neural networks with low-rank regularization.  arXiv preprint arXiv:1511.06067, 2015.
  * [93] Taiju Kanada, Masaki Onuki, and Yuichi Tanaka.  Low-rank sparse decomposition of graph adjacency matrices for extracting clean clusters.  In 2018 Asia-Pacific Signal and Information Processing Association Annual Summit and Conference (APSIPA ASC), pages 1153–1159. IEEE, 2018.
  * [94] Mikhail Belkin and Partha Niyogi.  Laplacian eigenmaps for dimensionality reduction and data representation.  Neural computation, 15(6):1373–1396, 2003.
  * [95] Xiaofei He and Partha Niyogi.  Locality preserving projections.  In Advances in neural information processing systems, pages 153–160, 2004.
  * [96] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  arXiv preprint arXiv:1611.07308, 2016.
  * [97] Yiyi Liao, Yue Wang, and Yong Liu.  Graph regularized auto-encoders for image representation.  IEEE Transactions on Image Processing, 26(6):2839–2852, 2016\.
  * [98] Arman Hasanzadeh, Ehsan Hajiramezanali, Nick Duffield, Krishna R Narayanan, Mingyuan Zhou, and Xiaoning Qian.  Semi-implicit graph variational auto-encoders.  arXiv preprint arXiv:1908.07078, 2019.
  * [99] Amin Salehi and Hasan Davulcu.  Graph attention auto-encoders.  arXiv preprint arXiv:1905.10715, 2019.
  * [100] Ting Zhang, Bang Liu, Di Niu, Kunfeng Lai, and Yu Xu.  Multiresolution graph attention networks for relevance matching.  In Proceedings of the 27th ACM International Conference on Information and Knowledge Management, pages 933–942, 2018.

##  Appendix A The Graph Convolutional Network

As mentioned in the Introduction, the Kipf and Welling Graph Convolutional
Network (GCN) [31] provides an innovative, expressive graph convolutional
network built in sequentially applied layers with local action. Here we give a
brief synopsis of their development.

The objective is to efficiently apply a graph filter
$g_{\boldsymbol{\theta}}={\operatorname{diag}}({\boldsymbol{\theta}})$ where
the parameter vector ${\boldsymbol{\theta}}$ has $N_{\text{nodes}}$ entries.
Convolution of the filter $g_{\boldsymbol{\theta}}$ and data $\mathsf{x}$ on a
graph can be expressed as [29]

$g_{\boldsymbol{\theta}}\ast\mathsf{x}=\mathsf{U}g_{\boldsymbol{\theta}}\
\mathsf{U}^{T}\mathsf{x}$ (20)

analogous to the classical convolution theorem. This formulation is connected
to the (normalized) graph Laplacian $\mathsf{L}$ and its spectral
representation

$\mathsf{L}=\mathsf{I}-\mathsf{D}^{-1/2}\mathsf{A}\mathsf{D}^{-1/2}=\mathsf{U}{\boldsymbol{\Lambda}}\mathsf{U}^{T}$
(21)

where $\mathsf{A}$ is the binary adjacency matrix, $\mathsf{D}$ is the
associated degree matrix, $\mathsf{U}$ is the matrix of eigenvectors and
${\boldsymbol{\Lambda}}$ is the diagonal matrix of eigenvalues. The
formulation for graph convolution in Eq. (20), in turn can be approximated by
an expansion of Chebyshev polynomials $T_{k}$ [30]

$g_{\boldsymbol{\theta}}\ast\mathsf{x}=\left(\mathsf{U}g_{\boldsymbol{\theta}}\mathsf{U}^{T}\right)\mathsf{x}\approx\sum_{k}\vartheta_{k=0}^{K}T_{k}(\tilde{\mathsf{L}})\,\mathsf{x}$
(22)

where $\tilde{\mathsf{L}}=\frac{2}{\lambda_{\text{max}}}\mathsf{L}-\mathsf{I}$
and $\lambda_{\text{max}}$ is the maximum eigenvalue of $\mathsf{L}$.

To this approximation Kipf and Welling make a number of additional
simplifications. First they approximate the maximum eigenvalue
$\lambda_{\text{max}}\approx 2$ so that
$\tilde{\mathsf{L}}=\mathsf{L}-\mathsf{I}$ i.e. $\tilde{\mathsf{L}}$ is the
graph Laplacian with added self-loops/interactions. Next, they truncate the
expansion in Eq. (22) at $K=1$ so that

$g_{\boldsymbol{\theta}}\ast\mathsf{x}\approx\vartheta_{0}\mathsf{I}\mathsf{x}-\vartheta_{1}\mathsf{D}^{-1/2}\mathsf{A}\mathsf{D}^{-1/2}\mathsf{x}$
(23)

This effectively reduces the number of free parameters in the filter from
$N_{\text{nodes}}$ to 2. This has the tremendous advantage of cheap and local
action. The expressiveness of a network built on these layers is controllable
by the GCNN depth (number of layers). Lastly they further collapse the number
of free parameters from 2 to 1 by setting $\vartheta_{1}=-\vartheta_{0}$