Gridap modalC0

[ioannisPApapadopoulos]

@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,n)

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