party.py

import bpy
import numpy as np

nx, ny, nz = 21, 21, 21
x0, y0, z0 = 10, 10, 10
r1, r2 = 6.0, 7.4

scale = 0.5

x, y, z = np.arange(nx), np.arange(ny), np.arange(nz)

rsq = (x[None, None, :]-x0)**2 + (y[None, :, None]-y0)**2 + (z[:, None, None]-z0)**2
sphere_1 = rsq < r1**2
sphere_2 = rsq < r2**2
shell = sphere_2 * (-sphere_1)    # I love python!

X, Y, Z = np.meshgrid(x, y, z)

ipoints = zip(
    X.flatten(), 
    Y.flatten(), 
    Z.flatten()
)

points = zip(
    scale*(X.flatten()-x0), 
    scale*(Y.flatten()-y0), 
    scale*(Z.flatten()-z0)
) 

points = [i for i in points]
ipoints = [i for i in ipoints]

points_shell = [thing for i, thing in enumerate(points) if shell[ipoints[i]]]

def create_repr_plane(obj_name, mesh_name):

    s = scale/2  
    verts = []
    v_add = verts.extend
    faces = []
    f_add = faces.extend

    def get_verts(x, y, z):
        return [
            (s+x, s+y, z-s),
            (s+x, -s+y, z-s),
            (-s+x, -s+y, z-s),
            (-s+x, s+y, z-s),
            (s+x, s+y, z+s),
            (s+x, -s+y, z+s),
            (-s+x, -s+y, z+s),
            (-s+x, s+y, z+s)
        ]

    def get_faces(o):
        return [
            (0+o, 1+o, 2+o, 3+o),
            (4+o, 7+o, 6+o, 5+o),
            (0+o, 4+o, 5+o, 1+o),
            (1+o, 5+o, 6+o, 2+o),
            (2+o, 6+o, 7+o, 3+o),
            (4+o, 0+o, 3+o, 7+o)
        ]

    # for each surface point
    total_idx = 0

    for point in points_shell:
        v_add(get_verts(*point))
        f_add(get_faces(total_idx))
        total_idx += 8

    profile_mesh = bpy.data.meshes.new(obj_name)
    profile_mesh.from_pydata(verts, [], faces)
    profile_mesh.update()

    profile_object = bpy.data.objects.new(obj_name, profile_mesh)

    bpy.context.scene.objects.link(profile_object)
    return profile_object


obj = create_repr_plane("dupli_object", "dupli_mesh")