1using Test, LibCEED, LinearAlgebra, StaticArrays 2 3include("buildmats.jl") 4 5showstr(x) = sprint(show, MIME("text/plain"), x) 6summarystr(x) = sprint(summary, x) 7getoutput(fname) = 8 chomp(read(joinpath(@__DIR__, "output", string(CeedScalar), fname), String)) 9 10function checkoutput(str, fname) 11 if str != getoutput(fname) 12 write(fname, str) 13 return false 14 end 15 return true 16end 17 18mutable struct CtxData 19 io::IOBuffer 20 x::Vector{Float64} 21end 22 23const run_dev_tests = !isrelease() || ("--run-dev-tests" in ARGS) 24 25if run_dev_tests 26 include("rundevtests.jl") 27end 28 29if !LibCEED.ceedversion_ge(LibCEED.minimum_libceed_version) && !run_dev_tests 30 @warn "Skipping tests because of incompatible libCEED versions." 31else 32 @testset "LibCEED Release Tests" begin 33 @testset "LibCEED" begin 34 @test ceedversion() isa VersionNumber 35 @test isrelease() isa Bool 36 @test isfile(get_libceed_path()) 37 end 38 39 @testset "Ceed" begin 40 res = "/cpu/self/ref/serial" 41 c = Ceed(res) 42 @test isdeterministic(c) 43 @test getresource(c) == res 44 @test !iscuda(c) 45 @test get_preferred_memtype(c) == MEM_HOST 46 @test_throws LibCEED.CeedError create_interior_qfunction(c, "") 47 @test showstr(c) == """ 48 Ceed 49 Ceed Resource: $res 50 Preferred MemType: host""" 51 end 52 53 @testset "Context" begin 54 c = Ceed() 55 data = zeros(CeedScalar, 3) 56 ctx = Context(c, data) 57 @test showstr(ctx) == """ 58 CeedQFunctionContext 59 Context Data Size: $(sizeof(data))""" 60 @test_throws Exception set_data!(ctx, MEM_HOST, OWN_POINTER, data) 61 end 62 63 @testset "CeedVector" begin 64 n = 10 65 c = Ceed() 66 v = CeedVector(c, n) 67 @test size(v) == (n,) 68 @test length(v) == n 69 @test axes(v) == (1:n,) 70 @test ndims(v) == 1 71 @test ndims(CeedVector) == 1 72 73 v[] = 0.0 74 @test @witharray(a = v, all(a .== 0.0)) 75 76 v1 = rand(CeedScalar, n) 77 v2 = CeedVector(c, v1) 78 @test @witharray_read(a = v2, mtype = MEM_HOST, a == v1) 79 @test Vector(v2) == v1 80 v[] = v1 81 for p ∈ [1, 2, Inf] 82 @test norm(v, p) ≈ norm(v1, p) 83 end 84 @test_throws Exception norm(v, 3) 85 @test witharray_read(sum, v) == sum(v1) 86 reciprocal!(v) 87 @test @witharray(a = v, mtype = MEM_HOST, all(a .== CeedScalar(1.0)./v1)) 88 89 witharray(x -> x .= 1.0, v) 90 @test @witharray(a = v, all(a .== 1.0)) 91 92 @test summarystr(v) == "$n-element CeedVector" 93 @test sprint(show, v) == @witharray_read(a = v, sprint(show, a)) 94 io = IOBuffer() 95 summary(io, v) 96 println(io, ":") 97 @witharray_read(a = v, Base.print_array(io, a)) 98 s1 = String(take!(io)) 99 @test showstr(v) == s1 100 101 setarray!(v, MEM_HOST, USE_POINTER, v1) 102 syncarray!(v, MEM_HOST) 103 @test @witharray_read(a = v, a == v1) 104 p = takearray!(v, MEM_HOST) 105 @test p == pointer(v1) 106 107 m = rand(CeedScalar, 10, 10) 108 vm = CeedVector(c, vec(m)) 109 @test @witharray_read(a = vm, size = size(m), a == m) 110 111 @test CeedVectorActive()[] == LibCEED.C.CEED_VECTOR_ACTIVE[] 112 @test CeedVectorNone()[] == LibCEED.C.CEED_VECTOR_NONE[] 113 114 w1 = rand(CeedScalar, n) 115 w2 = rand(CeedScalar, n) 116 w3 = rand(CeedScalar, n) 117 118 cv1 = CeedVector(c, w1) 119 cv2 = CeedVector(c, w2) 120 cv3 = CeedVector(c, w3) 121 122 alpha = rand(CeedScalar) 123 124 scale!(cv1, alpha) 125 w1 .*= alpha 126 @test @witharray_read(a = cv1, a == w1) 127 128 pointwisemult!(cv1, cv2, cv3) 129 w1 .= w2.*w3 130 @test @witharray_read(a = cv1, a == w1) 131 132 axpy!(alpha, cv2, cv1) 133 axpy!(alpha, w2, w1) 134 @test @witharray_read(a = cv1, a ≈ w1) 135 end 136 137 @testset "Basis" begin 138 c = Ceed() 139 dim = 3 140 ncomp = 1 141 p = 4 142 q = 6 143 b1 = create_tensor_h1_lagrange_basis(c, dim, ncomp, p, q, GAUSS_LOBATTO) 144 145 @test getdimension(b1) == 3 146 @test gettopology(b1) == HEX 147 @test getnumcomponents(b1) == ncomp 148 @test getnumnodes(b1) == p^dim 149 @test getnumnodes1d(b1) == p 150 @test getnumqpts(b1) == q^dim 151 @test getnumqpts1d(b1) == q 152 153 q1d, w1d = lobatto_quadrature(3, AbscissaAndWeights) 154 @test q1d ≈ CeedScalar[-1.0, 0.0, 1.0] 155 @test w1d ≈ CeedScalar[1/3, 4/3, 1/3] 156 157 q1d, w1d = gauss_quadrature(3) 158 @test q1d ≈ CeedScalar[-sqrt(3/5), 0.0, sqrt(3/5)] 159 @test w1d ≈ CeedScalar[5/9, 8/9, 5/9] 160 161 b1d = CeedScalar[1.0 0.0; 0.5 0.5; 0.0 1.0] 162 d1d = CeedScalar[-0.5 0.5; -0.5 0.5; -0.5 0.5] 163 q1d = CeedScalar[-1.0, 0.0, 1.0] 164 w1d = CeedScalar[1/3, 4/3, 1/3] 165 q, p = size(b1d) 166 d2d = zeros(CeedScalar, 2, q*q, p*p) 167 d2d[1, :, :] = kron(b1d, d1d) 168 d2d[2, :, :] = kron(d1d, b1d) 169 170 dim2 = 2 171 b2 = create_tensor_h1_basis(c, dim2, 1, p, q, b1d, d1d, q1d, w1d) 172 @test getinterp(b2) == kron(b1d, b1d) 173 @test getinterp1d(b2) == b1d 174 @test getgrad(b2) == d2d 175 @test getgrad1d(b2) == d1d 176 177 b3 = create_h1_basis( 178 c, 179 LINE, 180 1, 181 p, 182 q, 183 b1d, 184 reshape(d1d, 1, q, p), 185 reshape(q1d, 1, q), 186 w1d, 187 ) 188 @test getqref(b3) == reshape(q1d, 1, q) 189 @test getqweights(b3) == w1d 190 191 v = rand(CeedScalar, 2) 192 vq = apply(b3, v) 193 vd = apply(b3, v; emode=EVAL_GRAD) 194 @test vq ≈ b1d*v 195 @test vd ≈ d1d*v 196 end 197 198 @testset "Request" begin 199 @test RequestImmediate()[] == LibCEED.C.CEED_REQUEST_IMMEDIATE[] 200 @test RequestOrdered()[] == LibCEED.C.CEED_REQUEST_ORDERED[] 201 end 202 203 @testset "Misc" begin 204 for dim = 1:3 205 D = CeedDim(dim) 206 J = rand(CeedScalar, dim, dim) 207 @test det(J, D) ≈ det(J) 208 J = J + J' # make symmetric 209 @test setvoigt(SMatrix{dim,dim}(J)) == setvoigt(J, D) 210 @test getvoigt(setvoigt(J, D)) == J 211 V = zeros(CeedScalar, dim*(dim + 1)÷2) 212 setvoigt!(V, J, D) 213 @test V == setvoigt(J, D) 214 J2 = zeros(CeedScalar, dim, dim) 215 getvoigt!(J2, V, D) 216 @test J2 == J 217 end 218 end 219 220 @testset "Operator" begin 221 c = Ceed() 222 @interior_qf id = ( 223 c, 224 (input, :in, EVAL_INTERP), 225 (output, :out, EVAL_INTERP), 226 begin 227 output[] = input 228 end, 229 ) 230 b = create_tensor_h1_lagrange_basis(c, 3, 1, 3, 3, GAUSS_LOBATTO) 231 n = getnumnodes(b) 232 offsets = Vector{CeedInt}(0:n-1) 233 r = create_elem_restriction(c, 1, n, 1, 1, n, offsets) 234 op = Operator( 235 c; 236 qf=id, 237 fields=[ 238 (:input, r, b, CeedVectorActive()), 239 (:output, r, b, CeedVectorActive()), 240 ], 241 ) 242 243 v = rand(CeedScalar, n) 244 v1 = CeedVector(c, v) 245 v2 = CeedVector(c, n) 246 apply!(op, v1, v2) 247 @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 == a2)) 248 apply_add!(op, v1, v2) 249 @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 + a1 == a2)) 250 251 diag_vector = create_lvector(r) 252 LibCEED.assemble_diagonal!(op, diag_vector) 253 @test @witharray_read(a = diag_vector, a == ones(n)) 254 # TODO: change this test after bug-fix in libCEED 255 diag_vector[] = 0.0 256 LibCEED.assemble_add_diagonal!(op, diag_vector) 257 @test @witharray(a = diag_vector, a == fill(1.0, n)) 258 259 @test showstr(op) == """ 260 CeedOperator 261 1 elements with 27 quadrature points each 262 2 fields 263 1 input field: 264 Input field 0: 265 Name: "input" 266 Size: 1 267 EvalMode: interpolation 268 Active vector 269 1 output field: 270 Output field 0: 271 Name: "output" 272 Size: 1 273 EvalMode: interpolation 274 Active vector""" 275 end 276 277 @testset "ElemRestriction" begin 278 c = Ceed() 279 n = 10 280 offsets = Vector{CeedInt}([0:n-1; n-1:2*n-2]) 281 lsize = 2*n - 1 282 r = create_elem_restriction(c, 2, n, 1, lsize, lsize, offsets) 283 @test getcompstride(r) == lsize 284 @test getnumelements(r) == 2 285 @test getelementsize(r) == n 286 @test getlvectorsize(r) == lsize 287 @test getnumcomponents(r) == 1 288 @test length(create_lvector(r)) == lsize 289 @test length(create_evector(r)) == 2*n 290 lv, ev = create_vectors(r) 291 @test length(lv) == lsize 292 @test length(ev) == 2*n 293 mult = getmultiplicity(r) 294 mult2 = ones(lsize) 295 mult2[n] = 2 296 @test mult == mult2 297 rand_lv = rand(CeedScalar, lsize) 298 rand_ev = [rand_lv[1:n]; rand_lv[n:end]] 299 @test apply(r, rand_lv) == rand_ev 300 @test apply(r, rand_ev; tmode=TRANSPOSE) == rand_lv.*mult 301 @test showstr(r) == string( 302 "CeedElemRestriction from (19, 1) to 2 elements ", 303 "with 10 nodes each and component stride 19", 304 ) 305 306 strides = CeedInt[1, n, n] 307 rs = create_elem_restriction_strided(c, 1, n, 1, n, strides) 308 @test showstr(rs) == string( 309 "CeedElemRestriction from (10, 1) to 1 elements ", 310 "with 10 nodes each and strides [1, $n, $n]", 311 ) 312 313 @test ElemRestrictionNone()[] == LibCEED.C.CEED_ELEMRESTRICTION_NONE[] 314 end 315 316 @testset "QFunction" begin 317 c = Ceed() 318 @test showstr(create_interior_qfunction(c, "Poisson3DApply")) == """ 319 Gallery CeedQFunction - Poisson3DApply 320 2 input fields: 321 Input field 0: 322 Name: "du" 323 Size: 3 324 EvalMode: "gradient" 325 Input field 1: 326 Name: "qdata" 327 Size: 6 328 EvalMode: "none" 329 1 output field: 330 Output field 0: 331 Name: "dv" 332 Size: 3 333 EvalMode: "gradient\"""" 334 335 id = create_identity_qfunction(c, 1, EVAL_INTERP, EVAL_INTERP) 336 Q = 10 337 v = rand(CeedScalar, Q) 338 v1 = CeedVector(c, v) 339 v2 = CeedVector(c, Q) 340 apply!(id, Q, [v1], [v2]) 341 @test @witharray(a = v2, a == v) 342 343 @interior_qf id2 = (c, (a, :in, EVAL_INTERP), (b, :out, EVAL_INTERP), b .= a) 344 v2[] = 0.0 345 apply!(id2, Q, [v1], [v2]) 346 @test @witharray(a = v2, a == v) 347 348 ctxdata = CtxData(IOBuffer(), rand(CeedScalar, 3)) 349 ctx = Context(c, ctxdata) 350 dim = 3 351 @interior_qf qf = ( 352 c, 353 dim=dim, 354 ctxdata::CtxData, 355 (a, :in, EVAL_GRAD, dim), 356 (b, :in, EVAL_NONE), 357 (c, :out, EVAL_INTERP), 358 begin 359 c[] = b*sum(a) 360 show(ctxdata.io, MIME("text/plain"), ctxdata.x) 361 end, 362 ) 363 set_context!(qf, ctx) 364 in_sz, out_sz = LibCEED.get_field_sizes(qf) 365 @test in_sz == [dim, 1] 366 @test out_sz == [1] 367 v1 = rand(CeedScalar, dim) 368 v2 = rand(CeedScalar, 1) 369 cv1 = CeedVector(c, v1) 370 cv2 = CeedVector(c, v2) 371 cv3 = CeedVector(c, 1) 372 apply!(qf, 1, [cv1, cv2], [cv3]) 373 @test String(take!(ctxdata.io)) == showstr(ctxdata.x) 374 @test @witharray_read(v3 = cv3, v3[1] == v2[1]*sum(v1)) 375 376 @test QFunctionNone()[] == LibCEED.C.CEED_QFUNCTION_NONE[] 377 end 378 @testset "Basis" begin 379 @test BasisNone()[] == LibCEED.C.CEED_BASIS_NONE[] 380 381 c = Ceed() 382 dim = 3 383 ncomp = 1 384 p = 4 385 q = 6 386 b1 = create_tensor_h1_lagrange_basis(c, dim, ncomp, p, q, GAUSS_LOBATTO) 387 388 @test checkoutput(showstr(b1), "b1.out") 389 390 b1d = CeedScalar[1.0 0.0; 0.5 0.5; 0.0 1.0] 391 d1d = CeedScalar[-0.5 0.5; -0.5 0.5; -0.5 0.5] 392 q1d = CeedScalar[-1.0, 0.0, 1.0] 393 w1d = CeedScalar[1/3, 4/3, 1/3] 394 q, p = size(b1d) 395 d2d = zeros(CeedScalar, 2, q*q, p*p) 396 d2d[1, :, :] = kron(b1d, d1d) 397 d2d[2, :, :] = kron(d1d, b1d) 398 399 dim2 = 2 400 b2 = create_tensor_h1_basis(c, dim2, 1, p, q, b1d, d1d, q1d, w1d) 401 @test checkoutput(showstr(b2), "b2.out") 402 403 b3 = create_h1_basis( 404 c, 405 LINE, 406 1, 407 p, 408 q, 409 b1d, 410 reshape(d1d, 1, q, p), 411 reshape(q1d, 1, q), 412 w1d, 413 ) 414 @test checkoutput(showstr(b3), "b3.out") 415 416 dim = 2 417 ncomp = 1 418 p1 = 4 419 q1 = 4 420 qref1 = Array{Float64}(undef, dim, q1) 421 qweight1 = Array{Float64}(undef, q1) 422 interp1, div1 = build_mats_hdiv(qref1, qweight1) 423 b1 = create_hdiv_basis(c, QUAD, ncomp, p1, q1, interp1, div1, qref1, qweight1) 424 425 u1 = ones(Float64, p1) 426 v1 = apply(b1, u1) 427 428 for i = 1:q1 429 @test v1[i] ≈ -1.0 430 @test v1[q1+i] ≈ 1.0 431 end 432 433 p2 = 3 434 q2 = 4 435 qref2 = Array{Float64}(undef, dim, q2) 436 qweight2 = Array{Float64}(undef, q2) 437 interp2, curl2 = build_mats_hcurl(qref2, qweight2) 438 b2 = create_hcurl_basis( 439 c, 440 TRIANGLE, 441 ncomp, 442 p2, 443 q2, 444 interp2, 445 curl2, 446 qref2, 447 qweight2, 448 ) 449 450 u2 = [1.0, 2.0, 1.0] 451 v2 = apply(b2, u2) 452 453 for i = 1:q2 454 @test v2[i] ≈ 1.0 455 end 456 457 u2[1] = -1.0 458 u2[2] = 1.0 459 u2[3] = 2.0 460 v2 = apply(b2, u2) 461 462 for i = 1:q2 463 @test v2[q2+i] ≈ 1.0 464 end 465 end 466 467 @testset "ElemRestriction" begin 468 c = Ceed() 469 nelem = 3 470 elemsize = 2 471 offsets = Array{CeedInt}(undef, elemsize, nelem) 472 orients = Array{Bool}(undef, elemsize, nelem) 473 for i = 1:nelem 474 offsets[1, i] = i - 1 475 offsets[2, i] = i 476 # flip the dofs on element 1, 3, ... 477 orients[1, i] = (i - 1)%2 > 0 478 orients[2, i] = (i - 1)%2 > 0 479 end 480 r = create_elem_restriction_oriented( 481 c, 482 nelem, 483 elemsize, 484 1, 485 1, 486 nelem + 1, 487 offsets, 488 orients, 489 ) 490 491 lv = Vector{CeedScalar}(undef, nelem + 1) 492 for i = 1:nelem+1 493 lv[i] = 10 + i - 1 494 end 495 496 ev = apply(r, lv) 497 498 for i = 1:nelem 499 for j = 1:elemsize 500 k = j + elemsize*(i - 1) 501 @test 10 + k÷2 == ev[k]*(-1)^((i - 1)%2) 502 end 503 end 504 505 curlorients = Array{CeedInt8}(undef, 3*elemsize, nelem) 506 for i = 1:nelem 507 curlorients[1, i] = curlorients[6, i] = 0 508 if (i - 1)%2 > 0 509 # T = [0 -1] 510 # [-1 0] 511 curlorients[2, i] = 0 512 curlorients[3, i] = -1 513 curlorients[4, i] = -1 514 curlorients[5, i] = 0 515 else 516 # T = I 517 curlorients[2, i] = 1 518 curlorients[3, i] = 0 519 curlorients[4, i] = 0 520 curlorients[5, i] = 1 521 end 522 end 523 r = create_elem_restriction_curl_oriented( 524 c, 525 nelem, 526 elemsize, 527 1, 528 1, 529 nelem + 1, 530 offsets, 531 curlorients, 532 ) 533 534 ev = apply(r, lv) 535 536 for i = 1:nelem 537 for j = 1:elemsize 538 k = j + elemsize*(i - 1) 539 if (i - 1)%2 > 0 540 @test j == 2 || 10 + i == -ev[k] 541 @test j == 1 || 10 + i - 1 == -ev[k] 542 else 543 @test 10 + k÷2 == ev[k] 544 end 545 end 546 end 547 end 548 end 549end 550