Python Mesh.coordinates方法代码示例

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

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

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

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

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

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import coordinates [as 别名]
 def get_circle():
   filename = inspect.getframeinfo(inspect.currentframe()).filename
   home     = os.path.dirname(os.path.abspath(filename))
   mesh     = Mesh(home + '/test/circle.xml')
   mesh.coordinates()[:,2] /= 1000.0
   return mesh

# 需要导入模块: from dolfin import Mesh [as 别名]
# 或者: from dolfin.Mesh import coordinates [as 别名]
 def get_antarctica_coarse():
   filename = inspect.getframeinfo(inspect.currentframe()).filename
   home     = os.path.dirname(os.path.abspath(filename))
   mesh     = Mesh(home + '/antarctica/antarctica_50H_5l.xml')
   mesh.coordinates()[:,2] /= 1000.0
   return mesh