comparison src/SbpOperators/volumeops/laplace/laplace.jl @ 668:2d56a53a1646 feature/laplace_opset

Simplify construction of boundary-operator pairs in Laplace constructor
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 31 Jan 2021 22:19:53 +0100
parents f3a0d1f7d842
children 011863b3f24c
comparison
equal deleted inserted replaced
667:f3a0d1f7d842 668:2d56a53a1646
34 # Volume operators 34 # Volume operators
35 Δ = laplace(grid, D_inner_stecil, D_closure_stencils) 35 Δ = laplace(grid, D_inner_stecil, D_closure_stencils)
36 H = DiagonalQuadrature(grid, H_closure_stencils) 36 H = DiagonalQuadrature(grid, H_closure_stencils)
37 H⁻¹ = InverseDiagonalQuadrature(grid, H_closure_stencils) 37 H⁻¹ = InverseDiagonalQuadrature(grid, H_closure_stencils)
38 38
39 # Pair boundary operators and boundary quadratures with the boundary ids 39 # Pair operators with boundary ids
40 e_pairs = () 40 bids = boundary_identifiers(grid)
41 d_pairs = () 41 # Boundary operators
42 Hᵧ_pairs = () 42 e_pairs = ntuple(i -> Pair(bids[i],BoundaryRestriction(grid,e_closure_stencil,bids[i])),length(bids))
43 d_pairs = ntuple(i -> Pair(bids[i],NormalDerivative(grid,d_closure_stencil,bids[i])),length(bids))
44 # Boundary quadratures are constructed on the lower-dimensional grid defined
45 # by the coordinite directions orthogonal to that of the boundary.
43 dims = collect(1:dimension(grid)) 46 dims = collect(1:dimension(grid))
44 for id ∈ boundary_identifiers(grid) 47 orth_grids = ntuple(i -> restrict(grid,dims[dims .!= dim(bids[i])]),length(bids))
45 # Boundary operators 48 Hᵧ_pairs = ntuple(i -> Pair(bids[i],DiagonalQuadrature(orth_grids[i],H_closure_stencils)),length(bids))
46 e_pairs = (e_pairs...,Pair(id,BoundaryRestriction(grid,e_closure_stencil,id))) 49
47 d_pairs = (d_pairs...,Pair(id,NormalDerivative(grid,d_closure_stencil,id)))
48 # Boundary quadratures
49 # Construct these on the lower-dimensional grid in the
50 # coordinite directions orthogonal to dim(id)
51 orth_dims = dims[dims .!= dim(id)]
52 orth_grid = restrict(grid,orth_dims)
53 Hᵧ_pairs = (Hᵧ_pairs...,Pair(id,DiagonalQuadrature(orth_grid,H_closure_stencils)))
54 end
55 return Laplace(Δ, H, H⁻¹, Dict(e_pairs), Dict(d_pairs), Dict(Hᵧ_pairs)) 50 return Laplace(Δ, H, H⁻¹, Dict(e_pairs), Dict(d_pairs), Dict(Hᵧ_pairs))
56 end 51 end
57 52
58 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D) 53 LazyTensors.range_size(L::Laplace) = LazyTensors.range_size(L.D)
59 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D) 54 LazyTensors.domain_size(L::Laplace) = LazyTensors.domain_size(L.D)