Python dolfin.Mesh类代码示例

    def setupMeshes(cls, mesh, N, num_refine=0, randref=(1.0, 1.0)):
        """Create a set of N meshes based on provided mesh. Parameters
        num_refine>=0 and randref specify refinement
        adjustments. num_refine specifies the number of refinements
        per mesh, randref[0] specifies the probability that a given
        mesh is refined, and randref[1] specifies the probability that
        an element of the mesh is refined (if it is refined at all).
        """
        assert num_refine >= 0

        assert 0 < randref[0] <= 1.0
        assert 0 < randref[1] <= 1.0

        # create set of (refined) meshes
        meshes = list();
        for _ in range(N):
            m = Mesh(mesh)
            for _ in range(num_refine):
                if randref[0] == 1.0 and randref[1] == 1.0:
                    m = refine(m)
                elif random() <= randref[0]:
                    cell_markers = CellFunction("bool", m)
                    cell_markers.set_all(False)
                    cell_ids = range(m.num_cells())
                    shuffle(cell_ids)
                    num_ref_cells = int(ceil(m.num_cells() * randref[1]))
                    for cell_id in cell_ids[0:num_ref_cells]:
                        cell_markers[cell_id] = True
                    m = refine(m, cell_markers)
            meshes.append(m)
        return meshes

    def test_convert_diffpack_2d(self):

        from dolfin import Mesh, MPI, MeshFunction, mpi_comm_world

        fname = os.path.join(os.path.dirname(__file__), "data", "diffpack_tri")
        dfname = fname+".xml"

        # Read triangle file and convert to a dolfin xml mesh file
        meshconvert.diffpack2xml(fname+".grid", dfname)

        # Read in dolfin mesh and check number of cells and vertices
        mesh = Mesh(dfname)

        self.assertEqual(mesh.num_vertices(), 41)
        self.assertEqual(mesh.num_cells(), 64)
        self.assertEqual(len(mesh.domains().markers(2)), 64)

        mf_basename = dfname.replace(".xml", "_marker_%d.xml")
        for marker, num in [(1,10), (2,5), (3,5)]:

            mf_name = mf_basename % marker
            mf = MeshFunction("size_t", mesh, mf_name)
            self.assertEqual(sum(mf.array()==marker), num)
            os.unlink(mf_name)

        # Clean up
        os.unlink(dfname)

def vtk_ug_to_dolfin_mesh(ug):
    """
    Create a DOLFIN Mesh from a vtkUnstructuredGrid object
    """
    if not isinstance(ug, vtk.vtkUnstructuredGrid):
        raise TypeError("Expected a 'vtkUnstructuredGrid'")
    
    # Get mesh data
    num_cells = int(ug.GetNumberOfCells())
    num_vertices = int(ug.GetNumberOfPoints())
    
    # Get topological and geometrical dimensions
    cell = ug.GetCell(0)
    gdim = int(cell.GetCellDimension())
    cell_type = cell.GetCellType()                                                                                                                                          
    if cell_type not in [vtk.VTK_TETRA, vtk.VTK_TRIANGLE]:                                                                                                                  
        raise TypeError("DOLFIN only support meshes of triangles " + \
                        "and tetrahedrons.")
    
    tdim = 3 if cell_type == vtk.VTK_TETRA else 2
    
    # Create empty DOLFIN Mesh
    mesh = Mesh()
    editor = MeshEditor()
    editor.open(mesh, tdim, gdim)
    editor.init_cells(num_cells)
    editor.init_vertices(num_vertices)
    editor.close()
    
    # Assign the cell and vertex informations directly from the vtk data
    cells_array = array_handler.vtk2array(ug.GetCells().GetData())
    
    # Get the assumed fixed size of indices and create an index array
    cell_size = cell.GetPointIds().GetNumberOfIds()
    cellinds = np.arange(len(cells_array))
    
    # Each cell_ids_size:th data point need to be deleted from the
    # index array
    ind_delete = slice(0, len(cells_array), cell_size+1)
    
    # Check that the removed value all have the same value which should
    # be the size of the data
    if not np.all(cells_array[ind_delete]==cell_size):
        raise ValueError("Expected all cells to be of the same size")
    
    cellinds = np.delete(cellinds, ind_delete)
    
    # Get cell data from mesh and make it writeable (cell data is non
    # writeable by default) and update the values
    mesh_cells = mesh.cells()
    mesh_cells.flags.writeable = True
    mesh_cells[:] = np.reshape(cells_array[cellinds], \
                              (num_cells , cell_size))
    
    # Set coordinates from vtk data
    vertex_array = array_handler.vtk2array(ug.GetPoints().GetData())
    if vertex_array.shape[1] != gdim:
        vertex_array = vertex_array[:,:gdim]
    mesh.coordinates()[:] = vertex_array
    return mesh

    def test_convert_diffpack(self):
        from dolfin import Mesh, MPI, MeshFunction
        if MPI.num_processes() != 1:
            return
        fname = os.path.join("data", "diffpack_tet")
        dfname = fname+".xml"
        
        # Read triangle file and convert to a dolfin xml mesh file
        meshconvert.diffpack2xml(fname+".grid", dfname)

        # Read in dolfin mesh and check number of cells and vertices
        mesh = Mesh(dfname)
        self.assertEqual(mesh.num_vertices(), 27)
        self.assertEqual(mesh.num_cells(), 48)
        self.assertEqual(mesh.domains().markers(3).size(), 48)
        self.assertEqual(mesh.domains().markers(2).size(), 16)

        mf_basename = dfname.replace(".xml", "_marker_%d.xml")
        for marker, num in [(3, 9), (6, 9), (7, 3), (8, 1)]:

            mf_name = mf_basename % marker
            mf = MeshFunction("uint", mesh, mf_name)
            self.assertEqual(sum(mf.array()==marker), num)
            os.unlink(mf_name)
        
        # Clean up
        os.unlink(dfname)

def make_mesh(coordinates, cells, tdim, gdim, mesh=None):
    '''Mesh by MeshEditor from vertices and cells'''
    if mesh is None:
        mesh = Mesh()
        assert mesh.mpi_comm().size == 1

    module.fill_mesh(coordinates.flatten(), cells.flatten(), tdim, gdim, mesh)
    
    return mesh

def readmesh(mesh_file):
    ''' read HDF5 or DOLFIN XML mesh '''
    # TODO: exceptions, files exist?
    from dolfin import Mesh, MeshFunction, CellFunction, HDF5File, \
        FacetFunction
    tmp = mesh_file.split('/')[-1].split('.')
    mesh_type = tmp[-1]
    mesh_name = '.'.join(tmp[0:-1])
    if mesh_type == 'xml':
        mesh = Mesh(mesh_file)
        rank = mesh.mpi_comm().Get_rank()
        try:
            subdomains = MeshFunction("size_t", mesh,
                                      mesh_name+"_physical_region.xml")
        except:
	    if rank == 0:
	      print('no subdomain file found (%s)' %
		    (mesh_name+"_physical_region.xml"))
            subdomains = CellFunction("size_t", mesh)
        try:
            boundaries = MeshFunction("size_t", mesh,
                                      mesh_name+"_facet_region.xml")
        except:
	    if rank == 0:
	      print('no boundary file found (%s)' %
		    (mesh_name+"_physical_region.xml"))
            boundaries = FacetFunction("size_t", mesh)
    elif mesh_type == 'h5':
        mesh = Mesh()
        rank = mesh.mpi_comm().Get_rank()
        hdf = HDF5File(mesh.mpi_comm(), mesh_file, "r")
        hdf.read(mesh, "/mesh", False)
        subdomains = CellFunction("size_t", mesh)
        boundaries = FacetFunction("size_t", mesh)
        if hdf.has_dataset('subdomains'):
            hdf.read(subdomains, "/subdomains")
        else:
	    if rank == 0:
	      print('no <subdomains> datasets found in file %s' % mesh_file)
        if hdf.has_dataset('boundaries'):
            hdf.read(boundaries, "/boundaries")
        else:
	    if rank == 0:
	      print('no <boundaries> datasets found in file %s' % mesh_file)

    elif mesh_type in ['xdmf', 'xmf']:
        import sys
        sys.exit('XDMF not supported yet. Use HDF5 instead!')
    else:
        import sys
        sys.exit('mesh format not recognized. try XML (serial) or HDF5')

# NOTE http://fenicsproject.org/qa/5337/importing-marked-mesh-for-parallel-use
    # see file xml2xdmf.py
    return mesh, subdomains, boundaries

def cyclic3D(u):
    """ Symmetrize with respect to (xyz) cycle. """
    try:
        nrm = np.linalg.norm(u.vector())
        V = u.function_space()
        assert V.mesh().topology().dim() == 3
        mesh1 = Mesh(V.mesh())
        mesh1.coordinates()[:, :] = mesh1.coordinates()[:, [1, 2, 0]]
        W1 = FunctionSpace(mesh1, 'CG', 1)

        # testing if symmetric
        bc = DirichletBC(V, 1, DomainBoundary())
        test = Function(V)
        bc.apply(test.vector())
        test = interpolate(Function(W1, test.vector()), V)
        assert max(test.vector()) - min(test.vector()) < 1.1

        v1 = interpolate(Function(W1, u.vector()), V)

        mesh2 = Mesh(mesh1)
        mesh2.coordinates()[:, :] = mesh2.coordinates()[:, [1, 2, 0]]
        W2 = FunctionSpace(mesh2, 'CG', 1)
        v2 = interpolate(Function(W2, u.vector()), V)
        pr = project(u+v1+v2)
        assert np.linalg.norm(pr.vector())/nrm > 0.01
        return pr
    except:
        print "Cyclic symmetrization failed!"
        return u

def polyhedron_mesh(data):
    """
    Build polyhedral mesh. Must be strlike with respect to the origin.

    Input:
        data[0] - list of vertices
        data[1] - list of faces
        data[2] - optional other starlike point, instead of the origin
    """
    # TODO: Center of mass of the vertices as origin
    vertex_data = np.array(data[0], dtype='double')
    lowest = np.min(flatten(data[1]))
    face_data = [list(np.array(d) - lowest) for d in data[1]]
    numv = len(vertex_data)  # will be the index of the origin
    if len(data) > 2:
        origin = np.array(data[2], dtype='double')
    else:
        origin = [0.0, 0.0, 0.0]
    mesh = Mesh()
    editor = DynamicMeshEditor()
    editor.open(mesh, "tetrahedron", 3, 3, numv + 1, len(face_data))
    for i, vert in enumerate(vertex_data):
        editor.add_vertex(i, *vert)
    editor.add_vertex(numv, *origin)
    newv = numv + 1  # next vertex index
    newf = 0  # next face index
    for face in face_data:
        if len(face) == 3:
            # triangular face, no splitting
            editor.add_cell(newf, numv, *face)  # face + origin
            newf += 1
        else:
            # split face into triangles using center of mass
            # average face vertices to get the center
            vert = list(np.mean(vertex_data[np.array(face)], axis=0))
            editor.add_vertex(newv, *vert)  # new vertex: face center
            face.append(face[0])
            for i in zip(face[:-1], face[1:]):
                # pairs of vertices
                editor.add_cell(newf, numv, newv, *i)  # + face center + origin
                newf += 1
            newv += 1
    editor.close()
    mesh.order()
    return mesh

def convert_and_create_facet_mesh_function ( ifilename, ofilename ):
    # First convert the gmsh mesh
    meshconvert.convert2xml ( ifilename, ofilename )
    
    # Now load the created mesh and initialise the required connectivity information
    mesh = Mesh ( ofilename )
    mesh.order()
    
    File ( ofilename ) << mesh
    
    D = mesh.topology().dim()
    mesh.init(D-1, 0)
    
    # read the data from the gmsh file once again
    dim_count, vertices_used, tags = process_gmsh_elements( ifilename, D-1 )
    # Get the facet-node connectivity information (reshape as a row of node indices per facet)
    facets_as_nodes = mesh.topology()(D-1,0)().reshape ( mesh.num_facets(), D )
    
    # Create and initialise the mesh function
    facet_mark_function = MeshFunction ( 'uint', mesh, D-1 )
    facet_mark_function.set_all( 0 )
    
    # set the relevant values of the mesh function
    facets_to_check = range( mesh.num_facets() )
    for i in range(len(tags)):
        nodes = np.sort(np.array(vertices_used[2*i:(2*i+D)]))
        value  = tags[i][0]
        
        if value != 0:
            found = False
            for j in range(len(facets_to_check)):
                index = facets_to_check[j]
                if np.array_equal(facets_as_nodes[index,:], nodes):
                    found = True;
                    facets_to_check.pop(j)
                    # set the value of the mesh function
                    facet_mark_function[index] = value
                    break;
                
            if not found:
                raise Exception ( "The facet (%d) was not found to mark: %s" % (i, nodes) )
        
    # save the mesh function to file
    fname = os.path.splitext(ofilename)[0]
    mesh_function_file = File("%s_%s.xml" % (fname, "facet_regions"))
    
    mesh_function_file << facet_mark_function

def write_fenics_file(dim, ofilename):
    ofile  = File(ofilename + '.xml')
    mesh = Mesh()
    editor = MeshEditor()
    editor.open(mesh, dim, dim)
    editor.init_vertices(nodes.shape[1])
    editor.init_cells(len(cell_map))    
    for i in range(nodes.shape[1]):
        if dim == 2:
            editor.add_vertex(i, nodes[0, i], nodes[1, i])
        else:
            editor.add_vertex(i, nodes[0, i], nodes[1, i], nodes[2, i])
            
    for i in range(1, len(cell_map)+1):
        if dim == 2:
            editor.add_cell(i-1, cell_map[i][0]-1, cell_map[i][1]-1, cell_map[i][2]-1)
        else:
            editor.add_cell(i-1, cell_map[i][0]-1, cell_map[i][1]-1, cell_map[i][2]-1, cell_map[i][3]-1)
    
    mesh.order()
    mvc = mesh.domains().markers(dim-1)
    for zone, cells in boundary_cells.iteritems():
        for cell, nds in cells.iteritems():
            dolfin_cell = Cell(mesh, cell-1)
            nodes_of_cell = dolfin_cell.entities(0)
            #print cell
            #print nodes_of_cell
            nodes_of_face = nds - 1
            #print nodes_of_face
            for jj, ff in enumerate(facets(dolfin_cell)):
                facet_nodes = ff.entities(0)
                #print facet_nodes
                if all(map(lambda x: x in nodes_of_face, facet_nodes)):
                    local_index = jj
                    break
            mvc.set_value(cell-1, local_index, zone)
        
    ofile << mesh        
    from dolfin import plot
    plot(mesh, interactive=True)
    print 'Finished writing FEniCS mesh\n'

def rotational(u, n):
    """ Symmetrize with respect to n-fold symmetry. """
    # TODO: test one rotation only
    V = u.function_space()
    if V.mesh().topology().dim() > 2 or n < 2:
        return u
    mesh = V.mesh()
    sum = u
    nrm = np.linalg.norm(u.vector())
    rotation = np.array([[np.cos(2*np.pi/n), np.sin(2*np.pi/n)],
                         [-np.sin(2*np.pi/n), np.cos(2*np.pi/n)]])
    for i in range(1, n):
        mesh = Mesh(mesh)
        mesh.coordinates()[:, :] = np.dot(mesh.coordinates(), rotation)
        W = FunctionSpace(mesh, 'CG', 1)
        v = interpolate(Function(W, u.vector()), V)
        sum += v
    pr = project(sum)
    if np.linalg.norm(pr.vector())/nrm > 0.01:
        return pr
    else:
        return u

def symmetrize(u, d, sym):
    """ Symmetrize function u. """
    if len(d) == 3:
        # three dimensions -> cycle XYZ
        return cyclic3D(u)
    elif len(d) >= 4:
        # four dimensions -> rotations in 2D
        return rotational(u, d[-1])
    nrm = np.linalg.norm(u.vector())
    V = u.function_space()
    mesh = Mesh(V.mesh())

    # test if domain is symmetric using function equal 0 inside, 1 on boundary
    # extrapolation will force large values if not symmetric since the flipped
    # domain is different
    bc = DirichletBC(V, 1, DomainBoundary())
    test = Function(V)
    bc.apply(test.vector())

    if len(d) == 2:
        # two dimensions given: swap dimensions
        mesh.coordinates()[:, d] = mesh.coordinates()[:, d[::-1]]
    else:
        # one dimension given: reflect
        mesh.coordinates()[:, d[0]] *= -1
    # FIXME functionspace takes a long time to construct, maybe copy?
    W = FunctionSpace(mesh, 'CG', 1)
    try:
        # testing
        test = interpolate(Function(W, test.vector()), V)
        # max-min should be around 1 if domain was symmetric
        # may be slightly above due to boundary approximation
        assert max(test.vector()) - min(test.vector()) < 1.1

        v = interpolate(Function(W, u.vector()), V)
        if sym:
            # symmetric
            pr = project(u+v)
        else:
            # antisymmetric
            pr = project(u-v)
        # small solution norm most likely means that symmetrization gives
        # trivial function
        assert np.linalg.norm(pr.vector())/nrm > 0.01
        return pr
    except:
        # symmetrization failed for some reason
        print "Symmetrization " + str(d) + " failed!"
        return u

    def setUp(self, *args, **kwargs):
        self.mesh = Mesh()
        editor = MeshEditor()
        editor.open(self.mesh, 2, 2) # topo_dim = 2, geom dim = 2

        editor.init_vertices(6)
        editor.init_cells(2)

        vertex_0 = Vertex(self.mesh, 0)
        vertex_1 = Vertex(self.mesh, 1)
        vertex_2 = Vertex(self.mesh, 2)
        vertex_3 = Vertex(self.mesh, 3)

        vertex_4 = Vertex(self.mesh, 4)
        vertex_5 = Vertex(self.mesh, 5)

        editor.add_cell(0,1,2,3)
        editor.add_cell(1,0,2,3)

        editor.close()

class TestIonTag(unittest.TestCase):

    def setUp(self, *args, **kwargs):
        self.mesh = Mesh()
        editor = MeshEditor()
        editor.open(self.mesh, 2, 2) # topo_dim = 2, geom dim = 2

        editor.init_vertices(6)
        editor.init_cells(2)

        vertex_0 = Vertex(self.mesh, 0)
        vertex_1 = Vertex(self.mesh, 1)
        vertex_2 = Vertex(self.mesh, 2)
        vertex_3 = Vertex(self.mesh, 3)

        vertex_4 = Vertex(self.mesh, 4)
        vertex_5 = Vertex(self.mesh, 5)

        editor.add_cell(0,1,2,3)
        editor.add_cell(1,0,2,3)

        editor.close()


    def test_init(self):
        # test possible arguments and failure cases...
        #@todo implement the test once we have exception handling in the init method

        # pass in a string for values and see what happens
        values = 'i am just a string'

        t = IonTag('foo',3,'int', self.mesh)
        v = MeshEntity(self.mesh,0,1)

        try:
            t[v] = values
        except ValueError:
            pass
        else:
            raise AssertionError('It should have raised a value error!')


    def test_get_set_del(self):
        #Test the getter, setter and delete method

        values = [1,2,3]

        t = IonTag('foo',3,'int', self.mesh)

        for v in vertices(self.mesh):
            # test the setter
            t[v] = values

        # choose an entity in the mesh
        v = MeshEntity(self.mesh,0,1)

        # test the getter
        self.assertTrue((t[v] == values).all())

        #---------------------------------------------------------------------------------------
        # Check delete of a tag entry (for an entity)
        #---------------------------------------------------------------------------------------

        # choose an entity to delete
        entity_tuple = (v.dim(),v.index())

        # check that tag has the entity, v, in it
        self.assertTrue(t._entity_values.has_key(entity_tuple))

        # delete a tag entry for an entity
        del t[entity_tuple]

        # check that the tag no longer has the entity, v, in it
        self.assertFalse(t._entity_values.has_key(entity_tuple))

        #---------------------------------------------------------------------------------------
        # Add less number of values than the size defined in the tag object
        #---------------------------------------------------------------------------------------

        values = [1]
        t = IonTag('foo',3,'int', self.mesh)
        v = MeshEntity(self.mesh,0,1)

        #@todo check to see why self.assertRaises is not working for unittest:

#        with self.assertRaises(ValueError):
#            t[v] = values

        try:
            t[v] = values
        except ValueError:
            pass
        else:
            raise AssertionError('A Value Error should have been raised!')

        #---------------------------------------------------------------------------------------
        # Add more number of values that the size defined in the tag object
        #---------------------------------------------------------------------------------------

        values = [1,2,3,4]
#.........这里部分代码省略.........

 def get_greenland_detailed():
   filename = inspect.getframeinfo(inspect.currentframe()).filename
   home     = os.path.dirname(os.path.abspath(filename))
   mesh     = Mesh(home + '/greenland/greenland_detailed_mesh.xml')
   mesh.coordinates()[:,2] /= 100000.0
   return mesh