Mercurial > repos > public > sbplib

comparison +scheme/Beam2d.m @ 175:8f22829b69d0 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added and upgraded schemes for the beam equation in 1d and 2d.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 Feb 2016 16:21:47 +0100
parents cb2b12246b7e
children d095b5396103
comparison
equal deleted inserted replaced
174:513019e3d855 175:8f22829b69d0
1 classdef Beam2d < scheme.Scheme 1 classdef Beam2d < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 grid
4 N % Number of points total
5 h % Grid spacing
6 u,v % Grid
7 x,y % Values of x and y for each grid point
8 order % Order accuracy for the approximation 4 order % Order accuracy for the approximation
9 5
10 D % non-stabalized scheme operator 6 D % non-stabalized scheme operator
11 M % Derivative norm 7 M % Derivative norm
12 alpha 8 alpha
25 delt_x, delt_y 21 delt_x, delt_y
26 end 22 end
27 23
28 methods 24 methods
29 function obj = Beam2d(m,lim,order,alpha,opsGen) 25 function obj = Beam2d(m,lim,order,alpha,opsGen)
26 default_arg('alpha',1);
30 default_arg('opsGen',@sbp.Higher); 27 default_arg('opsGen',@sbp.Higher);
31 default_arg('a',1); 28
32 29 if ~isa(grid, 'Cartesian') || grid.D() ~= 2
33 if length(m) == 1 30 error('Grid must be 2d cartesian');
34 m = [m m];
35 end 31 end
36 32
37 m_x = m(1); 33 obj.grid = grid;
38 m_y = m(2); 34 obj.alpha = alpha;
39 35 obj.order = order;
40 xlim = lim{1}; 36
41 ylim = lim{2}; 37 m_x = grid.m(1);
42 38 m_y = grid.m(2);
43 [x, h_x] = util.get_grid(xlim{:},m_x); 39
44 [y, h_y] = util.get_grid(ylim{:},m_y); 40 h = grid.scaling();
41 h_x = h(1);
42 h_y = h(2);
45 43
46 ops_x = opsGen(m_x,h_x,order); 44 ops_x = opsGen(m_x,h_x,order);
47 ops_y = opsGen(m_y,h_y,order); 45 ops_y = opsGen(m_y,h_y,order);
48 46
49 I_x = speye(m_x); 47 I_x = speye(m_x);
50 I_y = speye(m_y); 48 I_y = speye(m_y);
51
52
53
54 49
55 D4_x = sparse(ops_x.derivatives.D4); 50 D4_x = sparse(ops_x.derivatives.D4);
56 H_x = sparse(ops_x.norms.H); 51 H_x = sparse(ops_x.norms.H);
57 Hi_x = sparse(ops_x.norms.HI); 52 Hi_x = sparse(ops_x.norms.HI);
58 e_l_x = sparse(ops_x.boundary.e_1); 53 e_l_x = sparse(ops_x.boundary.e_1);
103 obj.d3_w = kr(d3_l_x,I_y); 98 obj.d3_w = kr(d3_l_x,I_y);
104 obj.d3_e = kr(d3_r_x,I_y); 99 obj.d3_e = kr(d3_r_x,I_y);
105 obj.d3_s = kr(I_x,d3_l_y); 100 obj.d3_s = kr(I_x,d3_l_y);
106 obj.d3_n = kr(I_x,d3_r_y); 101 obj.d3_n = kr(I_x,d3_r_y);
107 102
108 obj.m = m;
109 obj.h = [h_x h_y];
110 obj.order = order;
111
112 obj.alpha = alpha;
113 obj.D = alpha*D4; 103 obj.D = alpha*D4;
114 obj.u = x;
115 obj.v = y;
116 obj.x = kr(x,ones(m_y,1));
117 obj.y = kr(ones(m_x,1),y);
118 104
119 obj.gamm_x = h_x*ops_x.borrowing.N.S2/2; 105 obj.gamm_x = h_x*ops_x.borrowing.N.S2/2;
120 obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2; 106 obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2;
121 107
122 obj.gamm_y = h_y*ops_y.borrowing.N.S2/2; 108 obj.gamm_y = h_y*ops_y.borrowing.N.S2/2;