comparison test/Grids/mapped_grid_test.jl @ 1777:819ab806960f feature/grids/manifolds

Merge feature/grids/curvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 15 Sep 2024 23:00:15 +0200
parents 03894fd7b132 ecec2b0eea0f
children 2b5f81e288f1
comparison
equal deleted inserted replaced
1750:e2abd72d7ce0 1777:819ab806960f
35 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) 35 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
36 mg = MappedGrid(lg, x̄, J) 36 mg = MappedGrid(lg, x̄, J)
37 37
38 @test mg isa Grid{SVector{2, Float64},2} 38 @test mg isa Grid{SVector{2, Float64},2}
39 @test jacobian(mg) isa Array{<:AbstractMatrix} 39 @test jacobian(mg) isa Array{<:AbstractMatrix}
40 @test logicalgrid(mg) isa Grid 40 @test logical_grid(mg) isa Grid
41 41
42 @test collect(mg) == x̄ 42 @test collect(mg) == x̄
43 @test jacobian(mg) == J 43 @test jacobian(mg) == J
44 @test logicalgrid(mg) == lg 44 @test logical_grid(mg) == lg
45 45
46 46
47 x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg) 47 x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
48 J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg) 48 J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
49 mg = MappedGrid(lg, x̄, J) 49 mg = MappedGrid(lg, x̄, J)
50 50
51 @test mg isa Grid{SVector{3, Float64},2} 51 @test mg isa Grid{SVector{3, Float64},2}
52 @test jacobian(mg) isa Array{<:AbstractMatrix} 52 @test jacobian(mg) isa Array{<:AbstractMatrix}
53 @test logicalgrid(mg) isa Grid 53 @test logical_grid(mg) isa Grid
54 54
55 @test collect(mg) == x̄ 55 @test collect(mg) == x̄
56 @test jacobian(mg) == J 56 @test jacobian(mg) == J
57 @test logicalgrid(mg) == lg 57 @test logical_grid(mg) == lg
58 58
59 sz1 = (10,11) 59 sz1 = (10,11)
60 sz2 = (10,12) 60 sz2 = (10,12)
61 @test_throws ArgumentError("Sizes must match") MappedGrid( 61 @test_throws ArgumentError("Sizes must match") MappedGrid(
62 equidistant_grid((0,0), (1,1), sz2...), 62 equidistant_grid((0,0), (1,1), sz2...),
245 0; 245 0;
246 1+ξ*(ξ-1); 246 1+ξ*(ξ-1);
247 ] 247 ]
248 248
249 function expected_bg(mg, bId, Jb) 249 function expected_bg(mg, bId, Jb)
250 lg = logicalgrid(mg) 250 lg = logical_grid(mg)
251 return MappedGrid( 251 return MappedGrid(
252 boundary_grid(lg, bId), 252 boundary_grid(lg, bId),
253 map(x̄, boundary_grid(lg, bId)), 253 map(x̄, boundary_grid(lg, bId)),
254 map(Jb, boundary_grid(lg, bId)), 254 map(Jb, boundary_grid(lg, bId)),
255 ) 255 )
277 x̄, J = _partially_curved_mapping() 277 x̄, J = _partially_curved_mapping()
278 mg = mapped_grid(x̄, J, 10, 11) 278 mg = mapped_grid(x̄, J, 10, 11)
279 @test mg isa MappedGrid{SVector{2,Float64}, 2} 279 @test mg isa MappedGrid{SVector{2,Float64}, 2}
280 280
281 lg = equidistant_grid((0,0), (1,1), 10, 11) 281 lg = equidistant_grid((0,0), (1,1), 10, 11)
282 @test logicalgrid(mg) == lg 282 @test logical_grid(mg) == lg
283 @test collect(mg) == map(x̄, lg) 283 @test collect(mg) == map(x̄, lg)
284 284
285 @test mapped_grid(lg, x̄, J) == mg 285 @test mapped_grid(lg, x̄, J) == mg
286 end
287
288 @testset "jacobian_determinant" begin
289 x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
290 J((ξ, η)) = @SMatrix[
291 η ξ;
292 1 2η;
293 ]
294
295 g = mapped_grid(x̄, J, 10, 11)
296 J = map(logicalgrid(g)) do (ξ,η)
297 2η^2 - ξ
298 end
299 @test jacobian_determinant(g) ≈ J
300
301
302 lg = equidistant_grid((0,0), (1,1), 11, 21)
303 x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
304 J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
305 mg = MappedGrid(lg, x̄, J)
306
307 @test_broken jacobian(mg) isa AbstractArray{2,Float64}
308 end 286 end
309 287
310 @testset "metric_tensor" begin 288 @testset "metric_tensor" begin
311 x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2] 289 x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
312 J((ξ, η)) = @SMatrix[ 290 J((ξ, η)) = @SMatrix[
313 η ξ; 291 η ξ;
314 1 2η; 292 1 2η;
315 ] 293 ]
316 294
317 g = mapped_grid(x̄, J, 10, 11) 295 g = mapped_grid(x̄, J, 10, 11)
318 G = map(logicalgrid(g)) do (ξ,η) 296 G = map(logical_grid(g)) do (ξ,η)
319 @SMatrix[ 297 @SMatrix[
320 1+η^2 ξ*η+2η; 298 1+η^2 ξ*η+2η;
321 ξ*η+2η ξ^2 + 4η^2; 299 ξ*η+2η ξ^2 + 4η^2;
322 ] 300 ]
323 end 301 end
330 1+ξ 0; 308 1+ξ 0;
331 ξ 1+η; 309 ξ 1+η;
332 ] 310 ]
333 311
334 g = mapped_grid(x̄, J, 10, 11) 312 g = mapped_grid(x̄, J, 10, 11)
335 G⁻¹ = map(logicalgrid(g)) do (ξ,η) 313 G⁻¹ = map(logical_grid(g)) do (ξ,η)
336 @SMatrix[ 314 @SMatrix[
337 (1+η)^2 -ξ*(1+η); 315 (1+η)^2 -ξ*(1+η);
338 -ξ*(1+η) (1+ξ)^2+ξ^2; 316 -ξ*(1+η) (1+ξ)^2+ξ^2;
339 ]/(((1+ξ)^2+ξ^2)*(1+η)^2 - ξ^2*(1+η)^2) 317 ]/(((1+ξ)^2+ξ^2)*(1+η)^2 - ξ^2*(1+η)^2)
340 318
382 g = mapped_grid(_partially_curved_mapping()...,10, 11) 360 g = mapped_grid(_partially_curved_mapping()...,10, 11)
383 361
384 @test normal(g, CartesianBoundary{1,LowerBoundary}()) == fill(@SVector[-1,0], 11) 362 @test normal(g, CartesianBoundary{1,LowerBoundary}()) == fill(@SVector[-1,0], 11)
385 @test normal(g, CartesianBoundary{1,UpperBoundary}()) == fill(@SVector[1,0], 11) 363 @test normal(g, CartesianBoundary{1,UpperBoundary}()) == fill(@SVector[1,0], 11)
386 @test normal(g, CartesianBoundary{2,LowerBoundary}()) == fill(@SVector[0,-1], 10) 364 @test normal(g, CartesianBoundary{2,LowerBoundary}()) == fill(@SVector[0,-1], 10)
387 @test normal(g, CartesianBoundary{2,UpperBoundary}()) ≈ map(boundary_grid(g,CartesianBoundary{2,UpperBoundary}())|>logicalgrid) do ξ̄ 365 @test normal(g, CartesianBoundary{2,UpperBoundary}()) ≈ map(boundary_grid(g,CartesianBoundary{2,UpperBoundary}())|>logical_grid) do ξ̄
388 α = 1-2ξ̄[1] 366 α = 1-2ξ̄[1]
389 @SVector[α,1]/√(α^2 + 1) 367 @SVector[α,1]/√(α^2 + 1)
390 end 368 end
391 369
392 g = mapped_grid(_fully_curved_mapping()...,5,4) 370 g = mapped_grid(_fully_curved_mapping()...,5,4)
393 371
394 unit(v) = v/norm(v) 372 unit(v) = v/norm(v)
395 @testset let bId = CartesianBoundary{1,LowerBoundary}() 373 @testset let bId = CartesianBoundary{1,LowerBoundary}()
396 lbg = boundary_grid(logicalgrid(g), bId) 374 lbg = boundary_grid(logical_grid(g), bId)
397 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 375 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
398 -unit(@SVector[1/2, η/3-1/6]) 376 -unit(@SVector[1/2, η/3-1/6])
399 end 377 end
400 end 378 end
401 379
402 @testset let bId = CartesianBoundary{1,UpperBoundary}() 380 @testset let bId = CartesianBoundary{1,UpperBoundary}()
403 lbg = boundary_grid(logicalgrid(g), bId) 381 lbg = boundary_grid(logical_grid(g), bId)
404 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 382 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
405 unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2)) 383 unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2))
406 end 384 end
407 end 385 end
408 386
409 @testset let bId = CartesianBoundary{2,LowerBoundary}() 387 @testset let bId = CartesianBoundary{2,LowerBoundary}()
410 lbg = boundary_grid(logicalgrid(g), bId) 388 lbg = boundary_grid(logical_grid(g), bId)
411 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 389 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
412 -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ)) 390 -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ))
413 end 391 end
414 end 392 end
415 393
416 @testset let bId = CartesianBoundary{2,UpperBoundary}() 394 @testset let bId = CartesianBoundary{2,UpperBoundary}()
417 lbg = boundary_grid(logicalgrid(g), bId) 395 lbg = boundary_grid(logical_grid(g), bId)
418 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 396 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
419 unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ)) 397 unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ))
420 end 398 end
421 end 399 end
422 end 400 end