Python mesh.Mesh类代码示例

    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)