1using Test, LibCEED, LinearAlgebra, StaticArrays 2 3include("buildmats.jl") 4 5function checkoutput(str, fname) 6 if str != getoutput(fname) 7 write(fname, str) 8 return false 9 end 10 return true 11end 12 13@testset "LibCEED Development Tests" begin 14 @testset "Basis" begin 15 @test BasisNone()[] == LibCEED.C.CEED_BASIS_NONE[] 16 17 c = Ceed() 18 dim = 3 19 ncomp = 1 20 p = 4 21 q = 6 22 b1 = create_tensor_h1_lagrange_basis(c, dim, ncomp, p, q, GAUSS_LOBATTO) 23 24 @test checkoutput(showstr(b1), "b1.out") 25 26 b1d = CeedScalar[1.0 0.0; 0.5 0.5; 0.0 1.0] 27 d1d = CeedScalar[-0.5 0.5; -0.5 0.5; -0.5 0.5] 28 q1d = CeedScalar[-1.0, 0.0, 1.0] 29 w1d = CeedScalar[1/3, 4/3, 1/3] 30 q, p = size(b1d) 31 d2d = zeros(CeedScalar, 2, q*q, p*p) 32 d2d[1, :, :] = kron(b1d, d1d) 33 d2d[2, :, :] = kron(d1d, b1d) 34 35 dim2 = 2 36 b2 = create_tensor_h1_basis(c, dim2, 1, p, q, b1d, d1d, q1d, w1d) 37 @test checkoutput(showstr(b2), "b2.out") 38 39 b3 = create_h1_basis( 40 c, 41 LINE, 42 1, 43 p, 44 q, 45 b1d, 46 reshape(d1d, 1, q, p), 47 reshape(q1d, 1, q), 48 w1d, 49 ) 50 @test checkoutput(showstr(b3), "b3.out") 51 52 dim = 2 53 ncomp = 1 54 p1 = 4 55 q1 = 4 56 qref1 = Array{Float64}(undef, dim, q1) 57 qweight1 = Array{Float64}(undef, q1) 58 interp1, div1 = build_mats_hdiv(qref1, qweight1) 59 b1 = create_hdiv_basis(c, QUAD, ncomp, p1, q1, interp1, div1, qref1, qweight1) 60 61 u1 = ones(Float64, p1) 62 v1 = apply(b1, u1) 63 64 for i = 1:q1 65 @test v1[i] ≈ -1.0 66 @test v1[q1+i] ≈ 1.0 67 end 68 69 p2 = 3 70 q2 = 4 71 qref2 = Array{Float64}(undef, dim, q2) 72 qweight2 = Array{Float64}(undef, q2) 73 interp2, curl2 = build_mats_hcurl(qref2, qweight2) 74 b2 = create_hcurl_basis(c, TRIANGLE, ncomp, p2, q2, interp2, curl2, qref2, qweight2) 75 76 u2 = [1.0, 2.0, 1.0] 77 v2 = apply(b2, u2) 78 79 for i = 1:q2 80 @test v2[i] ≈ 1.0 81 end 82 83 u2[1] = -1.0 84 u2[2] = 1.0 85 u2[3] = 2.0 86 v2 = apply(b2, u2) 87 88 for i = 1:q2 89 @test v2[q2+i] ≈ 1.0 90 end 91 end 92 93 @testset "ElemRestriction" begin 94 c = Ceed() 95 nelem = 3 96 elemsize = 2 97 offsets = Array{CeedInt}(undef, elemsize, nelem) 98 orients = Array{Bool}(undef, elemsize, nelem) 99 for i = 1:nelem 100 offsets[1, i] = i - 1 101 offsets[2, i] = i 102 # flip the dofs on element 1, 3, ... 103 orients[1, i] = (i - 1)%2 > 0 104 orients[2, i] = (i - 1)%2 > 0 105 end 106 r = create_elem_restriction_oriented( 107 c, 108 nelem, 109 elemsize, 110 1, 111 1, 112 nelem + 1, 113 offsets, 114 orients, 115 ) 116 117 lv = Vector{CeedScalar}(undef, nelem + 1) 118 for i = 1:nelem+1 119 lv[i] = 10 + i - 1 120 end 121 122 ev = apply(r, lv) 123 124 for i = 1:nelem 125 for j = 1:elemsize 126 k = j + elemsize*(i - 1) 127 @test 10 + k÷2 == ev[k]*(-1)^((i - 1)%2) 128 end 129 end 130 131 curlorients = Array{CeedInt8}(undef, 3*elemsize, nelem) 132 for i = 1:nelem 133 curlorients[1, i] = curlorients[6, i] = 0 134 if (i - 1)%2 > 0 135 # T = [0 -1] 136 # [-1 0] 137 curlorients[2, i] = 0 138 curlorients[3, i] = -1 139 curlorients[4, i] = -1 140 curlorients[5, i] = 0 141 else 142 # T = I 143 curlorients[2, i] = 1 144 curlorients[3, i] = 0 145 curlorients[4, i] = 0 146 curlorients[5, i] = 1 147 end 148 end 149 r = create_elem_restriction_curl_oriented( 150 c, 151 nelem, 152 elemsize, 153 1, 154 1, 155 nelem + 1, 156 offsets, 157 curlorients, 158 ) 159 160 ev = apply(r, lv) 161 162 for i = 1:nelem 163 for j = 1:elemsize 164 k = j + elemsize*(i - 1) 165 if (i - 1)%2 > 0 166 @test j == 2 || 10 + i == -ev[k] 167 @test j == 1 || 10 + i - 1 == -ev[k] 168 else 169 @test 10 + k÷2 == ev[k] 170 end 171 end 172 end 173 end 174end 175