Mercurial > repos > public > sbplib_julia
annotate diffOp.jl @ 76:81d9510cb2d0
Make Laplace take dimension as a parameter
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Thu, 24 Jan 2019 09:24:21 +0100 |
parents | d04569696918 |
children | 700a74c41b26 48079bd39969 |
rev | line source |
---|---|
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 abstract type DiffOp end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
2 |
33
e8d7137b3f07
Raname apply(DiffOp) to apply!(DiffOp)
Jonatan Werpers <jonatan@werpers.com>
parents:
31
diff
changeset
|
3 function apply!(D::DiffOp, u::AbstractVector, v::AbstractVector) |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
4 error("not implemented") |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
5 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
6 |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
7 function innerProduct(D::DiffOp, u::AbstractVector, v::AbstractVector)::Real |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
8 error("not implemented") |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
9 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
10 |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
11 function matrixRepresentation(D::DiffOp) |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
12 error("not implemented") |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
13 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
14 |
55
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
15 function boundaryCondition(D::DiffOp,b::Grid.BoundaryId,type)::(Closure, Penalty) |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
16 error("not implemented") |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
17 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
18 |
54
4300a3fbd818
switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents:
49
diff
changeset
|
19 function interface(Du::DiffOp, Dv::DiffOp, b::Grid.BoundaryId; type) |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
20 error("not implemented") |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
21 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
22 |
55
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
23 abstract type Closure end |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
24 |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
25 function apply(c::Closure, v::AbstractVector, i::Int) |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
26 error("not implemented") |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
27 end |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
28 |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
29 abstract type Penalty end |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
30 |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
31 function apply(c::Penalty, g, i::Int) |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
32 error("not implemented") |
c62ea0112d4d
Add abstract types for Closure and Penalty
Ylva Rydin <ylva.rydin@telia.com>
parents:
54
diff
changeset
|
33 end |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
34 |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
35 # Differential operator for a*d^2/dx^2 |
76
81d9510cb2d0
Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents:
73
diff
changeset
|
36 struct Laplace{D, T<:Real} <: DiffOp |
81d9510cb2d0
Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents:
73
diff
changeset
|
37 grid::Grid.EquidistantGrid{D,T} |
81d9510cb2d0
Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents:
73
diff
changeset
|
38 a::T |
73
d04569696918
Add types to Laplace operators
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
39 op::D2{Float64} |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
40 end |
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
41 |
6
cb8e50ca9e15
Add attempt att apply methods for Laplace
Jonatan Werpers <jonatan@werpers.com>
parents:
2
diff
changeset
|
42 # u = L*v |
76
81d9510cb2d0
Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents:
73
diff
changeset
|
43 function apply!(L::Laplace{1}, u::AbstractVector, v::AbstractVector) |
54
4300a3fbd818
switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents:
49
diff
changeset
|
44 h = Grid.spacings(L.grid)[1] |
43
ef060ab3b035
remove stride and remove some bugs
Ylva Rydin <ylva.rydin@telia.com>
parents:
34
diff
changeset
|
45 apply!(L.op, u, v, h) |
48
be53b58581a6
Fix bug with call by reference in Laplace1D apply
Jonatan Werpers <jonatan@werpers.com>
parents:
43
diff
changeset
|
46 u .= L.a * u |
9
aafc4bdbe675
Add h in stencil application
Jonatan Werpers <jonatan@werpers.com>
parents:
8
diff
changeset
|
47 return nothing |
2
43be32298ae2
Add function to get closure size
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
48 end |
49 | 49 |
50 # u = L*v | |
76
81d9510cb2d0
Make Laplace take dimension as a parameter
Jonatan Werpers <jonatan@werpers.com>
parents:
73
diff
changeset
|
51 function apply!(L::Laplace{2}, u::AbstractVector, v::AbstractVector) |
49 | 52 u .= 0*u |
54
4300a3fbd818
switch grid to Grid in diffOp
Ylva Rydin <ylva.rydin@telia.com>
parents:
49
diff
changeset
|
53 h = Grid.spacings(L.grid) |
49 | 54 |
55 li = LinearIndices(L.grid.numberOfPointsPerDim) | |
56 n_x, n_y = L.grid.numberOfPointsPerDim | |
57 | |
58 | |
59 # For each x | |
60 temp = zeros(eltype(u), n_y) | |
61 for i ∈ 1:n_x | |
62 | |
63 v_i = view(v, li[i,:]) | |
64 apply!(L.op, temp, v_i, h[2]) | |
65 | |
66 u[li[i,:]] += temp | |
67 end | |
68 | |
69 # For each y | |
70 temp = zeros(eltype(u), n_x) | |
71 for i ∈ 1:n_y | |
72 v_i = view(v, li[:,i]) | |
73 apply!(L.op, temp, v_i, h[1]) | |
74 | |
75 u[li[:,i]] += temp | |
76 end | |
77 | |
78 u .= L.a*u | |
79 | |
80 return nothing | |
81 end |