Mercurial > repos > public > sbplib

annotate +scheme/LaplaceCurvilinear.m @ 774:66eb4a2bbb72 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove default scaling of the system. The scaling doens't seem to help actual solutions. One example that fails in the flexural code. With large timesteps the solutions seems to blow up. One particular example is profilePresentation on the tdb_presentation_figures branch with k = 0.0005
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Jul 2018 15:42:52 -0700
parents 33b962620e24
children 07f8311374c6
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
450
8d455e49364f Copy Wave2dCurve to new scheme LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 449
diff changeset
1 classdef LaplaceCurvilinear < scheme.Scheme
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 m % Number of points in each direction, possibly a vector
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 h % Grid spacing
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
5
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
6 grid
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
7
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 order % Order accuracy for the approximation
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
552
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
10 a,b % Parameters of the operator
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
11
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
12
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
13 % Inner products and operators for physical coordinates
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
14 D % Laplace operator
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
15 H, Hi % Inner product
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
16 e_w, e_e, e_s, e_n
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
17 d_w, d_e, d_s, d_n % Normal derivatives at the boundary
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
18 H_w, H_e, H_s, H_n % Boundary inner products
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
19 Dx, Dy % Physical derivatives
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
20 M % Gradient inner product
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
21
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
22 % Metric coefficients
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 J, Ji
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 a11, a12, a22
552
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
25 x_u
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
26 x_v
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
27 y_u
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
28 y_v
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
552
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
30 % Inner product and operators for logical coordinates
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 H_u, H_v % Norms in the x and y directions
552
ffaf13533c27 Reorganize properties of the class
Jonatan Werpers <jonatan@werpers.com>
parents: 551
diff changeset
32 Hi_u, Hi_v
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 Hiu, Hiv
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 du_w, dv_w
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 du_e, dv_e
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 du_s, dv_s
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 du_n, dv_n
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 gamm_u, gamm_v
27
97a638f91fb8 Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents: 5
diff changeset
40 lambda
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 methods
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
44 % Implements a*div(b*grad(u)) as a SBP scheme
539
08b6281ba2a9 Add some todos
Jonatan Werpers <jonatan@werpers.com>
parents: 452
diff changeset
45 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
08b6281ba2a9 Add some todos
Jonatan Werpers <jonatan@werpers.com>
parents: 452
diff changeset
46
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
47 function obj = LaplaceCurvilinear(g ,order, a, b, opSet)
303
f18142c1530b Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents: 96
diff changeset
48 default_arg('opSet',@sbp.D2Variable);
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
49 default_arg('a', 1);
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
50 default_arg('b', 1);
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
51
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
52 if b ~=1
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
53 error('Not implemented yet')
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
54 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
56 assert(isa(g, 'grid.Curvilinear'))
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
58 m = g.size();
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 m_u = m(1);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 m_v = m(2);
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
61 m_tot = g.N();
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
63 h = g.scaling();
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
64 h_u = h(1);
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
65 h_v = h(2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
67
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
68 % 1D operators
303
f18142c1530b Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents: 96
diff changeset
69 ops_u = opSet(m_u, {0, 1}, order);
f18142c1530b Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents: 96
diff changeset
70 ops_v = opSet(m_v, {0, 1}, order);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 I_u = speye(m_u);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 I_v = speye(m_v);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
75 D1_u = ops_u.D1;
303
f18142c1530b Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents: 96
diff changeset
76 D2_u = ops_u.D2;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
77 H_u = ops_u.H;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
78 Hi_u = ops_u.HI;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
79 e_l_u = ops_u.e_l;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
80 e_r_u = ops_u.e_r;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
81 d1_l_u = ops_u.d1_l;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
82 d1_r_u = ops_u.d1_r;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
84 D1_v = ops_v.D1;
303
f18142c1530b Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents: 96
diff changeset
85 D2_v = ops_v.D2;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
86 H_v = ops_v.H;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
87 Hi_v = ops_v.HI;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
88 e_l_v = ops_v.e_l;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
89 e_r_v = ops_v.e_r;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
90 d1_l_v = ops_v.d1_l;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
91 d1_r_v = ops_v.d1_r;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
92
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
93
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
94 % Logical operators
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
95 Du = kr(D1_u,I_v);
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
96 Dv = kr(I_u,D1_v);
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
97 obj.Hu = kr(H_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
98 obj.Hv = kr(I_u,H_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
99 obj.Hiu = kr(Hi_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
100 obj.Hiv = kr(I_u,Hi_v);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
102 e_w = kr(e_l_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
103 e_e = kr(e_r_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
104 e_s = kr(I_u,e_l_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
105 e_n = kr(I_u,e_r_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
106 obj.du_w = kr(d1_l_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
107 obj.dv_w = (e_w'*Dv)';
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
108 obj.du_e = kr(d1_r_u,I_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
109 obj.dv_e = (e_e'*Dv)';
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
110 obj.du_s = (e_s'*Du)';
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
111 obj.dv_s = kr(I_u,d1_l_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
112 obj.du_n = (e_n'*Du)';
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
113 obj.dv_n = kr(I_u,d1_r_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
114
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
115
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
116 % Metric coefficients
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
117 coords = g.points();
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
118 x = coords(:,1);
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
119 y = coords(:,2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
120
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
121 x_u = Du*x;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
122 x_v = Dv*x;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
123 y_u = Du*y;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
124 y_v = Dv*y;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
125
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
126 J = x_u.*y_v - x_v.*y_u;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
127 a11 = 1./J .* (x_v.^2 + y_v.^2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
128 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
129 a22 = 1./J .* (x_u.^2 + y_u.^2);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
130 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
131
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
132 obj.x_u = x_u;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
133 obj.x_v = x_v;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
134 obj.y_u = y_u;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
135 obj.y_v = y_v;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
136
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
137
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
138 % Assemble full operators
551
c5a7a13c03dc Switch to using spdiag instead of spdiags
Jonatan Werpers <jonatan@werpers.com>
parents: 539
diff changeset
139 L_12 = spdiag(a12);
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
140 Duv = Du*L_12*Dv;
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
141 Dvu = Dv*L_12*Du;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
142
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
143 Duu = sparse(m_tot);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
144 Dvv = sparse(m_tot);
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
145 ind = grid.funcToMatrix(g, 1:m_tot);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
146
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
147 for i = 1:m_v
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
148 D = D2_u(a11(ind(:,i)));
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
149 p = ind(:,i);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
150 Duu(p,p) = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
151 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
152
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
153 for i = 1:m_u
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
154 D = D2_v(a22(ind(i,:)));
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
155 p = ind(i,:);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
156 Dvv(p,p) = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
157 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
158
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
159
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
160 % Physical operators
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
161 obj.J = spdiag(J);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
162 obj.Ji = spdiag(1./J);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
163
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
164 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
165 obj.H = obj.J*kr(H_u,H_v);
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
166 obj.Hi = obj.Ji*kr(Hi_u,Hi_v);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
167
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
168 obj.e_w = e_w;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
169 obj.e_e = e_e;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
170 obj.e_s = e_s;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
171 obj.e_n = e_n;
376
b0e95c81943e removed the white space
Ylva Rydin <ylva.rydin@telia.com>
parents: 371
diff changeset
172
557
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
173 %% normal derivatives
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
174 I_w = ind(1,:);
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
175 I_e = ind(end,:);
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
176 I_s = ind(:,1);
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
177 I_n = ind(:,end);
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
178
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
179 a11_w = spdiag(a11(I_w));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
180 a12_w = spdiag(a12(I_w));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
181 a11_e = spdiag(a11(I_e));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
182 a12_e = spdiag(a12(I_e));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
183 a22_s = spdiag(a22(I_s));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
184 a12_s = spdiag(a12(I_s));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
185 a22_n = spdiag(a22(I_n));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
186 a12_n = spdiag(a12(I_n));
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
187
565
f4b0d0e84305 Add missing square
Jonatan Werpers <jonatan@werpers.com>
parents: 564
diff changeset
188 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
f4b0d0e84305 Add missing square
Jonatan Werpers <jonatan@werpers.com>
parents: 564
diff changeset
189 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
f4b0d0e84305 Add missing square
Jonatan Werpers <jonatan@werpers.com>
parents: 564
diff changeset
190 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
f4b0d0e84305 Add missing square
Jonatan Werpers <jonatan@werpers.com>
parents: 564
diff changeset
191 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
562
11d8d6ccbcd7 Add boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 561
diff changeset
192
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
193 obj.d_w = -1*(spdiag(1./s_w)*(a11_w*obj.du_w' + a12_w*obj.dv_w'))';
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
194 obj.d_e = (spdiag(1./s_e)*(a11_e*obj.du_e' + a12_e*obj.dv_e'))';
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
195 obj.d_s = -1*(spdiag(1./s_s)*(a22_s*obj.dv_s' + a12_s*obj.du_s'))';
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
196 obj.d_n = (spdiag(1./s_n)*(a22_n*obj.dv_n' + a12_n*obj.du_n'))';
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
197
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
198 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
199 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
384
32df00102268 Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 360
diff changeset
200
562
11d8d6ccbcd7 Add boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 561
diff changeset
201 %% Boundary inner products
564
9bf49338f8e6 Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 563
diff changeset
202 obj.H_w = H_v*spdiag(s_w);
9bf49338f8e6 Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 563
diff changeset
203 obj.H_e = H_v*spdiag(s_e);
9bf49338f8e6 Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 563
diff changeset
204 obj.H_s = H_u*spdiag(s_s);
9bf49338f8e6 Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents: 563
diff changeset
205 obj.H_n = H_u*spdiag(s_n);
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
206
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
207 % Misc.
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
208 obj.m = m;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
209 obj.h = [h_u h_v];
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
210 obj.order = order;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
211 obj.grid = g;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
212
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
213 obj.a = a;
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
214 obj.b = b;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
215 obj.a11 = a11;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
216 obj.a12 = a12;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
217 obj.a22 = a22;
27
97a638f91fb8 Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents: 5
diff changeset
218 obj.lambda = lambda;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
219
389
42c89b5eedc0 Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents: 384
diff changeset
220 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
42c89b5eedc0 Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents: 384
diff changeset
221 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
222 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
223
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
224
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
225 % Closure functions return the opertors applied to the own doamin to close the boundary
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
226 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
227 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
228 % type is a string specifying the type of boundary condition if there are several.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
229 % data is a function returning the data that should be applied at the boundary.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
230 % neighbour_scheme is an instance of Scheme that should be interfaced to.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
231 % neighbour_boundary is a string specifying which boundary to interface to.
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
232 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
233 default_arg('type','neumann');
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
234 default_arg('parameter', []);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
235
567
33b962620e24 Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents: 566
diff changeset
236 [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
237 switch type
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
238 % Dirichlet boundary condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
239 case {'D','d','dirichlet'}
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
240 tuning = 1.2;
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
241 % tuning = 20.2;
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
242
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
243 b1 = gamm*obj.lambda./obj.a11.^2;
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
244 b2 = gamm*obj.lambda./obj.a22.^2;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
245
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
246 tau1 = tuning * spdiag(-1./b1 - 1./b2);
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
247 tau2 = 1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
248
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
249 tau = (tau1*e + tau2*d)*H_b;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
250
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
251 closure = obj.a*obj.Hi*tau*e';
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
252 penalty = -obj.a*obj.Hi*tau;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
253
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
254
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
255 % Neumann boundary condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
256 case {'N','n','neumann'}
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
257 tau1 = -1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
258 tau2 = 0;
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
259 tau = (tau1*e + tau2*d)*H_b;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
260
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
261 closure = obj.a*obj.Hi*tau*d';
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
262 penalty = -obj.a*obj.Hi*tau;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
263
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
264
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
265 % Unknown, boundary condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
266 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
267 error('No such boundary condition: type = %s',type);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
268 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
269 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
270
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
271 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
272 % u denotes the solution in the own domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
273 % v denotes the solution in the neighbour domain
5
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
274 tuning = 1.2;
27
97a638f91fb8 Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents: 5
diff changeset
275 % tuning = 20.2;
567
33b962620e24 Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents: 566
diff changeset
276 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
33b962620e24 Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents: 566
diff changeset
277 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
278
5
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
279 u = obj;
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
280 v = neighbour_scheme;
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
281
85
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
282 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
283 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
284 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
285 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
5
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
286
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
287 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
288 tau1 = tuning * spdiag(tau1);
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
289 tau2 = 1/2;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
290
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
291 sig1 = -1/2;
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
292 sig2 = 0;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
293
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
294 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
295 sig = (sig1*e_u + sig2*d_u)*H_b_u;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
296
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
297 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
298 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
299 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
300
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
301 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
302 % The right boundary is considered the positive boundary
85
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
303 %
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
304 % I -- the indecies of the boundary points in the grid matrix
567
33b962620e24 Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents: 566
diff changeset
305 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary)
85
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
306
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
307 % gridMatrix = zeros(obj.m(2),obj.m(1));
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
308 % gridMatrix(:) = 1:numel(gridMatrix);
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
309
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
310 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
85
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
311
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
312 switch boundary
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
313 case 'w'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
314 e = obj.e_w;
561
3a13916f8ff0 Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 560
diff changeset
315 d = obj.d_w;
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
316 H_b = obj.H_w;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
317 I = ind(1,:);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
318 case 'e'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
319 e = obj.e_e;
561
3a13916f8ff0 Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 560
diff changeset
320 d = obj.d_e;
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
321 H_b = obj.H_e;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
322 I = ind(end,:);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
323 case 's'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
324 e = obj.e_s;
561
3a13916f8ff0 Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 560
diff changeset
325 d = obj.d_s;
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
326 H_b = obj.H_s;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
327 I = ind(:,1)';
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
328 case 'n'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
329 e = obj.e_n;
561
3a13916f8ff0 Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents: 560
diff changeset
330 d = obj.d_n;
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
331 H_b = obj.H_n;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
332 I = ind(:,end)';
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
333 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
334 error('No such boundary: boundary = %s',boundary);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
335 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
336
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
337 switch boundary
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
338 case {'w','e'}
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
339 gamm = obj.gamm_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
340 case {'s','n'}
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
341 gamm = obj.gamm_v;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
342 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
343 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
344
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
345 function N = size(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
346 N = prod(obj.m);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
347 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
348 end
566
9c98a0526afc Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents: 565
diff changeset
349 end