def readMesh(self,ncoords,nelems,nplex,props,eltype,normals,sep,objtype='Mesh'):
        """Read a Mesh from a pyFormex geometry file.

        The following arrays are read from the file:
        - a coordinate array with `ncoords` points,
        - a connectivity array with `nelems` elements of plexitude `nplex`,
        - if present, a property number array for `nelems` elements.

        Returns the Mesh constructed from these data, or a subclass if
        an objtype is specified.
        """
        # Make sure to import the Mesh subclasses that can be read
        from plugins.trisurface import TriSurface

        ndim = 3
        x = readArray(self.fil,Float,(ncoords,ndim),sep=sep)
        e = readArray(self.fil,Int,(nelems,nplex),sep=sep)
        if props:
            p = readArray(self.fil,Int,(nelems,),sep=sep)
        else:
            p = None
        M = Mesh(x,e,p,eltype)
        if objtype != 'Mesh':
            try:
                clas = locals()[objtype]
            except:
                clas = globals()[objtype]
            M = clas(M)
        if normals:
            n = readArray(self.fil,Float,(nelems,nplex,ndim),sep=sep)
            M.normals = n
        return M

def readTetgen(fn):
    """Read and draw a tetgen file.

    This is an experimental function for the geometry import menu.
    """
    res = {}
    base,ext = os.path.splitext(fn)
    if ext == '.node':
        nodes = readNodeFile(fn)[0]
        res['tetgen'+ext] = nodes
    elif ext in [ '.ele', '.face' ]:
        nodes,nodenrs = readNodeFile(utils.changeExt(fn,'.node'))[:2]
        if ext == '.ele':
            elems = readEleFile(fn)[0]
        elif ext == '.face':
            elems = readFaceFile(fn)[0]
        if nodenrs.min() == 1 and nodenrs.max()==nodenrs.size:
            elems = elems-1
        M = Mesh(nodes,elems,eltype=elems.eltype)
        res['tetgen'+ext] = M

    elif ext == '.smesh':
        nodes,elems = readSmeshFile(fn)
        ML = [ Mesh(nodes,elems[e]) for e in elems ]
        res = dict([('Mesh-%s'%M.nplex(),M) for M in ML])

    elif ext == '.poly':
        nodes,elems = readPolyFile(fn)
        ML = [ Mesh(nodes,elems[e]) for e in elems ]
        res = dict([('Mesh-%s'%M.nplex(),M) for M in ML])

    return res

def areas(self):
    """area of elements

    For surface element the faces' area is returned.
    For volume elements the sum of the faces'areas is returned.

    """

    #In case of quadratic faces, the face's area should be 
    #the area inside the polygon of face vertices or 
    #the area of the equivalent linear face?

    ##this function would require some changes (here proposed inside the function as starting):
    ##create a _default_surfacetype to create quad8 instead of hex8 ?maybe also a _default_volumetype to create tet4 instead of quad4 ?

    def defaultSurfacetype(nplex):
        """Default face type for a surface mesh with given plexitude.

        For the most common cases of plexitudes, we define a default face
        type. The full list of default types can be found in
        mesh._default_facetype.
        """
        return _default_surfacetype.get(nplex,None)
    import geomtools
    nfacperel= len(self.eltype.faces[1])#nfaces per elem
    mf=Mesh(self.coords, self.getFaces())#mesh of all faces
    mf.eltype = elementType(defaultSurfacetype(mf.nplex()))
    ntriperfac= mf.select([0]).convert('tri3').nelems()#how many tri per face
    elfacarea = geomtools.areaNormals( mf.convert('tri3').toFormex()[:])[0].reshape(self.nelems(), nfacperel*ntriperfac)#elems'faces'areas
    return elfacarea.sum(axis=1)#elems'areas

    def crop(self, x_0, x_1, y_0, y_1, with_normals=False):
        """
        This function generates another mesh (not depthMesh) by cropping the organized set of veftices (a depth ROI)

        """
        w, h = self.depth.shape[1], self.depth.shape[0]
        if x_0 < 0 or y_0 < 0 or x_1 >= w or y_1 >= h or x_1 <=x_0 or y_1 <= y_0:
            return None
        mesh = Mesh()
        vcs_aux = self.vcs.reshape(self.depth.shape[0], self.depth.shape[1], 3)
        # Vertices
        mesh.vcs = vcs_aux[y_0:y_1, x_0:x_1, :].reshape(-1, 3)
        mesh.vcs_q = mesh.vcs.shape[0]
        # Normals
        if with_normals:
            vnorms_aux = self.vnorms.reshape(self.depth.shape[0], self.depth.shape[1], 3)
            mesh.vnorms = vnorms_aux[y_0:y_1, x_0:x_1, :].reshape(-1, 3)
        # Facets
        mesh.faces = DepthMesh.genFacets(x_1 - x_0, y_1 - y_0)
        mesh.faces_q = mesh.faces.shape[0]
        # texture mapping
        txcoord_aux = self.txcoord.reshape(self.depth.shape[0], self.depth.shape[1], 2)
        mesh.txcoord = txcoord_aux[y_0:y_1, x_0:x_1, :].reshape(-1, 2)
        mesh.texture = self.texture
        mesh.txwidth, mesh.txheight = self.txwidth, self.txheight
        return mesh

def gen_mesh_from_voxels_mc(voxels, voxelsize,
                            gmsh3d=False, scale_factor=0.25):
    import scipy.spatial as scsp

    tri = marching_cubes(voxels, voxelsize)

    nel, nnd, dim = tri.shape
    coors = tri.reshape((nel * nnd, dim))
    tree = scsp.ckdtree.cKDTree(coors)
    eps = nm.max(coors.max(axis=0) - coors.min(axis=0)) *1e-6
    dist, idx = tree.query(coors, k=24, distance_upper_bound=eps)

    uniq = set([])
    for ii in idx:
        ukey = ii[ii < tree.n]
        ukey.sort()
        uniq.add(tuple(ukey))

    ntri = nm.ones((nel * nnd,), dtype=nm.int32)
    nnod = len(uniq)
    ncoors = nm.zeros((nnod, 3), dtype=nm.float64)

    for ii, idxs in enumerate(uniq):
        ntri[nm.array(idxs)] = ii
        ncoors[ii] = coors[idxs[0]]

    mesh = Mesh.from_data('voxel_mc_data',
                          ncoors, nm.ones((nnod,), dtype=nm.int32),
                          {0: nm.ascontiguousarray(ntri.reshape((nel, nnd)))},
                          {0: nm.ones((nel,), dtype=nm.int32)},
                          {0: '%d_%d' % (2, 3)})

    if gmsh3d:
        from vtk2stl import vtk2stl
        import tempfile
        import os

        auxfile = os.path.join(tempfile.gettempdir(), 'dicom2fem_aux')
        vtk_fn = auxfile + '_surfmc.vtk'
        stl_fn = auxfile + '_surfmc.stl'
        geo_fn = auxfile + '_surf2vol.geo'
        mesh_fn = auxfile + '_volmv.mesh'
        mesh.write(vtk_fn)
        vtk2stl(vtk_fn, stl_fn)
        geofile = open(geo_fn, 'wt')
        geofile.write(gmsh3d_geo.replace('__INFILE__',
                                         stl_fn).replace('__SCFACTOR__',
                                                         str(scale_factor)))
        geofile.close()
        os.system('gmsh -3 -format mesh -o %s %s' % (mesh_fn, geo_fn))
        mesh = Mesh.from_file(mesh_fn)

    return mesh

    def createAnimatedModel(self, f11):
        wld = self.wld_container.wld_file_obj
        f10 = wld.getFragment(f11.fragRef)
        if f10.type != 0x10:
            print 'Model::createAnimatedModel() ERROR expected 0x10 fragment but got:', f10.type
            return
               
        # Lets initially try to only read all the mesh pieces and assemble the basic 
        # model. Once that is working we can start looking into animation
        
        # Loop over the parts of the model/skeleton: the entries list in the f10 fragment
        # define these
        root_mesh = None
        for i in range(0, f10.size1):
            
            if i > 0:
                f2d = wld.getFragment(f10.entries[i][3])    # entry[3] -> fragRef2
                # f2d.dump()
                f36 = wld.getFragment(f2d.fragRef)
                # f36.dump()

                m = Mesh(self.name+'_mesh_'+str(i))
                m.buildFromFragment(f36, self.wld_container, False)
                m.root.reparentTo(root_mesh.root)
            else: # the root node (index 0) does not have a mesh
                m = Mesh(self.name+'_mesh_'+str(i))  # empty dummy mesh
                root_mesh = m
                
            self.meshes.append(m)
            
            # get model part orientation data from 0x10->0x13->0x12 ref chain
            f13 = wld.getFragment(f10.entries[i][2])    # entry[2] -> fragRef1
            f12 = wld.getFragment(f13.fragRef)
            
            denom = float(f12.rotDenom)
            if denom != 0.0:
                rotx = f12.rotx/denom
                roty = f12.roty/denom
                rotz = f12.rotz/denom
                m.root.setHpr(rotx / 512.0 * 360.0, roty / 512.0 * 360.0, rotz / 512.0 * 360.0)
            
            
            denom = float(f12.shiftDenom)
            if denom != 0.0:
                shiftx = float(f12.shiftx)/denom
                shifty = float(f12.shifty)/denom
                shiftz = float(f12.shiftz)/denom
                # print shiftx, shifty, shiftz
                m.root.setPos(shiftx, shifty, shiftz)
            
        self.loaded = 1

def test_mesh():
    '''Basic test of the Mesh factory method and mesh extending'''
    ox, dx, nx = 0., 10., 5
    mesh = Mesh(type='uniform', ox=ox, lx=dx, nx=nx, block='B1')
    assert mesh.num_elem == nx
    assert mesh.boundary_nodes == [1, nx+1]
    assert np.allclose(mesh.boundary, [ox, dx])
    # verify extending the mesh
    dxb, nb = 4., 2
    mesh.extend(dxb, nb, block='B2')
    assert mesh.num_elem == nx + nb
    assert mesh.boundary_nodes == [1, nx+nb+1]
    assert np.allclose(mesh.boundary, [0., dx+dxb])
    assert len(mesh.nodes) == len(mesh.vertices)

    def load_ply(self, filename):
        Mesh.load_ply(self, filename)

        v_array=np.array(self.verts, dtype=np.float32)
        bbox=(np.min(v_array,0),  np.max(v_array,0) )

        self.v_array = v_array.astype(np.float32)
        self.bbox = bbox
        self.zoom=1.0/la.norm(bbox[1])

        self.tri_array = np.array(self.tris, dtype=np.uint32)
        self.n_array = self.vert_props['normal'].astype(np.float32)

        logging.debug( 'done matrix {}'.format(self.zoom) )
        return self.zoom

def gen_mesh_from_voxels_mc(voxels, voxelsize):
    import scipy.spatial as scsp

    tri = marching_cubes(voxels, voxelsize)
    
    nel, nnd, dim = tri.shape
    coors = tri.reshape((nel * nnd, dim))
    tree = scsp.ckdtree.cKDTree(coors)
    eps = nm.max(coors.max(axis=0) - coors.min(axis=0)) *1e-6
    dist, idx = tree.query(coors, k=24, distance_upper_bound=eps)

    uniq = set([])    
    for ii in idx:
        ukey = ii[ii < tree.n]
        ukey.sort()
        uniq.add(tuple(ukey))

    ntri = nm.ones((nel * nnd,), dtype=nm.int32)
    nnod = len(uniq)
    ncoors = nm.zeros((nnod, 3), dtype=nm.float64)

    for ii, idxs in enumerate(uniq):
        ntri[nm.array(idxs)] = ii
        ncoors[ii] = coors[idxs[0]]

    mesh = Mesh.from_data('voxel_mc_data',
                          ncoors, nm.ones((nnod,), dtype=nm.int32),
                          {0: nm.ascontiguousarray(ntri.reshape((nel, nnd)))},
                          {0: nm.ones((nel,), dtype=nm.int32)},
                          {0: '%d_%d' % (2, 3)})

    return mesh

    def from_conf(conf, init_fields=True, init_variables=True, init_equations=True, init_solvers=True):

        mesh = Mesh.from_file(conf.filename_mesh)

        eldesc_dir = op.join(install_dir, "eldesc")
        domain = Domain.from_mesh(mesh, eldesc_dir)
        domain.setup_groups()
        domain.fix_element_orientation()
        domain.setup_neighbour_lists()

        obj = ProblemDefinition(conf=conf, domain=domain, eldesc_dir=eldesc_dir)

        # Default output file trunk and format.
        obj.ofn_trunk = io.get_trunk(conf.filename_mesh)
        obj.output_format = "vtk"

        obj.set_regions(conf.regions, conf.materials, conf.funmod)

        if init_fields:
            obj.set_fields(conf.fields)

            if init_variables:
                obj.set_variables(conf.variables)

                if init_equations:
                    obj.set_equations(conf.equations)

        if init_solvers:
            obj.set_solvers(conf.solvers, conf.options)

        obj.ts = None

        return obj

    def save_state(
        self, filename, state=None, out=None, fill_value=None, post_process_hook=None, file_per_var=False, **kwargs
    ):
        extend = not file_per_var
        if (out is None) and (state is not None):
            out = self.state_to_output(state, fill_value=fill_value, extend=extend)
            if post_process_hook is not None:
                out = post_process_hook(out, self, state, extend=extend)

        float_format = get_default_attr(self.conf.options, "float_format", None)

        if file_per_var:
            import os.path as op

            meshes = {}
            for var in self.variables.iter_state():
                rname = var.field.region.name
                if meshes.has_key(rname):
                    mesh = meshes[rname]
                else:
                    mesh = Mesh.from_region(var.field.region, self.domain.mesh, localize=True)
                    meshes[rname] = mesh
                vout = {}
                for key, val in out.iteritems():
                    if val.var_name == var.name:
                        vout[key] = val
                base, suffix = op.splitext(filename)
                mesh.write(base + "_" + var.name + suffix, io="auto", out=vout, float_format=float_format, **kwargs)
        else:
            self.domain.mesh.write(filename, io="auto", out=out, float_format=float_format, **kwargs)

    def __init__(self, heightfield):
        mesh = Mesh()
       
        size = 32
        factor = 1.0
        vertices = []

        for z in xrange(size):
            z = float(z)/float(size-1)
            for x in xrange(size):
                x = float(x)/float(size-1)
                y = heightfield[x,z]
                vertices.append(mesh.vertex(x, y, z))

        for y in xrange(size-1):
            for x in xrange(size-1):
                v0 = vertices[(x+1) + (y+1)*size]
                v1 = vertices[(x+1) + (y+0)*size]
                v2 = vertices[(x+0) + (y+0)*size]
                v3 = vertices[(x+0) + (y+1)*size]

                mesh.face(v0, v1, v2)
                mesh.face(v3, v0, v2)

        splits = Splits(mesh, heightfield)
        while len(mesh.verts) < 21840:
        #while len(mesh.verts) < 3000:
            print len(mesh.faces), len(mesh.verts)
            splits.perform()

        mesh.save('mesh.bin')
        self.vbo = mesh.serialize()

	def draw( self ) :
		qm = tr.quaternion_matrix( self.Q[3:] )
		if self.drawstate & MESH :
			glPushMatrix()
			glMultTransposeMatrixf( qm )
			Mesh.draw( self )
			glPopMatrix()

		if self.drawstate & WIREFRAME :
			glPushMatrix()
			glMultTransposeMatrixf( qm )
			glDisable(GL_LIGHTING)
			glPolygonMode(GL_FRONT_AND_BACK, GL_LINE)
			glDisable(GL_CULL_FACE)
			Mesh.draw( self )
			glBegin(GL_LINES)
			glVertex3f(0,0,0)
			glVertex3f( self.x[-1,0] , self.x[-1,1] , self.x[-1,2] )
			glEnd()
			glEnable(GL_CULL_FACE)
			glPolygonMode(GL_FRONT_AND_BACK, GL_FILL)
			glEnable(GL_LIGHTING)
			glPopMatrix()

		if self.drawstate & TRACE :
			glDisable(GL_LIGHTING)
			glBegin(GL_POINTS)
			for p in self.trace : glVertex3f( *p[:3] )
			glEnd()
			glEnable(GL_LIGHTING)

		if self.drawstate & GRAVITY :
			glPushMatrix()
			glDisable(GL_LIGHTING)
			glTranslatef( 2 , 2 , 0 )
			glScalef(.1,.1,.1)
			glMultTransposeMatrixf( qm )
			glColor3f(1,.5,0)
			glBegin(GL_LINES)
			glVertex3f( 0 , 0 , 0 )
			glVertex3f( *self.G )
			glEnd()
			glEnable(GL_LIGHTING)
			glPopMatrix()

def test_linear_element():
    '''Test the implementation of the linear element with simple integrals'''
    dx = 5.
    num_elem = 5
    mesh = Mesh(type='uniform', ox=0., lx=dx, nx=num_elem)

    elem_num = 1
    connect = mesh.connectivity(elem_num)
    vertices = mesh.coordinates(connect)
    elem = LinearElement(elem_num, connect, vertices)

    # --- test some integrals
    ex = lambda x: x
    ex2 = lambda x: x ** 2

    # Integrate[phi]
    ans = elem.integrate()
    exact = integrate(N, (x, 0, 1))
    assert areclose(ans, exact)

    # Integrate[dphi]
    ans = elem.integrate(derivative=True)
    exact = integrate(dN, (x, 0, 1))
    assert areclose(ans, exact)

    # Integrate[x phi]
    ans = elem.integrate(ex)
    exact = integrate(x * N, (x, 0, 1))
    assert areclose(ans, exact)

    # Integrate[x dphi]
    ans = elem.integrate(ex, derivative=True)
    exact = integrate(x * dN, (x, 0, 1))
    assert areclose(ans, exact)

    # Integrate[x x phi]
    ans = elem.integrate(ex, ex)
    exact = integrate(x * x * N, (x, 0, 1))
    assert areclose(ans, exact)

    # Integrate[x x dphi]
    ans = elem.integrate(ex,ex,derivative=True)
    exact = integrate(x*x*dN,(x,0,1))
    assert areclose(ans, exact)

 def __init__(self, parent, color=0xFF0000, x=0, y=0):
     Mesh.__init__(self)
     self.parent = parent
     self.position.x = x
     self.position.y = y
     #self.matrix.translate(x, y, 0.0)
     self.color = color
     r, g, b = self.toRgb(self.color)
     v = [
         0.0, 0.0, r, g, b,
         0.0 + BLOCK_SIZE, 0.0, r, g, b,
         0.0 + BLOCK_SIZE, 0.0 + BLOCK_SIZE, r, g, b,
         0.0, 0.0 + BLOCK_SIZE, r, g, b,
     ]
     i = [
         0, 1, 2,
         2, 3, 0
     ]
     self.set_data(vertices=v, indices=i)

    def __init__(self, obj_file_name):
        self.default_shape_name = obj_file_name or 'initial_shape'
        self.default_surface_name = 'initial_surface'
        self.default_material_name = 'default'

        # Name of the shape/surface/material from the most recently parsed 'o'/'g'/'usemtl' line respectively
        self.curr_shape_name = self.default_shape_name
        self.curr_surf_name = self.default_surface_name
        self.curr_mtl_name = self.default_material_name

        # To keep of track of number of polygons parsed so far
        self.next_polygon_index = 0
        # A tuple of indices per vertex element
        self.indices = [ ]

        # Shortcut to be able to access the current surface in fewer characters and lookups
        self.curr_surf = Surface(0, self.default_material_name)
        # Dictionary of names -> shapes. Initialise to a default shape with a default surface
        self.shapes = { self.default_shape_name : Shape({ self.default_surface_name : self.curr_surf }) }
        Mesh.__init__(self)

def combine_meshes(mesh1, mesh2, ensure_continuity = False):
    """
    Combine two meshes into one disconnected mesh. This function
    relies on the user to make sure that nothing weird is going on.
    for example, the meshes probably should not intersect. I'm not
    sure what would happen if they do! Also, I assume that all
    the meshes are linear (linear mapping from real to reference space).

    Also, this function does not apply any mappings -- it assumes they have
    already been applied to the elements of the subordinate meshes.
    """
    vertices = copy.copy(mesh1.vertices)
    vertices.extend(mesh2.vertices)
    elements = copy.copy(mesh1.elements)
    elements.extend(mesh2.elements)

    result =  Mesh(vertices, elements)
    if ensure_continuity:
        result.condense_duplicate_vertices()
    return result

   def test_computeL1Error1D(self):
      # create mesh
      left = 2.5
      width = 1.5
      mesh = Mesh(3, width, x_start=left)

      # create "numerical solution" data
      numerical_solution = [(1.0,2.0),(1.6,2.2),(2.5,2.7)]

      # specify "exact solution" function
      def exact(x):
         return x**2 + 3.0

      # compute exact integral of difference
      exact_integral = 0.0
      for i in xrange(mesh.n_elems):
         # express local numerical solution as linear function y(x) = m*x + b
         el = mesh.getElement(i)
         xL = el.xl
         xR = el.xr
         yL = numerical_solution[i][0]
         yR = numerical_solution[i][1]
         dx = xR - xL
         dy = yR - yL
         m = dy / dx
         b = yL - xL*m

         # compute local integral of difference
         local_integral = (xR**3 - xL**3)/3.0 - 0.5*m*(xR**2 - xL**2)\
            + (3.0-b)*(xR - xL)
 
         # add to global integral of difference
         exact_integral += local_integral

      # compute numerical integral of difference
      numerical_integral = computeL1ErrorLD(mesh, numerical_solution, exact)

      # assert that numerical and exact integrals are approximately equal
      n_decimal_places = 14
      self.assertAlmostEqual(numerical_integral,exact_integral,n_decimal_places)

def close_loop_example():
    # one more example, originally not a closed loop curve
    F = Formex(pattern('11')).replic(2,1,1) + Formex(pattern('2')).replic(2,2,0)
    M = F.toMesh()
    draw(M,color='green')
    drawNumbers(M,color=red)
    drawNumbers(M.coords,color=blue)

    print "Original elements:",M.elems
    conn = connectivity.connectedLineElems(M.elems)
    if len(conn) > 1:
        message("This curve is not a closed circumference")
        return None
    
    sorted = conn[0]
    print "Sorted elements:",sorted

    showInfo('Click to continue')
    clear()
    M = Mesh(M.coords,sorted)
    drawNumbers(M)
    return M.toFormex()

	def set_block( self , s , d ) :
		aoy = m.atan2( s[2] , s[0] )
		aoz = m.atan2( s[1] , m.sqrt(s[0]**2+s[2]**2) )

		rot = tr.rotation_matrix( aoy , (0,1,0) )
		rot = np.dot( tr.rotation_matrix( -aoz , (0,0,1) ) , rot )
		rot = np.dot( tr.rotation_matrix( m.pi/2.0 , (0,0,1) ) , rot )

		v , n , t = self.gen_v( 1 , 1 , s )
		for x in range(v.shape[0]) :
			for y in range(v.shape[1]) :
				for z in range(v.shape[2]) :
					v[x,y,z] = np.dot(rot,v[x,y,z])
					n[x,y,z] = np.resize(np.dot(rot,np.resize(n[x,y,z],4)),3)
		Mesh.__init__( self , buffers = (v,n,t) )

		self.x = np.array( ((0,0,0,1),(s[0],0,0,1),(0,0,s[2],1),(s[0],0,s[2],1),(0,s[1],0,1),(s[0],s[1],0,1),(0,s[1],s[2],1),(s[0],s[1],s[2],1)) , np.float64 )
		for i in range(len(self.x)) : self.x[i] = np.dot(rot,self.x[i])
		self.r = np.resize( np.dot( rot , np.array((s[0],s[1],s[2],0) , np.float64 )/2.0 ) , 3 )
		self.m = np.array( [ d*s[0]*s[1]*s[2] / 8.0 ] * len(self.x) , np.float64 )
		self.M = self.calc_m( self.x , self.m )
		self.Mi = np.linalg.inv( self.M )