Periodic centroid

[eisenlohr]

As many DAMASK simulations utilize periodic grain structures, it is helpful to have a robust method to calculate the centroid of each of these grains. A plain averaging scheme that does not take the periodic boundaries into account will give wrong results, as is easily demonstrated for a one-dimensional example where a grain of interest is occupying grid points marked by x. Simple averaging yields C as center of gravity.

xxx.......C.....xxxx

The below Python function approximates the centroid of a point cloud (in arbitrary dimensions) considering a periodic tiling of space by a box of given size and origin.

import numpy as np

def centroid_periodic(points,
                      size=None,
                      origin=None):
    """
    Calculate approximate centroid of point cloud within periodic box of 'size' located at 'origin'.

    Parameters
    ----------
    points : numpy.ndarray, shape (:,d)
        Point cloud.
    size : numpy.ndarray, shape (d), optional
        Size of periodic domain.
        Defaults to unit cube of dimension d.
    origin : numpy.ndarray, shape (d), optional
        Origin of periodic domain.
        Defaults to zero.

    Returns
    -------
    centroid : numpy.ndarray, shape (d)

    Examples
    --------
    Periodic centroid of three data points along one dimension in domain [-1,1]

    >>> centroid_periodic([-1.1,-0.9,0.95],size=2,origin=-1)
    np.float64(0.9827856395454926)
    
    Periodic centroid of four data points around (0,0) in (1,1) periodic domain

    >>> centroid_periodic([[.1,.1],[.1,.9],[.9,.9],[.9,.1]],size=[1,1])
    array([0., 0.])
    
    """
    d = np.atleast_2d(points).shape[-1]
    size_ = np.ones(d) if size is None else np.asarray(size)
    origin_ = np.zeros(d) if origin is None else np.asarray(origin)
    points_ = 2*np.pi*(np.asarray(points)-origin_)/size_
    return origin_+size_*((np.atan2(np.average(np.sin(points_),axis=0),
                                    np.average(np.cos(points_),axis=0))%(2.*np.pi)/2./np.pi)%1.)

[MarDiehl]

for Grid objects, this is implemented as https://damask-multiphysics.org/documentation/processing_tools/pre-processing.html#damask.seeds.from_grid