Mimetic Operators Library Enhanced 4.0
Loading...
Searching...
No Matches
Public Member Functions | List of all members
Gradient Class Reference

#include <gradient.h>

Inheritance diagram for Gradient:
Inheritance graph
[legend]

Public Member Functions

 Gradient (u16 k, u32 m, Real dx)
 
 Gradient (u16 k, u32 m, u32 n, Real dx, Real dy)
 
 Gradient (u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)
 
vec getP ()
 

Detailed Description

Definition at line 7 of file gradient.h.

Constructor & Destructor Documentation

◆ Gradient() [1/3]

Gradient::Gradient ( u16 k,
u32 m,
Real dx )

Definition at line 4 of file gradient.cpp.

4 : 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}

◆ Gradient() [2/3]

Gradient::Gradient ( u16 k,
u32 m,
u32 n,
Real dx,
Real dy )

Definition at line 225 of file gradient.cpp.

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}
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

◆ Gradient() [3/3]

Gradient::Gradient ( u16 k,
u32 m,
u32 n,
u32 o,
Real dx,
Real dy,
Real dz )

Definition at line 253 of file gradient.cpp.

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}

Member Function Documentation

◆ getP()

vec Gradient::getP ( )

Definition at line 287 of file gradient.cpp.

288{
289 return P;
290}
The documentation for this class was generated from the following files: