Source code for clustering.citoperations

"""Green operator and cluster interaction tensors related operations.

This module includes several procedures required to compute the cluster
interaction tensors, namely (1) the frequency discretization of the spatial
domain and (2) the computation of the Green operator, as well as the assembly
of the global cluster interaction matrix.

Functions
---------
set_discrete_freqs
    Perform frequency discretization of the spatial domain.
gop_material_independent_terms
    Compute Green operator material independent terms in frequency domain.
assemble_cit
    Assemble global cluster interaction matrix.
"""
#
#                                                                       Modules
# =============================================================================
# Standard
import copy
import itertools as it
# Third-party
import numpy as np
import numpy.matlib
# Local
import tensor.matrixoperations as mop
#
#                                                          Authorship & Credits
# =============================================================================
__author__ = 'Bernardo Ferreira (bernardo_ferreira@brown.edu)'
__credits__ = ['Bernardo Ferreira', ]
__status__ = 'Stable'
# =============================================================================
#
# =============================================================================
#
#                                                          Discrete frequencies
# =============================================================================
def set_discrete_freqs(n_dim, rve_dims, n_voxels_dims):
    """Perform frequency discretization of the spatial domain.

    Perform frequency discretization by setting the spatial discrete
    frequencies (rad/m) for each dimension.

    *2D case*:

    .. math::

       \\boldsymbol{\\zeta}_{s_{1}, s_{2}} =
          \\left( \\dfrac{2\\pi}{(l_{\\mathrm{RVE}})_{1}} s_{1}, \\,
          \\dfrac{2\\pi}{(l_{\\mathrm{RVE}})_{2}} s_{2} \\right) \\, ,

    .. math::

       s_{i}=0, 1, \\dots, n_{i}-1 \\, , \\quad i=1,2 \\, .

    where :math:`\\boldsymbol{\\zeta}_{s_{1}, s_{2}} \\equiv
    \\boldsymbol{\\zeta}(s_{1}, s_{2})` denotes a sampling angular frequency,
    :math:`(l_{\\mathrm{RVE}})_{i}` is the RVE size in the :math:`i` th
    dimension, and :math:`n_{i}` is the number of voxels in the :math:`i` th
    dimension.

    ----

    *3D case*:

    .. math::

       \\boldsymbol{\\zeta}_{s_{1}, s_{2}, s_{3}} =
          \\left( \\dfrac{2\\pi}{(l_{\\mathrm{RVE}})_{1}} s_{1}, \\,
          \\dfrac{2\\pi}{(l_{\\mathrm{RVE}})_{2}} s_{2}, \\,
          \\dfrac{2\\pi}{(l_{\\mathrm{RVE}})_{3}} s_{3} \\right) \\, ,

    .. math::

       s_{i}=0, 1, \\dots, n_{i}-1 \\, , \\quad i=1,2,3 \\, .

    where :math:`\\boldsymbol{\\zeta}_{s_{1}, s_{2}, s_{3}} \\equiv
    \\boldsymbol{\\zeta}(s_{1}, s_{2}, s_{3})` denotes a sampling angular
    frequency, :math:`(l_{\\mathrm{RVE}})_{i}` is the RVE size in the
    :math:`i` th dimension, and :math:`n_{i}` is the number of voxels in the
    :math:`i` th dimension.

    ----

    Parameters
    ----------
    n_dim : int
        Problem number of spatial dimensions.
    rve_dims : list
        RVE size in each dimension.
    n_voxels_dims : list[int]
        Number of voxels in each dimension of the regular grid (spatial
        discretization of the RVE).

    Returns
    -------
    freq_dims : list[numpy.ndarray (1d)]
        List containing the sample frequencies (numpy.ndarray (1d) of shape
        (n_voxels_dim,)) for each dimension.
    """
    freqs_dims = list()
    for i in range(n_dim):
        # Set sampling spatial period
        sampling_period = rve_dims[i]/n_voxels_dims[i]
        # Set discrete frequencies
        freqs_dims.append(2*np.pi*np.fft.fftfreq(n_voxels_dims[i],
                                                 sampling_period))
    # Return
    return freqs_dims
#
#                                                                Green operator
# =============================================================================
def gop_material_independent_terms(strain_formulation, problem_type, rve_dims,
                                   n_voxels_dims):
    """Compute Green operator material independent terms in frequency domain.

    The Green operator is here assumed to be associated with an isotropic
    elastic reference material.

    The material independent terms of the Green operator are conveniently
    computed to perform an efficient update of the Green operator if the
    associated reference material elastic properties are updated by any means
    (e.g., self-consistent scheme).

    *Infinitesimal strains*:

    .. math::

       (\\breve{\\Phi}^{0}_{1})_{klij} (\\boldsymbol{\\zeta}) =
       \\dfrac{ \\delta_{ki} \\zeta_{j} \\zeta_{l}
              + \\delta_{kj} \\zeta_{i} \\zeta_{l}
              + \\delta_{li} \\zeta_{j} \\zeta_{k}
              + \\delta_{lj} \\zeta_{i} \\zeta_{k}}
              {||\\boldsymbol{\\zeta}||^{2}} \\, ,

    .. math::

        (\\breve{\\Phi}^{0}_{2})_{klij} (\\boldsymbol{\\zeta}) =
        - \\dfrac{\\zeta_{k}\\zeta_{l}\\zeta_{i}\\zeta_{j}}
                 {||\\boldsymbol{\\zeta}||^{4}} \\, ,

    .. math::

       i,j,k,l =1, \\dots, n_{\\text{dim}} \\, ,

    where :math:`(\\breve{\\Phi}^{0}_{1})_{klij}` and
    :math:`(\\breve{\\Phi}^{0}_{2})_{klij}` are the first and second
    fourth-order Green operator material independent terms, respectively,
    :math:`\\boldsymbol{\\zeta}` is the frequency wave vector,
    :math:`\\delta_{ij}` is the Kronecker delta, and :math:`n_{\\text{dim}}`
    is the number of spatial dimensions.

    ----

    *Finite strains*:

    .. math::

       (\\breve{\\Phi}^{0}_{1})_{klij} (\\boldsymbol{\\zeta}) =
       \\dfrac{ \\delta_{ki} \\zeta_{j} \\zeta_{l}}
              {||\\boldsymbol{\\zeta}||^{2}} \\, , \\qquad
        (\\breve{\\Phi}^{0}_{2})_{klij} (\\boldsymbol{\\zeta}) =
        - \\dfrac{\\zeta_{k}\\zeta_{l}\\zeta_{i}\\zeta_{j}}
                 {||\\boldsymbol{\\zeta}||^{4}} \\, ,

    .. math::

       i,j,k,l =1, \\dots, n_{\\text{dim}} \\, .

    where :math:`(\\breve{\\Phi}^{0}_{1})_{klij}` and
    :math:`(\\breve{\\Phi}^{0}_{2})_{klij}` are the first and second
    fourth-order Green operator material independent terms, respectively,
    :math:`\\boldsymbol{\\zeta}` is the frequency wave vector,
    :math:`\\delta_{ij}` is the Kronecker delta, and :math:`n_{\\text{dim}}`
    is the number of spatial dimensions.

    A detailed description of the computational implementation based on
    Hadamard (element-wise) operations can be found in Section 4.6 of
    Ferreira (2022) [#]_.

    .. [#] Ferreira, B.P. (2022). *Towards Data-driven Multi-scale
           Optimization of Thermoplastic Blends: Microstructural
           Generation, Constitutive Development and Clustering-based
           Reduced-Order Modeling.* PhD Thesis, University of Porto
           (see `here <https://repositorio-aberto.up.pt/handle/10216/
           146900?locale=en>`_)

    ----

    Parameters
    ----------
    strain_formulation: {'infinitesimal', 'finite'}
        Problem strain formulation.
    problem_type : int
        Problem type: 2D plane strain (1), 2D plane stress (2),
        2D axisymmetric (3) and 3D (4).
    rve_dims : list[float]
        RVE size in each dimension.
    n_voxels_dims : list[int]
        Number of voxels in each dimension of the regular grid (spatial
        discretization of the RVE).

    Returns
    -------
    gop_1_dft_vox : dict
        Regular grid shaped matrix (item, numpy.ndarray) containing each
        fourth-order matricial form component (key, str) of the first Green
        operator material independent term in the frequency domain (discrete
        Fourier transform).
    gop_2_dft_vox : dict
        Regular grid shaped matrix (item, numpy.ndarray) containing each
        fourth-order matricial form component (key, str) of the second Green
        operator material independent term in the frequency domain (discrete
        Fourier transform).
    gop_0_freq_dft_vox : dict
        Regular grid shaped matrix (item, numpy.ndarray) containing each
        fourth-order matricial form component (key, str) of the Green operator
        zero-frequency (material independent) term in the frequency domain
        (discrete Fourier transform).
    """
    # Get problem type parameters
    n_dim, comp_order_sym, comp_order_nsym = \
        mop.get_problem_type_parameters(problem_type)
    # Set strain/stress components order according to problem strain
    # formulation
    if strain_formulation == 'infinitesimal':
        comp_order = comp_order_sym
    elif strain_formulation == 'finite':
        comp_order = comp_order_nsym
    else:
        raise RuntimeError('Unknown problem strain formulation.')
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Set discrete frequencies (rad/m) for each dimension
    freqs_dims = set_discrete_freqs(n_dim, rve_dims, n_voxels_dims)
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Set Green operator matricial form components
    comps = list(it.product(comp_order, comp_order))
    # Set mapping between Green operator fourth-order tensor and matricial form
    # components
    fo_indexes = list()
    mf_indexes = list()
    for i in range(len(comp_order)**2):
        fo_indexes.append([int(x) - 1
                           for x in list(comps[i][0] + comps[i][1])])
        mf_indexes.append([x for x in [comp_order.index(comps[i][0]),
                                       comp_order.index(comps[i][1])]])
    # Set optimized variables
    var1 = [*np.meshgrid(*freqs_dims, indexing='ij')]
    var2 = dict()
    for fo_idx in fo_indexes:
        if str(fo_idx[1]) + str(fo_idx[3]) not in var2.keys():
            var2[str(fo_idx[1]) + str(fo_idx[3])] = \
                np.multiply(var1[fo_idx[1]], var1[fo_idx[3]])
        if ''.join([str(x) for x in fo_idx]) not in var2.keys():
            var2[''.join([str(x) for x in fo_idx])] = \
                np.multiply(np.multiply(var1[fo_idx[0]], var1[fo_idx[1]]),
                            np.multiply(var1[fo_idx[2]], var1[fo_idx[3]]))
        if strain_formulation == 'infinitesimal':
            if str(fo_idx[1]) + str(fo_idx[2]) not in var2.keys():
                var2[str(fo_idx[1]) + str(fo_idx[2])] = \
                    np.multiply(var1[fo_idx[1]], var1[fo_idx[2]])
            if str(fo_idx[0]) + str(fo_idx[2]) not in var2.keys():
                var2[str(fo_idx[0]) + str(fo_idx[2])] = \
                    np.multiply(var1[fo_idx[0]], var1[fo_idx[2]])
            if str(fo_idx[0]) + str(fo_idx[3]) not in var2.keys():
                var2[str(fo_idx[0]) + str(fo_idx[3])] = \
                    np.multiply(var1[fo_idx[0]], var1[fo_idx[3]])
    if n_dim == 2:
        var3 = np.sqrt(np.add(np.square(var1[0]), np.square(var1[1])))
    else:
        var3 = np.sqrt(np.add(np.add(np.square(var1[0]), np.square(var1[1])),
                              np.square(var1[2])))
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    # Initialize Green operator material independent terms
    gop_1_dft_vox = {''.join([str(x+1) for x in idx]):
                     np.zeros(tuple(n_voxels_dims)) for idx in fo_indexes}
    gop_2_dft_vox = {''.join([str(x+1) for x in idx]):
                     np.zeros(tuple(n_voxels_dims)) for idx in fo_indexes}
    gop_0_freq_dft_vox = {''.join([str(x+1) for x in idx]):
                          np.zeros(tuple(n_voxels_dims)) for idx in fo_indexes}
    # Compute Green operator matricial form components
    for i in range(len(mf_indexes)):
        # Get fourth-order tensor indexes
        fo_idx = fo_indexes[i]
        # Get Green operator component
        comp = ''.join([str(x+1) for x in fo_idx])
        # Set optimized variables
        var4 = [fo_idx[0] == fo_idx[2], fo_idx[0] == fo_idx[3],
                fo_idx[1] == fo_idx[3], fo_idx[1] == fo_idx[2]]
        var5 = [str(fo_idx[1]) + str(fo_idx[3]),
                str(fo_idx[1]) + str(fo_idx[2]),
                str(fo_idx[0]) + str(fo_idx[2]),
                str(fo_idx[0]) + str(fo_idx[3])]
        var4 = [fo_idx[0] == fo_idx[2], ]
        var5 = [str(fo_idx[1]) + str(fo_idx[3]), ]
        if strain_formulation == 'infinitesimal':
            var4 += [fo_idx[0] == fo_idx[3], fo_idx[1] == fo_idx[3],
                     fo_idx[1] == fo_idx[2]]
            var5 += [str(fo_idx[1]) + str(fo_idx[2]),
                     str(fo_idx[0]) + str(fo_idx[2]),
                     str(fo_idx[0]) + str(fo_idx[3])]
        # Compute first material independent term of Green operator
        first_term = np.zeros(tuple(n_voxels_dims))
        for j in range(len(var4)):
            if var4[j]:
                first_term = np.add(first_term, var2[var5[j]])
        first_term = np.divide(first_term, np.square(var3),
                               where=abs(var3) > 1e-10)
        gop_1_dft_vox[comp] = copy.copy(first_term)
        # Compute second material independent term of Green operator
        gop_2_dft_vox[comp] = -1.0*np.divide(
            var2[''.join([str(x) for x in fo_idx])],
            np.square(np.square(var3)), where=abs(var3) > 1e-10)
        # Compute Green operator zero-frequency term
        gop_0_freq_dft_vox[comp][tuple(n_dim*(0,))] = 1.0
    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    return gop_1_dft_vox, gop_2_dft_vox, gop_0_freq_dft_vox
#
#                                    Global cluster interaction matrix assembly
# =============================================================================
[docs]def assemble_cit(strain_formulation, problem_type, mat_prop_ref, material_phases, phase_n_clusters, phase_clusters, cit_1_mf, cit_2_mf, cit_0_freq_mf): """Assemble global cluster interaction matrix. Update the cluster interaction tensors by taking into account the material properties of the reference material and assemble them in the global cluster interaction matrix. The dependency on the material properties of the reference material stems from the definition of the Green operator as shown below. *Infinitesimal strains*: .. math:: \\breve{\\mathbf{\\Phi}}^{0} (\\boldsymbol{\\zeta}) = c_{1}(\\lambda^{0}, \\mu^{0}) \\, \\breve{\\mathbf{\\Phi}}^{0}_{1} + c_{2}(\\lambda^{0}, \\mu^{0}) \\, \\breve{\\mathbf{\\Phi}}^{0}_{2} \\, , .. math:: c_{1}(\\lambda^{0}, \\mu^{0}) = \\dfrac{1}{4 \\mu^{0}} \\, , \\qquad c_{2}(\\lambda^{0}, \\mu^{0}) = \\dfrac{\\lambda^{0} + \\mu^{0}}{\\mu^{0}(\\lambda^{0} + 2 \\mu^{0})} where :math:`\\breve{\\mathbf{\\Phi}}^{0}_{1}` and :math:`\\breve{\\mathbf{\\Phi}}^{0}_{2}` are the first and second fourth-order Green operator material independent terms, respectively, and :math:`(\\lambda^{0}, \\mu^{0})` are the reference elastic (isotropic) material Lamé parameters. ---- *Finite strains*: .. math:: \\breve{\\mathbf{\\Phi}}^{0} (\\boldsymbol{\\zeta}) = c_{1}(\\lambda^{0}, \\mu^{0}) \\, \\breve{\\mathbf{\\Phi}}^{0}_{1} + c_{2}(\\lambda^{0}, \\mu^{0}) \\, \\breve{\\mathbf{\\Phi}}^{0}_{2} \\, , .. math:: c_{1}(\\lambda^{0}, \\mu^{0}) = \\dfrac{1}{2 \\mu^{0}} \\, , \\qquad c_{2}(\\lambda^{0}, \\mu^{0}) = \\dfrac{\\lambda^{0}}{2 \\mu^{0}(\\lambda^{0} + 2 \\mu^{0})} where :math:`\\breve{\\mathbf{\\Phi}}^{0}_{1}` and :math:`\\breve{\\mathbf{\\Phi}}^{0}_{2}` are the first and second fourth-order Green operator material independent terms, respectively, and :math:`(\\lambda^{0}, \\mu^{0})` are the reference elastic (isotropic) material Lamé parameters. ---- Parameters ---------- strain_formulation: {'infinitesimal', 'finite'} Problem strain formulation. problem_type : int Problem type: 2D plane strain (1), 2D plane stress (2), 2D axisymmetric (3) and 3D (4). mat_prop_ref : dict Reference material properties. material_phases : list[str] RVE material phases labels (str). phase_n_clusters : dict Number of clusters (item, int) prescribed for each material phase (key, str). phase_clusters : dict Clusters labels (item, list[int]) associated to each material phase (key, str). cit_1_mf : dict Cluster interaction tensors associated with the first Green operator material independent term. Each tensor is stored in a dictionary (item, dict) for each pair of material phases (key, str), which in turn contains the corresponding matricial form (item, numpy.ndarray) associated to each pair of clusters (key, str). cit_2_mf : dict Cluster interaction tensors associated with the second Green operator material independent term. Each tensor is stored in a dictionary (item, dict) for each pair of material phases (key, str), which in turn contains the corresponding matricial form (item, numpy.ndarray) associated to each pair of clusters (key, str). cit_0_freq_mf : dict Cluster interaction tensors associated with the zero-frequency Green operator (material independent) term. Each tensor is stored in a dictionary (item, dict) for each pair of material phases (key, str), which in turn contains the corresponding matricial form (item, numpy.ndarray) associated to each pair of clusters (key, str). Returns ------- global_cit_mf : numpy.ndarray Global cluster interaction matrix. Assembly positions are assigned according to the order of material_phases (1st) and phase_clusters (2nd). """ # Get problem type parameters _, comp_order_sym, comp_order_nsym = \ mop.get_problem_type_parameters(problem_type) # Set strain/stress components order according to problem strain # formulation if strain_formulation == 'infinitesimal': comp_order = comp_order_sym elif strain_formulation == 'finite': comp_order = comp_order_nsym else: raise RuntimeError('Unknown problem strain formulation.') # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Get total number of clusters n_total_clusters = sum([phase_n_clusters[mat_phase] for mat_phase in material_phases]) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Get reference material Young modulus and Poisson ratio E_ref = mat_prop_ref['E'] v_ref = mat_prop_ref['v'] # Compute reference material Lamé parameters lam_ref = (E_ref*v_ref)/((1.0 + v_ref)*(1.0 - 2.0*v_ref)) miu_ref = E_ref/(2.0*(1.0 + v_ref)) # Compute Green operator's reference material coefficients if strain_formulation == 'infinitesimal': gop_factor_1 = 1.0/(4.0*miu_ref) gop_factor_2 = (lam_ref + miu_ref)/(miu_ref*(lam_ref + 2.0*miu_ref)) else: gop_factor_1 = 1.0/(2.0*miu_ref) gop_factor_2 = lam_ref/(2.0*miu_ref*(lam_ref + 2.0*miu_ref)) gop_factor_0_freq = numpy.matlib.repmat( np.zeros((len(comp_order), len(comp_order))), n_total_clusters, n_total_clusters) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Initialize global clustering interaction matrices dims = (n_total_clusters*len(comp_order), n_total_clusters*len(comp_order)) global_cit_1_mf = np.zeros(dims) global_cit_2_mf = np.zeros(dims) global_cit_0_freq_mf = np.zeros(dims) # Initialize column cluster index jclst = 0 # Loop over material phases for mat_phase_B in material_phases: # Loop over material phase B clusters for cluster_J in phase_clusters[mat_phase_B]: # Initialize row cluster index iclst = 0 # Loop over material phases for mat_phase_A in material_phases: # Set material phase pair mat_phase_pair = mat_phase_A + '_' + mat_phase_B # Loop over material phase A clusters for cluster_I in phase_clusters[mat_phase_A]: # Set cluster pair cluster_pair = str(cluster_I) + '_' + str(cluster_J) # Set assembling ranges i_init = iclst*len(comp_order) i_end = i_init + len(comp_order) j_init = jclst*len(comp_order) j_end = j_init + len(comp_order) # Assemble cluster interaction tensors global_cit_1_mf[i_init:i_end, j_init:j_end] = \ cit_1_mf[mat_phase_pair][cluster_pair] global_cit_2_mf[i_init:i_end, j_init:j_end] = \ cit_2_mf[mat_phase_pair][cluster_pair] global_cit_0_freq_mf[i_init:i_end, j_init:j_end] = \ cit_0_freq_mf[mat_phase_pair][cluster_pair] # Increment row cluster index iclst = iclst + 1 # Increment column cluster index jclst = jclst + 1 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # Assemble global cluster interaction matrix global_cit_mf = np.add( np.add(np.multiply(gop_factor_1, global_cit_1_mf), np.multiply(gop_factor_2, global_cit_2_mf)), np.multiply(gop_factor_0_freq, global_cit_0_freq_mf)) # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ return global_cit_mf