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