Mercurial > repos > public > sbplib

comparison +multiblock/+domain/Circle.m @ 588:42124009f940 feature/better_multiblock_defs

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add domain definition for circle.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 11 Sep 2017 13:50:29 +0200
parents
children 9be370486d36
comparison
equal deleted inserted replaced
587:25fdc7a625b6 588:42124009f940
1 classdef Circle < multiblock.DefCurvilinear
2 properties
3 r, c
4
5 hs
6 r_arc
7 omega
8 end
9
10 methods
11 function obj = Circle(r, c, hs)
12 default_arg('r', 1);
13 default_arg('c', [0; 0]);
14 default_arg('hs', 0.435);
15
16
17 % alpha = 0.75;
18 % hs = alpha*r/sqrt(2);
19
20 % Square should not be a square, it should be an arc. The arc radius
21 % is chosen so that the three angles of the meshes are all equal.
22 % This gives that the (half)arc opening angle of should be omega = pi/12
23 omega = pi/12;
24 r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); % = hs* 1/sin(omega)
25 c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0];
26
27 cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]);
28
29 c2 = cir(0);
30 c3 = cir(1);
31
32 s1 = [-hs; -hs];
33 s2 = [ hs; -hs];
34 s3 = [ hs; hs];
35 s4 = [-hs; hs];
36
37 sp2 = parametrization.Curve.line(s2,c2);
38 sp3 = parametrization.Curve.line(s3,c3);
39
40 Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]);
41 Se2 = Se1.rotate(c,pi/2);
42 Se3 = Se2.rotate(c,pi/2);
43 Se4 = Se3.rotate(c,pi/2);
44
45
46 S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1);
47
48 A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse);
49 B = A.rotate(c,1*pi/2).rotate_edges(-1);
50 C = A.rotate(c,2*pi/2).rotate_edges(-1);
51 D = A.rotate(c,3*pi/2).rotate_edges(0);
52
53 blocks = {S,A,B,C,D};
54 blocksNames = {'S','A','B','C','D'};
55
56 conn = cell(5,5);
57 conn{1,2} = {'e','w'};
58 conn{1,3} = {'n','s'};
59 conn{1,4} = {'w','s'};
60 conn{1,5} = {'s','w'};
61
62 conn{2,3} = {'n','e'};
63 conn{3,4} = {'w','e'};
64 conn{4,5} = {'w','s'};
65 conn{5,2} = {'n','s'};
66
67 boundaryGroups = struct();
68 boundaryGroups.E = multiblock.BoundaryGroup({2,'e'});
69 boundaryGroups.N = multiblock.BoundaryGroup({3,'n'});
70 boundaryGroups.W = multiblock.BoundaryGroup({4,'n'});
71 boundaryGroups.S = multiblock.BoundaryGroup({5,'e'});
72 boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}});
73
74 obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames);
75
76 obj.r = r;
77 obj.c = c;
78 obj.hs = hs;
79 obj.r_arc = r_arc;
80 obj.omega = omega;
81 end
82
83 function ms = getGridSizes(obj, m)
84 m_S = m;
85
86 % m_Radial
87 s = 2*obj.hs;
88 innerArc = obj.r_arc*obj.omega;
89 outerArc = obj.r*pi/2;
90 shortSpoke = obj.r-s/sqrt(2);
91 x = (1/(2-sqrt(3))-1)*obj.hs;
92 longSpoke = (obj.r+x)-obj.r_arc;
93 m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2);
94
95 ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]};
96 end
97 end
98 end