This repository was archived by the owner on Oct 18, 2022. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathRl_elem.m
More file actions
132 lines (120 loc) · 3.16 KB
/
Copy pathRl_elem.m
File metadata and controls
132 lines (120 loc) · 3.16 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
function [R, R_Q, R_U, R_UH] = Rl_elem(Q, U, uh, fd, lisb, risb, qd)
% function [R, R_Q, R_U, R_UH] = Rl_elem(Q, U, uh0, uh1, td, fd, sd, dx, qd)
%
% PURPOSE: Computes the bilinear form R((q,u,uh),m) : Qh x Uh x Mh x Mh -> Mh*
% restricted to an element and linearization wrt each input. This is:
%
% <[HH],m> = 0 forall m in Mh
%
% NOTE: On interior elements, computes a one-sided stabilized flux on both edges,
% on boundary edges, computes the flux H(q,uh)
%
% INPUTS:
% {Q,U} : basis coefficients for {grad(u),u} [nn{q,u}]
% uh : u|{0,1}, that is, the trace of u at the boundaries of the element
% fd : flux data [struct]
% lisb : true if left is boundary, else false [bool]
% risb : true if right is boundary, else false [bool]
% qd : quadrature data [struct]
%
% OUTPUTS:
% R : bilinear form [1]
% R_{Q,U,UH} : linearization of R wrt {q,u,uh}
%
nnq = size(qd.qPhi, 2);
nnu = size(qd.uPhi, 2);
q0 = qd.qPhi0 * Q;
q1 = qd.qPhi1 * Q;
u0 = qd.uPhi0 * U;
u1 = qd.uPhi1 * U;
n0 = -1.;
n1 = 1.;
[h0, h_q0, h_UH0] = flux(q0, uh(1), fd);
h_u0 = 0.;
% dot with normal
h0 = h0*n0;
h_q0 = h_q0*n0;
h_UH0 = h_UH0*n0;
[s0, s_u0, s_UH0] = flux_stab(u0, uh(1), n0, fd);
f0 = h0 + s0;
f_Q0 = h_q0*qd.qPhi0;
f_U0 = (h_u0+s_u0)*qd.uPhi0;
f_UH0 = h_UH0 + s_UH0;
[h1, h_q1, h_UH1] = flux(q1, uh(2), fd);
h_u1 = 0.;
% dot with normal
h1 = h1*n1;
h_q1 = h_q1*n1;
h_UH1 = h_UH1*n1;
[s1, s_u1, s_UH1] = flux_stab(u1, uh(2), n1, fd);
f1 = h1 + s1;
f_Q1 = h_q1*qd.qPhi1;
f_U1 = (h_u1+s_u1)*qd.uPhi1;
f_UH1 = h_UH1 + s_UH1;
% Rl = <[HH],m>
R(1,1) = f0;
R(2,1) = f1;
R_Q(1,:) = f_Q0;
R_Q(2,:) = f_Q1;
R_U(1,:) = f_U0;
R_U(2,:) = f_U1;
R_UH(1,1) = f_UH0;
R_UH(2,2) = f_UH1;
%{
%This is the version with the adjoint consistency fix, but since we
%aren't using a Riemann solver on the boundaries, this code is
%deprecated
if lisb
[h0, h_q0, h_u0, h_UH0] = flux_bc(q0, u0, uh(1), n0, fd);
if strcmp(scheme, 'dpg')
% adjoint consistency fix
taut = 0.5*(fd.a*n0+abs(fd.a*n0)) + fd.b/fd.vl;
c0 = -0.5*(fd.a*n0+abs(fd.a*n0))/taut;
else
c0 = 1.0;
end
else
[h0, h_q0, h_UH0] = flux(q0, uh(1), fd);
h_u0 = 0.;
% dot with normal
h0 = h0*n0;
h_q0 = h_q0*n0;
h_UH0 = h_UH0*n0;
if strcmp(scheme, 'dpg')
c0 = fd.c;
else
c0 = 1.0;
end
end
[s0, s_u0, s_UH0] = flux_stab(u0, uh(1), n0, fd);
f0 = h0 + c0*s0;
f_Q0 = h_q0*qd.qPhi0;
f_U0 = (h_u0+c0*s_u0)*qd.uPhi0;
f_UH0 = h_UH0 + c0*s_UH0;
if risb
[h1, h_q1, h_u1, h_UH1] = flux_bc(q1, u1, uh(2), n1, fd);
if strcmp(scheme, 'dpg')
taut = 0.5*(fd.a*n1+abs(fd.a*n1)) + fd.b/fd.vl;
c1 = -0.5*(fd.a*n1+abs(fd.a*n1))/taut;
else
c1 = 1.0;
end
else
[h1, h_q1, h_UH1] = flux(q1, uh(2), fd);
h_u1 = 0.;
% dot with normal
h1 = h1*n1;
h_q1 = h_q1*n1;
h_UH1 = h_UH1*n1;
if strcmp(scheme, 'dpg')
c1 = fd.c;
else
c1 = 1.0;
end
end
[s1, s_u1, s_UH1] = flux_stab(u1, uh(2), n1, fd);
f1 = h1 + c1*s1;
f_Q1 = h_q1*qd.qPhi1;
f_U1 = (h_u1+c1*s_u1)*qd.uPhi1;
f_UH1 = h_UH1 + c1*s_UH1;
%}