InertiaComputationsVoxelGrid Namespace Reference

Functions
def	BoundingBox (pebbles)
def	VoxelGridFromPebbles (pebbles, N)
def	ComputeVolCOMVoxelGrid (bbox, span, vox)
def	ShiftToVoxelCOMBBoxPebbles (com, pebbles, bbox)
def	ClumpTOIVoxelGrid (density, com, vox, bbox, span)
def	ComputeInertiaFromVoxelGrid (OPT, DATA)

Function Documentation

BoundingBox()

def InertiaComputationsVoxelGrid.BoundingBox

(

pebbles

)

def BoundingBox(pebbles):
    # returns CUBIC bounding box in the shape [x_min, x_max, y_min, y_max, z_min, z_max]
    # that encloses the clump. This bounding box is used for voxelization. Such an approach
    # is OK if the shapes are not too oblique.
 
    spans = 0.5 * np.array([np.max(pebbles[:, 0] + pebbles[:, 3]) - np.min(pebbles[:, 0] - pebbles[:, 3]),
                            np.max(pebbles[:, 1] + pebbles[:, 3]) - np.min(pebbles[:, 1] - pebbles[:, 3]),
                            np.max(pebbles[:, 2] + pebbles[:, 3]) - np.min(pebbles[:, 2] - pebbles[:, 3])])
    box_center = 0.5 * np.array([np.max(pebbles[:, 0] + pebbles[:, 3]) + np.min(pebbles[:, 0] - pebbles[:, 3]),
                                 np.max(pebbles[:, 1] + pebbles[:, 3]) + np.min(pebbles[:, 1] - pebbles[:, 3]),
                                 np.max(pebbles[:, 2] + pebbles[:, 3]) + np.min(pebbles[:, 2] - pebbles[:, 3])])
    max_span = np.max(spans)
 
    return np.array([box_center[0] - max_span,
                     box_center[0] + max_span,
                     box_center[1] - max_span,
                     box_center[1] + max_span,
                     box_center[2] - max_span,
                     box_center[2] + max_span]), max_span
 
@jit(nopython=True)

Referenced by VoxelGridFromPebbles().

ClumpTOIVoxelGrid()

def InertiaComputationsVoxelGrid.ClumpTOIVoxelGrid	(		density,
			com,
			vox,
			bbox,
			span
	)

def ClumpTOIVoxelGrid(density, com, vox, bbox, span):
    TOI = np.zeros(9).reshape(3,3)
    N = vox.shape[0]
    for i in range(N):
        for j in range(N):
            for k in range(N):
                if (vox[i,j,k]):
                    r = np.array([ bbox[0] + 2.*span *(1./ (2.*N) + i / N),
                               bbox[2] + 2.*span *(1./ (2.*N) + j / N),
                               bbox[4] + 2.*span *(1./ (2.*N) + k / N)]) - com
                    X = r[0]
                    Y = r[1]
                    Z = r[2]
                    TOI += np.array([[Y**2 + Z**2, -X*Y, -X*Z],
                                    [-X*Y, X**2 + Z**2, -Y*Z],
                                    [-X*Z, -Y*Z, X**2 + Y**2]])
    return density * (2*span/N)**3 * TOI
 

Referenced by ComputeInertiaFromVoxelGrid().

ComputeInertiaFromVoxelGrid()

def InertiaComputationsVoxelGrid.ComputeInertiaFromVoxelGrid	(		OPT,
			DATA
	)

def ComputeInertiaFromVoxelGrid(OPT, DATA):
    # take array of pebbles
    pebbles = DATA['pebbles']
    density = DATA['density']
 
    # Compute voxel grid
    N = OPT['voxNum']
    bbox, span, vox = VoxelGridFromPebbles(pebbles, N)
 
    # Compute mass, and center of mass
    vol, com = ComputeVolCOMVoxelGrid(bbox, span, vox)
 
    # Shift pebbles and span to make com an origin
    pebbles, bbox = ShiftToVoxelCOMBBoxPebbles(com, pebbles, bbox)
 
    # Compute tensor of inertia
    toi = ClumpTOIVoxelGrid(density, com, vox, bbox, span)
    if OPT['verbose']: print("Tensor of inertia of voxels: ", toi)
 
    # Compute principal directions
    v1, v2, v3 = ComputePrincipalDirections(toi)
    if OPT['verbose']: print("Principal directions (voxel grid): ", v1, v2, v3)
 
    # Rotate to principal directions
    if OPT['rotateToPD']:
        pebbles, toi, v1, v2, v3 = RotateToPDPebblesTOI(pebbles, toi, v1, v2, v3)
        if OPT['verbose']: print("Rotated principal directions: ", v1, v2, v3)
        if OPT['verbose']:
            print("Rotated toi: ")
            print(toi)
 
    # Save all the data
    DATA['pebbles'] = pebbles
    DATA['mass'] = vol * DATA['density']
    DATA['pd'] = v1, v2, v3
    DATA['toi'] = toi
    DATA['voxelGrid'] = vox, bbox, span
 
 
    return OPT, DATA

References ClumpTOIVoxelGrid(), InertiaComputationsPebbles.ComputePrincipalDirections(), ComputeVolCOMVoxelGrid(), Eigen::internal.print(), InertiaComputationsPebbles.RotateToPDPebblesTOI(), ShiftToVoxelCOMBBoxPebbles(), and VoxelGridFromPebbles().

Referenced by MClump.main().

ComputeVolCOMVoxelGrid()

def InertiaComputationsVoxelGrid.ComputeVolCOMVoxelGrid	(		bbox,
			span,
			vox
	)

def ComputeVolCOMVoxelGrid(bbox, span, vox):
    # This function defines the absolute position of center of gravity of all voxels
    N = vox.shape[0]
    vol = np.sum(vox) * (2. * span / N)**3
 
    com = np.array([0.,0.,0.])
    for i in range(N):
        for j in range(N):
            for k in range(N):
                com += vox[i,j,k] * np.array([ bbox[0] + 2.*span *(1./ (2.*N) + i / N),
                                               bbox[2] + 2.*span *(1./ (2.*N) + j / N),
                                               bbox[4] + 2.*span *(1./ (2.*N) + k / N)])
    return vol, com / np.sum(vox)
 
@jit(nopython=True)

Referenced by ComputeInertiaFromVoxelGrid().

ShiftToVoxelCOMBBoxPebbles()

def InertiaComputationsVoxelGrid.ShiftToVoxelCOMBBoxPebbles	(		com,
			pebbles,
			bbox
	)

def ShiftToVoxelCOMBBoxPebbles(com, pebbles, bbox):
    # Returns the coordinates of pebbles shifted such that the center of mass (COM) is at zero.
    for j in range(len(pebbles)):
        pebbles[j][:3] -= com
    for i in range(3):
        bbox[2 * i] -= com[i]
        bbox[2 * i + 1] -= com[i]
 
    return pebbles, bbox
 
@jit(nopython=True)

Referenced by ComputeInertiaFromVoxelGrid().

VoxelGridFromPebbles()

def InertiaComputationsVoxelGrid.VoxelGridFromPebbles	(		pebbles,
			N
	)

def VoxelGridFromPebbles(pebbles, N):
    # this function returns decomposition of cubic domain with void and material split in binary voxels
 
    vox = np.zeros(N ** 3).reshape(N, N, N)
    bbox, span = BoundingBox(pebbles)
 
    for i in range(N):
        for j in range(N):
            for k in range(N):
                for pebble in range(len(pebbles)):
                    pos_v = np.array([bbox[0] + 2. * span * (1. / (2. * N) + i / N),
                                      bbox[2] + 2. * span * (1. / (2. * N) + j / N),
                                      bbox[4] + 2. * span * (1. / (2. * N) + k / N)])
                    pos_c = pebbles[pebble, :3]
                    if (np.linalg.norm(pos_c - pos_v) < pebbles[pebble, 3]):
                        vox[i, j, k] = 1.
    return bbox, span, vox
 
@jit(nopython=True)

References BoundingBox().

Referenced by ComputeInertiaFromVoxelGrid().