Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1270:9a7cac4202d0 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
In CruvilinearAnisotropicUpwind: Make D1Upwind default opSet. Use D2Standard for metric derivatives. Bugfix global quadrature.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 31 May 2020 20:18:45 -0700
parents 02dfe3a56742
children 15865fbda16e
comparison
equal deleted inserted replaced
1269:5fa2f62a58a1 1270:9a7cac4202d0
62 62
63 % The coefficients can either be function handles or grid functions 63 % The coefficients can either be function handles or grid functions
64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization. 64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
65 function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag) 65 function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
66 default_arg('rho', @(x,y) 0*x+1); 66 default_arg('rho', @(x,y) 0*x+1);
67 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible}); 67 default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
68 default_arg('optFlag', false); 68 default_arg('optFlag', false);
69 dim = 2; 69 dim = 2;
70 70
71 C_default = cell(dim,dim,dim,dim); 71 C_default = cell(dim,dim,dim,dim);
72 for i = 1:dim 72 for i = 1:dim
103 103
104 m = g.size(); 104 m = g.size();
105 m_tot = g.N(); 105 m_tot = g.N();
106 106
107 % 1D operators 107 % 1D operators
108 opSetMetric = {@sbp.D2VariableCompatible, @sbp.D2VariableCompatible};
108 m_u = m(1); 109 m_u = m(1);
109 m_v = m(2); 110 m_v = m(2);
110 ops_u = opSet{1}(m_u, {0, 1}, order); 111 ops_u = opSetMetric{1}(m_u, {0, 1}, order);
111 ops_v = opSet{2}(m_v, {0, 1}, order); 112 ops_v = opSetMetric{2}(m_v, {0, 1}, order);
112 113
113 h_u = ops_u.h; 114 h_u = ops_u.h;
114 h_v = ops_v.h; 115 h_v = ops_v.h;
115 116
116 I_u = speye(m_u); 117 I_u = speye(m_u);
182 end 183 end
183 end 184 end
184 end 185 end
185 186
186 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1}); 187 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
187 refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, []); 188 refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, opSet);
188 189
189 %---- Set object properties ------ 190 %---- Set object properties ------
190 obj.RHO = spdiag(rho); 191 obj.RHO = spdiag(rho);
191 192
192 % Volume quadrature 193 % Volume quadrature
193 obj.J = spdiag(J); 194 obj.J = spdiag(J);
194 obj.Ji = spdiag(1./J); 195 obj.Ji = spdiag(1./J);
195 obj.H = obj.J*kr(H_u,H_v); 196 obj.H = obj.J*refObj.H;
196 197
197 % Boundary quadratures 198 % Boundary quadratures
198 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); 199 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
199 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); 200 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
200 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); 201 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);