Mercurial > repos > public > sbplib

annotate +scheme/LaplaceCurvilinear.m @ 558:54c775c3348a feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up use of the sign in boundary and interface routines
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Aug 2017 12:48:30 +0200
parents 2a856a589510
children 6132c52bf923
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
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
188 d_w = -1*(a11_w*obj.du_w' + a12_w*obj.dv_w')';
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
189 d_e = (a11_e*obj.du_e' + a12_e*obj.dv_e')';
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
190 d_s = -1*(a22_s*obj.dv_s' + a12_s*obj.du_s')';
2a856a589510 Add creation of physical boundary normal derivatives
Jonatan Werpers <jonatan@werpers.com>
parents: 554
diff changeset
191 d_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
192
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
193 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
194 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
195
553
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
196
63d5c07943dd Reorder and restructure constructor
Jonatan Werpers <jonatan@werpers.com>
parents: 552
diff changeset
197 % Misc.
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
198 obj.m = m;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
199 obj.h = [h_u h_v];
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
200 obj.order = order;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
201 obj.grid = g;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
202
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
203 obj.a = a;
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
204 obj.b = b;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
205 obj.a11 = a11;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
206 obj.a12 = a12;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
207 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
208 obj.lambda = lambda;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
209
389
42c89b5eedc0 Add borrowing constants for D2 operators in D4Variable
Jonatan Werpers <jonatan@werpers.com>
parents: 384
diff changeset
210 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
211 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
212 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
213
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
214
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
215 % 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
216 % 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
217 % 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
218 % 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
219 % 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
220 % neighbour_scheme is an instance of Scheme that should be interfaced to.
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
221 % 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
222 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
223 default_arg('type','neumann');
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
224 default_arg('parameter', []);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
225
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
226 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
227 switch type
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
228 % Dirichlet boundary condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
229 case {'D','d','dirichlet'}
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
230 % v denotes the solution in the neighbour domain
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
231 tuning = 1.2;
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
232 % tuning = 20.2;
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
233 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
234
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
235 A_n = spdiag(a_n);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
236 A_t = spdiag(a_t);
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
237
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
238 F = s*(A_n*d_n' + A_t*d_t')';
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
239
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
240 u = obj;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
241
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
242 b1 = gamm*u.lambda./u.a11.^2;
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
243 b2 = gamm*u.lambda./u.a22.^2;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
244
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
245 tau = -1./b1 - 1./b2;
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
246 tau = tuning * spdiag(tau);
360
447ceb41fb65 Fixed penalty parameter in Wave2Dcurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 347
diff changeset
247 sig1 = 1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
248
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
249 penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
250
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
251 closure = obj.Ji*obj.a * penalty_parameter_1*e';
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
252 penalty = -obj.Ji*obj.a * penalty_parameter_1;
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'}
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
257 A_n = spdiag(a_n);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
258 A_t = spdiag(a_t);
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
259 d = s*(A_n * d_n' + A_t*d_t')';
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
260
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
261 tau1 = -1;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
262 tau2 = 0;
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
263 tau = s*obj.a*obj.Ji*(tau1*e + tau2*d);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
264
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
265 closure = halfnorm_inv*tau*d';
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
266 penalty = -halfnorm_inv*tau;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
267
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
268 % Characteristic boundary condition
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
269 case {'characteristic', 'char', 'c'}
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
270 default_arg('parameter', 1);
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
271 beta = parameter;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
272
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
273 A_n = spdiag(a_n);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
274 A_t = spdiag(a_t);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
275 d = s*(A_n * d_n' + A_t*d_t')'; % outward facing normal derivative
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
276
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
277 tau = -obj.a * 1/beta*obj.Ji*e;
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
278
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
279 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
347
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
280 closure{2} = halfnorm_inv*tau*beta*d';
85c2fe06d551 Implemented characteristic form BC in Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 337
diff changeset
281 penalty = -halfnorm_inv*tau;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
282
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
283 % Unknown, boundary condition
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
284 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
285 error('No such boundary condition: type = %s',type);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
286 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
287 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
288
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
289 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
290 % u denotes the solution in the own domain
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
291 % 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
292 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
293 % tuning = 20.2;
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
294 [e_u, d_n_u, d_t_u, a_n_u, a_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
295 [e_v, d_n_v, d_t_v, a_n_v, a_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
296
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
297 A_n_u = spdiag(a_n_u);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
298 A_t_u = spdiag(a_t_u);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
299 A_n_v = spdiag(a_n_v);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
300 A_t_v = spdiag(a_t_v);
27
97a638f91fb8 Added function spdiag(). Fixed a bunch of bugs in the Wave2dCurve scheme.
Jonatan Werpers <jonatan@werpers.com>
parents: 5
diff changeset
301
558
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
302 F_u = s_u*(A_n_u * d_n_u' + A_t_u*d_t_u')';
54c775c3348a Clean up use of the sign in boundary and interface routines
Jonatan Werpers <jonatan@werpers.com>
parents: 557
diff changeset
303 F_v = s_v*(A_n_v * d_n_v' + A_t_v*d_t_v')';
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
304
5
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
305 u = obj;
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
306 v = neighbour_scheme;
5f6b0b6a012b First try at interface implementation in WaveCurve2s
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
307
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
308 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
309 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
310 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
311 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
312
82
df642750e8f7 Added Dirichelt BC to Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 55
diff changeset
313 tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
314 tau = tuning * spdiag(tau);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
315 sig1 = 1/2;
83
08fbc284718f Fixed bugs.
Jonatan Werpers <jonatan@werpers.com>
parents: 82
diff changeset
316 sig2 = -1/2;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
317
96
19d0c9325a3e Wave2dCurve: Removed some comments.
Jonatan Werpers <jonatan@werpers.com>
parents: 89
diff changeset
318 penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
319 penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
320
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
321
452
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
322 closure = obj.Ji*obj.a * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u');
56eb7c088bd4 Allow any coefficient in LaplaceCurvilinear
Jonatan Werpers <jonatan@werpers.com>
parents: 450
diff changeset
323 penalty = obj.Ji*obj.a * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v');
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
324 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
325
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
326 % 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
327 % 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
328 %
643bc513b8b8 Wave2DCurve: Fixed error in tau. Also made random guess at fixing theory from paper.
Jonatan Werpers <jonatan@werpers.com>
parents: 83
diff changeset
329 % I -- the indecies of the boundary points in the grid matrix
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
330 function [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = 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
331
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
332 % 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
333 % gridMatrix(:) = 1:numel(gridMatrix);
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
334
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
335 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
336
0
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'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
339 e = obj.e_w;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
340 d_n = obj.du_w;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
341 d_t = obj.dv_w;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
342 s = -1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
343
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
344 I = ind(1,:);
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
345 a_n = obj.a11(I);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
346 a_t = obj.a12(I);
384
32df00102268 Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 360
diff changeset
347 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
348 case 'e'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
349 e = obj.e_e;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
350 d_n = obj.du_e;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
351 d_t = obj.dv_e;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
352 s = 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
353
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
354 I = ind(end,:);
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
355 a_n = obj.a11(I);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
356 a_t = obj.a12(I);
384
32df00102268 Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 360
diff changeset
357 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
358 case 's'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
359 e = obj.e_s;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
360 d_n = obj.dv_s;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
361 d_t = obj.du_s;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
362 s = -1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
363
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
364 I = ind(:,1)';
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
365 a_n = obj.a22(I);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
366 a_t = obj.a12(I);
384
32df00102268 Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 360
diff changeset
367 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
368 case 'n'
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
369 e = obj.e_n;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
370 d_n = obj.dv_n;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
371 d_t = obj.du_n;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
372 s = 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
373
337
e070ebd94d9d Updated Wave2dCurve to use the new grid classes.
Jonatan Werpers <jonatan@werpers.com>
parents: 303
diff changeset
374 I = ind(:,end)';
554
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
375 a_n = obj.a22(I);
6e71d4f8b155 Better names for the metric coefficients
Jonatan Werpers <jonatan@werpers.com>
parents: 553
diff changeset
376 a_t = obj.a12(I);
384
32df00102268 Fix bug in Characteristic BC for Wave2dCurve.
Jonatan Werpers <jonatan@werpers.com>
parents: 360
diff changeset
377 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
378 otherwise
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
379 error('No such boundary: boundary = %s',boundary);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
380 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
381
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
382 switch boundary
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
383 case {'w','e'}
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
384 halfnorm_inv_n = obj.Hiu;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
385 halfnorm_inv_t = obj.Hiv;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
386 halfnorm_t = obj.Hv;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
387 gamm = obj.gamm_u;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
388 case {'s','n'}
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
389 halfnorm_inv_n = obj.Hiv;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
390 halfnorm_inv_t = obj.Hiu;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
391 halfnorm_t = obj.Hu;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
392 gamm = obj.gamm_v;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
393 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
394 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
395
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
396 function N = size(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
397 N = prod(obj.m);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
398 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
399
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
400
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
401 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
402 end