--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/boundaryOptimizedCurvilinear.m	Mon Feb 14 14:55:29 2022 +0100
@@ -0,0 +1,57 @@
+% Creates a curvilinear grid of dimension length(m).
+% over the logical domain xi_lim, eta_lim, ...
+% The grid point distribution is based on the boundary
+% optimized SBP operators, passed as arguments, order
+% and opt.
+% Examples:
+%   g = grid.boundaryOptimizedCurvilinear(mapping, [mx, my], xlim, ylim, order, opt)
+%   g = grid.boundaryOptimizedCurvilinear(mapping, [10, 15], {0,1}, {0,2}, 4) - defaults to 'accurate' stencils
+%   g = grid.boundaryOptimizedCurvilinear(mapping, [10, 15], {0,1}, {0,2}, 4, 'minimal')
+function g = boundaryOptimizedCurvilinear(mapping, m, varargin)
+    n = length(m);
+    % TODO: The input parameter check is also present in boundaryOptimizedGrid.m,
+    % Consider refactoring it.
+        
+    % Check that parameters matches dimensions
+    matchingParams = false;
+    if length(varargin) == n+1 % Minimal number of arguments
+            matchingParams = iscell([varargin{1:n}]) && ...
+                             isfloat([varargin{n+1}]);
+    elseif length(varargin) == n+2 % Stencil options supplied
+            matchingParams = iscell([varargin{1:n}]) && ...
+                             isfloat([varargin{n+1}]) && ...
+                             ischar([varargin{n+2}]);
+    end
+    assert(matchingParams,'grid:boundaryOptimizedCurvilinear:NonMatchingParameters','The number of parameters per dimensions do not match.');
+
+    % Check that stencil options are passed correctly (if supplied)
+    if length(varargin) == n+2 % Stencil options supplied
+        availabe_opts = ["Accurate","accurate","A","acc","Minimal","minimal","M","min"];
+        assert(any(varargin{n+2} == availabe_opts), ...
+            'grid:boundaryOptimizedCurvilinear:InvalidOption',"The operator option must be 'accurate' or 'minimal.'");
+    else %If not passed, populate varargin with default option 'accurate'
+        varargin(n+2) = {'accurate'};
+    end
+
+    % Specify generating function
+    switch varargin{n+2}
+        case {'Accurate','accurate','A','acc'}
+            gridgenerator = @sbp.grid.accurateBoundaryOptimizedGrid;
+        case {'Minimal','minimal','M','min'}
+            gridgenerator = @sbp.grid.minimalBoundaryOptimizedGrid;
+    end
+
+    X = {};
+    h = [];
+    for i = 1:n
+        try
+            [X{i},h(i)] = gridgenerator(varargin{i},m(i),varargin{n+1});
+        catch exception % Propagate any errors in the grid generation functions.
+            msgText = getReport(exception);
+            error('grid:boundaryOptimizedCurvilinear:InvalidParameter',msgText)
+        end
+    end
+
+    g = grid.Curvilinear(mapping, X{:});
+    g.logic.h = h;
+end
\ No newline at end of file