Mercurial > repos > public > sbplib
annotate +scheme/Beam2d.m @ 958:72cd29107a9a feature/poroelastic
Temporary changes in multiblock.DiffOp. Change traction operators in Elastic2dvariable to be true boundary operators. But adjoint FD conv test fails for dirichlet BC so need to debug!
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 05 Dec 2018 18:58:10 -0800 |
| parents | d095b5396103 |
| children | 459eeb99130f |
| rev | line source |
|---|---|
|
141
cb2b12246b7e
Fixed path to some superclasses.
Jonatan Werpers <jonatan@werpers.com>
parents:
125
diff
changeset
|
1 classdef Beam2d < scheme.Scheme |
| 0 | 2 properties |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
3 grid |
| 0 | 4 order % Order accuracy for the approximation |
| 5 | |
| 6 D % non-stabalized scheme operator | |
| 7 M % Derivative norm | |
| 8 alpha | |
| 9 | |
| 10 H % Discrete norm | |
| 11 Hi | |
| 12 H_x, H_y % Norms in the x and y directions | |
| 13 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. | |
| 14 Hi_x, Hi_y | |
| 15 Hix, Hiy | |
| 16 e_w, e_e, e_s, e_n | |
| 17 d1_w, d1_e, d1_s, d1_n | |
| 18 d2_w, d2_e, d2_s, d2_n | |
| 19 d3_w, d3_e, d3_s, d3_n | |
| 20 gamm_x, gamm_y | |
| 21 delt_x, delt_y | |
| 22 end | |
| 23 | |
| 24 methods | |
|
125
d52e5cdb6eff
Fixed some class name, file name errors.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
25 function obj = Beam2d(m,lim,order,alpha,opsGen) |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
26 default_arg('alpha',1); |
| 0 | 27 default_arg('opsGen',@sbp.Higher); |
| 28 | |
|
176
d095b5396103
Fixed some bugs in Beam schemes.
Jonatan Werpers <jonatan@werpers.com>
parents:
175
diff
changeset
|
29 if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 2 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
30 error('Grid must be 2d cartesian'); |
| 0 | 31 end |
| 32 | |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
33 obj.grid = grid; |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
34 obj.alpha = alpha; |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
35 obj.order = order; |
| 0 | 36 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
37 m_x = grid.m(1); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
38 m_y = grid.m(2); |
| 0 | 39 |
|
175
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
40 h = grid.scaling(); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
41 h_x = h(1); |
|
8f22829b69d0
Added and upgraded schemes for the beam equation in 1d and 2d.
Jonatan Werpers <jonatan@werpers.com>
parents:
141
diff
changeset
|
42 h_y = h(2); |
| 0 | 43 |
| 44 ops_x = opsGen(m_x,h_x,order); | |
| 45 ops_y = opsGen(m_y,h_y,order); | |
| 46 | |
| 47 I_x = speye(m_x); | |
| 48 I_y = speye(m_y); | |
| 49 | |
| 50 D4_x = sparse(ops_x.derivatives.D4); | |
| 51 H_x = sparse(ops_x.norms.H); | |
| 52 Hi_x = sparse(ops_x.norms.HI); | |
| 53 e_l_x = sparse(ops_x.boundary.e_1); | |
| 54 e_r_x = sparse(ops_x.boundary.e_m); | |
| 55 d1_l_x = sparse(ops_x.boundary.S_1); | |
| 56 d1_r_x = sparse(ops_x.boundary.S_m); | |
| 57 d2_l_x = sparse(ops_x.boundary.S2_1); | |
| 58 d2_r_x = sparse(ops_x.boundary.S2_m); | |
| 59 d3_l_x = sparse(ops_x.boundary.S3_1); | |
| 60 d3_r_x = sparse(ops_x.boundary.S3_m); | |
| 61 | |
| 62 D4_y = sparse(ops_y.derivatives.D4); | |
| 63 H_y = sparse(ops_y.norms.H); | |
| 64 Hi_y = sparse(ops_y.norms.HI); | |
| 65 e_l_y = sparse(ops_y.boundary.e_1); | |
| 66 e_r_y = sparse(ops_y.boundary.e_m); | |
| 67 d1_l_y = sparse(ops_y.boundary.S_1); | |
| 68 d1_r_y = sparse(ops_y.boundary.S_m); | |
| 69 d2_l_y = sparse(ops_y.boundary.S2_1); | |
| 70 d2_r_y = sparse(ops_y.boundary.S2_m); | |
| 71 d3_l_y = sparse(ops_y.boundary.S3_1); | |
| 72 d3_r_y = sparse(ops_y.boundary.S3_m); | |
| 73 | |
| 74 | |
| 75 D4 = kr(D4_x, I_y) + kr(I_x, D4_y); | |
| 76 | |
| 77 % Norms | |
| 78 obj.H = kr(H_x,H_y); | |
| 79 obj.Hx = kr(H_x,I_x); | |
| 80 obj.Hy = kr(I_x,H_y); | |
| 81 obj.Hix = kr(Hi_x,I_y); | |
| 82 obj.Hiy = kr(I_x,Hi_y); | |
| 83 obj.Hi = kr(Hi_x,Hi_y); | |
| 84 | |
| 85 % Boundary operators | |
| 86 obj.e_w = kr(e_l_x,I_y); | |
| 87 obj.e_e = kr(e_r_x,I_y); | |
| 88 obj.e_s = kr(I_x,e_l_y); | |
| 89 obj.e_n = kr(I_x,e_r_y); | |
| 90 obj.d1_w = kr(d1_l_x,I_y); | |
| 91 obj.d1_e = kr(d1_r_x,I_y); | |
| 92 obj.d1_s = kr(I_x,d1_l_y); | |
| 93 obj.d1_n = kr(I_x,d1_r_y); | |
| 94 obj.d2_w = kr(d2_l_x,I_y); | |
| 95 obj.d2_e = kr(d2_r_x,I_y); | |
| 96 obj.d2_s = kr(I_x,d2_l_y); | |
| 97 obj.d2_n = kr(I_x,d2_r_y); | |
| 98 obj.d3_w = kr(d3_l_x,I_y); | |
| 99 obj.d3_e = kr(d3_r_x,I_y); | |
| 100 obj.d3_s = kr(I_x,d3_l_y); | |
| 101 obj.d3_n = kr(I_x,d3_r_y); | |
| 102 | |
| 103 obj.D = alpha*D4; | |
| 104 | |
| 105 obj.gamm_x = h_x*ops_x.borrowing.N.S2/2; | |
| 106 obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2; | |
| 107 | |
| 108 obj.gamm_y = h_y*ops_y.borrowing.N.S2/2; | |
| 109 obj.delt_y = h_y^3*ops_y.borrowing.N.S3/2; | |
| 110 end | |
| 111 | |
| 112 | |
| 113 % Closure functions return the opertors applied to the own doamin to close the boundary | |
| 114 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
| 115 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
| 116 % type is a string specifying the type of boundary condition if there are several. | |
| 117 % data is a function returning the data that should be applied at the boundary. | |
| 118 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
| 119 % neighbour_boundary is a string specifying which boundary to interface to. | |
| 120 function [closure, penalty_e,penalty_d] = boundary_condition(obj,boundary,type,data) | |
| 121 default_arg('type','dn'); | |
| 122 default_arg('data',0); | |
| 123 | |
| 124 [e,d1,d2,d3,s,gamm,delt,halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 125 | |
| 126 switch type | |
| 127 % Dirichlet-neumann boundary condition | |
| 128 case {'dn'} | |
| 129 alpha = obj.alpha; | |
| 130 | |
| 131 % tau1 < -alpha^2/gamma | |
| 132 tuning = 1.1; | |
| 133 | |
| 134 tau1 = tuning * alpha/delt; | |
| 135 tau4 = s*alpha; | |
| 136 | |
| 137 sig2 = tuning * alpha/gamm; | |
| 138 sig3 = -s*alpha; | |
| 139 | |
| 140 tau = tau1*e+tau4*d3; | |
| 141 sig = sig2*d1+sig3*d2; | |
| 142 | |
| 143 closure = halfnorm_inv*(tau*e' + sig*d1'); | |
| 144 | |
| 145 pp_e = halfnorm_inv*tau; | |
| 146 pp_d = halfnorm_inv*sig; | |
| 147 switch class(data) | |
| 148 case 'double' | |
| 149 penalty_e = pp_e*data; | |
| 150 penalty_d = pp_d*data; | |
| 151 case 'function_handle' | |
| 152 penalty_e = @(t)pp_e*data(t); | |
| 153 penalty_d = @(t)pp_d*data(t); | |
| 154 otherwise | |
| 155 error('Wierd data argument!') | |
| 156 end | |
| 157 | |
| 158 % Unknown, boundary condition | |
| 159 otherwise | |
| 160 error('No such boundary condition: type = %s',type); | |
| 161 end | |
| 162 end | |
| 163 | |
| 164 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
| 165 % u denotes the solution in the own domain | |
| 166 % v denotes the solution in the neighbour domain | |
| 167 [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); | |
| 168 [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); | |
| 169 | |
| 170 tuning = 2; | |
| 171 | |
| 172 alpha_u = obj.alpha; | |
| 173 alpha_v = neighbour_scheme.alpha; | |
| 174 | |
| 175 tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; | |
| 176 % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; | |
| 177 tau4 = s_u*alpha_u/2; | |
| 178 | |
| 179 sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning; | |
| 180 sig3 = -s_u*alpha_u/2; | |
| 181 | |
| 182 phi2 = s_u*1/2; | |
| 183 | |
| 184 psi1 = -s_u*1/2; | |
| 185 | |
| 186 tau = tau1*e_u + tau4*d3_u; | |
| 187 sig = sig2*d1_u + sig3*d2_u ; | |
| 188 phi = phi2*d1_u ; | |
| 189 psi = psi1*e_u ; | |
| 190 | |
| 191 closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); | |
| 192 penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); | |
| 193 end | |
| 194 | |
| 195 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
| 196 % The right boundary is considered the positive boundary | |
| 197 function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) | |
| 198 switch boundary | |
| 199 case 'w' | |
| 200 e = obj.e_w; | |
| 201 d1 = obj.d1_w; | |
| 202 d2 = obj.d2_w; | |
| 203 d3 = obj.d3_w; | |
| 204 s = -1; | |
| 205 gamm = obj.gamm_x; | |
| 206 delt = obj.delt_x; | |
| 207 halfnorm_inv = obj.Hix; | |
| 208 case 'e' | |
| 209 e = obj.e_e; | |
| 210 d1 = obj.d1_e; | |
| 211 d2 = obj.d2_e; | |
| 212 d3 = obj.d3_e; | |
| 213 s = 1; | |
| 214 gamm = obj.gamm_x; | |
| 215 delt = obj.delt_x; | |
| 216 halfnorm_inv = obj.Hix; | |
| 217 case 's' | |
| 218 e = obj.e_s; | |
| 219 d1 = obj.d1_s; | |
| 220 d2 = obj.d2_s; | |
| 221 d3 = obj.d3_s; | |
| 222 s = -1; | |
| 223 gamm = obj.gamm_y; | |
| 224 delt = obj.delt_y; | |
| 225 halfnorm_inv = obj.Hiy; | |
| 226 case 'n' | |
| 227 e = obj.e_n; | |
| 228 d1 = obj.d1_n; | |
| 229 d2 = obj.d2_n; | |
| 230 d3 = obj.d3_n; | |
| 231 s = 1; | |
| 232 gamm = obj.gamm_y; | |
| 233 delt = obj.delt_y; | |
| 234 halfnorm_inv = obj.Hiy; | |
| 235 otherwise | |
| 236 error('No such boundary: boundary = %s',boundary); | |
| 237 end | |
| 238 end | |
| 239 | |
| 240 function N = size(obj) | |
| 241 N = prod(obj.m); | |
| 242 end | |
| 243 | |
| 244 end | |
| 245 end |