Python PETScIO.vecToArray方法代码示例

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def PicardToleranceDecouple(x,U,FSpaces,dim,NormType,iter):
    X = IO.vecToArray(x)
    uu = X[0:dim[0]]
    pp = X[dim[0]:dim[0]+dim[1]-1]
    bb = X[dim[0]+dim[1]-1:dim[0]+dim[1]+dim[2]-1]
    rr = X[dim[0]+dim[1]+dim[2]-1:]

    u = Function(FSpaces[0])
    u.vector()[:] = uu
    diffu = u.vector().array()

    p = Function(FSpaces[1])
    n = pp.shape
    p.vector()[:] = np.insert(pp,n,0)
    # ones = Function(FSpaces[1])
    # ones.vector()[:]=(0*ones.vector().array()+1)
    # pp = Function(FSpaces[1])
    # p.vector()[:] = p.vector().array()- assemble(p*dx)/assemble(ones*dx)

    b = Function(FSpaces[2])
    b.vector()[:] = bb
    diffb = b.vector().array()

    r = Function(FSpaces[3])
    print r.vector().array().shape
    print rr.shape
    r.vector()[:] = rr


    if (NormType == '2'):
        epsu = splin.norm(diffu)/sqrt(dim[0])
        epsp = splin.norm(p.vector().array())/sqrt(dim[1])
        epsb = splin.norm(diffb)/sqrt(dim[2])
        epsr = splin.norm(r.vector().array())/sqrt(dim[3])
    elif (NormType == 'inf'):
        epsu = splin.norm(diffu, ord=np.Inf)
        epsp = splin.norm(p.vector().array(),ord=np.inf)
        epsb = splin.norm(diffb, ord=np.Inf)
        epsr = splin.norm(r.vector().array(),ord=np.inf)
    else:
        print "NormType must be 2 or inf"
        quit()


    RHS = IO.vecToArray(U)
    u.vector()[:] = uu + RHS[0:dim[0]]
    p.vector()[:] = p.vector().array() + U.array[dim[0]:dim[0]+dim[1]]
    b.vector()[:] = bb + RHS[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]
    r.vector()[:] = rr + RHS[dim[0]+dim[1]+dim[2]:]




    print 'iter=%d:\n u-norm=%g   p-norm=%g  \n b-norm=%g   r-norm=%g' % (iter, epsu,epsp,epsb,epsr), '\n\n\n'
    return u,p,b,r,epsu+epsp+epsb+epsr

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def PicardTolerance(x,u_k,b_k,FSpaces,dim,NormType,iter):
    X = IO.vecToArray(x)
    uu = X[0:dim[0]]
    bb = X[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]

    u = Function(FSpaces[0])
    u.vector()[:] = u.vector()[:] + uu
    diffu = u.vector().array() - u_k.vector().array()

    b = Function(FSpaces[2])
    b.vector()[:] = b.vector()[:] + bb
    diffb = b.vector().array() - b_k.vector().array()
    if (NormType == '2'):
        epsu = splin.norm(diffu)/sqrt(dim[0])
        epsb = splin.norm(diffb)/sqrt(dim[0])
    elif (NormType == 'inf'):
        epsu = splin.norm(diffu, ord=np.Inf)
        epsb = splin.norm(diffb, ord=np.Inf)
    else:
        print "NormType must be 2 or inf"
        quit()

    print 'iter=%d: u-norm=%g   b-norm=%g ' % (iter, epsu,epsb)
    u_k.assign(u)
    b_k.assign(b)


    return u_k,b_k,epsu,epsb

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
        toc()
        ksp.solve(bb, x)

        time = toc()
        print time
        SolutionTime = SolutionTime +time
        print ksp.its
        outerit += ksp.its
        # r = bb.duplicate()
        # A.MUlt(x, r)
        # r.aypx(-1, bb)
        # rnorm = r.norm()
        # PETSc.Sys.Print('error norm = %g' % rnorm,comm=PETSc.COMM_WORLD)

        uu = IO.vecToArray(x)
        UU = uu[0:Vdim[xx-1][0]]
        # time = time+toc()
        u1 = Function(V)
        u1.vector()[:] = u1.vector()[:] + UU

        pp = uu[Vdim[xx-1][0]:]
        # time = time+toc()
        p1 = Function(Q)
        n = pp.shape

        p1.vector()[:] = p1.vector()[:] +  pp
        diff = u1.vector().array()
        eps = np.linalg.norm(diff, ord=np.Inf)

        print '\n\n\niter=%d: norm=%g' % (iter, eps)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
    # AvIt[xx-1] = np.ceil(outerit/iter)
    u = interpolate(ue,V)
    p = interpolate(pe,Q)

    ua = Function(V)
    ua.vector()[:] = u_k.vector().array()
    # nonlinear[xx-1] = assemble(inner((grad(ua)*ua),ua)*dx+(1/2)*div(ua)*inner(ua,ua)*dx- (1/2)*inner(ua,n)*inner(ua,ua)*ds)
    # parameters["form_compiler"]["quadrature_degree"] = 7
    VelocityE = VectorFunctionSpace(mesh,"CG",6)
    u = interpolate(ue,VelocityE)

    PressureE = FunctionSpace(mesh,"CG",5)

    Nv  = ua.vector().array().shape

    X = IO.vecToArray(r)
    xu = X[0:V.dim()]
    ua = Function(V)
    ua.vector()[:] = xu

    pp = X[V.dim():V.dim()+Q.dim()]


    n = pp.shape
    pa = Function(Q)
    pa.vector()[:] = pp

    pend = assemble(pa*dx)

    ones = Function(Q)
    ones.vector()[:]=(0*pp+1)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def PicardToleranceDecouple(x,U,FSpaces,dim,NormType,iter,SaddlePoint = "No"):
    X = IO.vecToArray(x)
    uu = X[0:dim[0]]

    if SaddlePoint == "Yes":
        bb = X[dim[0]:dim[0]+dim[1]]
        pp = X[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]
    else:
        pp = X[dim[0]:dim[0]+dim[1]]
        bb = X[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]

    rr = X[dim[0]+dim[1]+dim[2]:]

    u = Function(FSpaces[0])
    u.vector()[:] = uu
    u_ = assemble(inner(u,u)*dx)
    diffu = u.vector().array()

    # if SaddlePoint == "Yes":
    #     p = Function(FSpaces[2])
    #     p.vector()[:] = pp
    #     ones = Function(FSpaces[2])
    #     ones.vector()[:]=(0*ones.vector().array()+1)
    #     pp = Function(FSpaces[2])
    #     print ones
    #     pp.vector()[:] = p.vector().array()- assemble(p*dx)/assemble(ones*dx)
    #     p = pp.vector().array()
    #     b = Function(FSpaces[1])
    #     b.vector()[:] = bb
    #     diffb = b.vector().array()
    # else:
    print pp.shape
    p = Function(FSpaces[1])
    print FSpaces[1].dim()
    p.vector()[:] = pp
    p_ = assemble(p*p*dx)

    ones = Function(FSpaces[1])
    ones.vector()[:]=(0*ones.vector().array()+1)
    pp = Function(FSpaces[1])
    pp.vector()[:] = p.vector().array() - assemble(p*dx)/assemble(ones*dx)
    p_ = assemble(pp*pp*dx)
    p = pp.vector().array()
    b = Function(FSpaces[2])
    b.vector()[:] = bb
    b_ = assemble(inner(b,b)*dx)
    diffb = b.vector().array()

    r = Function(FSpaces[3])
    r.vector()[:] = rr
    r_ = assemble(r*r*dx)
    # print diffu
    if (NormType == '2'):
        epsu = splin.norm(diffu)/sqrt(dim[0])
        epsp = splin.norm(pp.vector().array())/sqrt(dim[1])
        epsb = splin.norm(diffb)/sqrt(dim[2])
        epsr = splin.norm(r.vector().array())/sqrt(dim[3])
    elif (NormType == 'inf'):
        epsu = splin.norm(diffu, ord=np.Inf)
        epsp = splin.norm(pp.vector().array(),ord=np.inf)
        epsb = splin.norm(diffb, ord=np.Inf)
        epsr = splin.norm(r.vector().array(),ord=np.inf)

    else:
        print "NormType must be 2 or inf"
        quit()
    # U.axpy(1,x)
    p = Function(FSpaces[1])
    RHS = IO.vecToArray(U+x)

    if SaddlePoint == "Yes":
        u.vector()[:] = RHS[0:dim[0]]
        p.vector()[:] = pp.vector().array()+U.array[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]
        b.vector()[:] = RHS[dim[0]:dim[0]+dim[1]]
        r.vector()[:] = RHS[dim[0]+dim[1]+dim[2]:]
    else:
        u.vector()[:] = RHS[0:dim[0]]
        p.vector()[:] = pp.vector().array()+U.array[dim[0]:dim[0]+dim[1]]
        b.vector()[:] = RHS[dim[0]+dim[1]:dim[0]+dim[1]+dim[2]]
        r.vector()[:] = RHS[dim[0]+dim[1]+dim[2]:]

    # print diffu.dot(diffu) + pp.dot(pp) + diffb.dot(diffb) + r.dot(r)
    # epsu = sqrt(u_)/sqrt(dim[0])
    # epsp = sqrt(p_)/sqrt(dim[1])
    # epsb = sqrt(b_)/sqrt(dim[2])
    # epsr = sqrt(r_)/sqrt(dim[3])
        # uOld = np.concatenate((diffu, pp.vector().array(), diffb, r.vector().array()), axis=0)
    # print np.linalg.norm(uOld)/sum(dim)

    print 'u-norm=%g   p-norm=%g  \n b-norm=%g   r-norm=%g' % (epsu,epsp,epsb,epsr), '\n\n\n'
    print 'u-norm=%g   p-norm=%g  \n b-norm=%g   r-norm=%g' % (sqrt(u_), sqrt(p_), sqrt(b_), sqrt(r_)), '\n\n\n'

    return u,p,b,r,epsu+epsp+epsb+epsr

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
            ksp.solve(bb, x)

            time = toc()
            print time
            SolutionTime = SolutionTime +time
            print ksp.its
            outerit += ksp.its
            del ksp,pc

            # r = bb.duplicate()
            # A.MUlt(x, r)
            # r.aypx(-1, bb)
            # rnorm = r.norm()
            # PETSc.Sys.Print('error norm = %g' % rnorm,comm=PETSc.COMM_WORLD)

            uu = IO.vecToArray(x)
            UU = uu[0:Vdim[xx-1][0]]
            # time = time+toc()
            u1 = Function(V)
            u1.vector()[:] = u1.vector()[:] + UU

            pp = uu[Vdim[xx-1][0]:]
            # time = time+toc()
            p1 = Function(Q)
            n = pp.shape

            p1.vector()[:] = p1.vector()[:] +  pp
            diff = u1.vector().array()
            eps = np.linalg.norm(diff, ord=np.Inf)

            print '\n\n\niter=%d: norm=%g' % (iter, eps)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
        time = toc()
        print time
        SolutionTime = SolutionTime +time
        print ":::::::::::::::::::::::::::::::::::::::::::::",ksp.its
        outerit += ksp.its



        # print (A*u-b).array
        # print (A*x-b).array
        print np.linalg.norm(np.abs((A*x).array)-np.abs((A*u).array))
        # ss
        # for line in reshist.values():
        #     print line
        uu = IO.vecToArray(u)
        UU = uu[0:Vdim[xx-1][0]]
        # UU[low_values_indices] = 0.0
        # print UU
        # time = time+toc()
        u1 = Function(V)
        u1.vector()[:] = UU - u_k.vector().array()


        # p1.vector()[:] = p1.vector()[:] +  pp
        # p1.vector()[:] += - assemble(p1*dx)/assemble(ones*dx)
        diff = u1.vector().array()
        # print p1.vector().array()
        eps = np.linalg.norm(diff, ord=np.Inf) #+np.linalg.norm(p1.vector().array(),ord=np.Inf)

        print '\n\n\niter=%d: norm=%g' % (iter, eps)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
 def Amatvec(v):
     x = IO.arrayToVec(v)
     B = x.duplicate()
     AA.mult(x,B)
     b = IO.vecToArray(B)
     return b

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def DirectErrors(x,mesh,FSpaces,ExactSolution,k,dim):

    Vdim = dim[0]
    Pdim = dim[1]
    Mdim = dim[2]
    Rdim = dim[3]

    VelocityE = VectorFunctionSpace(mesh,"CG",k+2)
    u = interpolate(ExactSolution[0],VelocityE)

    PressureE = FunctionSpace(mesh,"CG",k+1)


    MagneticE = FunctionSpace(mesh,"N1curl",k+2)
    b = interpolate(ExactSolution[2],MagneticE)

    LagrangeE = FunctionSpace(mesh,"CG",k+2)
    r = interpolate(ExactSolution[3],LagrangeE)

    X = IO.vecToArray(x)
    xu = X[0:Vdim]
    ua = Function(FSpaces[0])
    ua.vector()[:] = xu

    pp = X[Vdim:Vdim+Pdim-1]
    # xp[-1] = 0
    # pa = Function(Pressure)
    # pa.vector()[:] = xp

    n = pp.shape
    pp = np.insert(pp,n,0)
    pa = Function(FSpaces[1])
    pa.vector()[:] = pp

    pend = assemble(pa*dx)

    ones = Function(FSpaces[1])
    ones.vector()[:]=(0*pp+1)
    pp = Function(FSpaces[1])
    pp.vector()[:] = pa.vector().array()- assemble(pa*dx)/assemble(ones*dx)

    pInterp = interpolate(ExactSolution[1],PressureE)
    pe = Function(PressureE)
    pe.vector()[:] = pInterp.vector().array()
    const = - assemble(pe*dx)/assemble(ones*dx)
    pe.vector()[:] = pe.vector()[:]+const


    xb = X[Vdim+Pdim-1:Vdim+Pdim+Mdim-1]
    ba = Function(FSpaces[2])
    ba.vector()[:] = xb

    xr = X[Vdim+Pdim+Mdim-1:]
    ra = Function(FSpaces[3])
    ra.vector()[:] = xr

    ErrorU = Function(FSpaces[0])
    ErrorP = Function(FSpaces[1])
    ErrorB = Function(FSpaces[2])
    ErrorR = Function(FSpaces[3])

    ErrorU = u-ua
    ErrorP = pe-pp
    ErrorB = b-ba
    ErrorR = r-ra


    errL2u= sqrt(abs(assemble(inner(ErrorU, ErrorU)*dx)))
    errH1u= sqrt(abs(assemble(inner(grad(ErrorU), grad(ErrorU))*dx)))
    errL2p= sqrt(abs(assemble(inner(ErrorP, ErrorP)*dx)))
    errL2b= sqrt(abs(assemble(inner(ErrorB, ErrorB)*dx)))
    errCurlb = sqrt(abs(assemble(inner(curl(ErrorB), curl(ErrorB))*dx)))
    errL2r= sqrt(abs(assemble(inner(ErrorR, ErrorR)*dx)))
    errH1r= sqrt(abs(assemble(inner(grad(ErrorR), grad(ErrorR))*dx)))

    return errL2u, errH1u, errL2p, errL2b, errCurlb, errL2r, errH1r

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def Errors(x,mesh,FSpaces,ExactSolution,k,dim):
    Vdim = dim[0]
    Pdim = dim[1]
    Mdim = dim[2]
    Rdim = dim[3]
    # k +=2
    VelocityE = VectorFunctionSpace(mesh,"CG",k+4)
    u = interpolate(ExactSolution[0],VelocityE)

    PressureE = FunctionSpace(mesh,"CG",k+3)


    MagneticE = FunctionSpace(mesh,"N1curl",k+4)
    b = interpolate(ExactSolution[2],MagneticE)

    LagrangeE = FunctionSpace(mesh,"CG",k+4)
    r = interpolate(ExactSolution[3],LagrangeE)
    # parameters["form_compiler"]["quadrature_degree"] = 14
    X = IO.vecToArray(x)
    xu = X[0:Vdim]
    ua = Function(FSpaces[0])
    ua.vector()[:] = xu

    pp = X[Vdim:Vdim+Pdim-1]


    n = pp.shape
    pp = np.insert(pp,n,0)
    pa = Function(FSpaces[1])
    pa.vector()[:] = pp

    pend = assemble(pa*dx)

    ones = Function(FSpaces[1])
    ones.vector()[:]=(0*pp+1)
    pp = Function(FSpaces[1])
    pp.vector()[:] = pa.vector().array()- assemble(pa*dx)/assemble(ones*dx)

    pInterp = interpolate(ExactSolution[1],PressureE)
    pe = Function(PressureE)
    pe.vector()[:] = pInterp.vector().array()
    const = - assemble(pe*dx)/assemble(ones*dx)
    pe.vector()[:] = pe.vector()[:]+const


    xb = X[Vdim+Pdim-1:Vdim+Pdim+Mdim-1]
    ba = Function(FSpaces[2])
    ba.vector()[:] = xb

    xr = X[Vdim+Pdim+Mdim-1:]
    ra = Function(FSpaces[3])
    ra.vector()[:] = xr

    ErrorU = Function(FSpaces[0])
    ErrorP = Function(FSpaces[1])
    ErrorB = Function(FSpaces[2])
    ErrorR = Function(FSpaces[3])


    plot(ua)
    plot(pp)
    plot(ba)
    plot(ra)

    ErrorU = u-ua
    ErrorP = pe-pp
    ErrorB = b-ba
    ErrorR = r-ra

    # errL2u = errornorm(u,ua ,norm_type='L2',degree_rise=4)
    # errH1u = errornorm(u,ua ,norm_type='H10',degree_rise=4)
    # errL2p = errornorm(pe,pp ,norm_type='L2',degree_rise=4)
    # errL2b = errornorm(b,ba ,norm_type='L2',degree_rise=4)
    # errCurlb  =errornorm(b,ba ,norm_type='Hcurl0',degree_rise=4)
    # errL2r = errornorm(r,ra ,norm_type='L2',degree_rise=4)
    # errH1r = errornorm(r,ra ,norm_type='H10',degree_rise=4)
    errL2u= sqrt(abs(assemble(inner(ErrorU, ErrorU)*dx)))
    errH1u= sqrt(abs(assemble(inner(grad(ErrorU), grad(ErrorU))*dx)))
    errL2p= sqrt(abs(assemble(inner(ErrorP, ErrorP)*dx)))
    errL2b= sqrt(abs(assemble(inner(ErrorB, ErrorB)*dx)))
    errCurlb = sqrt(abs(assemble(inner(curl(ErrorB), curl(ErrorB))*dx)))
    errL2r= sqrt(abs(assemble(inner(ErrorR, ErrorR)*dx)))
    errH1r= sqrt(abs(assemble(inner(grad(ErrorR), grad(ErrorR))*dx)))

    return errL2u, errH1u, errL2p, errL2b, errCurlb, errL2r, errH1r

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def Stokes(V,Q,BC,f,mu,boundaries):

    parameters = CP.ParameterSetup()

    W = V*Q

    (u, p) = TrialFunctions(W)
    (v, q) = TestFunctions(W)

    def boundary(x, on_boundary):
        return on_boundary


    bcu1 = DirichletBC(W.sub(0),BC[0], boundaries,2)
    bcu2 = DirichletBC(W.sub(0),BC[1], boundaries,1)
    bcs = [bcu1,bcu2]
    u_k = Function(V)
    a11 = mu*inner(grad(v), grad(u))*dx
    a12 = div(v)*p*dx
    a21 = div(u)*q*dx
    L1  =  inner(v,f)*dx
    a = mu*a11-a12-a21
    i = p*q*dx

    tic()
    AA, bb = assemble_system(a, L1, bcs)
    # A = as_backend_type(AA).mat()
    A,b = CP.Assemble(AA,bb)
    print toc()
    b = bb.array()
    zeros = 0*b
    del bb
    bb = IO.arrayToVec(b)
    x = IO.arrayToVec(zeros)


    pp = inner(grad(v), grad(u))*dx + (1/mu)*p*q*dx
    PP, Pb = assemble_system(pp,L1,bcs)
    P = CP.Assemble(PP)


    u_is = PETSc.IS().createGeneral(range(V.dim()))
    p_is = PETSc.IS().createGeneral(range(V.dim(),V.dim()+Q.dim()))

    ksp = PETSc.KSP().create()
    ksp.setOperators(A,P)
    pc = ksp.getPC()
    pc.setType(pc.Type.FIELDSPLIT)
    fields = [ ("field1", u_is), ("field2", p_is)]
    pc.setFieldSplitIS(*fields)
    pc.setFieldSplitType(0)

    OptDB = PETSc.Options()

    OptDB["field_split_type"] = "multiplicative"

    OptDB["fieldsplit_field1_ksp_type"] = "preonly"
    OptDB["fieldsplit_field1_pc_type"] = "hypre"
    OptDB["fieldsplit_field2_ksp_type"] = "cg"
    OptDB["fieldsplit_field2_pc_type"] = "jacobi"
    OptDB["fieldsplit_field2_ksp_rtol"] = 1e-8



    ksp.setFromOptions()
    ksp.setTolerances(1e-8)


    ksp.solve(bb, x)


    X = IO.vecToArray(x)
    x = X[0:V.dim()]
    p =  X[V.dim():]
    # x =
    u = Function(V)
    u.vector()[:] = x

    p_k = Function(Q)
    n = p.shape
    p_k.vector()[:] = p


    u_k.assign(u)


    return u_k, p_k

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
        time = toc()
        print time
        SolutionTime = SolutionTime +time
        print ":::::::::::::::::::::::::::::::::::::::::::::",ksp.its
        outerit += ksp.its



        # print (A*u-b).array
        # print (A*x-b).array
        print np.linalg.norm(np.abs((A*x).array)-np.abs((A*u).array))
        # ss
        # for line in reshist.values():
        #     print line
        uu = IO.vecToArray(u)
        UU = uu[0:Vdim[xx-1][0]]
        # UU[low_values_indices] = 0.0
        # print UU
        # time = time+toc()
        u1 = Function(V)
        u1.vector()[:] = u1.vector()[:] + UU

        pp = uu[V.dim():V.dim()+Q.dim()]
        # print
        # time = time+toc()
        p1 = Function(Q)
        # n = pp.shape
        # pend = assemble(pa*dx)

# 需要导入模块: import PETScIO [as 别名]
# 或者: from PETScIO import vecToArray [as 别名]
def Maxwell(u,V,Q,BC,f,params,HiptmairMatrices,Hiptmairtol,InitialTol,Neumann=None):
    # parameters['linear_algebra_backend'] = 'uBLAS'
    dim = f.shape()[0]
    W = V*Q

    (b, r) = TrialFunctions(W)
    (c,s) = TestFunctions(W)

    def boundary(x, on_boundary):
        return on_boundary

    bcb = DirichletBC(W.sub(0),BC[0], boundary)
    bcr = DirichletBC(W.sub(1),BC[1], boundary)
    if params[0] == 0:
        a11 = params[1]*inner(curl(c),curl(b))*dx
    else:
        if dim == 2:
            a11 = params[1]*params[0]*inner(curl(c),curl(b))*dx
        elif dim == 3:
            a11 = params[1]*params[0]*inner(curl(c),curl(b))*dx

    a12 = inner(c,grad(r))*dx
    a21 = inner(b,grad(s))*dx
    Lmaxwell  = inner(c, f)*dx
    maxwell = a11+a12+a21



    bcs =  [bcb,bcr]
    tic()
    AA, bb = assemble_system(maxwell, Lmaxwell, bcs)
    A,bb = CP.Assemble(AA,bb)
    del AA
    x = bb.duplicate()

    u_is = PETSc.IS().createGeneral(range(V.dim()))
    p_is = PETSc.IS().createGeneral(range(V.dim(),V.dim()+Q.dim()))




    ksp = PETSc.KSP().create()
    ksp.setTolerances(InitialTol)
    pc = ksp.getPC()
    pc.setType(PETSc.PC.Type.PYTHON)
    # ksp.setType('preonly')
    # pc.setType('lu')
    # [G, P, kspVector, kspScalar, kspCGScalar, diag]
    reshist = {}
    def monitor(ksp, its, rnorm):
        print rnorm
        reshist[its] = rnorm


    # ksp.setMonitor(monitor)
    if V.__str__().find("N1curl2") == -1:
        ksp.setOperators(A)
        pc.setPythonContext(MP.Hiptmair(W, HiptmairMatrices[3], HiptmairMatrices[4], HiptmairMatrices[2], HiptmairMatrices[0], HiptmairMatrices[1], HiptmairMatrices[6],Hiptmairtol))
    else:
        p = params[1]*params[0]*inner(curl(c),curl(b))*dx+inner(c,b)*dx + inner(grad(r),grad(s))*dx
        P = assemble(p)
        for bc in bcs:
            bc.apply(P)
        P = CP.Assemble(P)
        ksp.setOperators(A,P)
        pc.setType(PETSc.PC.Type.LU)
#        pc.setPythonContext(MP.Direct(W))

    scale = bb.norm()
    bb = bb/scale
    start_time = time.time()
    ksp.solve(bb, x)
    print ("{:40}").format("Maxwell solve, time: "), " ==>  ",("{:4f}").format(time.time() - start_time),("{:9}").format("   Its: "), ("{:4}").format(ksp.its),  ("{:9}").format("   time: "), ("{:4}").format(time.strftime('%X %x %Z')[0:5])

    x = x*scale


    ksp.destroy()
    X = IO.vecToArray(x)
    x = X[0:V.dim()]
    ua = Function(V)
    ua.vector()[:] = x


    p =  X[V.dim():]
    pa = Function(Q)
    pa.vector()[:] = p
    del ksp,pc
    return ua, pa