Python LinSystem.set_pss方法代码示例

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_matrix():
    set_verbose(False)

    mesh = Mesh()
    mesh.load(domain_mesh)
    mesh.refine_element(0)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)
    space.assign_dofs()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

    # assemble the stiffness matrix and solve the system
    sln = Solution()
    sys.assemble()
    A = sys.get_matrix()

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_ScalarView_mpl_unknown():
    mesh = Mesh()
    mesh.load(domain_mesh)
    mesh.refine_element(0)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)
    space.assign_dofs()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

    # assemble the stiffness matrix and solve the system
    sys.assemble()
    A = sys.get_matrix()
    b = sys.get_rhs()
    from scipy.sparse.linalg import cg

    x, res = cg(A, b)
    sln = Solution()
    sln.set_fe_solution(space, pss, x)

    view = ScalarView("Solution")

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_03():
    set_verbose(False)

    mesh = Mesh()
    mesh.load(domain_mesh)
    mesh.refine_element(0)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)
    from hermes2d.examples.c03 import set_bc

    set_bc(space)
    space.assign_dofs()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

    # assemble the stiffness matrix and solve the system
    sln = Solution()
    sys.assemble()
    sys.solve_system(sln)
    assert abs(sln.l2_norm() - 0.25493) < 1e-4
    assert abs(sln.h1_norm() - 0.89534) < 1e-4

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_08():
    from hermes2d.examples.c08 import set_bc, set_forms

    set_verbose(False)

    mesh = Mesh()
    mesh.load(cylinder_mesh)
    #mesh.refine_element(0)
    #mesh.refine_all_elements()
    mesh.refine_towards_boundary(5, 3)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    xvel = H1Space(mesh, shapeset)
    yvel = H1Space(mesh, shapeset)
    press = H1Space(mesh, shapeset)
    xvel.set_uniform_order(2)
    yvel.set_uniform_order(2)
    press.set_uniform_order(1)

    set_bc(xvel, yvel, press)

    ndofs = 0
    ndofs += xvel.assign_dofs(ndofs)
    ndofs += yvel.assign_dofs(ndofs)
    ndofs += press.assign_dofs(ndofs)

    xprev = Solution()
    yprev = Solution()

    xprev.set_zero(mesh)
    yprev.set_zero(mesh)

    # initialize the discrete problem
    wf = WeakForm(3)
    set_forms(wf, xprev, yprev)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(xvel, yvel, press)
    sys.set_pss(pss)
    #dp.set_external_fns(xprev, yprev)

    # visualize the solution

    EPS_LOW = 0.0014

    for i in range(3):
        psln = Solution()
        sys.assemble()
        sys.solve_system(xprev, yprev, psln)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_07():
    from hermes2d.examples.c07 import set_bc, set_forms

    set_verbose(False)

    mesh = Mesh()
    mesh.load(sample_mesh)
    #mesh.refine_element(0)
    #mesh.refine_all_elements()
    #mesh.refine_towards_boundary(5, 3)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    xdisp = H1Space(mesh, shapeset)
    ydisp = H1Space(mesh, shapeset)
    xdisp.set_uniform_order(8)
    ydisp.set_uniform_order(8)

    set_bc(xdisp, ydisp)

    ndofs = xdisp.assign_dofs(0)
    ndofs += ydisp.assign_dofs(ndofs)

    # initialize the discrete problem
    wf = WeakForm(2)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(xdisp, ydisp)
    sys.set_pss(pss)

    xsln = Solution()
    ysln = Solution()
    old_flag = set_warn_integration(False)
    sys.assemble()
    set_warn_integration(old_flag)
    sys.solve_system(xsln, ysln)

    E = float(200e9)
    nu = 0.3
    stress = VonMisesFilter(xsln, ysln, E / (2*(1 + nu)),
            (E * nu) / ((1 + nu) * (1 - 2*nu)))

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_02():
    set_verbose(False)

    mesh = Mesh()
    mesh.load(domain_mesh)
    mesh.refine_element(0)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)
    space.assign_dofs()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_04():
    from hermes2d.examples.c04 import set_bc

    set_verbose(False)

    mesh = Mesh()
    mesh.load(domain_mesh)
    # mesh.refine_element(0)
    # mesh.refine_all_elements()
    mesh.refine_towards_boundary(5, 3)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)

    set_bc(space)

    space.assign_dofs()

    xprev = Solution()
    yprev = Solution()

    # initialize the discrete problem
    wf = WeakForm()
    set_forms(wf, -4)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

    # assemble the stiffness matrix and solve the system
    sys.assemble()
    sln = Solution()
    sys.solve_system(sln)
    assert abs(sln.l2_norm() - 1.22729) < 1e-4
    assert abs(sln.h1_norm() - 2.90006) < 1e-4

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_05():
    from hermes2d.examples.c05 import set_bc
    from hermes2d.examples.c05 import set_forms as set_forms_surf

    set_verbose(False)

    mesh = Mesh()
    mesh.load(domain_mesh)
    mesh.refine_towards_vertex(3, 12)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(4)

    set_bc(space)

    space.assign_dofs()

    xprev = Solution()
    yprev = Solution()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf, -1)
    set_forms_surf(wf)

    sln = Solution()
    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)
    sys.assemble()
    sys.solve_system(sln)
    assert abs(sln.l2_norm() - 0.535833) < 1e-4
    assert abs(sln.h1_norm() - 1.332908) < 1e-4

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_07():
    from hermes2d.examples.c07 import set_bc, set_forms

    set_verbose(False)

    P_INIT = 2  # Initial polynomial degree of all mesh elements.

    mesh = Mesh()
    mesh.load(get_07_mesh())

    # Initialize the shapeset and the cache
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create finite element space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(P_INIT)
    set_bc(space)

    # Enumerate basis functions
    space.assign_dofs()

    # weak formulation
    wf = WeakForm(1)
    set_forms(wf)

    # matrix solver
    solver = DummySolver()

    # Solve the problem
    sln = Solution()
    ls = LinSystem(wf, solver)
    ls.set_spaces(space)
    ls.set_pss(pss)
    ls.assemble()
    ls.solve_system(sln)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def poisson_solver(mesh_tuple):
    """
    Poisson solver.

    mesh_tuple ... a tuple of (nodes, elements, boundary, nurbs)
    """
    set_verbose(False)
    mesh = Mesh()
    mesh.create(*mesh_tuple)
    mesh.refine_element(0)
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # create an H1 space
    space = H1Space(mesh, shapeset)
    space.set_uniform_order(5)
    space.assign_dofs()

    # initialize the discrete problem
    wf = WeakForm(1)
    set_forms(wf)

    solver = DummySolver()
    sys = LinSystem(wf, solver)
    sys.set_spaces(space)
    sys.set_pss(pss)

    # assemble the stiffness matrix and solve the system
    sys.assemble()
    A = sys.get_matrix()
    b = sys.get_rhs()
    from scipy.sparse.linalg import cg
    x, res = cg(A, b)
    sln = Solution()
    sln.set_fe_solution(space, pss, x)
    return sln

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]

#.........这里部分代码省略.........
    ISO_ONLY = False  # Isotropic refinement flag (concerns quadrilateral elements only).
    # ISO_ONLY = false ... anisotropic refinement of quad elements
    # is allowed (default),
    # ISO_ONLY = true ... only isotropic refinements of quad elements
    # are allowed.
    MESH_REGULARITY = -1  # Maximum allowed level of hanging nodes:
    # MESH_REGULARITY = -1 ... arbitrary level hangning nodes (default),
    # MESH_REGULARITY = 1 ... at most one-level hanging nodes,
    # MESH_REGULARITY = 2 ... at most two-level hanging nodes, etc.
    # Note that regular meshes are not supported, this is due to
    # their notoriously bad performance.
    MAX_ORDER = 10  # Maximum allowed element degree
    ERR_STOP = 0.5  # Stopping criterion for adaptivity (rel. error tolerance between the
    # fine mesh and coarse mesh solution in percent).
    NDOF_STOP = 40000  # Adaptivity process stops when the number of degrees of freedom grows over
    # this limit. This is mainly to prevent h-adaptivity to go on forever.

    # Problem constants
    E = 200e9  # Young modulus for steel: 200 GPa
    nu = 0.3  # Poisson ratio
    lamda = (E * nu) / ((1 + nu) * (1 - 2 * nu))
    mu = E / (2 * (1 + nu))

    # Load the mesh
    xmesh = Mesh()
    ymesh = Mesh()
    xmesh.load(get_bracket_mesh())

    # initial mesh refinements
    xmesh.refine_element(1)
    xmesh.refine_element(4)

    # Create initial mesh for the vertical displacement component,
    # identical to the mesh for the horizontal displacement
    # (bracket.mesh becomes a master mesh)
    ymesh.copy(xmesh)

    # Initialize the shapeset and the cache
    shapeset = H1Shapeset()
    xpss = PrecalcShapeset(shapeset)
    ypss = PrecalcShapeset(shapeset)

    # Create the x displacement space
    xdisp = H1Space(xmesh, shapeset)
    set_bc(xdisp)
    xdisp.set_uniform_order(P_INIT)

    # Create the x displacement space
    ydisp = H1Space(ymesh, shapeset)
    set_bc(ydisp)
    ydisp.set_uniform_order(P_INIT)

    # Enumerate basis functions
    ndofs = xdisp.assign_dofs()
    ydisp.assign_dofs(ndofs)

    # Initialize the weak formulation
    wf = WeakForm(2)
    set_wf_forms(wf)

    # Matrix solver
    solver = DummySolver()

    # adaptivity loop
    it = 1
    done = False
    cpu = 0.0

    x_sln_coarse = Solution()
    y_sln_coarse = Solution()

    x_sln_fine = Solution()
    y_sln_fine = Solution()

    # Calculating the number of degrees of freedom
    ndofs = xdisp.assign_dofs()
    ndofs += ydisp.assign_dofs(ndofs)

    # Solve the coarse mesh problem
    ls = LinSystem(wf, solver)
    ls.set_spaces(xdisp, ydisp)
    ls.set_pss(xpss, ypss)
    ls.assemble()
    ls.solve_system(x_sln_coarse, y_sln_coarse)

    # View the solution -- this can be slow; for illustration only
    stress_coarse = VonMisesFilter(x_sln_coarse, y_sln_coarse, mu, lamda)

    # Solve the fine mesh problem
    rs = RefSystem(ls)
    rs.assemble()
    rs.solve_system(x_sln_fine, y_sln_fine)

    # Calculate element errors and total error estimate
    hp = H1OrthoHP(xdisp, ydisp)
    set_hp_forms(hp)
    err_est = hp.calc_error_2(x_sln_coarse, y_sln_coarse, x_sln_fine, y_sln_fine) * 100

    # Show the fine solution - this is the final result
    stress_fine = VonMisesFilter(x_sln_fine, y_sln_fine, mu, lamda)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
def test_example_09():
    from hermes2d.examples.c09 import set_bc, temp_ext, set_forms

    # The following parameters can be played with:
    P_INIT = 1  # polynomial degree of elements
    INIT_REF_NUM = 4  # number of initial uniform refinements
    TAU = 300.0  # time step in seconds

    # Problem constants
    T_INIT = 10  # temperature of the ground (also initial temperature)
    FINAL_TIME = 86400  # length of time interval (24 hours) in seconds

    # Global variable
    TIME = 0

    # Load the mesh
    mesh = Mesh()
    mesh.load(get_cathedral_mesh())

    # for i in range(INIT_REF_NUM):
    #    mesh.refine_all_elements()
    # mesh.refine_towards_boundary(2, 5)

    # Set up shapeset
    shapeset = H1Shapeset()
    pss = PrecalcShapeset(shapeset)

    # Set up spaces
    space = H1Space(mesh, shapeset)
    set_bc(space)
    space.set_uniform_order(P_INIT)

    # Enumerate basis functions
    space.assign_dofs()

    # Set initial condition
    tsln = Solution()
    tsln.set_const(mesh, T_INIT)

    # Weak formulation
    wf = WeakForm(1)
    set_forms(wf, tsln)

    # Matrix solver
    solver = DummySolver()

    # Linear system
    ls = LinSystem(wf, solver)
    ls.set_spaces(space)
    ls.set_pss(pss)

    # Visualisation
    sview = ScalarView("Temperature", 0, 0, 450, 600)
    # title = "Time %s, exterior temperature %s" % (TIME, temp_ext(TIME))
    # Tview.set_min_max_range(0,20);
    # Tview.set_title(title);
    # Tview.fix_scale_width(3);

    # Time stepping
    nsteps = int(FINAL_TIME / TAU + 0.5)
    rhsonly = False

    # Assemble and solve
    ls.assemble()
    rhsonly = True
    ls.solve_system(tsln)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
yomview = MeshView()
while(not done):

    print ("\n---- Adaptivity step %d ---------------------------------------------\n" % it)
    it += 1

    # Calculating the number of degrees of freedom
    ndofs = xdisp.assign_dofs()
    ndofs += ydisp.assign_dofs(ndofs)

    print("xdof=%d, ydof=%d\n" % (xdisp.get_num_dofs(), ydisp.get_num_dofs()) )

    # Solve the coarse mesh problem
    ls = LinSystem(wf, solver)
    ls.set_spaces(xdisp, ydisp)
    ls.set_pss(xpss, ypss)
    ls.assemble()
    ls.solve_system(x_sln_coarse, y_sln_coarse, lib="scipy")

    # View the solution -- this can be slow; for illustration only
    stress_coarse = VonMisesFilter(x_sln_coarse, y_sln_coarse, mu, lamda)
    #sview.set_min_max_range(0, 3e4)
    sview.show(stress_coarse, lib='mayavi')
    #xoview.show(xdisp, lib='mayavi')
    #yoview.show(ydisp, lib='mayavi')
    xomview.show(xmesh, space=xdisp, lib="mpl", method="orders", notebook=False)
    yomview.show(ymesh, space=ydisp, lib="mpl", method="orders", notebook=False)

    # Solve the fine mesh problem
    rs = RefSystem(ls)
    rs.assemble()

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]
done = False
cpu = 0.0

sln_coarse = Solution()
sln_fine = Solution()

while(not done):

    print("\n---- Adaptivity step %d ---------------------------------------------\n" % it)
    it += 1

    # Solve the coarse mesh problem
    ls = LinSystem(wf, solver)
    ls.set_spaces(space)
    ls.set_pss(pss)
    ls.assemble()
    ls.solve_system(sln_coarse)

    # View the solution
    sview.show(sln_coarse, lib='mayavi')

    # View the mesh
    mview = MeshView("Example 7", 100, 100, 500, 500)
    mview.show(mesh, lib="mpl", method="orders", notebook=False)

    # Solve the fine mesh problem
    rs = RefSystem(ls)
    rs.assemble()
    rs.solve_system(sln_fine)

# 需要导入模块: from hermes2d import LinSystem [as 别名]
# 或者: from hermes2d.LinSystem import set_pss [as 别名]

#.........这里部分代码省略.........
    wf2 = WeakForm(1)
    # this is induced by set_verbose():
    #dp2.set_quiet(not verbose)
    set_forms7(wf2)

    solver = DummySolver()

    w = 320
    h = 320
    views = [ScalarView("", i*w, 0, w, h) for i in range(4)]
    viewsm = [ScalarView("", i*w, h, w, h) for i in range(4)]
    viewse = [ScalarView("", i*w, 2*h, w, h) for i in range(4)]
    #for v in viewse:
    #    v.set_min_max_range(0, 10**-4)
    ord = OrderView("Polynomial Orders", 0, 2*h, w, h)

    rs = None

    precision = 30.0

    if verbose_level >= 1:
        print "Problem initialized. Starting calculation."

    for it in range(iter):
        if verbose_level >= 1:
            print "-"*80
            print "Starting iteration %d." % it
        if report:
            iteration["n"] = it

        #mesh.save("refined2.mesh")
        sys1 = LinSystem(wf1, solver)
        sys1.set_spaces(space)
        sys1.set_pss(pss)
        sys2 = LinSystem(wf2, solver)
        sys2.set_spaces(space)
        sys2.set_pss(pss)

        if verbose_level >= 1:
            print "Assembling the matrices A, B."
        sys1.assemble()
        sys2.assemble()
        if verbose_level == 2:
            print "converting matrices A, B"
        A = sys1.get_matrix()
        B = sys2.get_matrix()
        if verbose_level >= 1:
            n = A.shape[0]
            print "Solving the problem Ax=EBx  (%d x %d)." % (n, n)
        if report:
            n = A.shape[0]
            iteration["DOF"] = n
        if report:
            t = clock()
        eigs, sols = solve(A, B, n_eigs, verbose_level == 2)
        if report:
            t = clock() - t
            iteration["cpu_solve"] = t
            iteration["eigenvalues"] = array(eigs)
            #h5eigs.append(sols)
        if verbose_level >= 1:
            print "   \-Done."
            print_eigs(eigs, E_exact)
        s = []

        n = sols.shape[1]