Mercurial > repos > public > sbplib
annotate +scheme/LaplaceCurvilinear.m @ 743:f4595f14d696 feature/utux2D
Change schemes to work for special coefficients.
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Mon, 21 May 2018 23:11:34 -0700 |
parents | 07f8311374c6 |
children | 14b093a344eb |
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 | 2 properties |
3 m % Number of points in each direction, possibly a vector | |
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 | 8 order % Order accuracy for the approximation |
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 | 23 J, Ji |
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 | 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 | 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 | 33 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. |
34 Hiu, Hiv | |
35 du_w, dv_w | |
36 du_e, dv_e | |
37 du_s, dv_s | |
38 du_n, dv_n | |
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 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
41 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
42 interpolation_type |
0 | 43 end |
44 | |
45 methods | |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
46 % Implements a*div(b*grad(u)) as a SBP scheme |
539 | 47 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) |
48 | |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
49 function obj = LaplaceCurvilinear(g ,order, a, b, opSet, interpolation_type) |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
50 default_arg('opSet',@sbp.D2Variable); |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
51 default_arg('a', 1); |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
52 default_arg('b', 1); |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
53 default_arg('interpolation_type','AWW'); |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
54 |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
55 if b ~=1 |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
56 error('Not implemented yet') |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
57 end |
0 | 58 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
59 % assert(isa(g, 'grid.Curvilinear')) |
743
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
60 if isa(a, 'function_handle') |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
61 a = grid.evalOn(g, a); |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
62 a = spdiag(a); |
f4595f14d696
Change schemes to work for special coefficients.
Martin Almquist <malmquist@stanford.edu>
parents:
720
diff
changeset
|
63 end |
0 | 64 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
65 m = g.size(); |
0 | 66 m_u = m(1); |
67 m_v = m(2); | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
68 m_tot = g.N(); |
0 | 69 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
70 h = g.scaling(); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
71 h_u = h(1); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
72 h_v = h(2); |
0 | 73 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
74 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
75 % 1D operators |
303
f18142c1530b
Fixed Wave2dCurve for new operators. Some cleanup in d2_variable_4.
Jonatan Werpers <jonatan@werpers.com>
parents:
96
diff
changeset
|
76 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
|
77 ops_v = opSet(m_v, {0, 1}, order); |
0 | 78 |
79 I_u = speye(m_u); | |
80 I_v = speye(m_v); | |
81 | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
82 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
|
83 D2_u = ops_u.D2; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
84 H_u = ops_u.H; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
85 Hi_u = ops_u.HI; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
86 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
|
87 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
|
88 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
|
89 d1_r_u = ops_u.d1_r; |
0 | 90 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
91 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
|
92 D2_v = ops_v.D2; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
93 H_v = ops_v.H; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
94 Hi_v = ops_v.HI; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
95 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
|
96 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
|
97 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
|
98 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
|
99 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
100 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
101 % Logical operators |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
102 Du = kr(D1_u,I_v); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
103 Dv = kr(I_u,D1_v); |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
104 obj.Hu = kr(H_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
105 obj.Hv = kr(I_u,H_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
106 obj.Hiu = kr(Hi_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
107 obj.Hiv = kr(I_u,Hi_v); |
0 | 108 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
109 e_w = kr(e_l_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
110 e_e = kr(e_r_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
111 e_s = kr(I_u,e_l_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
112 e_n = kr(I_u,e_r_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
113 obj.du_w = kr(d1_l_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
114 obj.dv_w = (e_w'*Dv)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
115 obj.du_e = kr(d1_r_u,I_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
116 obj.dv_e = (e_e'*Dv)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
117 obj.du_s = (e_s'*Du)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
118 obj.dv_s = kr(I_u,d1_l_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
119 obj.du_n = (e_n'*Du)'; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
120 obj.dv_n = kr(I_u,d1_r_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
121 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
122 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
123 % Metric coefficients |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
124 coords = g.points(); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
125 x = coords(:,1); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
126 y = coords(:,2); |
0 | 127 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
128 x_u = Du*x; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
129 x_v = Dv*x; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
130 y_u = Du*y; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
131 y_v = Dv*y; |
0 | 132 |
133 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
|
134 a11 = 1./J .* (x_v.^2 + y_v.^2); |
0 | 135 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); |
136 a22 = 1./J .* (x_u.^2 + y_u.^2); | |
137 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); | |
138 | |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
139 obj.x_u = x_u; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
140 obj.x_v = x_v; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
141 obj.y_u = y_u; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
142 obj.y_v = y_v; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
143 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
144 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
145 % Assemble full operators |
551
c5a7a13c03dc
Switch to using spdiag instead of spdiags
Jonatan Werpers <jonatan@werpers.com>
parents:
539
diff
changeset
|
146 L_12 = spdiag(a12); |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
147 Duv = Du*L_12*Dv; |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
148 Dvu = Dv*L_12*Du; |
0 | 149 |
150 Duu = sparse(m_tot); | |
151 Dvv = sparse(m_tot); | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
152 ind = grid.funcToMatrix(g, 1:m_tot); |
0 | 153 |
154 for i = 1:m_v | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
155 D = D2_u(a11(ind(:,i))); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
156 p = ind(:,i); |
0 | 157 Duu(p,p) = D; |
158 end | |
159 | |
160 for i = 1:m_u | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
161 D = D2_v(a22(ind(i,:))); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
162 p = ind(i,:); |
0 | 163 Dvv(p,p) = D; |
164 end | |
165 | |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
166 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
167 % Physical operators |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
168 obj.J = spdiag(J); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
169 obj.Ji = spdiag(1./J); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
170 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
171 obj.D = obj.Ji*a*(Duu + Duv + Dvu + Dvv); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
172 obj.H = obj.J*kr(H_u,H_v); |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
173 obj.Hi = obj.Ji*kr(Hi_u,Hi_v); |
0 | 174 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
175 obj.e_w = e_w; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
176 obj.e_e = e_e; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
177 obj.e_s = e_s; |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
178 obj.e_n = e_n; |
376 | 179 |
557
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
180 %% normal derivatives |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
181 I_w = ind(1,:); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
182 I_e = ind(end,:); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
183 I_s = ind(:,1); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
184 I_n = ind(:,end); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
185 |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
186 a11_w = spdiag(a11(I_w)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
187 a12_w = spdiag(a12(I_w)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
188 a11_e = spdiag(a11(I_e)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
189 a12_e = spdiag(a12(I_e)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
190 a22_s = spdiag(a22(I_s)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
191 a12_s = spdiag(a12(I_s)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
192 a22_n = spdiag(a22(I_n)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
193 a12_n = spdiag(a12(I_n)); |
2a856a589510
Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents:
554
diff
changeset
|
194 |
565 | 195 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); |
196 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); | |
197 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); | |
198 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
|
199 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
200 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
|
201 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
|
202 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
|
203 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
|
204 |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
205 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
|
206 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
|
207 |
562
11d8d6ccbcd7
Add boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
561
diff
changeset
|
208 %% Boundary inner products |
564
9bf49338f8e6
Fix bug in creation of boundary inner products
Jonatan Werpers <jonatan@werpers.com>
parents:
563
diff
changeset
|
209 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
|
210 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
|
211 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
|
212 obj.H_n = H_u*spdiag(s_n); |
553
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
213 |
63d5c07943dd
Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents:
552
diff
changeset
|
214 % Misc. |
0 | 215 obj.m = m; |
216 obj.h = [h_u h_v]; | |
217 obj.order = order; | |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
218 obj.grid = g; |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
219 obj.interpolation_type = interpolation_type; |
0 | 220 |
452
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
221 obj.a = a; |
56eb7c088bd4
Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents:
450
diff
changeset
|
222 obj.b = b; |
0 | 223 obj.a11 = a11; |
224 obj.a12 = a12; | |
225 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
|
226 obj.lambda = lambda; |
0 | 227 |
389
42c89b5eedc0
Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents:
384
diff
changeset
|
228 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
|
229 obj.gamm_v = h_v*ops_v.borrowing.M.d1; |
0 | 230 end |
231 | |
232 | |
233 % Closure functions return the opertors applied to the own doamin to close the boundary | |
234 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. | |
235 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. | |
236 % type is a string specifying the type of boundary condition if there are several. | |
237 % data is a function returning the data that should be applied at the boundary. | |
238 % neighbour_scheme is an instance of Scheme that should be interfaced to. | |
239 % 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
|
240 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) |
0 | 241 default_arg('type','neumann'); |
347
85c2fe06d551
Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
337
diff
changeset
|
242 default_arg('parameter', []); |
0 | 243 |
567
33b962620e24
Remove unneeded sign information from get_boundary_ops
Jonatan Werpers <jonatan@werpers.com>
parents:
566
diff
changeset
|
244 [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); |
0 | 245 switch type |
246 % Dirichlet boundary condition | |
247 case {'D','d','dirichlet'} | |
82
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
248 tuning = 1.2; |
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
249 % tuning = 20.2; |
df642750e8f7
Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents:
55
diff
changeset
|
250 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
251 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
|
252 b2 = gamm*obj.lambda./obj.a22.^2; |
0 | 253 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
254 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
|
255 tau2 = 1; |
0 | 256 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
257 tau = (tau1*e + tau2*d)*H_b; |
0 | 258 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
259 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
|
260 penalty = -obj.a*obj.Hi*tau; |
0 | 261 |
262 | |
263 % Neumann boundary condition | |
264 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
|
265 tau1 = -1; |
0 | 266 tau2 = 0; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
267 tau = (tau1*e + tau2*d)*H_b; |
0 | 268 |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
269 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
|
270 penalty = -obj.a*obj.Hi*tau; |
0 | 271 |
272 | |
273 % Unknown, boundary condition | |
274 otherwise | |
275 error('No such boundary condition: type = %s',type); | |
276 end | |
277 end | |
278 | |
279 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
280 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
281 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
282 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
283 Hi = obj.Hi; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
284 a = obj.a; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
285 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
286 m_u = length(H_b_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
287 m_v = length(H_b_v); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
288 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
289 grid_ratio = m_u/m_v; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
290 if grid_ratio ~= 1 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
291 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
292 [ms, index] = sort([m_u, m_v]); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
293 orders = [obj.order, neighbour_scheme.order]; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
294 orders = orders(index); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
295 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
296 switch obj.interpolation_type |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
297 case 'MC' |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
298 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
299 if grid_ratio < 1 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
300 I_v2u_good = interpOpSet.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
301 I_v2u_bad = interpOpSet.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
302 I_u2v_good = interpOpSet.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
303 I_u2v_bad = interpOpSet.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
304 elseif grid_ratio > 1 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
305 I_v2u_good = interpOpSet.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
306 I_v2u_bad = interpOpSet.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
307 I_u2v_good = interpOpSet.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
308 I_u2v_bad = interpOpSet.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
309 end |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
310 case 'AWW' |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
311 %String 'C2F' indicates that ICF2 is more accurate. |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
312 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
313 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
314 if grid_ratio < 1 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
315 % Local is coarser than neighbour |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
316 I_v2u_good = interpOpSetF2C.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
317 I_v2u_bad = interpOpSetC2F.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
318 I_u2v_good = interpOpSetC2F.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
319 I_u2v_bad = interpOpSetF2C.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
320 elseif grid_ratio > 1 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
321 % Local is finer than neighbour |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
322 I_v2u_good = interpOpSetC2F.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
323 I_v2u_bad = interpOpSetF2C.IC2F; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
324 I_u2v_good = interpOpSetF2C.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
325 I_u2v_bad = interpOpSetC2F.IF2C; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
326 end |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
327 otherwise |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
328 error(['Interpolation type ' obj.interpolation_type ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
329 ' is not available.' ]); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
330 end |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
331 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
332 else |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
333 % No interpolation required |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
334 I_v2u_good = speye(m_u,m_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
335 I_v2u_bad = speye(m_u,m_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
336 I_u2v_good = speye(m_u,m_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
337 I_u2v_bad = speye(m_u,m_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
338 end |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
339 |
0 | 340 % u denotes the solution in the own domain |
341 % 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
|
342 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
|
343 % tuning = 20.2; |
5
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
344 u = obj; |
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
345 v = neighbour_scheme; |
5f6b0b6a012b
First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
346 |
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
|
347 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
|
348 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
|
349 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
|
350 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
|
351 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
352 tau_u = -1./(4*b1_u) -1./(4*b2_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
353 tau_v = -1./(4*b1_v) -1./(4*b2_v); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
354 |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
355 tau_u = tuning * spdiag(tau_u); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
356 tau_v = tuning * spdiag(tau_v); |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
357 beta_u = tau_v; |
0 | 358 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
359 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
360 a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
361 a*1/2*Hi*d_u*H_b_u*e_u' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
362 -a*1/2*Hi*e_u*H_b_u*d_u'; |
0 | 363 |
720
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
364 penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
365 -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
366 -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ... |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
367 -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v'; |
07f8311374c6
Add interpolation to LaplaceCurveilinear.
Martin Almquist <malmquist@stanford.edu>
parents:
567
diff
changeset
|
368 |
0 | 369 |
370 end | |
371 | |
372 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. | |
373 % 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
|
374 % |
643bc513b8b8
Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents:
83
diff
changeset
|
375 % 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
|
376 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
|
377 |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
378 % 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
|
379 % gridMatrix(:) = 1:numel(gridMatrix); |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
380 |
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
381 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
|
382 |
0 | 383 switch boundary |
384 case 'w' | |
385 e = obj.e_w; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
386 d = obj.d_w; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
387 H_b = obj.H_w; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
388 I = ind(1,:); |
0 | 389 case 'e' |
390 e = obj.e_e; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
391 d = obj.d_e; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
392 H_b = obj.H_e; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
393 I = ind(end,:); |
0 | 394 case 's' |
395 e = obj.e_s; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
396 d = obj.d_s; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
397 H_b = obj.H_s; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
398 I = ind(:,1)'; |
0 | 399 case 'n' |
400 e = obj.e_n; | |
561
3a13916f8ff0
Switch to using boundary ops for normal derivative
Jonatan Werpers <jonatan@werpers.com>
parents:
560
diff
changeset
|
401 d = obj.d_n; |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
402 H_b = obj.H_n; |
337
e070ebd94d9d
Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents:
303
diff
changeset
|
403 I = ind(:,end)'; |
0 | 404 otherwise |
405 error('No such boundary: boundary = %s',boundary); | |
406 end | |
407 | |
408 switch boundary | |
409 case {'w','e'} | |
410 gamm = obj.gamm_u; | |
411 case {'s','n'} | |
412 gamm = obj.gamm_v; | |
413 end | |
414 end | |
415 | |
416 function N = size(obj) | |
417 N = prod(obj.m); | |
418 end | |
419 end | |
566
9c98a0526afc
Switch implementation of boundary and interface to SBP notation
Jonatan Werpers <jonatan@werpers.com>
parents:
565
diff
changeset
|
420 end |