diff -r 6d5953fc090e -r 0aff87f6fb2c +multiblock/+domain/Circle.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Circle.m	Fri Apr 20 12:06:38 2018 -0700
@@ -0,0 +1,98 @@
+classdef Circle < multiblock.DefCurvilinear
+    properties
+        r, c
+
+        hs
+        r_arc
+        omega
+    end
+
+    methods
+        function obj = Circle(r, c, hs)
+            default_arg('r', 1);
+            default_arg('c', [0; 0]);
+            default_arg('hs', 0.435);
+
+
+            % alpha = 0.75;
+            % hs = alpha*r/sqrt(2);
+
+            % Square should not be a square, it should be an arc. The arc radius
+            % is chosen so that the three angles of the meshes are all equal.
+            % This gives that the (half)arc opening angle of should be omega = pi/12
+            omega = pi/12;
+            r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); %  = hs* 1/sin(omega)
+            c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0];
+
+            cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]);
+
+            c2 = cir(0);
+            c3 = cir(1);
+
+            s1 = [-hs; -hs];
+            s2 = [ hs; -hs];
+            s3 = [ hs;  hs];
+            s4 = [-hs;  hs];
+
+            sp2 = parametrization.Curve.line(s2,c2);
+            sp3 = parametrization.Curve.line(s3,c3);
+
+            Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]);
+            Se2 = Se1.rotate(c,pi/2);
+            Se3 = Se2.rotate(c,pi/2);
+            Se4 = Se3.rotate(c,pi/2);
+
+
+            S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1);
+
+            A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse);
+            B = A.rotate(c,1*pi/2).rotate_edges(-1);
+            C = A.rotate(c,2*pi/2).rotate_edges(-1);
+            D = A.rotate(c,3*pi/2).rotate_edges(0);
+
+            blocks = {S,A,B,C,D};
+            blocksNames = {'S','A','B','C','D'};
+
+            conn = cell(5,5);
+            conn{1,2} = {'e','w'};
+            conn{1,3} = {'n','s'};
+            conn{1,4} = {'w','s'};
+            conn{1,5} = {'s','w'};
+
+            conn{2,3} = {'n','e'};
+            conn{3,4} = {'w','e'};
+            conn{4,5} = {'w','s'};
+            conn{5,2} = {'n','s'};
+
+            boundaryGroups = struct();
+            boundaryGroups.E = multiblock.BoundaryGroup({2,'e'});
+            boundaryGroups.N = multiblock.BoundaryGroup({3,'n'});
+            boundaryGroups.W = multiblock.BoundaryGroup({4,'n'});
+            boundaryGroups.S = multiblock.BoundaryGroup({5,'e'});
+            boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}});
+
+            obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames);
+
+            obj.r     = r;
+            obj.c     = c;
+            obj.hs    = hs;
+            obj.r_arc = r_arc;
+            obj.omega = omega;
+        end
+
+        function ms = getGridSizes(obj, m)
+            m_S = m;
+
+            % m_Radial
+            s = 2*obj.hs;
+            innerArc = obj.r_arc*obj.omega;
+            outerArc = obj.r*pi/2;
+            shortSpoke = obj.r-s/sqrt(2);
+            x = (1/(2-sqrt(3))-1)*obj.hs;
+            longSpoke =  (obj.r+x)-obj.r_arc;
+            m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2);
+
+            ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]};
+        end
+    end
+end

diff -r 6d5953fc090e -r 0aff87f6fb2c +multiblock/+domain/Rectangle.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Rectangle.m	Fri Apr 20 12:06:38 2018 -0700
@@ -0,0 +1,188 @@
+classdef Rectangle < multiblock.Definition
+    properties
+
+    blockTi % Transfinite interpolation objects used for plotting
+    xlims
+    ylims
+    blockNames % Cell array of block labels
+    nBlocks
+    connections % Cell array specifying connections between blocks
+    boundaryGroups % Structure of boundaryGroups
+
+    end
+
+
+    methods
+        % Creates a divided rectangle
+        % x and y are vectors of boundary and interface positions.
+        % blockNames: cell array of labels. The id is default.
+        function obj = Rectangle(x,y,blockNames)
+            default_arg('blockNames',[]);
+
+            n = length(y)-1; % number of blocks in the y direction.
+            m = length(x)-1; % number of blocks in the x direction.
+            N = n*m; % number of blocks
+
+            if ~issorted(x)
+                error('The elements of x seem to be in the wrong order');
+            end
+            if ~issorted(flip(y))
+                error('The elements of y seem to be in the wrong order');
+            end
+
+            % Dimensions of blocks and number of points
+            blockTi = cell(N,1);
+            xlims = cell(N,1);
+            ylims = cell(N,1);
+            for i = 1:n
+                for j = 1:m
+                    p1 = [x(j), y(i+1)];
+                    p2 = [x(j+1), y(i)];
+                    I = flat_index(m,j,i);
+                    blockTi{I} = parametrization.Ti.rectangle(p1,p2);
+                    xlims{I} = {x(j), x(j+1)};
+                    ylims{I} = {y(i+1), y(i)};
+                end
+            end
+
+            % Interface couplings
+            conn = cell(N,N);
+            for i = 1:n
+                for j = 1:m
+                    I = flat_index(m,j,i);
+                    if i < n
+                        J = flat_index(m,j,i+1);
+                        conn{I,J} = {'s','n'};
+                    end
+
+                    if j < m
+                        J = flat_index(m,j+1,i);
+                        conn{I,J} = {'e','w'};
+                    end
+                end
+            end
+
+            % Block names (id number as default)
+            if isempty(blockNames)
+                obj.blockNames = cell(1, N);
+                for i = 1:N
+                    obj.blockNames{i} = sprintf('%d', i);
+                end
+            else
+                assert(length(blockNames) == N);
+                obj.blockNames = blockNames;
+            end
+            nBlocks = N;
+
+            % Boundary groups
+            boundaryGroups = struct();
+            nx = m;
+            ny = n;
+            E = cell(1,ny);
+            W = cell(1,ny);
+            S = cell(1,nx);
+            N = cell(1,nx);
+            for i = 1:ny
+                E_id = flat_index(m,nx,i);
+                W_id = flat_index(m,1,i);
+                E{i} = {E_id,'e'};
+                W{i} = {W_id,'w'};
+            end
+            for j = 1:nx
+                S_id = flat_index(m,j,ny);
+                N_id = flat_index(m,j,1);
+                S{j} = {S_id,'s'};
+                N{j} = {N_id,'n'};
+            end  
+            boundaryGroups.E = multiblock.BoundaryGroup(E);
+            boundaryGroups.W = multiblock.BoundaryGroup(W);
+            boundaryGroups.S = multiblock.BoundaryGroup(S);
+            boundaryGroups.N = multiblock.BoundaryGroup(N);
+            boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]);
+            boundaryGroups.WS = multiblock.BoundaryGroup([W,S]);
+            boundaryGroups.WN = multiblock.BoundaryGroup([W,N]);
+            boundaryGroups.ES = multiblock.BoundaryGroup([E,S]);
+            boundaryGroups.EN = multiblock.BoundaryGroup([E,N]);
+
+            obj.connections = conn;
+            obj.nBlocks = nBlocks;
+            obj.boundaryGroups = boundaryGroups;
+            obj.blockTi = blockTi;
+            obj.xlims = xlims;
+            obj.ylims = ylims;
+
+        end
+
+
+        % Returns a multiblock.Grid given some parameters
+        % ms: cell array of [mx, my] vectors
+        % For same [mx, my] in every block, just input one vector.
+        function g = getGrid(obj, ms, varargin)
+
+            default_arg('ms',[21,21])
+
+            % Extend ms if input is a single vector
+            if (numel(ms) == 2) && ~iscell(ms)
+                m = ms;
+                ms = cell(1,obj.nBlocks);
+                for i = 1:obj.nBlocks
+                    ms{i} = m;
+                end
+            end
+
+            grids = cell(1, obj.nBlocks);
+            for i = 1:obj.nBlocks
+                grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i});
+            end
+
+            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
+        end
+
+        % label is the type of label used for plotting,
+        % default is block name, 'id' show the index for each block.
+        function show(obj, label, gridLines, varargin)
+            default_arg('label', 'name')
+            default_arg('gridLines', false);
+
+            if isempty('label') && ~gridLines
+                for i = 1:obj.nBlocks
+                    obj.blockTi{i}.show(2,2);
+                end
+                axis equal
+                return
+            end
+
+            if gridLines
+                m = 10;
+                for i = 1:obj.nBlocks
+                    obj.blockTi{i}.show(m,m);
+                end
+            end
+
+
+            switch label
+                case 'name'
+                    labels = obj.blockNames;
+                case 'id'
+                    labels = {};
+                    for i = 1:obj.nBlocks
+                        labels{i} = num2str(i);
+                    end
+                otherwise
+                    axis equal
+                    return
+            end
+
+            for i = 1:obj.nBlocks
+                parametrization.Ti.label(obj.blockTi{i}, labels{i});
+            end
+
+            axis equal
+        end
+
+        % Returns the grid size of each block in a cell array
+        % The input parameters are determined by the subclass
+        function ms = getGridSizes(obj, varargin)
+        end
+    end
+end

diff -r 6d5953fc090e -r 0aff87f6fb2c +multiblock/Def.m
--- a/+multiblock/Def.m	Fri Apr 20 11:32:04 2018 -0700
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,101 +0,0 @@
-classdef Def
-    properties
-        nBlocks
-        blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation
-        blockNames
-        connections % Cell array specifying connections between blocks
-        boundaryGroups % Structure of boundaryGroups
-    end
-
-    methods
-        % Defines a multiblock setup for transfinite interpolation blocks
-        % TODO: How to bring in plotting of points?
-        function obj = Def(blockMaps, connections, boundaryGroups, blockNames)
-            default_arg('boundaryGroups', struct());
-            default_arg('blockNames',{});
-
-            nBlocks = length(blockMaps);
-
-            obj.nBlocks = nBlocks;
-
-            obj.blockMaps = blockMaps;
-
-            assert(all(size(connections) == [nBlocks, nBlocks]));
-            obj.connections = connections;
-
-
-            if isempty(blockNames)
-                obj.blockNames = cell(1, nBlocks);
-                for i = 1:length(blockMaps)
-                    obj.blockNames{i} = sprintf('%d', i);
-                end
-            else
-                assert(length(blockNames) == nBlocks);
-                obj.blockNames = blockNames;
-            end
-
-            obj.boundaryGroups = boundaryGroups;
-        end
-
-        function g = getGrid(obj, varargin)
-            ms = obj.getGridSizes(varargin{:});
-
-            grids = cell(1, obj.nBlocks);
-            for i = 1:obj.nBlocks
-                grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i});
-            end
-
-            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
-        end
-
-        function show(obj, label, gridLines, varargin)
-            default_arg('label', 'name')
-            default_arg('gridLines', false);
-
-            if isempty('label') && ~gridLines
-                for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(2,2);
-                end
-                axis equal
-                return
-            end
-
-            if gridLines
-                ms = obj.getGridSizes(varargin{:});
-                for i = 1:obj.nBlocks
-                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
-                end
-            end
-
-
-            switch label
-                case 'name'
-                    labels = obj.blockNames;
-                case 'id'
-                    labels = {};
-                    for i = 1:obj.nBlocks
-                        labels{i} = num2str(i);
-                    end
-                otherwise
-                    axis equal
-                    return
-            end
-
-            for i = 1:obj.nBlocks
-                parametrization.Ti.label(obj.blockMaps{i}, labels{i});
-            end
-
-            axis equal
-        end
-    end
-
-    methods (Abstract)
-        % Returns the grid size of each block in a cell array
-        % The input parameters are determined by the subclass
-        ms = getGridSizes(obj, varargin)
-        % end
-    end
-
-end
-
-

diff -r 6d5953fc090e -r 0aff87f6fb2c +multiblock/DefCurvilinear.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/DefCurvilinear.m	Fri Apr 20 12:06:38 2018 -0700
@@ -0,0 +1,101 @@
+classdef DefCurvilinear < multiblock.Definition
+    properties
+        nBlocks
+        blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation
+        blockNames
+        connections % Cell array specifying connections between blocks
+        boundaryGroups % Structure of boundaryGroups
+    end
+
+    methods
+        % Defines a multiblock setup for transfinite interpolation blocks
+        % TODO: How to bring in plotting of points?
+        function obj = DefCurvilinear(blockMaps, connections, boundaryGroups, blockNames)
+            default_arg('boundaryGroups', struct());
+            default_arg('blockNames',{});
+
+            nBlocks = length(blockMaps);
+
+            obj.nBlocks = nBlocks;
+
+            obj.blockMaps = blockMaps;
+
+            assert(all(size(connections) == [nBlocks, nBlocks]));
+            obj.connections = connections;
+
+
+            if isempty(blockNames)
+                obj.blockNames = cell(1, nBlocks);
+                for i = 1:length(blockMaps)
+                    obj.blockNames{i} = sprintf('%d', i);
+                end
+            else
+                assert(length(blockNames) == nBlocks);
+                obj.blockNames = blockNames;
+            end
+
+            obj.boundaryGroups = boundaryGroups;
+        end
+
+        function g = getGrid(obj, varargin)
+            ms = obj.getGridSizes(varargin{:});
+
+            grids = cell(1, obj.nBlocks);
+            for i = 1:obj.nBlocks
+                grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i});
+            end
+
+            g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
+        end
+
+        function show(obj, label, gridLines, varargin)
+            default_arg('label', 'name')
+            default_arg('gridLines', false);
+
+            if isempty('label') && ~gridLines
+                for i = 1:obj.nBlocks
+                    obj.blockMaps{i}.show(2,2);
+                end
+                axis equal
+                return
+            end
+
+            if gridLines
+                ms = obj.getGridSizes(varargin{:});
+                for i = 1:obj.nBlocks
+                    obj.blockMaps{i}.show(ms{i}(1),ms{i}(2));
+                end
+            end
+
+
+            switch label
+                case 'name'
+                    labels = obj.blockNames;
+                case 'id'
+                    labels = {};
+                    for i = 1:obj.nBlocks
+                        labels{i} = num2str(i);
+                    end
+                otherwise
+                    axis equal
+                    return
+            end
+
+            for i = 1:obj.nBlocks
+                parametrization.Ti.label(obj.blockMaps{i}, labels{i});
+            end
+
+            axis equal
+        end
+    end
+
+    methods (Abstract)
+        % Returns the grid size of each block in a cell array
+        % The input parameters are determined by the subclass
+        ms = getGridSizes(obj, varargin)
+        % end
+    end
+
+end
+
+

diff -r 6d5953fc090e -r 0aff87f6fb2c +multiblock/Definition.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/Definition.m	Fri Apr 20 12:06:38 2018 -0700
@@ -0,0 +1,11 @@
+classdef Definition
+    methods (Abstract)
+
+        % Returns a multiblock.Grid given some parameters
+        g = getGrid(obj, varargin)
+
+        % label is the type of label used for plotting,
+        % default is block name, 'id' show the index for each block.
+        show(obj, label, gridLines, varargin)
+    end
+end