Mercurial > repos > public > sbplib

comparison +scheme/Beam.m @ 175:8f22829b69d0 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added and upgraded schemes for the beam equation in 1d and 2d.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 Feb 2016 16:21:47 +0100
parents
children d095b5396103
comparison
equal deleted inserted replaced
174:513019e3d855 175:8f22829b69d0
1 classdef Beam < scheme.Scheme
2 properties
3 order % Order accuracy for the approximation
4 grid
5
6 D % non-stabalized scheme operator
7 alpha
8
9 H % Discrete norm
10 Hi
11
12 e_l, e_r
13 d1_l, d1_r
14 d2_l, d2_r
15 d3_l, d3_r
16 gamm
17 delt
18 end
19
20 methods
21 function obj = Beam(grid, order, alpha, opsGen)
22 default_arg('alpha', 1);
23 default_arg('opsGen', @sbp.Higher);
24
25 if ~isa(grid, 'Cartesian') || grid.D() ~= 1
26 error('Grid must be 1d cartesian');
27 end
28
29 obj.grid = grid;
30 obj.order = order;
31 obj.alpha = alpha;
32
33 m = grid.m;
34 h = grid.spacing();
35
36 ops = opsGen(m, h, order);
37
38 I = speye(m);
39
40 D4 = sparse(ops.derivatives.D4);
41 obj.H = sparse(ops.norms.H);
42 obj.Hi = sparse(ops.norms.HI);
43 obj.e_l = sparse(ops.boundary.e_1);
44 obj.e_r = sparse(ops.boundary.e_m);
45 obj.d1_l = sparse(ops.boundary.S_1);
46 obj.d1_r = sparse(ops.boundary.S_m);
47 obj.d2_l = sparse(ops.boundary.S2_1);
48 obj.d2_r = sparse(ops.boundary.S2_m);
49 obj.d3_l = sparse(ops.boundary.S3_1);
50 obj.d3_r = sparse(ops.boundary.S3_m);
51
52 obj.D = alpha*D4;
53
54 obj.gamm = h*ops.borrowing.N.S2/2;
55 obj.delt = h^3*ops.borrowing.N.S3/2;
56 end
57
58
59 % Closure functions return the opertors applied to the own doamin to close the boundary
60 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
61 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
62 % type is a string specifying the type of boundary condition if there are several.
63 % neighbour_scheme is an instance of Scheme that should be interfaced to.
64 % neighbour_boundary is a string specifying which boundary to interface to.
65 function [closure, penalty_e, penalty_d] = boundary_condition(obj,boundary,type)
66 default_arg('type','dn');
67
68 [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary);
69 gamm = obj.gamm;
70 delt = obj.delt;
71
72 switch type
73 case {'dn'} % Dirichlet-neumann boundary condition
74 alpha = obj.alpha;
75
76 % tau1 < -alpha^2/gamma
77 tuning = 1.1;
78
79 tau1 = tuning * alpha/delt;
80 tau4 = s*alpha;
81
82 sig2 = tuning * alpha/gamm;
83 sig3 = -s*alpha;
84
85 tau = tau1*e+tau4*d3;
86 sig = sig2*d1+sig3*d2;
87
88 closure = obj.Hi*(tau*e' + sig*d1');
89
90 penalty_e = obj.Hi*tau;
91 penalty_d = obj.Hi*sig;
92 otherwise % Unknown, boundary condition
93 error('No such boundary condition: type = %s',type);
94 end
95 end
96
97 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
98 % u denotes the solution in the own domain
99 % v denotes the solution in the neighbour domain
100 [e_u,d1_u,d2_u,d3_u,s_u] = obj.get_boundary_ops(boundary);
101 [e_v,d1_v,d2_v,d3_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
102
103 gamm_u = obj.gamm;
104 delt_u = obj.delt;
105
106 gamm_v = neighbour_scheme.gamm;
107 delt_v = neighbour_scheme.delt;
108
109 tuning = 2;
110
111 alpha_u = obj.alpha;
112 alpha_v = neighbour_scheme.alpha;
113
114 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning;
115 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning;
116 tau4 = s_u*alpha_u/2;
117
118 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning;
119 sig3 = -s_u*alpha_u/2;
120
121 phi2 = s_u*1/2;
122
123 psi1 = -s_u*1/2;
124
125 tau = tau1*e_u + tau4*d3_u;
126 sig = sig2*d1_u + sig3*d2_u ;
127 phi = phi2*d1_u ;
128 psi = psi1*e_u ;
129
130 closure = obj.Hi*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u');
131 penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v');
132 end
133
134 % Returns the boundary ops and sign for the boundary specified by the string boundary.
135 % The right boundary is considered the positive boundary
136 function [e, d1, d2, d3, s] = get_boundary_ops(obj,boundary)
137 switch boundary
138 case 'l'
139 e = obj.e_l;
140 d1 = obj.d1_l;
141 d2 = obj.d2_l;
142 d3 = obj.d3_l;
143 s = -1;
144 case 'r'
145 e = obj.e_e;
146 d1 = obj.d1_e;
147 d2 = obj.d2_e;
148 d3 = obj.d3_e;
149 s = 1;
150 otherwise
151 error('No such boundary: boundary = %s',boundary);
152 end
153 end
154
155 function N = size(obj)
156 N = prod(obj.m);
157 end
158
159 end
160 end