Mimetic Operators Library Enhanced 4.0
Loading...
Searching...
No Matches
gradient.cpp
Go to the documentation of this file.
1#include "gradient.h"
2
3// 1-D Constructor
4Gradient::Gradient(u16 k, u32 m, Real dx) : sp_mat(m+1, m+2)
5{
6 assert(!(k%2));
7 assert(k > 1 && k < 9);
8 assert(m >= 2*k);
9
10 switch (k) {
11 case 2:
12 // A
13 at(0, 0) = -8.0/3.0;
14 at(0, 1) = 3.0;
15 at(0, 2) = -1.0/3.0;
16 // A'
17 at(m, m+1) = 8.0/3.0;
18 at(m, m) = -3.0;
19 at(m, m-1) = 1.0/3.0;
20 // Middle
21 for (u32 i = 1; i < m; i++) {
22 at(i, i) = -1.0;
23 at(i, i+1) = 1.0;
24 }
25 // Weights
26 P << 3.0/8.0 << 9.0/8.0 << 1.0 << 9.0/8.0 << 3.0/8.0;
27 break;
28 case 4:
29 // A
30 at(0, 0) = -352.0/105.0;
31 at(0, 1) = 35.0/8.0;
32 at(0, 2) = -35.0/24.0;
33 at(0, 3) = 21.0/40.0;
34 at(0, 4) = -5.0/56.0;
35 at(1, 0) = 16.0/105.0;
36 at(1, 1) = -31.0/24.0;
37 at(1, 2) = 29.0/24.0;
38 at(1, 3) = -3.0/40.0;
39 at(1, 4) = 1.0/168.0;
40 // A'
41 at(m, m+1) = 352.0/105.0;
42 at(m, m) = -35.0/8.0;
43 at(m, m-1) = 35.0/24.0;
44 at(m, m-2) = -21.0/40.0;
45 at(m, m-3) = 5.0/56.0;
46 at(m-1, m+1) = -16.0/105.0;
47 at(m-1, m) = 31.0/24.0;
48 at(m-1, m-1) = -29.0/24.0;
49 at(m-1, m-2) = 3.0/40.0;
50 at(m-1, m-3) = -1.0/168.0;
51 // Middle
52 for (u32 i = 2; i < m-1; i++) {
53 at(i, i-1) = 1.0/24.0;
54 at(i, i) = -9.0/8.0;
55 at(i, i+1) = 9.0/8.0;
56 at(i, i+2) = -1.0/24.0;
57 }
58 // Weights
59 P << 1606.0/4535.0 << 941.0/766.0 << 1384.0/1541.0 << 1371.0/1346.0
60 << 701.0/700.0 << 1371.0/1346.0 << 1384.0/1541.0 << 941.0/766.0
61 << 1606.0/4535.0;
62 break;
63 case 6:
64 // A
65 at(0, 0) = -13016.0/3465.0;
66 at(0, 1) = 693.0/128.0;
67 at(0, 2) = -385.0/128.0;
68 at(0, 3) = 693.0/320.0;
69 at(0, 4) = -495.0/448.0;
70 at(0, 5) = 385.0/1152.0;
71 at(0, 6) = -63.0/1408.0;
72 at(1, 0) = 496.0/3465.0;
73 at(1, 1) = -811.0/640.0;
74 at(1, 2) = 449.0/384.0;
75 at(1, 3) = -29.0/960.0;
76 at(1, 4) = -11.0/448.0;
77 at(1, 5) = 13.0/1152.0;
78 at(1, 6) = -37.0/21120.0;
79 at(2, 0) = -8.0/385.0;
80 at(2, 1) = 179.0/1920.0;
81 at(2, 2) = -153.0/128.0;
82 at(2, 3) = 381.0/320.0;
83 at(2, 4) = -101.0/1344.0;
84 at(2, 5) = 1.0/128.0;
85 at(2, 6) = -3.0/7040.0;
86 // A'
87 at(m, m+1) = 13016.0/3465.0;
88 at(m, m) = -693.0/128.0;
89 at(m, m-1) = 385.0/128.0;
90 at(m, m-2) = -693.0/320.0;
91 at(m, m-3) = 495.0/448.0;
92 at(m, m-4) = -385.0/1152.0;
93 at(m, m-5) = 63.0/1408.0;
94 at(m-1, m+1) = -496.0/3465.0;
95 at(m-1, m) = 811.0/640.0;
96 at(m-1, m-1) = -449.0/384.0;
97 at(m-1, m-2) = 29.0/960.0;
98 at(m-1, m-3) = 11.0/448.0;
99 at(m-1, m-4) = -13.0/1152.0;
100 at(m-1, m-5) = 37.0/21120.0;
101 at(m-2, m+1) = 8.0/385.0;
102 at(m-2, m) = -179.0/1920.0;
103 at(m-2, m-1) = 153.0/128.0;
104 at(m-2, m-2) = -381.0/320.0;
105 at(m-2, m-3) = 101.0/1344.0;
106 at(m-2, m-4) = -1.0/128.0;
107 at(m-2, m-5) = 3.0/7040.0;
108 // Middle
109 for (u32 i = 3; i < m-2; i++) {
110 at(i, i-2) = -3.0/640.0;
111 at(i, i-1) = 25.0/384.0;
112 at(i, i) = -75.0/64.0;
113 at(i, i+1) = 75.0/64.0;
114 at(i, i+2) = -25.0/384.0;
115 at(i, i+3) = 3.0/640.0;
116 }
117 // Weights
118 P << 420249.0/1331069.0 << 2590978.0/1863105.0 << 882762.0/1402249.0
119 << 1677712.0/1359311.0 << 239985.0/261097.0 << 664189.0/657734.0
120 << 756049.0/754729.0 << 664189.0/657734.0 << 239985.0/261097.0
121 << 1677712.0/1359311.0 << 882762.0/1402249.0 << 2590978.0/1863105.0
122 << 420249.0/1331069.0;
123 break;
124 case 8:
125 // A
126 at(0, 0) = -4856215.0/1200963.0;
127 at(0, 1) = 45858154.0/7297397.0;
128 at(0, 2) = -23409299.0/4789435.0;
129 at(0, 3) = 3799178.0/719717.0;
130 at(0, 4) = -4892189.0/1089890.0;
131 at(0, 5) = 1789111.0/658879.0;
132 at(0, 6) = -1406819.0/1289899.0;
133 at(0, 7) = 1154863.0/4436807.0;
134 at(0, 8) = -2936602.0/105142673.0;
135 at(1, 0) = 86048.0/675675.0;
136 at(1, 1) = -131093.0/107520.0;
137 at(1, 2) = 5503131.0/5166017.0;
138 at(1, 3) = 305249.0/2136437.0;
139 at(1, 4) = -1763845.0/8250973.0;
140 at(1, 5) = 1562032.0/10745723.0;
141 at(1, 6) = -270419.0/4422611.0;
142 at(1, 7) = 2983.0/199680.0;
143 at(1, 8) = -2621.0/1612800.0;
144 at(2, 0) = -3776.0/225225.0;
145 at(2, 1) = 8707.0/107520.0;
146 at(2, 2) = -17947.0/15360.0;
147 at(2, 3) = 29319.0/25600.0;
148 at(2, 4) = -533.0/21504.0;
149 at(2, 5) = -263.0/9216.0;
150 at(2, 6) = 903.0/56320.0;
151 at(2, 7) = -283.0/66560.0;
152 at(2, 8) = 257.0/537600.0;
153 at(3, 0) = 32.0/9009.0;
154 at(3, 1) = -543.0/35840.0;
155 at(3, 2) = 265.0/3072.0;
156 at(3, 3) = -1233.0/1024.0;
157 at(3, 4) = 8625.0/7168.0;
158 at(3, 5) = -775.0/9216.0;
159 at(3, 6) = 639.0/56320.0;
160 at(3, 7) = -15.0/13312.0;
161 at(3, 8) = 1.0/21504.0;
162 // A'
163 at(m, m+1) = 4856215.0/1200963.0;
164 at(m, m) = -45858154.0/7297397.0;
165 at(m, m-1) = 23409299.0/4789435.0;
166 at(m, m-2) = -3799178.0/719717.0;
167 at(m, m-3) = 4892189.0/1089890.0;
168 at(m, m-4) = -1789111.0/658879.0;
169 at(m, m-5) = 1406819.0/1289899.0;
170 at(m, m-6) = -1154863.0/4436807.0;
171 at(m, m-7) = 2936602.0/105142673.0;
172 at(m-1, m+1) = -86048.0/675675.0;
173 at(m-1, m) = 131093.0/107520.0;
174 at(m-1, m-1) = -5503131.0/5166017.0;
175 at(m-1, m-2) = -305249.0/2136437.0;
176 at(m-1, m-3) = 1763845.0/8250973.0;
177 at(m-1, m-4) = -1562032.0/10745723.0;
178 at(m-1, m-5) = 270419.0/4422611.0;
179 at(m-1, m-6) = -2983.0/199680.0;
180 at(m-1, m-7) = 2621.0/1612800.0;
181 at(m-2, m+1) = 3776.0/225225.0;
182 at(m-2, m) = -8707.0/107520.0;
183 at(m-2, m-1) = 17947.0/15360.0;
184 at(m-2, m-2) = -29319.0/25600.0;
185 at(m-2, m-3) = 533.0/21504.0;
186 at(m-2, m-4) = 263.0/9216.0;
187 at(m-2, m-5) = -903.0/56320.0;
188 at(m-2, m-6) = 283.0/66560.0;
189 at(m-2, m-7) = -257.0/537600.0;
190 at(m-3, m+1) = -32.0/9009.0;
191 at(m-3, m) = 543.0/35840.0;
192 at(m-3, m-1) = -265.0/3072.0;
193 at(m-3, m-2) = 1233.0/1024.0;
194 at(m-3, m-3) = -8625.0/7168.0;
195 at(m-3, m-4) = 775.0/9216.0;
196 at(m-3, m-5) = -639.0/56320.0;
197 at(m-3, m-6) = 15.0/13312.0;
198 at(m-3, m-7) = -1.0/21504.0;
199 // Middle
200 for (u32 i = 4; i < m-3; i++) {
201 at(i, i-3) = 5.0/7168.0;
202 at(i, i-2) = -49.0/5120.0;
203 at(i, i-1) = 245.0/3072.0;
204 at(i, i) = -1225.0/1024.0;
205 at(i, i+1) = 1225.0/1024.0;
206 at(i, i+2) = -245.0/3072.0;
207 at(i, i+3) = 49.0/5120.0;
208 at(i, i+4) = -5.0/7168.0;
209 }
210 // Weights
211 P << 267425.0/904736.0 << 2307435.0/1517812.0 << 847667.0/3066027.0
212 << 4050911.0/2301238.0 << 498943.0/1084999.0 << 211042.0/170117.0
213 << 2065895.0/2191686.0 << 1262499.0/1258052.0 << 1314891.0/1312727.0
214 << 1262499.0/1258052.0 << 2065895.0/2191686.0 << 211042.0/170117.0
215 << 498943.0/1084999.0 << 4050911.0/2301238.0 << 847667.0/3066027.0
216 << 2307435.0/1517812.0 << 267425.0/904736.0;
217 break;
218 }
219
220 // Scaling
221 *this /= dx;
222}
223
224// 2-D Constructor
225Gradient::Gradient(u16 k, u32 m, u32 n, Real dx, Real dy)
226{
227 Gradient Gx(k, m, dx);
228 Gradient Gy(k, n, dy);
229
230 sp_mat Im = speye(m+2, m+2);
231 sp_mat In = speye(n+2, n+2);
232
233 Im.shed_row(0);
234 Im.shed_row(m);
235 In.shed_row(0);
236 In.shed_row(n);
237
238 sp_mat G1 = Utils::spkron(In, Gx);
239 sp_mat G2 = Utils::spkron(Gy, Im);
240
241 // Dimensions = 2*m*n+m+n, (m+2)*(n+2)
242 if (m != n)
243 *this = Utils::spjoin_cols(G1, G2);
244 else {
245 sp_mat A1(2, 1);
246 sp_mat A2(2, 1);
247 A1(0, 0) = A2(1, 0) = 1.0;
248 *this = Utils::spkron(A1, G1) + Utils::spkron(A2, G2);
249 }
250}
251
252// 3-D Constructor
253Gradient::Gradient(u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)
254{
255 Gradient Gx(k, m, dx);
256 Gradient Gy(k, n, dy);
257 Gradient Gz(k, o, dz);
258
259 sp_mat Im = speye(m+2, m+2);
260 sp_mat In = speye(n+2, n+2);
261 sp_mat Io = speye(o+2, o+2);
262
263 Im.shed_row(0);
264 Im.shed_row(m);
265 In.shed_row(0);
266 In.shed_row(n);
267 Io.shed_row(0);
268 Io.shed_row(o);
269
270 sp_mat G1 = Utils::spkron(Utils::spkron(Io, In), Gx);
271 sp_mat G2 = Utils::spkron(Utils::spkron(Io, Gy), Im);
272 sp_mat G3 = Utils::spkron(Utils::spkron(Gz, In), Im);
273
274 // Dimensions = 3*m*n*o+m*n+m*o+n*o, (m+2)*(n+2)*(o+2)
275 if ((m != n) || (n != o))
276 *this = Utils::spjoin_cols(Utils::spjoin_cols(G1, G2), G3);
277 else {
278 sp_mat A1(3, 1);
279 sp_mat A2(3, 1);
280 sp_mat A3(3, 1);
281 A1(0, 0) = A2(1, 0) = A3(2, 0) = 1.0;
282 *this = Utils::spkron(A1, G1) + Utils::spkron(A2, G2) + Utils::spkron(A3, G3);
283 }
284}
285
286// Returns weights
287vec Gradient::getP()
288{
289 return P;
290}
Gradient::getP
vec getP()
Definition gradient.cpp:287
Gradient::Gradient
Gradient(u16 k, u32 m, Real dx)
Definition gradient.cpp:4
Utils::spjoin_cols
static sp_mat spjoin_cols(const sp_mat &A, const sp_mat &B)
Definition utils.cpp:120
Utils::spkron
static sp_mat spkron(const sp_mat &A, const sp_mat &B)
Definition utils.cpp:53
Real
double Real
Definition utils.h:8