Add Julia API tests

julia/LibCEED.jl/src/Basis.jl

@@ -126,7 +126,7 @@ Create a non tensor-product basis for $H^1$ discretizations
              at quadrature points.
 - `grad`:    Array of size `(dim, nqpts, nnodes)` expressing derivatives of nodal basis
              functions at quadrature points.
-- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+- `qref`:    Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the
              reference element $[-1, 1]$.
 - `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
              element.
@@ -145,12 +145,13 @@ function create_h1_basis(
     dim = getdimension(topo)
     @assert size(interp) == (nqpts, nnodes)
     @assert size(grad) == (dim, nqpts, nnodes)
-    @assert length(qref) == nqpts
+    @assert size(qref) == (dim, nqpts)
     @assert length(qweight) == nqpts
 
     # Convert from Julia matrices and tensors (column-major) to row-major format
     interp_rowmajor = collect(interp')
     grad_rowmajor = permutedims(grad, [3, 2, 1])
+    qref_rowmajor = collect(qref')
 
     ref = Ref{C.CeedBasis}()
     C.CeedBasisCreateH1(
@@ -161,7 +162,7 @@ function create_h1_basis(
         nqpts,
         interp_rowmajor,
         grad_rowmajor,
-        qref,
+        qref_rowmajor,
         qweight,
         ref,
     )
@@ -183,7 +184,7 @@ Create a non tensor-product basis for H(div) discretizations
              at quadrature points.
 - `div`:     Array of size `(nqpts, nnodes)` expressing divergence of basis functions at
              quadrature points.
-- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+- `qref`:    Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the
              reference element $[-1, 1]$.
 - `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
              element.
@@ -202,12 +203,13 @@ function create_hdiv_basis(
     dim = getdimension(topo)
     @assert size(interp) == (dim, nqpts, nnodes)
     @assert size(div) == (nqpts, nnodes)
-    @assert length(qref) == nqpts
+    @assert size(qref) == (dim, nqpts)
     @assert length(qweight) == nqpts
 
     # Convert from Julia matrices and tensors (column-major) to row-major format
     interp_rowmajor = permutedims(interp, [3, 2, 1])
     div_rowmajor = collect(div')
+    qref_rowmajor = collect(qref')
 
     ref = Ref{C.CeedBasis}()
     C.CeedBasisCreateHdiv(
@@ -218,7 +220,7 @@ function create_hdiv_basis(
         nqpts,
         interp_rowmajor,
         div_rowmajor,
-        qref,
+        qref_rowmajor,
         qweight,
         ref,
     )
@@ -240,7 +242,7 @@ Create a non tensor-product basis for H(curl) discretizations
              at quadrature points.
 - `curl`:    Matrix of size `(curlcomp, nqpts, nnodes)`, `curlcomp = 1 if dim < 3 else dim`)
              matrix expressing curl of basis functions at quadrature points.
-- `qref`:    Array of length `nqpts` holding the locations of quadrature points on the
+- `qref`:    Matrix of size `(dim, nqpts)` holding the locations of quadrature points on the
              reference element $[-1, 1]$.
 - `qweight`: Array of length `nqpts` holding the quadrature weights on the reference
              element.
@@ -260,12 +262,13 @@ function create_hcurl_basis(
     curlcomp = dim < 3 ? 1 : dim
     @assert size(interp) == (dim, nqpts, nnodes)
     @assert size(curl) == (curlcomp, nqpts, nnodes)
-    @assert length(qref) == nqpts
+    @assert size(qref) == (dim, nqpts)
     @assert length(qweight) == nqpts
 
     # Convert from Julia matrices and tensors (column-major) to row-major format
     interp_rowmajor = permutedims(interp, [3, 2, 1])
     curl_rowmajor = permutedims(curl, [3, 2, 1])
+    qref_rowmajor = collect(qref')
 
     ref = Ref{C.CeedBasis}()
     C.CeedBasisCreateHcurl(
@@ -276,7 +279,7 @@ function create_hcurl_basis(
         nqpts,
         interp_rowmajor,
         curl_rowmajor,
-        qref,
+        qref_rowmajor,
         qweight,
         ref,
     )
@@ -328,10 +331,9 @@ function apply(b::Basis, u::AbstractVector; nelem=1, tmode=NOTRANSPOSE, emode=EV
 
     u_vec = CeedVector(c, u)
 
-    len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : getnumqpts(b)
-    if emode == EVAL_GRAD
-        len_v *= getdimension(b)
-    end
+    qcomp = Ref{CeedInt}()
+    C.CeedBasisGetNumQuadratureComponents(b[], emode, qcomp)
+    len_v = (tmode == TRANSPOSE) ? getnumnodes(b) : qcomp[]*getnumqpts(b)
 
     v_vec = CeedVector(c, len_v)
 
@@ -437,11 +439,8 @@ function getqref(b::Basis)
     ref = Ref{Ptr{CeedScalar}}()
     C.CeedBasisGetQRef(b[], ref)
     copy(
-        unsafe_wrap(
-            Array,
-            ref[],
-            istensor[] ? getnumqpts1d(b) : (getnumqpts(b)*getdimension(b)),
-        ),
+        istensor[] ? unsafe_wrap(Array, ref[], getnumqpts1d(b)) :
+        unsafe_wrap(Array, ref[], (getnumqpts(b), getdimension(b)))',
     )
 end
 

julia/LibCEED.jl/src/LibCEED.jl

@@ -79,6 +79,8 @@ export @interior_qf,
     create_elem_restriction_strided,
     create_evector,
     create_h1_basis,
+    create_hdiv_basis,
+    create_hcurl_basis,
     create_identity_qfunction,
     create_interior_qfunction,
     create_lvector,

julia/LibCEED.jl/test/buildmats.jl

@@ -0,0 +1,78 @@
+function build_mats_hdiv(qref, qweight, ::Type{T}=Float64) where {T}
+    P, Q, dim = 4, 4, 2
+    interp = Array{T}(undef, dim, Q, P)
+    div = Array{T}(undef, Q, P)
+
+    qref[1, 1] = -1.0/sqrt(3.0)
+    qref[1, 2] = qref[1, 1]
+    qref[1, 3] = qref[1, 1]
+    qref[1, 4] = -qref[1, 1]
+    qref[2, 1] = -qref[1, 1]
+    qref[2, 2] = -qref[1, 1]
+    qref[2, 3] = qref[1, 1]
+    qref[2, 4] = qref[1, 1]
+    qweight[1] = 1.0
+    qweight[2] = 1.0
+    qweight[3] = 1.0
+    qweight[4] = 1.0
+
+    # Loop over quadrature points
+    for i = 1:Q
+        x1 = qref[1, i]
+        x2 = qref[2, i]
+        # Interp
+        interp[1, i, 1] = 0.0
+        interp[2, i, 1] = 1.0 - x2
+        interp[1, i, 2] = x1 - 1.0
+        interp[2, i, 2] = 0.0
+        interp[1, i, 3] = -x1
+        interp[2, i, 3] = 0.0
+        interp[1, i, 4] = 0.0
+        interp[2, i, 4] = x2
+        # Div
+        div[i, 1] = -1.0
+        div[i, 2] = 1.0
+        div[i, 3] = -1.0
+        div[i, 4] = 1.0
+    end
+
+    return interp, div
+end
+
+function build_mats_hcurl(qref, qweight, ::Type{T}=Float64) where {T}
+    P, Q, dim = 3, 4, 2
+    interp = Array{T}(undef, dim, Q, P)
+    curl = Array{T}(undef, 1, Q, P)
+
+    qref[1, 1] = 0.2
+    qref[1, 2] = 0.6
+    qref[1, 3] = 1.0/3.0
+    qref[1, 4] = 0.2
+    qref[2, 1] = 0.2
+    qref[2, 2] = 0.2
+    qref[2, 3] = 1.0/3.0
+    qref[2, 4] = 0.6
+    qweight[1] = 25.0/96.0
+    qweight[2] = 25.0/96.0
+    qweight[3] = -27.0/96.0
+    qweight[4] = 25.0/96.0
+
+    # Loop over quadrature points
+    for i = 1:Q
+        x1 = qref[1, i]
+        x2 = qref[2, i]
+        # Interp
+        interp[1, i, 1] = -x2
+        interp[2, i, 1] = x1
+        interp[1, i, 2] = x2
+        interp[2, i, 2] = 1.0 - x1
+        interp[1, i, 3] = 1.0 - x2
+        interp[2, i, 3] = x1
+        # Curl
+        curl[1, i, 1] = 2.0
+        curl[1, i, 2] = -2.0
+        curl[1, i, 3] = -2.0
+    end
+
+    return interp, curl
+end

julia/LibCEED.jl/test/rundevtests.jl

@@ -1,3 +1,48 @@
 using Test, LibCEED, LinearAlgebra, StaticArrays
 
-@testset "LibCEED Development Tests" begin end
+include("buildmats.jl")
+
+@testset "LibCEED Development Tests" begin
+    @testset "Basis" begin
+        c = Ceed()
+        dim = 2
+        ncomp = 1
+        p1 = 4
+        q1 = 4
+        qref1 = Array{Float64}(undef, dim, q1)
+        qweight1 = Array{Float64}(undef, q1)
+        interp1, div1 = build_mats_hdiv(qref1, qweight1)
+        b1 = create_hdiv_basis(c, QUAD, ncomp, p1, q1, interp1, div1, qref1, qweight1)
+
+        u1 = ones(Float64, p1)
+        v1 = apply(b1, u1)
+
+        for i = 1:q1
+            @test v1[i] ≈ -1.0
+            @test v1[q1+i] ≈ 1.0
+        end
+
+        p2 = 3
+        q2 = 4
+        qref2 = Array{Float64}(undef, dim, q2)
+        qweight2 = Array{Float64}(undef, q2)
+        interp2, curl2 = build_mats_hcurl(qref2, qweight2)
+        b2 = create_hcurl_basis(c, TRIANGLE, ncomp, p2, q2, interp2, curl2, qref2, qweight2)
+
+        u2 = [1.0, 2.0, 1.0]
+        v2 = apply(b2, u2)
+
+        for i = 1:q2
+            @test v2[i] ≈ 1.0
+        end
+
+        u2[1] = -1.0
+        u2[2] = 1.0
+        u2[3] = 2.0
+        v2 = apply(b2, u2)
+
+        for i = 1:q2
+            @test v2[q2+i] ≈ 1.0
+        end
+    end
+end

julia/LibCEED.jl/test/runtests.jl

@@ -174,7 +174,17 @@ else
             @test getgrad1d(b2) == d1d
             @test checkoutput(showstr(b2), "b2.out")
 
-            b3 = create_h1_basis(c, LINE, 1, p, q, b1d, reshape(d1d, 1, q, p), q1d, w1d)
+            b3 = create_h1_basis(
+                c,
+                LINE,
+                1,
+                p,
+                q,
+                b1d,
+                reshape(d1d, 1, q, p),
+                reshape(q1d, 1, q),
+                w1d,
+            )
             @test getqref(b3) == q1d
             @test getqweights(b3) == w1d
             @test checkoutput(showstr(b3), "b3.out")