gistfile1.txt

from dolfin import *
import numpy

def project_assign_rho(w, rho_init):
  W = w.function_space()
  rho, u, p, s, lmbd_u, lmbd_p, lmbd_r = split(w)
  vec = as_vector((rho_init, u[0], u[1], p, s, lmbd_u[0], lmbd_u[1], lmbd_p, lmbd_r))
  w_proj = project(vec, W)
  w.assign(w_proj)  


def newton_solve(F, u, bcs=None):
  # Project to feasible solution
  m = u.split(deepcopy=True)[0]
  eps = 0.01
  m.vector()[:] = numpy.maximum(m.vector().array(), eps) 
  m.vector()[:] = numpy.minimum(m.vector().array(), 1 - eps) 
  #assign(u.sub(0), m)  # Assign subspace, requires trunk dolfin
  project_assign_rho(u, m)


  i = 0 
  maxiter = 10
  du = Function(u.function_space())
  mvec = u.split(deepcopy=True)[0].vector()
  print "initial m: (min, max) == (%s, %s)" % (mvec.min(), mvec.max())
  #plot(u.split(deepcopy=True)[0], interactive=True)

  while True:
    J = assemble(derivative(F, u))
    du.vector().zero()
    b = assemble(-F)

    oldu = Function(u.function_space())
    oldu.assign(u)

    if bcs is not None:
      [bc.apply(J, b, u.vector()) for bc in bcs]

    print "%d: residual norm %s" % (i, b.norm("l2"))

    solve(J, du.vector(), b, "lu")

    mvec = u.split(deepcopy=True)[0].vector()

    print "%d: update norm %s" % (i, du.vector().norm("l2"))
    #plot(du[0], interactive=True)
    u.assign(u + 0.1 * du)
    plot(u.split(deepcopy=True)[0], interactive=True)

    # Project to feasible solution
    m = u.split(deepcopy=True)[0]
    eps = 0.01
    m.vector()[:] = numpy.maximum(m.vector().array(), eps) 
    m.vector()[:] = numpy.minimum(m.vector().array(), 1 - eps) 
    #assign(u.sub(0), m)  # Assign subspace, requires trunk dolfin
    project_assign_rho(u, m)

    mvec = u.split(deepcopy=True)[0].vector()
    print "updated m: (min, max) == (%s, %s)" % (mvec.min(), mvec.max())
    print "integral of updated m: ", assemble(u[0]*dx)
    i = i + 1

    plot(u.split(deepcopy=True)[2], title="Pressure")

    if (oldu.vector() - u.vector()).norm("l2") < 1.0e-9: break
    if i > maxiter: break 

    oldu.assign(u)

def picard_solve(F, u, u_new, bcs=None):
  i = 0 
  while True:
    solve(F == 0, u_new)
    print "%d: update norm %s" % (i, (u_new.vector() - u.vector()).norm("l2"))

    u.assign(u_new)
    i = i + 1