Mercurial > repos > public > sbplib_julia

comparison DiffOps/src/laplace.jl @ 248:05e7bbe0af97 boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge refactoring of laplace
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 27 Jun 2019 13:36:04 +0200
parents ed29ee13e92e d9e262cb2e8d
children 863a98d9e798
comparison
equal deleted inserted replaced
247:ed29ee13e92e 248:05e7bbe0af97
36 36
37 37
38 """ 38 """
39 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 39 BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
40 40
41 Implements the boundary operator `e` as a TensorMapping 41 Implements the boundary operator `e` as a TensorMapping
42 """ 42 """
43 struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1} 43 struct BoundaryValue{T,N,M,K} <: TensorMapping{T,2,1}
44 op::D2{T,N,M,K} 44 op::D2{T,N,M,K}
45 grid::EquidistantGrid 45 grid::EquidistantGrid
46 bId::CartesianBoundary 46 bId::CartesianBoundary
47 end 47 end
48 export BoundaryValue 48 export BoundaryValue
49 49
50 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? 50 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
51 # Can we give special treatment to TensorMappings that go to a higher dim? 51 # Can we give special treatment to TensorMappings that go to a higher dim?
52 LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid) 52 LazyTensors.range_size(e::BoundaryValue{T}, domain_size::NTuple{1,Integer}) where T = size(e.grid)
53 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],) 53 LazyTensors.domain_size(e::BoundaryValue{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
54 54
55 function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int}) 55 function LazyTensors.apply(e::BoundaryValue, v::AbstractArray, I::NTuple{2,Int})
56 i = I[dim(e.bId)] 56 i = I[dim(e.bId)]
57 j = I[3-dim(e.bId)] 57 j = I[3-dim(e.bId)]
58 N_i = size(e.grid)[dim(e.bId)] 58 N_i = size(e.grid)[dim(e.bId)]
59 59 return apply_e(e.op, v[j], N_i, i, region(e.bId))
60 if region(e.bId) == Lower
61 # NOTE: closures are indexed by offset. Fix this D:<
62 return e.op.eClosure[i-1]*v[j]
63 elseif region(e.bId) == Upper
64 return e.op.eClosure[N_i-i]*v[j]
65 end
66 end 60 end
67 61
68 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int}) 62 function LazyTensors.apply_transpose(e::BoundaryValue, v::AbstractArray, I::NTuple{1,Int})
69 u = selectdim(v,3-dim(e.bId),I[1]) 63 u = selectdim(v,3-dim(e.bId),I[1])
70 return apply_e(e.op, u, region(e.bId)) 64 return apply_e_T(e.op, u, region(e.bId))
71 end 65 end
72 66
73 67
74 68
75 """ 69 """
76 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 70 NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
77 71
78 Implements the boundary operator `d` as a TensorMapping 72 Implements the boundary operator `d` as a TensorMapping
79 """ 73 """
80 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1} 74 struct NormalDerivative{T,N,M,K} <: TensorMapping{T,2,1}
81 op::D2{T,N,M,K} 75 op::D2{T,N,M,K}
82 grid::EquidistantGrid 76 grid::EquidistantGrid
83 bId::CartesianBoundary 77 bId::CartesianBoundary
84 end 78 end
85 export NormalDerivative 79 export NormalDerivative
86 80
87 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue? 81 # TODO: This is obviouly strange. Is domain_size just discarded? Is there a way to avoid storing grid in BoundaryValue?
88 # Can we give special treatment to TensorMappings that go to a higher dim? 82 # Can we give special treatment to TensorMappings that go to a higher dim?
90 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],) 84 LazyTensors.domain_size(e::NormalDerivative{T}, range_size::NTuple{2,Integer}) where T = (range_size[3-dim(e.bId)],)
91 85
92 # Not correct abstraction level 86 # Not correct abstraction level
93 # TODO: Not type stable D:< 87 # TODO: Not type stable D:<
94 function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int}) 88 function LazyTensors.apply(d::NormalDerivative, v::AbstractArray, I::NTuple{2,Int})
95 i = I[dim(d.bId)] 89 i = I[dim(d.bId)]
96 j = I[3-dim(d.bId)] 90 j = I[3-dim(d.bId)]
97 N_i = size(d.grid)[dim(d.bId)] 91 N_i = size(d.grid)[dim(d.bId)]
98 92 h_inv = d.grid.inverse_spacing[dim(d.bId)]
99 if region(d.bId) == Lower 93 return apply_d(d.op, h_inv, v[j], N_i, i, region(d.bId))
100 # Note, closures are indexed by offset. Fix this D:<
101 return d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[i-1]*v[j]
102 elseif region(d.bId) == Upper
103 return -d.grid.inverse_spacing[dim(d.bId)]*d.op.dClosure[N_i-i]*v[j]
104 end
105 end 94 end
106 95
107 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int}) 96 function LazyTensors.apply_transpose(d::NormalDerivative, v::AbstractArray, I::NTuple{1,Int})
108 u = selectdim(v,3-dim(d.bId),I[1]) 97 u = selectdim(v,3-dim(d.bId),I[1])
109 return apply_d(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId)) 98 return apply_d_T(d.op, d.grid.inverse_spacing[dim(d.bId)], u, region(d.bId))
110 end 99 end
111 100
112 101
113 102
114 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end 103 struct Neumann{Bid<:BoundaryIdentifier} <: BoundaryCondition end
136 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g) 125 return -L.Hi*(tau/h*e + d)*Hᵧ*(e'*v - g)
137 # Need to handle scalar multiplication and addition of TensorMapping 126 # Need to handle scalar multiplication and addition of TensorMapping
138 end 127 end
139 128
140 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D 129 # function apply(s::MyWaveEq{D}, v::AbstractArray{T,D}, i::CartesianIndex{D}) where D
141 # return apply(s.L, v, i) + 130 # return apply(s.L, v, i) +
142 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) + 131 # sat(s.L, Dirichlet{CartesianBoundary{1,Lower}}(s.tau), v, s.g_w, i) +
143 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) + 132 # sat(s.L, Dirichlet{CartesianBoundary{1,Upper}}(s.tau), v, s.g_e, i) +
144 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) + 133 # sat(s.L, Dirichlet{CartesianBoundary{2,Lower}}(s.tau), v, s.g_s, i) +
145 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i) 134 # sat(s.L, Dirichlet{CartesianBoundary{2,Upper}}(s.tau), v, s.g_n, i)
146 # end 135 # end