--- a/+multiblock/DefCurvilinear.m	Mon Feb 14 11:14:46 2022 +0100
+++ b/+multiblock/DefCurvilinear.m	Mon Feb 14 14:55:29 2022 +0100
@@ -37,13 +37,39 @@
             obj.boundaryGroups = boundaryGroups;
         end
 
-        function g = getGrid(obj, varargin)
-            ms = obj.getGridSizes(varargin{:});
+        % Returns a multiblock.Grid given some parameters
+        % ms: cell array of [mx, my] vectors
+        % Currently defaults to an equidistant curvilinear grid if varargin is empty.
+        % If varargin is non-empty, the first argument should supply the grid type, followed by
+        % additional arguments required to construct the grid.
+        % Grid types:
+        %          'equidist' - equidistant curvilinear grid
+        %                       Additional argumets: none
+        %          'boundaryopt' - boundary optimized grid based on boundary
+        %                          optimized SBP operators
+        %                          Additional arguments: order, stencil option
+        function g = getGrid(obj, ms, varargin)
+            % If a scalar is passed, defer to getGridSizes implemented by subclass
+            % TODO: This forces the interface of subclasses.
+            % Should ms be included in varargin? Figure out bow to do it properly
+            if length(ms) == 1
+                ms = obj.getGridSizes(ms);
+            end
 
-            grids = cell(1, obj.nBlocks);
-            for i = 1:obj.nBlocks
-                grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i});
-            end
+            if isempty(varargin) || strcmp(varargin{1},'equidist')
+                gridgenerator = @(blockMap,m) grid.equidistantCurvilinear(blockMap, m);
+            elseif strcmp(varargin{1},'boundaryopt')
+                 order = varargin{2};
+                 stenciloption = varargin{3};
+                 gridgenerator = @(blockMap,m) grid.boundaryOptimizedCurvilinear(blockMap,m,{0,1},{0,1},...
+                     order,stenciloption);
+             else
+                 error('No grid type supplied!');
+             end
+             grids = cell(1, obj.nBlocks);
+             for i = 1:obj.nBlocks
+                grids{i} = gridgenerator(obj.blockMaps{i}.S, ms{i});
+             end
 
             g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
         end