Mercurial > repos > public > sbplib

changeset 1191:881afc40a3d2 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change to use the interpolation operators in sbplib
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 29 Jul 2019 16:43:08 +0200
parents 79618b58b0a0
children b3c47a716d57
files +rv/+diffops/constructDiffOps.m
diffstat 1 files changed, 21 insertions(+), 19 deletions(-) [+]
line wrap: on
line diff
--- a/+rv/+diffops/constructDiffOps.m	Mon Jul 08 15:12:19 2019 +0200
+++ b/+rv/+diffops/constructDiffOps.m	Mon Jul 29 16:43:08 2019 +0200
@@ -8,7 +8,7 @@
             diffOpsStruct = diffOpsMultiStage(varargin{:});
         case {'mg1','multi-grid1'}
             diffOpsStruct = diffOpsMultiGrid1(varargin{:});
-        case {'mg2','multi-grid2'}
+        case {'mg','mg2','multi-grid2'} % Default multigrid diffops
             diffOpsStruct = diffOpsMultiGrid2(varargin{:});
     end
 end
@@ -85,30 +85,30 @@
                      'D_t', D_t,...
                      'penalties_scheme', {penalties_scheme});
 end
-%% Multigrid functions:
-% TODO: Implement properly.
-function v_c = coarse_v(g, g_c, v)
-    V = grid.funcToMatrix(g,v);
-    v_c = reshape(V(1:2:end,1:2:end)',g_c.N(),1);
-end
 
 % Note: Only works for equidistant grids.
-% TODO: Change from using matlabs interp to using proper interpolation operators..
 function D_c = coarserGridDiffOp(scheme, g, schemeOrder, schemeParams, opSet, BCs)
     m = g.m();
     lim = g.lim();
     m_c = (m-1)/2 + 1;
     switch g.D()
-    case 1
-        g_c = grid.equidistant(m_c, lim{1});
-        D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
-        x = g.x{1};
-        x_c = x(1:2:end);
-        D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
-    case 2
-        g_c = grid.equidistant(m_c, lim{1}, lim{2});
-        D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
-        D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1);
+        case 1
+            interpOps = sbp.InterpOpsOP(m_c, m, schemeOrder, schemeOrder);
+            I = interpOps.Iu2v.good;
+            g_c = grid.equidistant(m_c, lim{1});
+            D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
+            D_c = @(v) I*D_c(v(1:2:end));
+        case 2
+            interpOps_x = sbp.InterpOpsOP(m_c(1), m(1), schemeOrder, schemeOrder);
+            I_x = interpOps_x.Iu2v.good;
+            interpOps_y = sbp.InterpOpsOP(m_c(2), m(2), schemeOrder, schemeOrder);
+            I_y = interpOps_y.Iu2v.good;
+            I = kron(I_x,I_y);
+            g_c = grid.equidistant(m_c, lim{1}, lim{2});
+            ind = grid.funcToMatrix(g, 1:g.N());
+            ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
+            D_c = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
+            D_c = @(v) I*D_c(v(ind_c));
     end
 end
 
@@ -126,9 +126,11 @@
         D_c = @(v) interp1(x_c,D_c(v(1:2:end)),x,'spline');
     case 2
         g_c = grid.equidistant(m_c, lim{1}, lim{2});
+        ind = grid.funcToMatrix(g, 1:g.N());
+        ind_c = reshape(ind(1:2:end,1:2:end)',g_c.N(),1);
         [D_c, closures] = rv.diffops.constructFluxDiffOps(scheme, g_c, schemeOrder, schemeParams, opSet, BCs);
         D_c = rv.diffops.addClosuresToDiffOp(D_c, closures);
-        D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(coarse_v(g,g_c,v))),'spline')',g.N(),1);
+        D_c = @(v) reshape(interp2(grid.funcToMatrix(g_c,D_c(v(ind_c))),'spline')',g.N(),1);
     end
 end