--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+sbp/D2Nonequidistant.m	Sun Feb 13 19:32:34 2022 +0100
@@ -0,0 +1,97 @@
+classdef D2Nonequidistant < sbp.OpSet
+    % Implements the boundary optimized variable coefficient
+    % second derivative.
+    %
+    % The boundary closure uses the first and last rows of the
+    % boundary-optimized D1 operator, i.e. the operators are
+    % fully compatible.
+    properties
+        H % Norm matrix
+        HI % H^-1
+        D1 % SBP operator approximating first derivative
+        D2 % SBP operator approximating first derivative
+        DI % Dissipation operator
+        e_l % Left boundary operator
+        e_r % Right boundary operator
+        d1_l % Left boundary first derivative
+        d1_r % Right boundary first derivative
+        m % Number of grid points.
+        h % Step size
+        x % grid
+        borrowing % Struct with borrowing limits for different norm matrices
+        options % Struct holding options used to create the operator
+    end
+
+    methods
+
+        function obj = D2Nonequidistant(m, lim, order, options)
+            % m - number of gridpoints
+            % lim - cell array holding the limits of the domain
+            % order - order of the operator
+            % options - struct holding options used to construct the operator
+            %           struct.stencil_width:   {'minimal', 'nonminimal', 'wide'}
+            %               minimal: minimal compatible stencil width (default)
+            %               nonminimal: a few additional stencil points compared to minimal
+            %               wide: wide stencil obtained by applying D1 twice
+            %           struct.AD: {'op', 'upwind'}
+            %               'op': order-preserving AD (preserving interior stencil order) (default)
+            %               'upwind': upwind AD (order-1 upwind interior stencil)
+            %           struct.variable_coeffs: {true, false}
+            %               true: obj.D2 is a function handle D2(c) returning a matrix 
+            %                     for coefficient vector c (default)
+            %               false: obj.D2 is a matrix.
+            default_arg('options', struct);
+            default_field(options,'stencil_width','minimal');
+            default_field(options,'AD','op');
+            default_field(options,'variable_coeffs',true);
+            [x, h] = sbp.grid.accurateBoundaryOptimizedGrid(lim, m, order);
+            switch order
+                case 4
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_4(m, h, options);
+                case 6
+                    % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_6(m,h,2);
+                    % obj.D2 = @(c)D2;
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_6(m, h, options);
+                case 8
+                    % [obj.H, obj.HI, obj.D1, D2, M, Q, obj.e_l, obj.e_r, obj.d1_l, obj.d1_r] = sbp.implementations.d2_noneq_8(m,h,1);
+                    % obj.D2 = @(c)D2;
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_8(m, h, options);
+                case 10
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_10(m, h, options);
+                case 12
+                    [obj.H, obj.HI, obj.D1, obj.D2, obj.DI] = sbp.implementations.d2_noneq_variable_12(m, h, options);
+                otherwise
+                    error('Invalid operator order %d.', order);
+            end
+
+            if ~options.variable_coeffs
+                obj.D2 = obj.D2(ones(m,1));
+            end
+
+            % Boundary operators
+            obj.e_l = sparse(m, 1); obj.e_l(1) = 1;
+            obj.e_r = sparse(m, 1); obj.e_r(m) = 1;
+            obj.d1_l = (obj.e_l' * obj.D1)';
+            obj.d1_r = (obj.e_r' * obj.D1)';
+
+            % Borrowing coefficients
+            obj.borrowing.H11 = obj.H(1, 1) / h; % First element in H/h,
+            obj.borrowing.M.d1 = obj.H(1, 1) / h; % First element in H/h is borrowing also for M
+            obj.borrowing.R.delta_D = inf;
+
+            % grid data
+            obj.x = x;
+            obj.h = h;
+            obj.m = m;
+
+            % misc
+            obj.options = options;
+        end
+
+        function str = string(obj)
+            str = [class(obj) '_' num2str(obj.order)];
+        end
+
+    end
+
+end