Mercurial > repos > public > sbplib

diff +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
line wrap: on
line diff
--- a/+scheme/Elastic2dCurvilinear.m	Mon May 07 14:35:54 2018 -0700
+++ b/+scheme/Elastic2dCurvilinear.m	Mon May 07 15:55:35 2018 -0700
@@ -138,10 +138,18 @@
             x = coords(:,1);
             y = coords(:,2);
 
-            x_xi = obj.D1{1}*x;
-            x_eta = obj.D1{2}*x;
-            y_xi = obj.D1{1}*y;
-            y_eta = obj.D1{2}*y;
+            % Use non-periodic difference operators for metric even if opSet is periodic.
+            xmax = max(ops{1}.x);
+            ymax = max(ops{2}.x);
+            opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
+            opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
+            D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
+            D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); 
+
+            x_xi = D1Metric{1}*x;
+            x_eta = D1Metric{2}*x;
+            y_xi = D1Metric{1}*y;
+            y_eta = D1Metric{2}*y;
 
             J = x_xi.*y_eta - x_eta.*y_xi;
 
@@ -156,7 +164,6 @@
             beta{1} = sqrt(x_eta.^2 + y_eta.^2);
             beta{2} = sqrt(x_xi.^2 + y_xi.^2);
 
-
             J = spdiag(J);
             Ji = inv(J);
             for i = 1:dim