Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/DiffOps.jl @ 231:fbabfd4e8f20

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in boundary_conditions
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 26 Jun 2019 15:07:47 +0200
parents cd60382f392b
children a5fdc00d5070
comparison
equal deleted inserted replaced
144:ce56727e4232 231:fbabfd4e8f20
1 module DiffOps
2
3 using RegionIndices
4 using SbpOperators
5 using Grids
6
7 """
8 DiffOp
9
10 Supertype of differential operator discretisations.
11 The action of the DiffOp is defined in the method
12 apply(D::DiffOp, v::AbstractVector, I...)
13 """
14 abstract type DiffOp end
15
16 function apply end
17
18 function matrixRepresentation(D::DiffOp)
19 error("not implemented")
20 end
21
22 abstract type DiffOpCartesian{Dim} <: DiffOp end
23
24 # DiffOp must have a grid of dimension Dim!!!
25 function apply!(D::DiffOpCartesian{Dim}, u::AbstractArray{T,Dim}, v::AbstractArray{T,Dim}) where {T,Dim}
26 for I ∈ eachindex(D.grid)
27 u[I] = apply(D, v, I)
28 end
29
30 return nothing
31 end
32 export apply!
33
34 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
35 apply_region!(D, u, v, Lower, Lower)
36 apply_region!(D, u, v, Lower, Interior)
37 apply_region!(D, u, v, Lower, Upper)
38 apply_region!(D, u, v, Interior, Lower)
39 apply_region!(D, u, v, Interior, Interior)
40 apply_region!(D, u, v, Interior, Upper)
41 apply_region!(D, u, v, Upper, Lower)
42 apply_region!(D, u, v, Upper, Interior)
43 apply_region!(D, u, v, Upper, Upper)
44 return nothing
45 end
46
47 # Maybe this should be split according to b3fbef345810 after all?! Seems like it makes performance more predictable
48 function apply_region!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
49 for I ∈ regionindices(D.grid.size, closureSize(D.op), (r1,r2))
50 @inbounds indextuple = (Index{r1}(I[1]), Index{r2}(I[2]))
51 @inbounds u[I] = apply(D, v, indextuple)
52 end
53 return nothing
54 end
55 export apply_region!
56
57 function apply_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}) where T
58 apply_region_tiled!(D, u, v, Lower, Lower)
59 apply_region_tiled!(D, u, v, Lower, Interior)
60 apply_region_tiled!(D, u, v, Lower, Upper)
61 apply_region_tiled!(D, u, v, Interior, Lower)
62 apply_region_tiled!(D, u, v, Interior, Interior)
63 apply_region_tiled!(D, u, v, Interior, Upper)
64 apply_region_tiled!(D, u, v, Upper, Lower)
65 apply_region_tiled!(D, u, v, Upper, Interior)
66 apply_region_tiled!(D, u, v, Upper, Upper)
67 return nothing
68 end
69
70 using TiledIteration
71 function apply_region_tiled!(D::DiffOpCartesian{2}, u::AbstractArray{T,2}, v::AbstractArray{T,2}, r1::Type{<:Region}, r2::Type{<:Region}) where T
72 ri = regionindices(D.grid.size, closureSize(D.op), (r1,r2))
73 # TODO: Pass Tilesize to function
74 for tileaxs ∈ TileIterator(axes(ri), padded_tilesize(T, (5,5), 2))
75 for j ∈ tileaxs[2], i ∈ tileaxs[1]
76 I = ri[i,j]
77 u[I] = apply(D, v, (Index{r1}(I[1]), Index{r2}(I[2])))
78 end
79 end
80 return nothing
81 end
82 export apply_region_tiled!
83
84 function apply(D::DiffOp, v::AbstractVector)::AbstractVector
85 u = zeros(eltype(v), size(v))
86 apply!(D,v,u)
87 return u
88 end
89
90 export apply
91
92 """
93 A BoundaryCondition should implement the method
94 sat(::DiffOp, v::AbstractArray, data::AbstractArray, ...)
95 """
96 abstract type BoundaryCondition end
97
98
99 include("laplace.jl")
100 export Laplace
101
102
103 end # module