changeset 588:42124009f940 feature/better_multiblock_defs

Add domain definition for circle.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 11 Sep 2017 13:50:29 +0200
parents 25fdc7a625b6
children 53d422efaabd
files +multiblock/+domain/Circle.m
diffstat 1 files changed, 98 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
diff -r 25fdc7a625b6 -r 42124009f940 +multiblock/+domain/Circle.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+multiblock/+domain/Circle.m	Mon Sep 11 13:50:29 2017 +0200
@@ -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