def create_mesh():

    nx = 3 #how many ellipsoids in x direction
    x0 = 3. #spacing of ellipsoids in x
    ny = 3 #how many ellipsoids in y direction
    y0 = 3. #spacing of ellipsoids in x
    nz = 3 #how many ellipsoids in y direction
    z0 = 3. #spacing of ellipsoids in x
    rx,ry,rz = 0.8,1.2,1#radii of ellipsoids

    #create list 'objects' with ellipsoids
    objects = []
    for i in range(nx):
        for j in range(ny):
            for k in range(nz):
                objects.append( nmesh.ellipsoid([rx,ry,rz], [("shift",[i*x0,j*y0,k*z0])]))

    #bounding box
    bbox = [[-2,-2,-2],[(nx-1)*x0+2,(ny-1)*y0+2,(nz-1)*z0+2]]

    #create the mesh
    mesh = nmesh.mesh(objects=objects,a0=0.75,bounding_box=bbox,mesh_bounding_box=True)
    #save plot to file
    mesh.save("test.nmesh",directory='.')

    return mesh

    # MayaVi captures the screen to do this, so ensure there are no
    # windows on top of MayaVi. (hence the unattended mode should be
    # set up to savefig() but not to draw() 

    pylab.draw()
#    pylab.savefig('hist_'+str( globals()['call_counter']).zfill(4)+'.png')
    nmesh.visual.export_visualisation(globals()['v'],'mesh_'+str( globals()['call_counter']).zfill(4)+'.png')

    # increment the call counter
    globals()['call_counter'] += 1

    # to run attended, comment in the following 2 lines

# create the mesh and visualise it every n steps
N = 1000
mesh = nmesh.mesh(objects = [box,cone], a0=rod, bounding_box=bbox, callback=(my_function, N))

# show the finished mesh
mesh_info = mesh.tolists()

# to do this by writing to disk, uncomment the following 2 lines
globals()['v'], in2circ = nmesh.visual.solid_in2circ(mesh_info, myv=globals()['v'])
globals()['v'].renwin.z_plus_view()

# if we have written to disk first (ie the above line is being used)
# then we can save the VTK file by uncommenting the following line
shutil.copy('in2circ_solid.vtk', 'mesh_'+str(globals()['call_counter']).zfill(4)+'.vtk')

# save an image from MayaVi of the finished mesh
nmesh.visual.export_visualisation(globals()['v'],'mesh_'+str(call_counter).zfill(4)+'.png')

# (1) It is somewhat problematic/annoying that the first script does
# not really make clear when a name is a MWE name and when it is a
# field/vector name. Corrected that for this example.

import os,time,sys,math
import nmesh

execfile("../../interface/nsim/linalg_machine.py")

objects=[nmesh.difference(nmesh.ellipsoid([3.0,3.0],[("shift",[0,0])]),
                          [nmesh.ellipsoid([2.5,2.5],[("shift",[0,0])])])
         ]

mesh = nmesh.mesh(objects=objects,
                  a0=0.25,
                  bounding_box=[[-5.0,-5.0],[5.0,5.0]],
                  cache_name="geometry.mesh",
                  )

raw_mesh=mesh.raw_mesh

if ocaml.petsc_is_mpi():
    print "*** PARALLEL EXECUTION ***"
    nr_nodes=ocaml.petsc_mpi_nr_nodes()
    nr_points=ocaml.mesh_nr_points(raw_mesh)
    z=nr_points/nr_nodes
    distrib = [int(round(z*(i+1)))-int(round(z*i)) for i in range(0,nr_nodes)]
    slack=nr_points-reduce(lambda x,y:x+y,distrib)
    distrib[0] = distrib[0] + slack
    print "*** RAW MESH %s *** DISTRIB %s ***" %(repr(raw_mesh),repr(distrib))
    ocaml.mesh_set_vertex_distribution(raw_mesh,distrib)

"""This script demonstrates the function 'separate_parts()' compounded with
   'outer_skin()'.
   In the first dataset displayed a triangle and square are present. In the
   second dataset the square has been removed. 

   Author: James Kenny  Last modified: $Date$
"""
import nmesh
import nmesh.visual as viz


# define a simple mesh and submit request to mesher
box=nmesh.box( [0.0,0.0], [1.0,1.0] )
cone = nmesh.conic([3.0,0.0],1.0,[3.0,4.0],0.0)
bbox = [[-1.,-1.],[7.,6.]]
mesh = nmesh.mesh(objects = [box,cone],a0=0.4, bounding_box=bbox)
mesh_info=mesh.tolists()




# create mesh_info lists for parts 1&2 combined and for part 2 alone
[mesh_info1, mesh_info2] = viz.separate_parts(mesh_info, listOfParts=[[1,2],[2]])


# generate a VTK dataset for parts 1&2 combined
vtkData, points, simplices,indices, icradii, ccradii = \
         viz.mesh2vtk(mesh_info1, VTKonly=False)
in2circ=viz.findRatios(icradii, ccradii, factor=2)
vtkData=viz.append2vtk(vtkData, in2circ, "in2circ")
vtkData=viz.append2vtk(vtkData, indices, "part indices")

    stack5=[]
    stack6=[]

    iteration = [1,5,10,50,200,500,1000,2000]

    #mesh object at different maximum number of iterations N
    for N in iteration:
        timing.start()
        mesh = nmesh.mesh(objects=[box, ellipsoid], bounding_box=bbox,\
                          a0=0.3, mesh_bounding_box=True,

                          fixed_points = fix,

                          max_steps = N,

                          density = density,

                          neigh_force_scale = 0.0,

                          shape_force_scale = shape_elem,

                          volume_force_scale = 0.0
                          )
        timing.finish()
        time.append(timing.seconds())
        #iteration.append(N)

        #save the mesh as a .ps file in temp dir
        nmesh.visual.plot2d_ps(mesh,"fig_mesh_shape_%d_iter%06d.ps"\
                               % (shape_elem,N) ) 

                     nmesh.difference(nmesh.ellipsoid([3.0,3.0],
                                                      transform=[("shift",[2.5,0.0])]),
                                      [nmesh.ellipsoid([1.5,1.5],
                                                       transform=[("shift",[2.5,0.0])])])])

boxed_rings=nmesh.intersect([rings,nmesh.box([-8.0,-2.5],[8.0,2.5])])

N = 100
density = "density=1.;"

the_mesh = nmesh.mesh(objects = [boxed_rings],
                      cache_name="rings-mesh",
                      a0=0.3,
                      bounding_box=[[-10.0,-3.5],[10.0,3.5]],
                      neigh_force_scale = 1.,
                      density = density,
                      initial_settling_steps = 50,
                      max_relaxation = 4,
                      # callback=(my_function, N),
                      # max_steps=677
                      max_steps=200
                      )

nfem.set_default_mesh(the_mesh)

##### Making the elements... #####

empty_element=ocaml.empty_element;

# conductivity (scalar)
element_sigma     = nfem.make_element("sigma",[]);
element_drho_by_dt= nfem.make_element("drho_by_dt",[]);

# temporarily commented out by turning this into a string...
"""
the_mesh = nm.mesh(([-5.0,-5.0,-5.0],[5.0,5.0,5.0]), # bounding box
                   objects= \
                    [nm.difference( \
                      nm.intersect([nm.conic([0.0,0.0,-2.0],7.0,
                                             [0.0,0.0,2.0],7.0),
                                    nm.box([-4.0,-4.0,-4.0],[4.0,4.0,4.0])]))],
                   a0=1.2,
                   max_steps=500,
                   cache_name="ddiffop_mesh")
"""

mesh_obj = nm.mesh(([-5.0,-5.0,-5.0],[5.0,5.0,5.0]), # bounding box
                   objects= [nm.conic([0.0,0.0,-2.0],4.0,[0.0,0.0,2.0],4.0)],
                   a0=1.2,
                   max_steps=500,
                   cache_name="ddiffop_mesh")

the_mesh = mesh_obj.raw_mesh

print "MESH: ",the_mesh,type(the_mesh)


elem_T = ocaml.make_element("T",[],3,1)
elem_sigma = ocaml.make_element("sigma",[],3,1)
elem_j_q = ocaml.make_element("j_q",[3],3,1)

mwe_T=ocaml.make_mwe("mwe_T",the_mesh,[(1,elem_T)],[])
mwe_sigma=ocaml.make_mwe("mwe_sigma",the_mesh,[(1,elem_sigma)],[])
mwe_j_q=ocaml.make_mwe("mwe_j_q",the_mesh,[(1,elem_j_q)],[])

import nmesh

ellipsoid = nmesh.ellipsoid([0.75,1.25,1])

# create mesh
bbox = [[-0.75,-1.25,-1],[0.75,1.25,1]]

mesh = nmesh.mesh(objects = [ellipsoid], bounding_box=bbox,a0=0.5)

#create 3d-plot of surfaces and export eps
vis = nmesh.visual.show_bodies_mayavi(mesh)
nmesh.visual.export_visualisation(vis,"simple3d.eps")

#save mesh also as nmesh file
mesh.save('simple3d.nmesh')

bbox = [[-1.25,-0.75],[1.25,0.75]]

#Define all parameters, for example in string:
myconf = """
[nmesh-2D]
a0                     : 1.0
shape_force_scale      : 0.1
volume_force_scale     : 0.0
neigh_force_scale      : 1.0
irrel_elem_force_scale : 0.0
thresh_add             : 0.6
thresh_del             : 1.6
initial_settling_steps : 200
sliver_correction      : 1.0
max_relaxation         : 3.0
topology_threshold     : 0.2
max_relaxation         : 3.0
topology_threshold     : 0.2
time_step_scale        : 0.1
tolerated_rel_move     : 0.002
max_steps              : 2000
"""

#Then create MeshingParameters object 
mp = nmesh.MeshingParameters(string=myconf)

#and use MeshingParameter object when computing the mesh:
mesh = nmesh.mesh(objects = [ellipsoid], bounding_box=bbox,a0=0.5,\
                  meshing_parameters=mp)

# For now, we use a very very simple mesh...

double_disc_2d = nmesh.union([nmesh.ellipsoid([3.0,3.0], transform=[("shift",[-1.0,0.0])]),
                              nmesh.ellipsoid([3.0,3.0], transform=[("shift",[1.0,0.0])])])

double_disc_3d = nmesh.union([nmesh.conic([-1.0,0.0,thickness2d*0.5],3.0,[-1.0,0.0,-thickness2d*0.5],3.0),
                              nmesh.conic([ 1.0,0.0,thickness2d*0.5],3.0,[ 1.0,0.0,-thickness2d*0.5],3.0)])

density = "density=1.;"


mesh_2d = nmesh.mesh(objects = [double_disc_2d],
                     cache_name="double-disc-2d",
                     a0=0.4,
                     bounding_box=[[-5.0,-5.0],[5.0,5.0]],
                     neigh_force_scale = 1.,
                     density = density,
                     initial_settling_steps = 50,
                     max_relaxation = 4,
                     max_steps=500
                     )

mesh_3d = nmesh.mesh(objects = [double_disc_3d],
                     cache_name="double-disc-3d",
                     a0=0.4,
                     bounding_box=[[-5.0,-5.0,-0.3],[5.0,5.0,0.3]],
                     neigh_force_scale = 1.,
                     density = density,
                     initial_settling_steps = 50,
                     max_relaxation = 4,
                     max_steps=500
                     )

import nmesh

large = nmesh.box( [-4,-4], [4,4] )
small = nmesh.box( [-2,-2], [2,2] )
diff = nmesh.difference(large,[small])

bbox=[[-4,-4],[4,4]]
mesh = nmesh.mesh(bounding_box=bbox, objects = [diff] )
nmesh.visual.plot2d_ps( mesh, "fixedpoints_faulty.ps")


#can provide some 'fixed points' to avoid round corners
fixed_points = [[-2,-2],[2,-2],[-2,2],[2,2]]
mesh = nmesh.mesh(bounding_box=bbox, objects = [diff], fixed_points=fixed_points )
nmesh.visual.plot2d_ps( mesh, "fixedpoints.ps")

# Using the discretized differential operator machinery
# from within python...
#

# NOTE: Mesher (version from 30-04-2007) has difficulty here, cannot
# generate the mesh. Reverting to the mesher version from 01-04-2007 works.

import nmag2 as nmag, nmesh as nm, sys, math, time
import nfem
import nfem.visual


mesh_obj = nm.mesh(([-6.0,-6.0,-6.0],[6.0,6.0,6.0]), # bounding box
                   objects= [nm.conic([0.0,0.0,-4.0],2.0,[0.0,0.0,4.0],2.0),
                             nm.conic([0.0,0.0,-5.0],2.0,[0.0,0.0,-4.0],2.0),
                             nm.conic([0.0,0.0,5.0],2.0,[0.0,0.0,4.0],2.0)
                             ],
                   a0=0.6,
                   max_steps=100,
                   cache_name="ddiffop_mesh_rod")

the_mesh = mesh_obj.raw_mesh

print "MESH: ",the_mesh,type(the_mesh)


elem_scalar = ocaml.make_element("S",[],3,1)
elem_vector = ocaml.make_element("V",[3],3,1)

mwe_scalar=ocaml.make_mwe("mwe_S",the_mesh,[(1,elem_scalar),(2,elem_scalar),(3,elem_scalar)],[])
mwe_vector=ocaml.make_mwe("mwe_V",the_mesh,[(1,elem_vector),(2,elem_vector),(3,elem_vector)],[])

import nmesh


# create square
box = nmesh.box([0.0,0.0],[1.0,1.0])

# create cone
cone = nmesh.conic([3.0,0.0],1.0,[3.0,4.0],0.0)

rod= 0.4
bbox = [[-1.,-1.],[7.,6.]]

# create mesh and save it in a file every 50 steps
# during its construction; with the show flag it
# is possible also to visualise it. 
mesh_ex = nmesh.mesh(objects = [box,cone], a0=rod, bounding_box=bbox)

#visualisation missing


# plot mesh
nmesh.visual.plot2d_ps(mesh_ex,"visual_timeseries2d_mesh.ps")

import nmesh

a=2.3

P1 = [0,0,0] #one corner of box
P2 = [a,a,a]#other corner
box = nmesh.box( P1,P2, use_fixed_corner_points=True)

bbox = [[0,0,0],[a,a,a]]

mesh = nmesh.mesh(bounding_box=bbox, objects=[box])

mv=nmesh.visual.show_bodies_mayavi( mesh )

nmesh.visual.export_visualisation(mv,'box_with_fixed_corner_points.eps')

    # the following two commands save images from pylab and MayaVi
    # MayaVi captures the screen to do this, so ensure there are no
    # windows on top of MayaVi. (hence the unattended mode should be
    # set up to savefig() but not to draw() 

    pylab.draw()
    raw_input()

C_bottom = [0,0,0]  #center of spiral
R_spiral = 3        #radius of spiral
C_top = [0,0,10]    #top of spiral
R_circle = 3        #radius of max circle along the spiral


helix = nmesh.helix( C_bottom, R_spiral, C_top, R_circle )

bbox = [[-7,-7,-1],[7,7,11]]

N = 5
rod = 0.5
# create mesh of three objects and bounding box
mesh_ex = nmesh.mesh(objects = [helix],
                     a0=rod,
                     bounding_box=bbox,

                     #callback= (my_function,N),
                     max_steps=1000
                     )

mesh_ex.save("helix.nmesh")

import nmesh

# create a number of objects
one = nmesh.ellipsoid([3.0,3.0])
two = nmesh.ellipsoid([3.0,3.0],transform=[("shift",[7,0])])
three=nmesh.box( [-4.0,-6], [10,-4] )

bbox = [[-5.,-8.],[11.,5.]]

# create mesh of three objects and bounding box
mesh_ex = nmesh.mesh(objects = [one,two,three], bounding_box=bbox,
                     mesh_bounding_box=True)

# plot mesh
nmesh.visual.plot2d_ps(mesh_ex,"multiobjects.ps")

import nfem, nmesh, math, sys

nfem.set_default_dimension(3)
nfem.set_default_order(1)


bar = nmesh.box([-5.0,-2.0,-2.0],[5.0,2.0,2.0])


the_mesh = nmesh.mesh(objects = [bar],
                      cache_name="bar-mesh3",
                      a0=0.7,
                      bounding_box=[[-6.0,-3.0,-3.0],[6.0,3.0,3.0]],
                      neigh_force_scale = 1.,
                      # density = density,
                      initial_settling_steps = 50,
                      max_relaxation = 4,
                      # callback=(my_function, N),
                      # max_steps=677
                      max_steps=500
                      )

nfem.set_default_mesh(the_mesh)

element_M = nfem.make_element("M",[3])
element_H = nfem.make_element("H",[3])

mwe_M      = nfem.make_mwe("mwe_M", [(1,element_M)])
mwe_H      = nfem.make_mwe("mwe_H", [(1,element_H)])

diffop_v_laplace = nfem.diffop("-<d/dxj H(k) || d/dxj M(k)>, j:3, k:3")

    x, y = node
    
    for i in range(nr_layers+1):
        layer_heigth = -thickness/2. + layer_thick*i

        mobile_pts.append([x, y, layer_heigth])
        

bbox = [[-edge,-edge,-thickness/2.],[edge,edge,thickness/2.]]


# box
box = nmesh.box([-edge,-edge,-thickness/2.],[edge,edge,thickness/2.])

density = """density = 1.;"""
# create mesh of three objects and bounding box
mesh_ex = nmesh.mesh(objects = [box],
                     cache_name = "",
                     a0=rod,
                     mobile_points = mobile_pts,
                     bounding_box=bbox,
                     density = density,
                     neigh_force_scale = 0.0,
                     shape_force_scale = 0.0,
                     max_steps = 1
                     )

mesh_ex.save("box-2d-to-3d.nmesh")

nfem.set_default_dimension(3)
nfem.set_default_order(1)

##### Creating the mesh #####

# For now, we use a very very simple mesh...

ball = nmesh.ellipsoid([3.0,3.0,3.0])
density = "density=1.;"


the_mesh = nmesh.mesh(objects = [ball],
                      cache_name="bem-ball",
                      a0=1.0,
                      bounding_box=[[-5.0,-5.0,-5.0],[5.0,5.0,5.0]],
                      neigh_force_scale = 1.,
                      density = density,
                      initial_settling_steps = 50,
                      max_relaxation = 4,
                      max_steps=500
                      )

nfem.set_default_mesh(the_mesh)

##### Making the elements... #####

element_M     = nfem.make_element("M",[3]);
element_H     = nfem.make_element("H",[3]);
element_rho_M = nfem.make_element("rho_M",[]);
element_phi_M = nfem.make_element("phi_M",[]);

mwe_M      = nfem.make_mwe("mwe_M",     [(1,element_M)])

##    pylab.plot([px],[py],"bo",markersize = 10)
##    print p
##pylab.show()
##raw_input()

N = 100

# create mesh of three objects and bounding box
mesh_ex = nmesh.mesh(objects = [rings_array],#[outer,voltage_probes[0], voltage_probes[1], rings_array],
                     a0=rod,
                     bounding_box=bbox,
                     fixed_points = fix_pts,
                    
                     neigh_force_scale = 1.,
                     shape_force_scale = 0.1,

                     density = density,
                     thresh_del = 1.4,
                     
                     initial_settling_steps = 100,
                     max_relaxation = 4,
                     
                     #callback=(my_function, N),
                     max_steps=1000
                     )
#mesh_loaded = nmesh.load("9rings_alone.nmesh")
mesh_ex.save("9rings.nmesh")
nmesh.visual.plot2d_ps(mesh_ex,"9rings.ps")