Python Function.eval_cell方法代码示例

# 需要导入模块: from function import Function [as 别名]
# 或者: from function.Function import eval_cell [as 别名]
def solve(n_cells, degree=3, with_plot=False):
    # Problem
    x = Symbol('x')
    vvvv = -0.0625*x**3/pi**3 + 0.0625*x/pi**3 + sin(2*pi*x)/(16*pi**4)
    u = -vvvv.diff(x, 2)
    f = sin(2*pi*x)

    # As Expr
    u = Expression(u)
    f = Expression(f)
    vvvv = Expression(vvvv)

    # Space
    element = HermiteElement(degree)

    mesh = IntervalMesh(a=-1, b=1, n_cells=n_cells)
    V = FunctionSpace(mesh, element)
    bc = DirichletBC(V, vvvv)

    # Need mass matrix to intefrate the rhs
    M = assemble_matrix(V, 'mass', get_geom_tensor=None, timer=0)
    # NOTE We cannot you apply the alpha transform idea because the functions
    # are mapped with this selective weight on 2nd, 3rd functions. So some rows
    # of alpha would have to be multiplied by weights which are cell specific.
    # And then on top of this there would be a dx = J*dy term. Better just to
    # use the qudrature representations
    # Mpoly_matrix = leg.mass_matrix(degree)
    # M_ = assemble_matrix(V, Mpoly_matrix, Mget_geom_tensor, timer=0)
    
    # Stiffness matrix for Laplacian
    A = assemble_matrix(V, 'bending', get_geom_tensor=None, timer=0)
    # NOTE the above
    # Apoly_matrix = leg.stiffness_matrix(degree)
    # A_ = assemble_matrix(V, Apoly_matrix, Aget_geom_tensor, timer=0)
  
    # Interpolant of source
    fV = V.interpolate(f)
    # Integrate in L2 to get the vector
    b = M.dot(fV.vector)
    
    # Apply boundary conditions
    bc.apply(A, b, True)
    x = spsolve(A, b)

    print '>>>>', np.linalg.norm(x - V.interpolate(vvvv).vector)

    ALaplace = assemble_matrix(V, 'stiffness', get_geom_tensor=None, timer=0)

    c = ALaplace.dot(x)
    y = spsolve(M, c)
    uh = Function(V, y)

    # This is a (slow) way of plotting the high order
    if with_plot:
        fig = plt.figure()
        ax = fig.gca()
        uV = V.interpolate(u)
        
        for cell in Cells(mesh):
            a, b = cell.vertices[0, 0], cell.vertices[1, 0]
            x = np.linspace(a, b, 100)

            y = uh.eval_cell(x, cell)
            ax.plot(x, y, color=random.choice(['b', 'g', 'm', 'c']))
            
            y = uV.eval_cell(x, cell)
            ax.plot(x, y, color='r')

            y = u.eval_cell(x, cell)
            ax.plot(x, y, color='k')

        plt.show()

    # Error norm
    fine_degree = degree + 3
    element = HermiteElement(fine_degree)
    
    V_fine = FunctionSpace(mesh, element)
    # Interpolate exact solution to fine
    u_fine = V_fine.interpolate(u)
    # Interpolate approx solution fine
    uh_fine = V_fine.interpolate(uh)

    # Difference vector
    e = u_fine.vector - uh_fine.vector

    # Integrate the error
    A_fine = assemble_matrix(V_fine, 'mass', get_geom_tensor=None, timer=0)
    e = sqrt(np.sum(e*A_fine.dot(e)))
    # Mesh size
    hmin = mesh.hmin()

    # Add the cond number
    kappa = np.linalg.cond(A.toarray())

    return hmin, e, kappa, A.shape[0]

# 需要导入模块: from function import Function [as 别名]
# 或者: from function.Function import eval_cell [as 别名]
def solve(n_cells, degree=3, with_plot=False):
    # Problem
    w = 3 * np.pi
    x = Symbol("x")
    u = sin(w * x)
    f = -u.diff(x, 2)

    # As Expr
    u = Expression(u)
    f = Expression(f)

    # Space
    # element = HermiteElement(degree)
    poly_set = leg.basis_functions(degree)
    dof_set = chebyshev_points(degree)
    element = LagrangeElement(poly_set, dof_set)

    mesh = IntervalMesh(a=-1, b=1, n_cells=n_cells)
    V = FunctionSpace(mesh, element)
    bc = DirichletBC(V, u)

    # Need mass matrix to intefrate the rhs
    M = assemble_matrix(V, "mass", get_geom_tensor=None, timer=0)
    # NOTE We cannot you apply the alpha transform idea because the functions
    # are mapped with this selective weight on 2nd, 3rd functions. So some rows
    # of alpha would have to be multiplied by weights which are cell specific.
    # And then on top of this there would be a dx = J*dy term. Better just to
    # use the qudrature representations
    # Mpoly_matrix = leg.mass_matrix(degree)
    # M_ = assemble_matrix(V, Mpoly_matrix, Mget_geom_tensor, timer=0)

    # Stiffness matrix for Laplacian
    A = assemble_matrix(V, "stiffness", get_geom_tensor=None, timer=0)
    # NOTE the above
    # Apoly_matrix = leg.stiffness_matrix(degree)
    # A_ = assemble_matrix(V, Apoly_matrix, Aget_geom_tensor, timer=0)

    # Interpolant of source
    fV = V.interpolate(f)
    # Integrate in L2 to get the vector
    b = M.dot(fV.vector)

    # Apply boundary conditions
    bc.apply(A, b, True)
    x = spsolve(A, b)

    # As function
    uh = Function(V, x)

    # This is a (slow) way of plotting the high order
    if with_plot:
        fig = plt.figure()
        ax = fig.gca()
        uV = V.interpolate(u)

        for cell in Cells(mesh):
            a, b = cell.vertices[0, 0], cell.vertices[1, 0]
            x = np.linspace(a, b, 100)

            y = uh.eval_cell(x, cell)
            ax.plot(x, y, color=random.choice(["b", "g", "m", "c"]))

            y = uV.eval_cell(x, cell)
            ax.plot(x, y, color="r")

            y = u.eval_cell(x, cell)
            ax.plot(x, y, color="k")

        plt.show()

    # Error norm in CG high order
    fine_degree = degree + 3
    poly_set = leg.basis_functions(fine_degree)
    dof_set = chebyshev_points(fine_degree)
    element = LagrangeElement(poly_set, dof_set)

    V_fine = FunctionSpace(mesh, element)
    # Interpolate exact solution to fine
    u_fine = V_fine.interpolate(u)
    # Interpolate approx solution fine
    uh_fine = V_fine.interpolate(uh)

    # Difference vector
    e = u_fine.vector - uh_fine.vector

    # L2
    if False:
        Apoly_matrix = leg.mass_matrix(fine_degree)
        get_geom_tensor = lambda cell: 1.0 / cell.Jac

    # Need matrix for integration of H10 norm
    else:
        Apoly_matrix = leg.stiffness_matrix(fine_degree)
        get_geom_tensor = lambda cell: cell.Jac

    A_fine = assemble_matrix(V_fine, Apoly_matrix, get_geom_tensor, timer=0)

    # Integrate the error
    e = sqrt(np.sum(e * A_fine.dot(e)))
    # Mesh size
#.........这里部分代码省略.........