add PostponedTopology, BasisArgument

[joostvanzwieten]

Problem

Scripts with hierarchical topologies are require residuals and namespaces to be redefined for each new topology. A typical script looks like:

topo = ...
basis = topo.basis(...)
lhs = ...
while ...:
  ns = function.Namespace()
  ns.basis = basis
  ns.u = '?lhs_n basis_n'
  res = topo.integral('basis_n,i u_,i d:x' @ ns, degree=...)
  lhs = solver.newton('lhs', res, lhs0=lhs, ...)
  if ...:
    break
  topo = topo.refined_by(...)
  ns.newbasis = topo.basis(...)
  ns.newu = '?newlhs_n newbasis_n'
  proj = topo.integral('(newu - u)^2 d:x' @ ns, degree=...)
  lhs = solver.optimize('newlhs', proj)

Proposal

Add PostponedTopology as a postponed Topology and BasisArgument as postponed Basis. The PostponedTopology is identified by one or more roots. Example:

X, Y = Root('X', 1), Root('Y', 1)
pptopo = PostponedTopology(X, Y)

The PostponedTopology supports boundaries, interfaces and groups, returning PostponedTopology objects, and integrals, returning PostponedIntegral objects. Upon evaluation of a PostponedIntegral a Topology object with the same roots as the PostponedTopology object should be specified. Example:

res = pptopo.integral(..., degree=1)
res.eval(LineTopo(X)*LineTopo(Y))

The function.BasisArgument(name) object is a function.Array with unspecified coefficients and dofs, to be specified upon evaluation by assigning a function.Basis object to name via evalargs, similar to function.Argument. The namespace gains notation for creating a BasisArgument: #name_n is notation for BasisArgument('name'). Example:

ns = Namespace()
ns.p = '#basis_n ?p_n'
ns.u_i = '#basis_n ?u_ni'

Hierarchical example:

ns = Namespace()
ns.u = '#basis_n ?lhs_n'
ns.prevu = '#prevbasis_n ?lhs_n'
pptopo = PostponedTopology(X, Y)
proj = pptopo.integral('(u - uprev)^2 d:x', degree=...)
res = pptopo.integral('#basis_n,i u_,i d:x', degree=...)
topo = LineTopo(X)*LineTopo(Y)
basis = topo.basis(...)
lhs = ...
while ...:
  lhs = newton('lhs', res, topo, lhs0=lhs, evalargs=dict(basis=basis)).solve()
  if ...:
    break
  topo = topo.refined_by(...)
  prevbasis, basis = basis, topo.basis(...)
  lhs = optimize('lhs', proj, evalargs=dict(basis=basis, prevbasis=prevbasis, prevlhs=lhs))

This paves the way for a general purpose goal adaptivity helper function that performs the solves, refinements and projections based on an initial topology, a residual, a goal and a callback for generating the required bases.