module MultiGridBarrier

Module MultiGridBarrier solves convex optimization problems in function spaces, for example, solving p-Laplace problems. We recommend to start with the functions fem1d_solve(), fem2d_solve(), spectral1d_solve(), spectral2d_solve(). These functions are sufficient to solve p-Laplace problems in 1d or 2d, using finite or spectral elements.

For more general use, the user will need to familiarize themselves with the basic ideas of convex optimization.

• Overview of convex optimization in function spaces by MultiGrid Barrier method.

The general idea is to build a multigrid hierarchy, represented by an AMG object, and barrier for a convex set, represented by a Barrier object, and then solve a convex optimization problem using the amgb() solver.

To generate the multigrid hierarchy represented by the AMG object, use either fem1d(), fem2d(), spectral1d() or spectral2d() functions. These constructors will assemble suitable AMG objects for either FEM or spectral discretizations, in 1d or 2d. One should think of these four constructors as being specialized in constructing some specific function spaces. A user can use the amg() constructor directly if custom function spaces are required, but this is more difficult.

We now describe the barrier function.

Assume that $\Omega \subset \mathbb{R}^d$ is some open set. Consider the example of the p-Laplace problem on $\Omega$. Let $f(x)$ be a "forcing" (a function) on $\Omega$, and $1 \leq p < \infty$. One wishes to solve the minimization problem

\[\begin{equation} \inf_u \int_{\Omega} fu + \|\nabla u\|_2^p \, dx. \tag{1} \end{equation}\]

Generally speaking, $u$ will range in some function space, e.g. a space of differentiable functions satisfying homogeneous Dirichlet conditions. Under some conditions, minimizing (1) is equivalent to solving the p-Laplace PDE:

\[\nabla \cdot (\|\nabla u\|_2^{p-2}\nabla u) = {1 \over p} f.\]

We introduct the "slack function" $s(x)$ and replace (1) with the following equivalent problem:

\[\begin{equation} \inf_{s(x) \geq \|\nabla u(x)\|_2^p} \int_{\Omega} fu + s \, dx. \tag{2} \end{equation}\]

Define the convex set $\mathcal{Q} = \{ (u(x),q(x),s(x)) \; : \; s(x) \geq \|q(x)\|_2^p \}$, and

\[z = \begin{bmatrix} u \\ s \end{bmatrix}, \qquad c^T = [f,0,1], \qquad Dz = \begin{bmatrix} u \\ \nabla u \\ s \end{bmatrix}.\]

Then, (2) can be rewritten as

\[\begin{equation} \inf_{Dz \in \mathcal{Q}} \int_{\Omega} c^T(x)Dz(x) \, dx. \tag{3} \end{equation}\]

Recall that a barrier for $\mathcal{Q}$ is a convex function $\mathcal{F}$ on $\mathcal{Q}$ such that $\mathcal{F} < \infty$ in the interior of $\mathcal{Q}$ and $\mathcal{F} = \infty$ on the boundary of $\mathcal{Q}$. A barrier for the p-Laplace problem is:

\[\mathcal{F}(u,q,s) = \int_{\Omega} -\log(s^{2 \over p} - \|q\|_2^2) - 2\log s \, dx = \int_{\Omega} F(Dz(x)) \, dx.\]

The central path $z^*(t)$ minimizes, for each fixed $t>0$, the quantity

\[\int_{\Omega} tc^TDz + F(Dz) \, dx.\]

As $t \to \infty$, $z^*(t)$ forms a minimizing sequence (or filter) for (3). We think of the function $c(x)$ as the "functional" that we seek to minimize.

The Convex{T} type describes various convex sets (denoted $Q$ above) by way of functions barrier(), cobarrier() and slack(). barrier is indeed a barrier for $Q$, cobarrier() is a barrier for a related feasibility problems, and $slack()$ is used in solving the feasibility problem. Convex{T} objects can be created using the various convex_...() constructors, e.g. convex_Euclidian_power() for the p-Laplace problem.

Once one has AMG and Convex objects, and a suitable "functional" c, one uses the amgb() function to solve the optimization problem by the MultiGrid Barrier method, a variant of the barrier method (or interior point method) that is quasi-optimal for sufficiently regular problems.

@kwdef struct AMG{T,M,Geometry}
    ...
end

Objects of this type should probably be assembled by the constructor amg().

A multigrid with L level. Denote by l between 1 and L, a grid level. Fields are:

• x::Matrix{T} the vertices of the fine grid.
• w::Vector{T} corresponding quadrature weights.
• R_fine::Array{M,1} an array of L matrices. The columns of R_fine[l] are basis functions for the function space on grid level l, interpolated to the fine grid.
• R_coarse::Array{M,1} an array of L matrices. The columns of R_coarse[l] are basis functions for the function space on grid level l. Unlike R_fine[l], these basis functions are on grid level l, not interpolated to the fine grid.
• D::Array{M,2} an array of differential operators. For example, if the barrier parameters are to be u,ux,s, with ux the derivative of u, then D[l,:] = [I,Dx,I], where Dx is a numerical differentiation operator on grid level l.
• refine_u::Array{M,1} an array of L grid refinement matrices. If x[l] has n[l] vertices, then refine_u[l] is n[l+1] by n[l].
• coarsen_u::Array{M,1} an array of L grid coarsening matrices. coarsen_u[l] is n[l] by n[l+1].
• refine_z::Array{M,1} an array of L grid refining matrices for the "state vector" z. For example, if z contains the state functions u and s, then there are k=2 state functions, and refine_z[l] is k*n[l+1] by k*n[l].
• coarsen_z::Array{M,1} an array of L grid coarsening matrices for the "state vector" z. coarsen_z[l] is k*n[l] by k*n[l+1].

These various matrices must satisfy a wide variety of algebraic relations. For this reason, it is recommended to use the constructor amg().

Barrier

A type for holding barrier functions. Fields are:

f0::Function
f1::Function
f2::Function

f0 is the barrier function itself, while f1 is its gradient and f2 is the Hessian.

struct Convex
    barrier::Function
    cobarrier::Function
    slack::Function
end

The Convex data structure represents a convex domain $Q$ implicitly by way of three functions. The barrier function is a barrier for $Q$. cobarrier is a barrier for the feasibility subproblem, and slack is a function that initializes a valid slack value for the feasibility subproblem. The various convex_ functions can be used to generate various convex domains.

These function are called as follows: barrier(x,y). x is a vertex in a grid, as per the AMG object. y is some vector. For each fixed x variable, y -> barrier(x,y) defines a barrier for a convex set in y.

abstract type FEM1D end

abstract type FEM2D end

abstract type SPECTRAL1D end

abstract type SPECTRAL2D end

function amg(::Type{Geometry};
    x::Matrix{T},
    w::Vector{T},
    state_variables::Matrix{Symbol},
    D::Matrix{Symbol},
    subspaces::Dict{Symbol,Vector{M}},
    operators::Dict{Symbol,M},
    refine::Vector{M},
    coarsen::Vector{M},
    full_space=:full,
    id_operator=:id,
    feasibility_slack=:feasibility_slack,
    generate_feasibility=true) where {T,M,Geometry}

Construct an AMG object for use with the amgb solver. In many cases, this constructor is not called directly by the user. For 1d and 2d finite elements, use the fem1d() or fem2d(). For 1d and 2d spectral elements, use spectral1d() or spectral2d(). You use amg() directly if you are implementing your own function spaces.

The AMG object shall represent all L grid levels of the multigrid hierarchy. Parameters are:

• x: the vertices of the fine grid.
• w: the quadrature weights for the fine grid.
• state_variables: a matrix of symbols. The first column indicates the names of the state vectors or functions, and the second column indicates the names of the corresponding subspaces. A typical example is: state_variables = [:u :dirichlet; :s :full]. This would define the solution as being functions named u(x) and s(x). The u function would lie in the space :dirichlet, presumably consisting of functions with homogeneous Dirichlet conditions. The s function would lie in the space :full, presumably being the full function space, without boundary conditions.
• D: a matrix of symbols. The first column indicates the names of various state variables, and the second column indicates the corresponding differentiation operator(s). For example: D = [:u :id ; :u :dx ; :s :id]. This would indicate that the barrier should be called as F(x,y) with y = [u,ux,s], where ux denotes the derivative of u with respect to the space variable x.
• subspaces: a Dict mapping each subspace symbol to an array of L matrices, e.g. for each l, subspaces[:dirichlet][l] is a matrix whose columns span the homogeneous Dirichlet subspace of grid level l.
• operators: a Dict mapping each differential operator symbol to a matrix, e.g. operators[:id] is an identity matrix, while operators[:dx] is a numerical differentiation matrix, on the fine grid level L.
• refine: an array of length L of matrices. For each l, refine[l] interpolates from grid level l to grid level l+1. refine[L] should be the identity, and coarsen[l]*refine[l] should be the identity.
• coarsen: an array of length L of matrices. For each l, coarsen[l] interpolates or projects from grid level l+1 to grid level l. coarsen[L] should be the identity.
• generate_feasibility: if true, amg() returns a pair M of AMG objects. M[1] is an AMG object to be used for the main optimization problem, while M[2] is an AMG object for the preliminary feasibility sub problem. In this case, amg() also needs to be provided with the following additional information: feasibility_slack is the name of a special slack variable that must be unique to the feasibility subproblem (default: :feasibility_slack); full_space is the name of the "full" vector space (i.e. no boundary conditions, default: :full); and id_operator is the name of the identity operator (default: :id).

amg_construct(::Type{T},::Type{FEM1D};rest...) where {T} = fem1d(T;rest...)

amg_construct(::Type{T},::Type{FEM2D};rest...) where {T} = fem2d(T;rest...)

amg_construct(::Type{T},::Type{SPECTRAL1D};rest...) where {T} = spectral1d(T;rest...)

amg_construct(::Type{T},::Type{SPECTRAL2D},L,n,K) where {T} = spectral2d(T,n=n,L=L)

amg_dim(::Type{FEM1D}) = 1

amg_dim(::Type{FEM2D}) = 2

amg_dim(::Type{SPECTRAL1D}) = 1

amg_dim(::Type{SPECTRAL2D}) = 2

amg_plot(M::AMG{T,Mat,FEM1D}, z::Vector{T}) where {T,Mat} = plot(M.x[end],z)

function amg_plot(M::AMG{T, Mat,FEM2D}, z::Array{T}) where {T,Mat}

Plot a piecewise quadratic (plus cubic "bubble") solution z on the given mesh. Note that the solution is drawn as (linear) triangles, even though the underlying solution is piecewise cubic. To obtain a more accurate depiction, especially when the mesh is coarse, it would be preferable to apply a few levels of additional subdivision, so as to capture the curve of the quadratic basis functions.

function amg_plot(M::AMG{T,Mat,SPECTRAL1D},y;x=Array(-1:T(0.01):1),rest...) where {T,Mat}

Plot a solution using pyplot.

• M: a mesh.
• x: x values where the solution should be evaluated and plotted.
• y: the solution, to be interpolated at the given x values via spectral1d_interp.
• rest... parameters are passed directly to pyplot.plot.

function amg_plot(M::AMG{T,Mat,SPECTRAL2D},z::Array{T,1};x=-1:T(0.01):1,y=-1:T(0.01):1,rest...) where {T,Mat}

Plot a 2d solution.

• M a 2d mesh.
• x, y should be ranges like -1:0.01:1.
• z the solution to plot.

function amg_solve(::Type{T}=Float64; 
    L::Integer=2, n=nothing,
    method=FEM1D,
    K = nothing,
    state_variables::Matrix{Symbol} = [:u :dirichlet
                       :s :full],
    dim::Integer = amg_dim(method),
    D::Matrix{Symbol} = default_D[dim],
    M = amg_construct(T,method,L=L,n=n,K=K,state_variables=state_variables,D=D),
    p::T = T(1.0),
    g::Union{Function,Matrix{T}} = default_g(T)[dim],
    f::Union{Function,Matrix{T}} = default_f(T)[dim],
    Q::Convex{T} = convex_Euclidian_power(T,idx=2:dim+2,p=x->p),
    show=true,         
    return_details=false, rest...) where {T}

A simplified interface for module MultiGridBarrier to "quickly get started". To solve a p-Laplace problem, do: amg_solve().

Different behaviors can be obtained by supplying various optional keyword arguments, as follows.

• L=2: the number of times to subdivide the base mesh.
• The n parameter is only used in spectral methods, in which case, if n is an integer, then L is disregarded. n is the number of quadrature nodes along each axis.
• method=FEM1D: this must be either FEM1D, FEM2D, SPECTRAL1D or SPECTRAL2D. This parameter is used twice: once to choose the constructor for the M parameter, and again to plot the solution if show is true. If show is false and if M is constructed "manually", not by its default value, then the method parameter is ignored.
• K: In most cases, this is nothing, but in the fem2d case, K is the initial mesh.
• state_variables = [:u :dirichlet ; :s :full]: the names of the components of the solution vector z.
• dim = size(M[1].x[end],2), the dimension of the problem, should be 1 or 2. This is only used in the default values of the g, f, Q, D parameters, and is ignored if these parameters do not use default values.
• D: the differential operators, see below for defaults.
• M: a mesh obtained by one of the constructors fem1d, fem2d, spectral1d or spectral2d, corresponding to the method parameter.
• x = M[1].x: a matrix, same number of rows as M[1].x. This matrix will be passed, row by row, to the barrier function, as the x parameter.
• p = T(1.0): the parameter of the p-Laplace operator. This is only relevant if the default value is used for the Q parameter, and is ignored otherwise.
• g: the "boundary conditions" function. See below for defaults.
• f: the "forcing" or "cost functional" to be minimized. See below for defaults.
• Q: the convex domain to which the solution should belong. Defaults to convex_Euclidian_power(T,idx=2:dim+2,p=x->p), which corresponds to p-Laplace problems. Change this to solve other variational problems.
• show=true: if true, plot the solution.
• return_details=false: if false, return a Vector{T} of the solution. If true, returned a named tuple with some more details about the solution process.
• rest...: any further keyword arguments are passed on to amgb.

The default values for the parameters f, g, D are as follows

function amgb(M::Tuple{AMG{T,Mat,Geometry},AMG{T,Mat,Geometry}},
          f::Union{Function,Matrix{T}}, 
          g::Union{Function,Matrix{T}}, 
          Q::Convex;
          x::Matrix{T} = M[1].x,
          t=T(0.1),
          t_feasibility=t,
          verbose=true,
          return_details=false,
          rest...) where {T,Mat,Geometry}

A thin wrapper around amgb_core(). Parameters are:

• M: obtained from the amg constructor, a pair of AMG structures. M[1] is the main problem while M[2] is the feasibility problem.
• f: the functional to minimize.
• g: the "boundary conditions".
• Q: a Convex domain for the convex optimization problem.
• rest...: any further named arguments are passed on to amgb_core.

The initial z0 guess, and the cost functional c0, are computed as follows:

m = size(M[1].x,1)
for k=1:m
    z0[k,:] .= g(M[1].x[k,:])
    c0[k,:] .= f(M[1].x[k,:])
end

By default, the return value z is an m×n matrix, where n is the number of state_variables, see either fem1d(), fem2d(), spectral1d() or spectral2d(). If return_details=true then the return value is a named tuple with fields z, SOL_feasibility and SOL_main; the latter two fields are named tuples with detailed information regarding the various solves.

function amgb_core(B::Barrier,
    M::AMG{T,Mat,Geometry},
    x::Matrix{T},
    z::Array{T,1},
    c::Array{T,2};
    tol=(eps(T)),
    t=T(0.1),
    maxit=10000,
    kappa=T(10.0),
    early_stop=z->false,
    verbose=true) where {T,Mat,Geometry}

The "Algebraic MultiGrid Barrier" method.

• B a Barrier object.
• M an AMG object.
• x a matrix with the same number of rows as M.x. This is passed as the x parameter of the barrier. Often, x = M.x.
• z a starting point for the minimization, which should be admissible, i.e. B.f0(z)<∞.
• c an objective functional to minimize. Concretely, we minimize the integral of c.*(D*z), as computed by the finest quadrature in M, subject to B.f0(z)<∞. Here, D is the differential operator provided in M.

Optional parameters:

• t: the initial value of t
• tol: we stop when 1/t<tol.
• maxit: the maximum number of t steps.
• kappa: the initial size of the t-step. Stepsize adaptation is used in the AMGB algorithm, where the t-step size may be made smaller or large, but it will never exceed the initial size provided here.
• verbose: set to true to see a progress bar.
• early_stop: if early_stop(z) is true then the minimization is stopped early. This is used when solving the preliminary feasibility problem.

Return value is a named tuple SOL with the following fields:

• SOL.converged is true if convergence was obtained, else it is false.
• SOL.z the computed solution.

Further SOL fields contain various statistics about the solve process.

function barrier(F;
    F1=(x,y)->ForwardDiff.gradient(z->F(x,z),y),
    F2=(x,y)->ForwardDiff.hessian(z->F(x,z),y))::Barrier

Constructor for barriers.

• F is the actual barrier function. It should take parameters (x,y).
• F1 is the gradient of F with respect to y.
• F2 is the Hessian of F with respect to y.

By default, F1 and F2 are automatically generated by the module ForwardDiff.

A more specific description of the Barrier object is as follows. The function Barrier.f0 has parameters:

function Barrier.f0(z,x,w,c,R,D,z0)

Here, R is a matrix and D is an array of matrices; x is a matrix of quadrature nodes with weights w, and c is a matrix describing the functional we seek to minimize. The value of Barrier.f0 is given by:

        p = length(w)
        n = length(D)
        Rz = z0+R*z
        Dz = hcat([D[k]*Rz for k=1:n]...)
        y = [F(x[k,:],Dz[k,:]) for k=1:p]
        dot(w,y)+sum([dot(w.*c[:,k],Dz[:,k]) for k=1:n])

Thus, Barrier.f0 can be regarded as a quadrature approximation of the integral

\[\int_{\Omega} \left(\sum_{k=1}^nc_k(x)v_k(x)\right) + F(x,v_1(x),\ldots,v_n(x)) \, dx \text{ where } v_k = D_k(z_0 + Rz).\]

Functions Barrier.f1 and Barrier.f2 are the gradient and Hessian, respectively, of Barrier.f0, with respect to the z parameter. If the underlying matrices are sparse, then sparse arithmetic is used for Barrier.f2.

function convex_Euclidian_power(::Type{T}=Float64;idx=Colon(),A::Function=(x)->I,b::Function=(x)->T(0),p::Function=x->T(2)) where {T}

Generate a Convex object corresponding to the convex set defined by $z[end] \geq \|z[1:end-1]\|_2^p$ where $z = A(x)*y[idx] .+ b(x)$.

function convex_linear(::Type{T}=Float64;idx=Colon(),A::Function=(x)->I,b::Function=(x)->T(0)) where {T}

Generate a Convex structure corresponding to the convex domain A(x,k)*y[idx] .+ b(x,k) ≤ 0.

function fem1d(::Type{T}=Float64; L::Int=4,
                state_variables = [:u :dirichlet
                                   :s :full],
                D = [:u :id
                     :u :dx
                     :s :id],
                generate_feasibility=true) where {T}

Construct an AMG object for a 1d piecewise linear finite element grid. The interval is [-1,1]. Parameters are:

• L: divide the interval into 2^L subintervals (L for Levels).
• state_variables: the "state vector" consists of functions, by default this is u(x) and s(x), on the finite element grid.
• D: the set of differential operator. The barrier function F will eventually be called with the parameters F(x,Dz), where z is the state vector. By default, this results in F(x,u,ux,s), where ux is the derivative of u.
• generate_feasibility: if true, returns a pair M of AMG objects. M[1] is the AMG object for the main problem, and M[2] is for the feasibility subproblem.

function fem1d_interp(x::Vector{T},
                  y::Vector{T},
                  t::T) where{T}

Interpolate a 1d piecewise linear function at the given t value. If u(xi) is the piecewise linear function such that u(x[k])=y[k] then this function returns u(t).

function fem1d_interp(x::Vector{T},
                  y::Vector{T},
                  t::Vector{T}) where{T}

Returns [fem1d_interp(x,y,t[k]) for k=1:length(t)].

fem1d_solve(::Type{T}=Float64;rest...) where {T} = amg_solve(T;method=FEM1D,rest...)

function fem2d(::Type{T}, L::Int, K::Matrix{T};
                state_variables = [:u :dirichlet
                                   :s :full],
                D = [:u :id
                     :u :dx
                     :u :dy
                     :s :id],
                generate_feasibility=true) where {T}

Construct an AMG object for a 2d finite element grid on the domain K with piecewise quadratic elements. Parameters are:

• K: a triangular mesh. If there are n triangles, then K should be a 3n by 2 matrix of vertices. The first column of K represents x coordinates, the second column represents y coordinates.
• L: divide the interval into 2^L subintervals (L for Levels).
• state_variables: the "state vector" consists of functions, by default this is u(x) and s(x), on the finite element grid.
• D: the set of differential operator. The barrier function F will eventually be called with the parameters F(x,y,Dz), where z is the state vector. By default, this results in F(x,y,u,ux,uy,s), where (ux,uy) is the gradient of u.
• generate_feasibility: if true, returns a pair M of AMG objects. M[1] is the AMG object for the main problem, and M[2] is for the feasibility subproblem.

fem2d_solve(::Type{T}=Float64;rest...) where {T} = amg_solve(T;method=FEM2D,rest...)

function illinois(f,a::T,b::T;fa=f(a),fb=f(b),maxit=10000) where {T}

Find a root of f between a and b using the Illinois algorithm. If f(a)*f(b)>=0, returns b.

function newton(::Type{Mat},
                   F0::Function,
                   F1::Function,
                   F2::Function,
                   x::Array{T,1};
                   maxit=10000,
                   theta=T(0.1),
                   beta=T(0.1),
                   tol=eps(T)) where {T,Mat}

Damped Newton iteration for minimizing a function.

• F0 the objective function
• F1 and F2 are the gradient and Hessian of F0, respectively.
• x the starting point of the minimization procedure.

The Hessian F2 return value should be of type Mat.

The optional parameters are:

• maxit, the iteration aborts with a failure message if convergence is not achieved within maxit iterations.
• tol is used as a stopping criterion. We stop when the decrement in the objective is sufficiently small.

function parabolic_plot(method,M::AMG{T, Mat,Geometry}, U::Matrix{T}; 
    interval=200, embed_limit=200.0,
    printer=(animation)->display("text/html", animation.to_html5_video(embed_limit=embed_limit))) where {T,Mat,Geometry}

Animate the solution of the parabolic problem.

function parabolic_solve(::Type{T}=Float64;
    method = FEM2D,
    state_variables = [:u  :dirichlet
                       :s1 :full
                       :s2 :full],
    dim = amg_dim(method),
    f1 = x->T(0.5),
    f_default = default_f_parabolic[dim],
    p = T(1),
    h = T(0.2),
    f = (t,x)->f_default(h*f1(x)-x[1+dim],T(0.5),h/p),
    g = default_g_parabolic[dim],
    D = default_D_parabolic[dim],
    L = 2,
    t0 = T(0),
    t1 = T(1),
    M = amg_construct(T,method,L=L,D=D,state_variables=state_variables),
    Q = (convex_Euclidian_power(;idx=[1,2+dim],p=x->T(2)) 
        ∩ convex_Euclidian_power(;idx=vcat(2:1+dim,3+dim),p=x->p)),
    verbose = true,
    show = true,
    interval = 200,
    printer=(animation)->display("text/html", animation.to_html5_video(embed_limit=200.0)),
    rest...) where {T}

Solves a parabolic (i.e. time-dependent) p-Laplace problem of the form:

\[u_t - \nabla \cdot (\|\nabla u\|_2^{p-2}\nabla u) = -f_1.\]

We use the implicit Euler scheme $u_t \approx (u_{k+1}-u_k)/h$ to arrive at:

\[u_{k+1} - h\nabla \cdot (\|\nabla u_{k+1}\|_2^{p-2}\nabla u_{k+1}) = u_k-hf_1.\]

According to the calculus of variation, we look for a weak solution minimizing

\[J(u) = \int_{\Omega}{1 \over 2} u^2 + h {1 \over p} \|\nabla u\|_2^p + (hf_1-u_k)u \, dx\]

We introduce the slack functions $s_1 \geq u^2$ and $s_2 \geq \|\nabla u\|_2^p$ and minimize instead

\[\int_{\Omega} {1 \over 2}s_1 + {h \over p} s_2 + (hf_1-u_k)u \, dx.\]

The canonical form is:

\[z = \begin{bmatrix} u \\ s_1 \\ s_2 \end{bmatrix} \qquad f^T = \left[hf_1-u_k,0,0,{1 \over 2},{h \over p}\right] \qquad Dz = \begin{bmatrix} u \\ u_x \\ u_y \\ s_1 \\ s_2 \end{bmatrix} \qquad g = \begin{bmatrix} g_1 \\ 0 \\ 0 \end{bmatrix}.\]

Here, $g_1$ encodes boundary conditions for $u$. Then we minimize:

\[\int_{\Omega} f^TDz\]

The named arguments rest... are passed verbatim to amg_solve.

function spectral1d(::Type{T}=Float64; n::Integer=5,
                state_variables = [:u :dirichlet
                                   :s :full],
                D = [:u :id
                     :u :dx
                     :s :id],
                generate_feasibility=true) where {T}

Construct an AlgebraicMultiGridBarrier.AMG object for a 1d spectral grid of polynomials of degree n-1. See also fem1d for a description of the parameters state_variables and D.

function spectral1d_interp(MM::AMG{T,Mat,SPECTRAL1D}, y::Array{T,1},x) where {T,Mat}

A function to interpolate a solution y at some point(s) x.

• MM the mesh of the solution.
• y the solution.
• x point(s) at which the solution should be evaluated.

spectral1d_solve(::Type{T}=Float64;rest...) where {T} = amg_solve(T;method=SPECTRAL1D,rest...)

function spectral2d(::Type{T}=Float64; n=5::Integer,
                state_variables = [:u :dirichlet
                                   :s :full],
                D = [:u :id
                     :u :dx
                     :u :dy
                     :s :id],
                generate_feasibility=true) where {T}

Construct an AMG object for a 2d spectral grid of degree n-1. See also fem2d for a description of state_variables and D.

function spectral2d_interp(MM::AMG{T,Mat,SPECTRAL2D},z::Array{T,1},x::Array{T,2}) where {T,Mat}

Interpolate a solution z at point(s) x, given the mesh MM. See also spectral1d_interp.

spectral2d_solve(::Type{T}=Float64;rest...) where {T} = amg_solve(T;method=SPECTRAL2D,rest...)