Python Mesh.init方法代码示例

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import init [as 别名]
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

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import init [as 别名]
# Rossby radius.
LR=c/params["f"]

class InitialConditions(Expression):
    def __init__(self):
        pass
    def eval(self, values, X):
        r=(X[0]**2+X[1]**2)**0.5
        if r>0.0001:
            values[0]=-0.05*c*exp((r-r0)/LR)*X[0]/r*X[1]/r
            values[1]= 0.05*c*exp((r-r0)/LR)*X[0]/r*X[0]/r
            values[2]= 0.05*exp((r-r0)/LR)*X[0]/r
        else:
            values[0]=0.
            values[1]=0.
            values[2]=0.
    def value_shape(self):
        return (3,)

try:
    mesh=Mesh("basin.xml")
except RuntimeError:
    import sys
    import os.path

    mesh=Mesh(os.path.dirname(sys.argv[0]) + os.path.sep + "basin.xml")

mesh.order()
mesh.init()

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import init [as 别名]

#.........这里部分代码省略.........
        # TODO: Output a warning or an error here
        me = None

    while state != 10:

        # Read next line
        line = ifile.readline()
        if not line: break

        # Skip comments
        if line[0] == '#':
            continue

        # Remove newline
        if line[-1] == "\n":
            line = line[:-1]

        if state == 0:
            if line == "$MeshFormat":
                state = 1
        elif state == 1:
            (version, file_type, data_size) = line.split()
            state = 2
        elif state == 2:
            if line == "$EndMeshFormat":
                state = 3
        elif state == 3:
            if line == "$Nodes":
                state = 4
        elif state == 4:
            num_vertices = len(vertex_dict)
            handler.start_vertices(num_vertices)
            if process_facets:
                mesh_editor.init_vertices ( num_vertices )
            state = 5
        elif state == 5:
            (node_no, x, y, z) = line.split()
            node_no = int(node_no)
            x,y,z = [float(xx) for xx in (x,y,z)]
            if vertex_dict.has_key(node_no):
                node_no = vertex_dict[node_no]
            else:
                continue
            nodelist[int(node_no)] = num_vertices_read
            handler.add_vertex(num_vertices_read, [x, y, z])
            if process_facets:
                if highest_dim == 1:
                    coords = numpy.array([x])
                elif highest_dim == 2:
                    coords = numpy.array([x, y])
                elif highest_dim == 3:
                    coords = numpy.array([x, y, z])
                mesh_editor.add_vertex(num_vertices_read, coords)

            num_vertices_read +=1

            if num_vertices == num_vertices_read:
                handler.end_vertices()
                state = 6
        elif state == 6:
            if line == "$EndNodes":
                state = 7
        elif state == 7:
            if line == "$Elements":
                state = 8
        elif state == 8:

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import init [as 别名]
    def test_convert_triangle(self): # Disabled because it fails, see FIXME below
        # test no. 1
        from dolfin import Mesh, MPI
        if MPI.num_processes() != 1:
            return
        fname = os.path.join("data", "triangle")
        dfname = fname+".xml"
        
        # Read triangle file and convert to a dolfin xml mesh file
        meshconvert.triangle2xml(fname, dfname)

        # Read in dolfin mesh and check number of cells and vertices
        mesh = Mesh(dfname)
        self.assertEqual(mesh.num_vertices(), 96)
        self.assertEqual(mesh.num_cells(), 159)

        # Clean up
        os.unlink(dfname)


        # test no. 2
        from dolfin import MPI, Mesh, MeshFunction, \
                           edges, Edge, faces, Face, \
                           SubsetIterator, facets, CellFunction
        if MPI.num_processes() != 1:
            return
        fname = os.path.join("data", "test_Triangle_3")
        dfname = fname+".xml"
        dfname0 = fname+".attr0.xml"

        # Read triangle file and convert to a dolfin xml mesh file
        meshconvert.triangle2xml(fname, dfname)

        # Read in dolfin mesh and check number of cells and vertices
        mesh = Mesh(dfname)
        mesh.init()
        mfun = MeshFunction('double', mesh, dfname0)
        self.assertEqual(mesh.num_vertices(), 58)
        self.assertEqual(mesh.num_cells(), 58)

        # Create a size_t CellFunction and assign the values based on the
        # converted Meshfunction
        cf = CellFunction("size_t", mesh)
        cf.array()[mfun.array()==10.0] = 0
        cf.array()[mfun.array()==-10.0] = 1

        # Meassure total area of cells with 1 and 2 marker
        add = lambda x, y : x+y
        area0 = reduce(add, (Face(mesh, cell.index()).area() \
                             for cell in SubsetIterator(cf, 0)), 0.0)
        area1 = reduce(add, (Face(mesh, cell.index()).area() \
                             for cell in SubsetIterator(cf, 1)), 0.0)
        total_area = reduce(add, (face.area() for face in faces(mesh)), 0.0)

        # Check that all cells in the two domains are either above or below y=0
        self.assertTrue(all(cell.midpoint().y()<0 for cell in SubsetIterator(cf, 0)))
        self.assertTrue(all(cell.midpoint().y()>0 for cell in SubsetIterator(cf, 1)))
        
        # Check that the areas add up
        self.assertAlmostEqual(area0+area1, total_area)
        
        # Measure the edge length of the two edge domains
        edge_markers = mesh.domains().facet_domains()
        self.assertTrue(edge_markers is not None)
        length0 = reduce(add, (Edge(mesh, e.index()).length() \
                            for e in SubsetIterator(edge_markers, 0)), 0.0)
        length1 = reduce(add, (Edge(mesh, e.index()).length() \
                            for e in SubsetIterator(edge_markers, 1)), 0.0)
        
        # Total length of all edges and total length of boundary edges
        total_length = reduce(add, (e.length() for e in edges(mesh)), 0.0)
        boundary_length = reduce(add, (Edge(mesh, f.index()).length() \
                          for f in facets(mesh) if f.exterior()), 0.0)
        
        # Check that the edges add up
        self.assertAlmostEqual(length0+length1, total_length)
        self.assertAlmostEqual(length1, boundary_length)

        # Clean up
        os.unlink(dfname)
        os.unlink(dfname0)