Following on from my issue in ClassicalOrthogonalPolynomials..

My goal is to apply a tensor-product plan_transform(ContinuousPolynomial{0}(r), (Block(p), Block(p)) to multiple vectors simultaneously. Focusing on the 1D case first, it seems to me that the infrastructure is not completely there yet, right?

From the tests I can see that the following works:

r, p = range(0,1,3), 3
P = ContinuousPolynomial{0}(r)
x, F = plan_grid_transform(P, (Block(p),2), 1)
c = F*[exp.(x) ;;; cos.(x)]

but then F \ c errors.

Do you think the goal is attainable @dlfivefifty ?

Hi @dlfivefifty. Do you have a plan to implement support for a different truncation degree on each cell? In particular for the plan_transform of ContinuousPolynomial{0}.

I feel since the transform on each cell are independent it should not be so difficult but I am not sure whether one expects things to be really slow if the transforms on each cell have a different degree.

@dlfivefifty Does the @assert in https://github.com/JuliaApproximation/PiecewiseOrthogonalPolynomials.jl/blob/5c9cc9448b25c62b5b3b5301f12262dca864be29/src/arrowhead.jl#L281C5-L281C5 need to be there?

The matrices for the disk do not fit the exact format for reversecholesky currently implemented in PiecewisePolys and I error at this assert.

@dlfivefifty I think there is a bug in Derivative(x) * ContinuousPolynomial{1}(r). E.g.

n = 5; r = range(0, 1; length=n+1) # solve on unit interval with n cells

C¹ = ContinuousPolynomial{1}(r) # hat functions and (1-x^2) * P_k^(1,1)(x)
C⁰ = ContinuousPolynomial{0}(r)

x = axes(C,1)
D = C⁰ \ (Derivative(x) * C¹)

# D[1,1] should be -5 not -0.5

I think M = BlockVcat(Hcat(Ones{T}(N) / 2, -Ones{T}(N) / 2), H) in line 241 of PiecewiseOrthogonalPolynomials.jl should be

M = BlockVcat(Hcat(Ones{T}(N) * (N-1), -Ones{T}(N) * (N-1)), H)

@ioannisPApapadopoulos @KarsKnook

I've been playing with QL instead of reverse Cholesky for non-SPD operators. It suffers from fill in in the top row of blocks so the complexity would be $O(np+n^3)$.

BUT if we did strong form (so non-symmetric matrices) I believe everything apart from the first block column will be diagonal. Then QL will be optimal complexity.

Note what we need is ContinuousPolynomial{-1} which would be P^(1,1) instead of Legendre but also with delta functions between the elements.

@DanielVandH @dlfivefifty I thought the following info might be useful:

Gridap has an implementation of ContinuousPolynomial{1} called modalC0 although the implementation isn't aware of the sparsity pattern of the mass/weak Laplacian matrices and I assume their transforms/quadrature schemes are not (quasi)optimal. It has also been implemented in 2D and 3D:

1D:

using Gridap, SparseArrays, Plots

domain = (0,1) # unit interval domain
n = (5,) # 5 cells
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
p = 6
reffe_u = ReferenceFE(modalC0,Float64, p)

Vu = TestFESpace(model,reffe_u,labels=labels,dirichlet_tags="boundary",conformity=:H1)
Uu = TrialFESpace(Vu, 0.0)

Ω = Triangulation(model)
dΩ = Measure(Ω,p+3)

a(uh, v) =∫( ∇(uh) ⋅ ∇(v)) * dΩ
jac(uh, duh, v) = ∫( ∇(duh) ⋅ ∇(v)) * dΩ
op = FEOperator(a, jac, Uu, Vu)

u0 = ones(Float64,num_free_dofs(Vu))
uh = FEFunction(Uu,u0);
J  = Gridap.Algebra.jacobian(op, uh) # J is the weak Laplacian matrix

# Reveal actual sparsity pattern
J[abs.(J) .< 1e-10] .=0 # reveal
J = sparse(Matrix(J))

Plots.spy(J)

To look at the 2D version, change the domain and n to

domain = (0,1,0,1) # unit square domain
n = (5,5) # 25 cells

I think the only thing that's not trivial is the A block. The D blocks should just be changing step(r) to diff(r)

Duh! 🙈

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!