--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/Cartesian.m	Thu Nov 15 16:36:21 2018 -0800
@@ -0,0 +1,197 @@
+classdef Cartesian < grid.Structured
+    properties
+        n % Number of points in the grid
+        d % Number of dimensions
+        m % Number of points in each direction
+        x % Cell array of vectors with node placement for each dimension.
+        h % Spacing/Scaling
+        lim % Cell array of left and right boundaries for each dimension.
+    end
+
+    % General d dimensional grid with n points
+    methods
+        % Creates a cartesian grid given vectors conatining the coordinates
+        % in each direction
+        function obj = Cartesian(varargin)
+            obj.d = length(varargin);
+
+            for i = 1:obj.d
+                assert(isnumeric(varargin{i}), 'Coordinate inputs must be vectors.')
+
+                obj.x{i} = varargin{i};
+                obj.m(i) = length(varargin{i});
+            end
+
+            obj.n = prod(obj.m);
+            if obj.n == 0
+                error('grid:Cartesian:EmptyGrid','Input parameter gives an empty grid.')
+            end
+
+            obj.h = [];
+
+            obj.lim = cell(1,obj.d);
+            for i = 1:obj.d
+                obj.lim{i} = {obj.x{i}(1), obj.x{i}(end)};
+            end
+        end
+        % n returns the number of points in the grid
+        function o = N(obj)
+            o = obj.n;
+        end
+
+        % d returns the spatial dimension of the grid
+        function o = D(obj)
+            o = obj.d;
+        end
+
+        function m = size(obj)
+            m = obj.m;
+        end
+
+        % points returns a n x d matrix containing the coordianets for all points.
+        % points are ordered according to the kronecker product with X*Y*Z
+        function X = points(obj)
+            X = zeros(obj.n, obj.d);
+
+            for i = 1:obj.d
+                if iscolumn(obj.x{i})
+                    c = obj.x{i};
+                else
+                    c = obj.x{i}';
+                end
+
+                m_before = prod(obj.m(1:i-1));
+                m_after = prod(obj.m(i+1:end));
+
+                X(:,i) = kr(ones(m_before,1),c,ones(m_after,1));
+            end
+        end
+
+        % matrices returns a cell array with coordinates in matrix form.
+        % For 2d case these will have to be transposed to work with plotting routines.
+        function X = matrices(obj)
+
+            if obj.d == 1 % There is no 1d matrix data type in matlab, handle special case
+                X{1} = reshape(obj.x{1}, [obj.m 1]);
+                return
+            end
+
+            X = cell(1,obj.d);
+            for i = 1:obj.d
+                s = ones(1,obj.d);
+                s(i) = obj.m(i);
+
+                t = reshape(obj.x{i},s);
+
+                s = obj.m;
+                s(i) = 1;
+                X{i} = repmat(t,s);
+            end
+        end
+
+        function h = scaling(obj)
+            if isempty(obj.h)
+                error('grid:Cartesian:NoScalingSet', 'No scaling set')
+            end
+
+            h = obj.h;
+        end
+
+        % Restricts the grid function gf on obj to the subgrid g.
+        % Only works for even multiples
+        function gf = restrictFunc(obj, gf, g)
+            m1 = obj.m;
+            m2 = g.m;
+
+            % Check the input
+            if prod(m1) ~= numel(gf)
+                error('grid:Cartesian:restrictFunc:NonMatchingFunctionSize', 'The grid function has to few or too many points.');
+            end
+
+            if ~all(mod(m1-1,m2-1) == 0)
+                error('grid:Cartesian:restrictFunc:NonMatchingGrids', 'Only integer downsamplings are allowed');
+            end
+
+            % Calculate stride for each dimension
+            stride = (m1-1)./(m2-1);
+
+            % Create downsampling indecies
+            I = {};
+            for i = 1:length(m1)
+                I{i} = 1:stride(i):m1(i);
+            end
+
+            gf = reshapeRowMaj(gf, m1);
+            gf = gf(I{:});
+            gf = reshapeRowMaj(gf, prod(m2));
+        end
+
+        % Projects the grid function gf on obj to the grid g.
+        function gf = projectFunc(obj, gf, g)
+            error('grid:Cartesian:NotImplemented')
+        end
+
+        % Return the names of all boundaries in this grid.
+        function bs = getBoundaryNames(obj)
+            switch obj.D()
+                case 1
+                    bs = {'l', 'r'};
+                case 2
+                    bs = {'w', 'e', 's', 'n'};
+                case 3
+                    bs = {'w', 'e', 's', 'n', 'd', 'u'};
+                otherwise
+                    error('not implemented');
+            end
+        end
+
+        % Return coordinates for the given boundary
+        function X = getBoundary(obj, name)
+            % In what dimension is the boundary?
+            switch name
+                case {'l', 'r', 'w', 'e'}
+                    D = 1;
+                case {'s', 'n'}
+                    D = 2;
+                case {'d', 'u'}
+                    D = 3;
+                otherwise
+                    error('not implemented');
+            end
+
+            % At what index is the boundary?
+            switch name
+                case {'l', 'w', 's', 'd'}
+                    index = 1;
+                case {'r', 'e', 'n', 'u'}
+                    index = obj.m(D);
+                otherwise
+                    error('not implemented');
+            end
+
+
+
+            I = cell(1, obj.d);
+            for i = 1:obj.d
+                if i == D
+                    I{i} = index;
+                else
+                    I{i} = ':';
+                end
+            end
+
+            % Calculate size of result:
+            m = obj.m;
+            m(D) = [];
+            N = prod(m);
+
+            X = zeros(N, obj.d);
+
+            coordMat = obj.matrices();
+            for i = 1:length(coordMat)
+                Xtemp = coordMat{i}(I{:});
+                X(:,i) = reshapeRowMaj(Xtemp, [N,1]);
+            end
+        end
+    end
+end
\ No newline at end of file