Mercurial > repos > public > sbplib

comparison +scheme/Laplace1D.m @ 896:09c5fbc783d3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename and mordernize scheme.Wave to scheme.Laplace1d. Not fully converted
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 22 Nov 2018 07:58:11 +0100
parents +scheme/Wave.m@cb2b12246b7e
children bd79326ebcd0
comparison
equal deleted inserted replaced
894:f30eafd6d4dc 896:09c5fbc783d3
1 classdef Laplace1D < scheme.Scheme
2 properties
3 g
4 order % Order accuracy for the approximation
5
6 D % non-stabalized scheme operator
7 H % Discrete norm
8 M % Derivative norm
9 a
10
11 D2
12 Hi
13 e_l
14 e_r
15 d_l
16 d_r
17 gamm
18 end
19
20 methods
21 function obj = Laplace1D(g, order, a)
22 default_arg('a', 1);
23
24 assertType(g, 'grid.Cartesian');
25
26 ops = sbp.Ordinary(g.size(), g.h, order);
27
28 obj.D2 = sparse(ops.derivatives.D2);
29 obj.H = sparse(ops.norms.H);
30 obj.Hi = sparse(ops.norms.HI);
31 obj.M = sparse(ops.norms.M);
32 obj.e_l = sparse(ops.boundary.e_1);
33 obj.e_r = sparse(ops.boundary.e_m);
34 obj.d_l = sparse(ops.boundary.S_1);
35 obj.d_r = sparse(ops.boundary.S_m);
36
37
38 obj.g = g;
39 obj.order = order;
40
41 obj.a = a;
42 obj.D = a*obj.D2;
43
44 obj.gamm = h*ops.borrowing.M.S;
45 end
46
47
48 % Closure functions return the opertors applied to the own doamin to close the boundary
49 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
50 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
51 % type is a string specifying the type of boundary condition if there are several.
52 % data is a function returning the data that should be applied at the boundary.
53 % neighbour_scheme is an instance of Scheme that should be interfaced to.
54 % neighbour_boundary is a string specifying which boundary to interface to.
55 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
56 default_arg('type','neumann');
57 default_arg('data',0);
58
59 [e,d,s] = obj.get_boundary_ops(boundary);
60
61 switch type
62 % Dirichlet boundary condition
63 case {'D','dirichlet'}
64 alpha = obj.alpha;
65
66 % tau1 < -alpha^2/gamma
67 tuning = 1.1;
68 tau1 = -tuning*alpha/obj.gamm;
69 tau2 = s*alpha;
70
71 p = tau1*e + tau2*d;
72
73 closure = obj.Hi*p*e';
74
75 pp = obj.Hi*p;
76 switch class(data)
77 case 'double'
78 penalty = pp*data;
79 case 'function_handle'
80 penalty = @(t)pp*data(t);
81 otherwise
82 error('Wierd data argument!')
83 end
84
85
86 % Neumann boundary condition
87 case {'N','neumann'}
88 alpha = obj.alpha;
89 tau1 = -s*alpha;
90 tau2 = 0;
91 tau = tau1*e + tau2*d;
92
93 closure = obj.Hi*tau*d';
94
95 pp = obj.Hi*tau;
96 switch class(data)
97 case 'double'
98 penalty = pp*data;
99 case 'function_handle'
100 penalty = @(t)pp*data(t);
101 otherwise
102 error('Wierd data argument!')
103 end
104
105 % Unknown, boundary condition
106 otherwise
107 error('No such boundary condition: type = %s',type);
108 end
109 end
110
111 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
112 % u denotes the solution in the own domain
113 % v denotes the solution in the neighbour domain
114 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
115 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
116
117 tuning = 1.1;
118
119 alpha_u = obj.alpha;
120 alpha_v = neighbour_scheme.alpha;
121
122 gamm_u = obj.gamm;
123 gamm_v = neighbour_scheme.gamm;
124
125 % tau1 < -(alpha_u/gamm_u + alpha_v/gamm_v)
126
127 tau1 = -(alpha_u/gamm_u + alpha_v/gamm_v) * tuning;
128 tau2 = s_u*1/2*alpha_u;
129 sig1 = s_u*(-1/2);
130 sig2 = 0;
131
132 tau = tau1*e_u + tau2*d_u;
133 sig = sig1*e_u + sig2*d_u;
134
135 closure = obj.Hi*( tau*e_u' + sig*alpha_u*d_u');
136 penalty = obj.Hi*(-tau*e_v' - sig*alpha_v*d_v');
137 end
138
139 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
140 % The right boundary is considered the positive boundary
141 function [e,d,s] = get_boundary_ops(obj,boundary)
142 switch boundary
143 case 'l'
144 e = obj.e_l;
145 d = obj.d_l;
146 s = -1;
147 case 'r'
148 e = obj.e_r;
149 d = obj.d_r;
150 s = 1;
151 otherwise
152 error('No such boundary: boundary = %s',boundary);
153 end
154 end
155
156 function N = size(obj)
157 N = obj.m;
158 end
159
160 end
161
162 methods(Static)
163 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
164 % and bound_v of scheme schm_v.
165 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
166 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
167 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
168 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
169 end
170 end
171 end