Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinear.m @ 740:f4e2a6a2df08 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make computation of metric derivatives work for periodic.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 07 May 2018 15:55:35 -0700
parents 8efc04e97da4
children 08f3ffe63f48
comparison
equal deleted inserted replaced
739:8efc04e97da4 740:f4e2a6a2df08
136 % -- Metric coefficients ---- 136 % -- Metric coefficients ----
137 coords = g.points(); 137 coords = g.points();
138 x = coords(:,1); 138 x = coords(:,1);
139 y = coords(:,2); 139 y = coords(:,2);
140 140
141 x_xi = obj.D1{1}*x; 141 % Use non-periodic difference operators for metric even if opSet is periodic.
142 x_eta = obj.D1{2}*x; 142 xmax = max(ops{1}.x);
143 y_xi = obj.D1{1}*y; 143 ymax = max(ops{2}.x);
144 y_eta = obj.D1{2}*y; 144 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
145 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
146 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
147 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
148
149 x_xi = D1Metric{1}*x;
150 x_eta = D1Metric{2}*x;
151 y_xi = D1Metric{1}*y;
152 y_eta = D1Metric{2}*y;
145 153
146 J = x_xi.*y_eta - x_eta.*y_xi; 154 J = x_xi.*y_eta - x_eta.*y_xi;
147 155
148 b = cell(dim,dim); 156 b = cell(dim,dim);
149 b{1,1} = y_eta./J; 157 b{1,1} = y_eta./J;
153 161
154 % Scale factors for boundary integrals 162 % Scale factors for boundary integrals
155 beta = cell(dim,1); 163 beta = cell(dim,1);
156 beta{1} = sqrt(x_eta.^2 + y_eta.^2); 164 beta{1} = sqrt(x_eta.^2 + y_eta.^2);
157 beta{2} = sqrt(x_xi.^2 + y_xi.^2); 165 beta{2} = sqrt(x_xi.^2 + y_xi.^2);
158
159 166
160 J = spdiag(J); 167 J = spdiag(J);
161 Ji = inv(J); 168 Ji = inv(J);
162 for i = 1:dim 169 for i = 1:dim
163 beta{i} = spdiag(beta{i}); 170 beta{i} = spdiag(beta{i});