Python IterOperations类代码示例

def foo():
    m = 6
    mm = 5


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,3*(mm-1)))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))
    KappaSave = np.zeros((mm-1,1))
    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'
    qq = -1
    MU[0]= 1e0
    kappa = 0.01
    qq = -1
    for yy in xrange(1,mm):
        kappa = kappa*10
        KappaSave[yy-1] = kappa
        IterTypes = ['Full','MD','CD']
        for kk in range(len(IterTypes)):
            qq += 1
            for xx in xrange(1,m):
                print xx
                level[xx-1] = xx+ 2
                nn = 2**(level[xx-1])



                # Create mesh and define function space
                nn = int(nn)
                NN[xx-1] = nn/2
                # parameters["form_compiler"]["quadrature_degree"] = 6
                # parameters = CP.ParameterSetup()
                mesh = UnitSquareMesh(nn,nn)

                order = 1
                parameters['reorder_dofs_serial'] = False
                Velocity = VectorFunctionSpace(mesh, "CG", order+1)
                Pressure = FunctionSpace(mesh, "CG", order)
                Magnetic = FunctionSpace(mesh, "N1curl", order)
                Lagrange = FunctionSpace(mesh, "CG", order)
                W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
                # W = Velocity*Pressure*Magnetic*Lagrange
                Velocitydim[xx-1] = Velocity.dim()
                Pressuredim[xx-1] = Pressure.dim()
                Magneticdim[xx-1] = Magnetic.dim()
                Lagrangedim[xx-1] = Lagrange.dim()
                Wdim[xx-1] = W.dim()
                print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
                dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


                def boundary(x, on_boundary):
                    return on_boundary

                u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD2D(4,1)


                bcu = DirichletBC(W.sub(0),u0, boundary)
                bcb = DirichletBC(W.sub(2),b0, boundary)
                bcr = DirichletBC(W.sub(3),r0, boundary)

                # bc = [u0,p0,b0,r0]
                bcs = [bcu,bcb,bcr]
                FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


#.........这里部分代码省略.........

    MO.PrintStr("Setting up MHD initial guess",5,"+","\n\n","\n\n")
    u_k,p_k,b_k,r_k = common.InitialGuess(FSpaces,[u0,p0,b0,r0],[F_NS,F_M],params,HiptmairMatrices,1e-6,Neumann=Expression(("0","0")),options ="New", FS = "DG")
    #plot(p_k, interactive = True) 
    b_t = TrialFunction(Velocity)
    c_t = TestFunction(Velocity)
    #print assemble(inner(b,c)*dx).array().shape
    #print mat
    #ShiftedMass = assemble(inner(mat*b,c)*dx)
    #as_vector([inner(b,c)[0]*b_k[0],inner(b,c)[1]*(-b_k[1])])

    ones = Function(Pressure)
    ones.vector()[:]=(0*ones.vector().array()+1)
    # pConst = - assemble(p_k*dx)/assemble(ones*dx)
    p_k.vector()[:] += - assemble(p_k*dx)/assemble(ones*dx)
    x = Iter.u_prev(u_k,b_k,p_k,r_k)

    KSPlinearfluids, MatrixLinearFluids = PrecondSetup.FluidLinearSetup(Pressure, MU)
    kspFp, Fp = PrecondSetup.FluidNonLinearSetup(Pressure, MU, u_k)
    #plot(b_k)

    ns,maxwell,CoupleTerm,Lmaxwell,Lns = forms.MHD2D(mesh, W,F_M,F_NS, u_k,b_k,params,IterType,"DG", SaddlePoint = "Yes")
    RHSform = forms.PicardRHS(mesh, W, u_k, p_k, b_k, r_k, params,"DG",SaddlePoint = "Yes")


    bcu = DirichletBC(W.sub(0),Expression(("0.0","0.0","0.0")), boundary)
    bcb = DirichletBC(W.sub(1),Expression(("0.0","0.0","0.0")), boundary)
    bcr = DirichletBC(W.sub(3),Expression(("0.0")), boundary)
    bcs = [bcu,bcb,bcr]

    eps = 1.0           # error measure ||u-u_k||

def foo():
    m = 7


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))
    
    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'
    kappa = 0.01
    for yy in xrange(1,5):
        kappa = kappa*10
        for xx in xrange(1,m):
            print xx
            level[xx-1] = xx+ 2
            nn = 2**(level[xx-1])



            # Create mesh and define function space
            nn = int(nn)
            NN[xx-1] = nn/2
            # parameters["form_compiler"]["quadrature_degree"] = 6
            # parameters = CP.ParameterSetup()
            mesh = UnitSquareMesh(nn,nn)

            order = 1
            parameters['reorder_dofs_serial'] = False
            Velocity = VectorFunctionSpace(mesh, "CG", order)
            Pressure = FunctionSpace(mesh, "DG", order-1)
            Magnetic = FunctionSpace(mesh, "N1curl", order)
            Lagrange = FunctionSpace(mesh, "CG", order)
            W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
            # W = Velocity*Pressure*Magnetic*Lagrange
            Velocitydim[xx-1] = Velocity.dim()
            Pressuredim[xx-1] = Pressure.dim()
            Magneticdim[xx-1] = Magnetic.dim()
            Lagrangedim[xx-1] = Lagrange.dim()
            Wdim[xx-1] = W.dim()
            print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
            dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


            def boundary(x, on_boundary):
                return on_boundary

            u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD2D(4,1)


            bcu = DirichletBC(W.sub(0),u0, boundary)
            bcb = DirichletBC(W.sub(2),b0, boundary)
            bcr = DirichletBC(W.sub(3),r0, boundary)

            # bc = [u0,p0,b0,r0]
            bcs = [bcu,bcb,bcr]
            FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


            (u, b, p, r) = TrialFunctions(W)
            (v, c, q, s) = TestFunctions(W)
            #kappa = 1.0
            MU = 1.0
            Mu_m =1e1
            IterType = 'Combined'

            F_NS = -MU*Laplacian+Advection+gradPres-kappa*NS_Couple
#.........这里部分代码省略.........

def foo():
    m = 6


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))
    
    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'

    MU[0]= 1e0
    for xx in xrange(1,m):
        print xx
        level[xx-1] = xx+ 0
        nn = 2**(level[xx-1])



        # Create mesh and define function space
        nn = int(nn)
        NN[xx-1] = nn/2
        # parameters["form_compiler"]["quadrature_degree"] = 6
        # parameters = CP.ParameterSetup()
        mesh = UnitCubeMesh(nn,nn,nn)

        order = 2
        parameters['reorder_dofs_serial'] = False
        Velocity = VectorFunctionSpace(mesh, "CG", order)
        Pressure = FunctionSpace(mesh, "CG", order-1)
        Magnetic = FunctionSpace(mesh, "N1curl", order)
        Lagrange = FunctionSpace(mesh, "CG", order)
        W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
        # W = Velocity*Pressure*Magnetic*Lagrange
        Velocitydim[xx-1] = Velocity.dim()
        Pressuredim[xx-1] = Pressure.dim()
        Magneticdim[xx-1] = Magnetic.dim()
        Lagrangedim[xx-1] = Lagrange.dim()
        Wdim[xx-1] = W.dim()
        print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
        dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


        def boundary(x, on_boundary):
            return on_boundary

        u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD3D(4,1)


        bcu = DirichletBC(Velocity,u0, boundary)
        bcb = DirichletBC(Magnetic,b0, boundary)
        bcr = DirichletBC(Lagrange,r0, boundary)

        # bc = [u0,p0,b0,r0]
        bcs = [bcu,bcb,bcr]
        FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


        (u, b, p, r) = TrialFunctions(W)
        (v, c, q, s) = TestFunctions(W)
        kappa = 1.0
        Mu_m =10.0
        MU = 1.0/1
        IterType = 'Full'
        Split = "No"
        Saddle = "No"
        Stokes = "No"
#.........这里部分代码省略.........

def foo():
    m = 2


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))
    
    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'

    MU[0]= 1e0
    for xx in xrange(1,m):
        print xx
        level[xx-1] = xx+ 10
        nn = 2**(level[xx-1])



        # Create mesh and define function space
        nn = int(nn)
        NN[xx-1] = nn/2
        # parameters["form_compiler"]["quadrature_degree"] = 6
        # parameters = CP.ParameterSetup()
        mesh = UnitSquareMesh(nn,nn)

        order = 1
        parameters['reorder_dofs_serial'] = False
        Velocity = VectorFunctionSpace(mesh, "CG", order)
        Pressure = FunctionSpace(mesh, "DG", order-1)
        Magnetic = FunctionSpace(mesh, "N1curl", order)
        Lagrange = FunctionSpace(mesh, "CG", order)
        W = MixedFunctionSpace([Velocity,Magnetic, Pressure, Lagrange])
        # W = Velocity*Pressure*Magnetic*Lagrange
        Velocitydim[xx-1] = Velocity.dim()
        Pressuredim[xx-1] = Pressure.dim()
        Magneticdim[xx-1] = Magnetic.dim()
        Lagrangedim[xx-1] = Lagrange.dim()
        Wdim[xx-1] = W.dim()
        print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
        dim = [Velocity.dim(), Magnetic.dim(), Pressure.dim(), Lagrange.dim()]


        def boundary(x, on_boundary):
            return on_boundary

        u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD2D(4,1)


        bcu = DirichletBC(W.sub(0),u0, boundary)
        bcb = DirichletBC(W.sub(1),b0, boundary)
        bcr = DirichletBC(W.sub(3),r0, boundary)

        # bc = [u0,p0,b0,r0]
        bcs = [bcu,bcb,bcr]
        FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


        (u, b, p, r) = TrialFunctions(W)
        (v, c, q, s) = TestFunctions(W)
        kappa = 1.0
        Mu_m =1e1
        MU = 1.0/1
        IterType = 'Full'

        F_NS = -MU*Laplacian+Advection+gradPres-kappa*NS_Couple
        if kappa == 0:
#.........这里部分代码省略.........

def foo():
    m = 6


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))

    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'

    MU[0]= 1e0
    for xx in xrange(1,m):
        print xx
        level[xx-1] = xx + 0
        nn = 2**(level[xx-1])



        # Create mesh and define function space
        nn = int(nn)
        NN[xx-1] = nn/2
        # parameters["form_compiler"]["quadrature_degree"] = 6
        # parameters = CP.ParameterSetup()
        mesh = UnitSquareMesh(nn,nn)

        order = 2
        parameters['reorder_dofs_serial'] = False
        Velocity = VectorFunctionSpace(mesh, "CG", order)
        Pressure = FunctionSpace(mesh, "CG", order-1)
        Magnetic = FunctionSpace(mesh, "N1curl", order-1)
        Lagrange = FunctionSpace(mesh, "CG", order-1)
        W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
        # W = Velocity*Pressure*Magnetic*Lagrange
        Velocitydim[xx-1] = Velocity.dim()
        Pressuredim[xx-1] = Pressure.dim()
        Magneticdim[xx-1] = Magnetic.dim()
        Lagrangedim[xx-1] = Lagrange.dim()
        Wdim[xx-1] = W.dim()
        print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
        dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


        def boundary(x, on_boundary):
            return on_boundary

        u0, p0,b0, r0, Laplacian, Advection, gradPres,CurlCurl, gradR, NS_Couple, M_Couple = ExactSol.MHD2D(4,1, mesh)


        bcu = DirichletBC(Velocity,u0, boundary)
        bcb = DirichletBC(Magnetic,b0, boundary)
        bcr = DirichletBC(Lagrange,r0, boundary)

        # bc = [u0,p0,b0,r0]
        bcs = [bcu,bcb,bcr]
        FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


        (u, b, p, r) = TrialFunctions(W)
        (v, c, q, s) = TestFunctions(W)
        kappa = 10.0
        Mu_m =10.0
        MU = 1.0/1
        IterType = 'Full'
        Split = "No"
        Saddle = "No"
        Stokes = "No"
#.........这里部分代码省略.........

def foo():
    m = 2


    errL2u =np.zeros((m-1,1))
    errH1u =np.zeros((m-1,1))
    errL2p =np.zeros((m-1,1))
    errL2b =np.zeros((m-1,1))
    errCurlb =np.zeros((m-1,1))
    errL2r =np.zeros((m-1,1))
    errH1r =np.zeros((m-1,1))



    l2uorder =  np.zeros((m-1,1))
    H1uorder =np.zeros((m-1,1))
    l2porder =  np.zeros((m-1,1))
    l2border =  np.zeros((m-1,1))
    Curlborder =np.zeros((m-1,1))
    l2rorder =  np.zeros((m-1,1))
    H1rorder = np.zeros((m-1,1))

    NN = np.zeros((m-1,1))
    DoF = np.zeros((m-1,1))
    Velocitydim = np.zeros((m-1,1))
    Magneticdim = np.zeros((m-1,1))
    Pressuredim = np.zeros((m-1,1))
    Lagrangedim = np.zeros((m-1,1))
    Wdim = np.zeros((m-1,1))
    iterations = np.zeros((m-1,1))
    SolTime = np.zeros((m-1,1))
    udiv = np.zeros((m-1,1))
    MU = np.zeros((m-1,1))
    level = np.zeros((m-1,1))
    NSave = np.zeros((m-1,1))
    Mave = np.zeros((m-1,1))
    TotalTime = np.zeros((m-1,1))

    nn = 2

    dim = 2
    ShowResultPlots = 'yes'
    split = 'Linear'

    MU[0]= 1e0
    for xx in xrange(1,m):
        # print xx
        level[xx-1] = xx + 6
        nn = 2**(level[xx-1])



        # Create mesh and define function space
        #nn = int(nn)
        #NN[xx-1] = nn/2
        # parameters["form_compiler"]["quadrature_degree"] = 6
        # parameters = CP.ParameterSetup()
        #mesh = UnitSquareMesh(nn,nn)

        N = int(sys.argv[1])
        n = int(2**N)
        mesh = UnitCubeMesh(n,n,n)
        order = 2
        parameters['reorder_dofs_serial'] = False
        Velocity = VectorFunctionSpace(mesh, "CG", order)
        Pressure = FunctionSpace(mesh, "CG", order-1)
        Magnetic = FunctionSpace(mesh, "N1curl", order-1)
        Lagrange = FunctionSpace(mesh, "CG", order-1)
        W = MixedFunctionSpace([Velocity, Pressure, Magnetic,Lagrange])
        # W = Velocity*Pressure*Magnetic*Lagrange
        Velocitydim[xx-1] = Velocity.dim()
        Pressuredim[xx-1] = Pressure.dim()
        Magneticdim[xx-1] = Magnetic.dim()
        Lagrangedim[xx-1] = Lagrange.dim()
        Wdim[xx-1] = W.dim()
        print W.dim()
        #print "\n\nW:  ",Wdim[xx-1],"Velocity:  ",Velocitydim[xx-1],"Pressure:  ",Pressuredim[xx-1],"Magnetic:  ",Magneticdim[xx-1],"Lagrange:  ",Lagrangedim[xx-1],"\n\n"
        dim = [Velocity.dim(), Pressure.dim(), Magnetic.dim(), Lagrange.dim()]


        # def boundary(x, on_boundary):
        #     return on_boundary




        # bc = [u0,p0,b0,r0]
        FSpaces = [Velocity,Pressure,Magnetic,Lagrange]


        (u, b, p, r) = TrialFunctions(W)
        (v, c, q, s) = TestFunctions(W)
        kappa = 10.0
        Mu_m =10.0
        MU = 1.0/1
        IterType = 'Full'
        Split = "No"
        Saddle = "No"
        Stokes = "No"
        SetupType = 'python-class'
#.........这里部分代码省略.........