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

# This example handles both demag and exchange.

import os,time,sys,math
import nmesh

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

#ocaml.init_hlib("/home/fangohr/build/HLib-1.3/Library/.libs/libhmatrix-1.3.so")
ocaml.init_hlib("/home/tf/HLib-1.3/Library/.libs/libhmatrix-1.3.so")

objects=[nmesh.ellipsoid([3.0,3.0,3.0])
         ]


mesh = nmesh.mesh(objects=objects,
                  a0=1.0,
                  bounding_box=[[-5.0,-5.0,-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))

import nmesh

cigar = nmesh.ellipsoid( [4,2] )

bbox = [[-5,-5],[5,5]]

mesh = nmesh.mesh(objects=[cigar], bounding_box=bbox, a0=0.5,
		  mesh_bounding_box=True)

nmesh.visual.plot2d_ps( mesh, "tutorial3.ps")

nfem.set_default_dimension(2)
nfem.set_default_order(1)

# Simulation parameters
sigma0 = 1.0

print
"""
** Example: meshing half ring and compute current density for contacts at either side.
This is for possible collaboration with Uni Hamburg. **
"""

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

ring = nmesh.difference(
    nmesh.ellipsoid([4.0, 4.0]), [nmesh.ellipsoid([2.5, 2.5])])

halfring = nmesh.intersect([ring, nmesh.box([-4.0, -0.0], [4.0, 4.0])])

the_mesh = nmesh.mesh(
    objects=[halfring],
    cache_name="halfring",
    a0=0.2,
    bounding_box=[[-4.0, -4], [4.0, 4.0]],
    neigh_force_scale=1.,
    initial_settling_steps=50,
    max_relaxation=4,
    max_steps=400)

nfem.set_default_mesh(the_mesh)

import nmesh


nx = 2  #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
rx,ry = 1,1.5 #radii of ellipsoids
angle = 45    #orientation of ellipsoids

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

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

#create the mesh
mesh = nmesh.mesh(objects=objects,a0=0.5,bounding_box=bbox)

#Create post script plot of mesh
nmesh.visual.plot2d_ps(mesh,"ellipsoid_array.ps")

#save plot to file
mesh.save('ellipsoid_array.nmesh')

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')

import nmesh

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

bbox = [[-1, -1.5, -1.5, -1], [1, 1.5, 1.5, 1]]

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

mesh.save("simple4d.nmesh")
print "should have saved mesh now"

# create gnuplot plot for manual (in other file not relevant here)


import os, os.path, sys, time

filename = "run_simple4d/simple4d.nmesh"
if not os.path.exists(filename):
    print "You need to run simple4d.py first to compute %s" % filename
    sys.exit(1)


def gnuplot_2d_points(points, filename):

    print "Starting to write %s" % filename,
    print time.asctime()

    """Given a list of pairs like
      points = [ [x0,y0], [x1,y1], [x2,y2], ..., [xN,yN]]

      and a filename, this will create a postscriptfile of name

                          J=SI(13.0e-12,"J/m"),
                          #Ms=1.,                #Matteo, very strange: if I use the SI units above,
                          #J=13.,                #then I get a wierd magsim-brain failure!
                          anisotropy_order=2,
                          anisotropy=test_anisotropy_energy
                          )

sim=nmag.SimulationContext("sphere")

# sim.timestepper_tuning_params=[1e-6,1e-6,2,300] # These are the defaults...
sim.timestepper_tuning_params=[1e-6,1e-6,4,300]

sim.set_magnetization([0.0,1.0,0.0])
# ^ just to show we can place set_magnetization() where we want...

sim.defregion("Py", nm.ellipsoid([4.0,4.0,5.0]), mag_mat=mat_Py)

def initial_magnetization(coords,mag_type):
    return [0.8*math.sin(coords[0]/1.0)*1e6,
            0.8*math.cos(coords[0]/1.0)*1e6,
            0.6*1e6
            ]

sim.generate_mesh(([-30.0,-30.0,-50.0],[30.0,30.0,50.0]), # bounding box
                  a0=2.0,
                  max_steps=400,
                  cache_name="sphere3"
                  )

#
#sim.set_magnetization(initial_magnetization)

"""Compounds outer_skin() with order_mesh() to obtain a cross-section
   of the outer skin of an ellipsoid.

   Author: James Kenny           Last modified: $Date$
"""
import nmesh

ellipsoid = nmesh.ellipsoid([1.5,2.5,2])

# create mesh
bbox = [[-1.5,-2.5,-2],[1.5,2.5,2]]

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

# now extract the outer layer of cells and order the mesh in the y-axis
# to see an outer skin with only elements containing at least 2 surface
# points, the user should uncomment the end of the 2nd following line
mesh_info = mesh.tolists()
mesh_info = nmesh.visual.outer_skin(mesh_info)#, condition='>=2')


# visualise with in2circ as a solid and order it in the y-axis
v = nmesh.visual.solid_in2circ(mesh_info, order=1)

    

import nmesh

ell = nmesh.ellipsoid([3,2])            # create ellipsoid
cone = nmesh.conic([-3,0],2,[3,0],0)    # create cone

inters = nmesh.intersect([ell, cone])   # create intersection of objects 

bbox = [[-5.,-4.],[5.,4.]]
mesh_ex = nmesh.mesh(objects = [inters], a0=0.4, bounding_box=bbox)

nmesh.visual.plot2d_ps(mesh_ex,"intersection.ps")

import nmesh

radius = 4
sphere = nmesh.ellipsoid([radius,radius])

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

# periodicity on x-axis
mesh_ex = nmesh.mesh(objects = [sphere], bounding_box=bbox, 
		     mesh_bounding_box=True,periodic=[True,False])

nmesh.visual.plot2d_ps(mesh_ex,"periodic.ps")


mesh_ex.save('periodic.nmesh')

import nmesh

big = nmesh.ellipsoid([4.0,3.0])    # create a big ellipsoid
small = nmesh.ellipsoid([3.0,2.0])  # small ellipsoid

diff = nmesh.difference(big,[small])# create difference of ellipsoids 

bbox = [[-5.,-4.],[5.,4.]]
mesh_ex = nmesh.mesh(objects = [diff], a0=0.4, bounding_box=bbox)

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

import nmesh

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)
#save plot to file
mesh.save("ellipsoid_array3d.nmesh")

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

h_total_Funny[2] += 0.0;
""")

mag.set_default_material(PermAlloy)


mag.set_intensive_parameters(["T","p","H_x","H_y","H_z"])

#mag.defregion("Ball 1",nm.ellipsoid([3.0,3.0,3.0],transform=[("shift",[-3.0,0.0,0.0])]))

print "OK 1"
sys.stdout.flush()

#mag.defregion("Ball 2",nm.ellipsoid([3.0,3.0,3.0],transform=[("shift",[3.0,0.0,0.0])]))

mag.defregion("Ball 2",nm.ellipsoid([2.0,2.0,2.0],transform=[("shift",[3.0,0.0,0.0])]))

# Note: clearly, we DO need a better way to specify geometries. Ideally, I would like to be
# able to write instead:
#
# mag.defregion("Ball 1",nm.shifted([-3,0,0],nm.sphere(3)))
# mag.defregion("Ball 2",nm.shifted([ 3,0,0],nm.sphere(3)))
#
# or alternatively:
#
# sphere = nm.sphere(3)
# mag.defregion("Ball 1",nm.shifted([-3,0,0],sphere))
# mag.defregion("Ball 2",nm.shifted([ 3,0,0],sphere))


mag.set_meshing_parameters(cache_name="two-balls")

import nmesh


def scale_list_of_lists( vec_data, factor ):
    """silly function -- just to change the vector data somehow.

    In real life, the simulation would provide this data, obviously."""
    v2=[]
    for i in range(len(vec_data)):
        v2.append(vec_data[i] + Numeric.array([-1,0,0])*factor)

    return v2


#create mesh
ellipsoid = nmesh.ellipsoid([1,1,0.5])
bbox = [[-1,-1,-1],[1,1,1]]
mesh = nmesh.mesh(objects=[ellipsoid],bounding_box=bbox,a0=0.5,cache_name="visual_plot_vectorfield")

#visualise
import nmesh.visual
meshinfo=mesh.tolists()

#create some vector field defind on vertices (=nodes)
vec_data = []
points=meshinfo[0][2]
origin = Numeric.array([0,0,0])
for point in points:
    vec_data.append(Numeric.array(point) - origin) 

#plot the vector field with changing data

rings_box = nmesh.box([-16.95, -14.95],
                      [16.95,14.95])

bbox = [[-half_bbox_dimx,-half_bbox_dimy],[half_bbox_dimx,half_bbox_dimy]]

fix_pts = []
rings = []
voltage_probes = []

for x in range(-1,2):
    for y in range(-1,2):
        
        # external ring
        R = nmesh.ellipsoid(
            [outer_radius, outer_radius],
            transform=[("shift",[x*centres_distance, y*centres_distance])]
            )
        #internal ring
        r = nmesh.ellipsoid(
            [inner_radius, inner_radius],
            transform=[("shift",[x*centres_distance, y*centres_distance])]
            )
        # take the difference
        rings.append(nmesh.difference(R,[r]))

# cut the rings at top and bottom
union_rings = nmesh.union(rings)
rings_array = nmesh.intersect([union_rings,rings_box])


# semiaxis of a rhombus on the diagonals of the array 

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")

                          )

mat_Fe = nmag.MagMaterial("Fe",
                          Ms=SI(1.7e6,"A/m"),
                          J=SI(2.07e-11,"J/m")
                          )

sim=nmag.SimulationContext("sphere")

# sim.timestepper_tuning_params=[1e-6,1e-6,2,300] # These are the defaults...
sim.timestepper_tuning_params=[1e-6,1e-6,4,300]

sim.set_magnetization([0.0,1.0,0.0])
# ^ just to show we can place set_magnetization() where we want...

sim.defregion("Py", nmesh.ellipsoid([2.0,2.0,2.0],transform=[('shift',[5,0,0])]), mag_mat=mat_Py)
sim.defregion("Fe", nmesh.ellipsoid([2.0,2.0,2.0]), mag_mat=mat_Fe)


def initial_magnetization(coords,mag_type):
    if mag_type == "m_Fe":
        return [1,0,0]
    else:
        return [0.8*math.sin(coords[0]/3.0)*1e6,
                0.8*math.cos(coords[0]/3.0)*1e6,
                0.6*1e6
                ]

sim.generate_mesh(([-2.5,-2.5,-2.5],[7.5,2.5,2.5]), # bounding box
                  a0=1.0,
                  max_steps=1200,

"""
** Example: meshing four rings & solving the laplace equation
   for a space-dependent resistivity that depends on outer parameters

   (Presumably, we will encounter bugs at our first try)
**
"""

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

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

rings = nmesh.union(
    [
        nmesh.difference(
            nmesh.ellipsoid([3.0, 3.0], transform=[("shift", [-2.5, 0.0])]),
            [nmesh.ellipsoid([1.0, 1.0], transform=[("shift", [-2.5, 0.0])])],
        ),
        nmesh.difference(
            nmesh.ellipsoid([3.0, 3.0], transform=[("shift", [2.5, 0.0])]),
            [nmesh.ellipsoid([1.0, 1.0], 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(

#nmag.set_log_level('debug')

intensive_param_by_name={"H_x":0.1,"H_y":0.0,"H_z":0.0}
# very slightly pulling in x-direction

mat_Py = nmag.MagMaterial("Py",
                          Ms=SI(1e6,"A/m"),
                          exchange_coupling=SI(13.0e-12,"J/m"),
                          )

sim=nmag.SimulationContext("sphere")

# sim.timestepper_tuning_params=[1e-6,1e-6,2,300] # These are the defaults...
sim.timestepper_tuning_params=[1e-6,1e-6,2,300]

sim.defregion("Py", nm.ellipsoid([3.0,3.0,3.0]), mag_mat=mat_Py)

def initial_magnetization(coords,mag_type):
    return [0.8*math.sin(coords[0]/1.0)*1e6,
            0.8*math.cos(coords[0]/1.0)*1e6,
            0.6*1e6
            ]

sim.generate_mesh(([-5.0,-5.0,-5.0],[5.0,5.0,5.0]), # bounding box
                  a0=0.7,
                  max_steps=450,
                  cache_name="sphere"
                  )

#sim.set_magnetization([1.0,0.0,0.0])
#sim.set_magnetization([0.0,1.0,0.0])