Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1338:da61892884a4 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for using periodic grids.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:29:05 +0200
parents df8c71b80c33
children
comparison
equal deleted inserted replaced
1337:bf2554f1825d 1338:da61892884a4
150 150
151 % Metric coefficients 151 % Metric coefficients
152 coords = g.points(); 152 coords = g.points();
153 x = coords(:,1); 153 x = coords(:,1);
154 y = coords(:,2); 154 y = coords(:,2);
155 155
156 x_u = Du*x; 156 % When computing the metric derivatives, we can't use periodic operators
157 x_v = Dv*x; 157 if isequal(opSet{1},@sbp.D2VariablePeriodic)
158 y_u = Du*y; 158 ops_metric = sbp.D2Standard(m_u, {0, 1-1/m_u}, order);
159 y_v = Dv*y; 159 Du_metric = kr(ops_metric.D1,I_v);
160 else
161 Du_metric = Du;
162 end
163
164 if isequal(opSet{2},@sbp.D2VariablePeriodic)
165 ops_metric = sbp.D2Standard(m_v, {0, 1-1/m_v}, order);
166 Dv_metric = kr(I_u,ops_metric.D1);
167 else
168 Dv_metric = Dv;
169 end
170
171 x_u = Du_metric*x;
172 x_v = Dv_metric*x;
173 y_u = Du_metric*y;
174 y_v = Dv_metric*y;
160 175
161 J = x_u.*y_v - x_v.*y_u; 176 J = x_u.*y_v - x_v.*y_u;
162 177
163 K = cell(dim, dim); 178 K = cell(dim, dim);
164 K{1,1} = y_v./J; 179 K{1,1} = y_v./J;
166 K{2,1} = -x_v./J; 181 K{2,1} = -x_v./J;
167 K{2,2} = x_u./J; 182 K{2,2} = x_u./J;
168 obj.K = K; 183 obj.K = K;
169 184
170 % Physical derivatives 185 % Physical derivatives
186 % TODO: Fix for periodic operators
171 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; 187 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
172 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; 188 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
173 189
174 % Wrap around Aniosotropic Cartesian 190 % Wrap around Aniosotropic Cartesian
175 rho_tilde = J.*rho; 191 rho_tilde = J.*rho;