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

Mimetic Gradient operator. More...

#include <gradient.h>

Inheritance diagram for Gradient:
Inheritance graph
[legend]

Public Member Functions

 Gradient (u16 k, u32 m, Real dx)
 1-D Mimetic Gradient Constructor
 
 Gradient (u16 k, u32 m, u32 n, Real dx, Real dy)
 2-D Mimetic Gradient Constructor
 
 Gradient (u16 k, u32 m, u32 n, u32 o, Real dx, Real dy, Real dz)
 3-D Mimetic Gradient Constructor
 
vec getP ()
 Returns the weights used in the Mimeitc Gradient Operators.
 

Detailed Description

Mimetic Gradient operator.

Definition at line 27 of file gradient.h.

Constructor & Destructor Documentation

◆ Gradient() [1/3]

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

1-D Mimetic Gradient Constructor

Parameters
kOrder of accuracy
mNumber of cells
dxSpacing between cells

Definition at line 20 of file gradient.cpp.

20 : sp_mat(m + 1, m + 2) {
21 assert(!(k % 2));
22 assert(k > 1 && k < 9);
23 assert(m >= 2 * k);
24
25 switch (k) {
26 case 2:
27 // A
28 at(0, 0) = -8.0 / 3.0;
29 at(0, 1) = 3.0;
30 at(0, 2) = -1.0 / 3.0;
31 // A'
32 at(m, m + 1) = 8.0 / 3.0;
33 at(m, m) = -3.0;
34 at(m, m - 1) = 1.0 / 3.0;
35 // Middle
36 for (u32 i = 1; i < m; i++) {
37 at(i, i) = -1.0;
38 at(i, i + 1) = 1.0;
39 }
40 // Weights
41 P = { 3.0 / 8.0 , 9.0 / 8.0 , 1.0 , 9.0 / 8.0 , 3.0 / 8.0 };
42 break;
43 case 4:
44 // A
45 at(0, 0) = -352.0 / 105.0;
46 at(0, 1) = 35.0 / 8.0;
47 at(0, 2) = -35.0 / 24.0;
48 at(0, 3) = 21.0 / 40.0;
49 at(0, 4) = -5.0 / 56.0;
50 at(1, 0) = 16.0 / 105.0;
51 at(1, 1) = -31.0 / 24.0;
52 at(1, 2) = 29.0 / 24.0;
53 at(1, 3) = -3.0 / 40.0;
54 at(1, 4) = 1.0 / 168.0;
55 // A'
56 at(m, m + 1) = 352.0 / 105.0;
57 at(m, m) = -35.0 / 8.0;
58 at(m, m - 1) = 35.0 / 24.0;
59 at(m, m - 2) = -21.0 / 40.0;
60 at(m, m - 3) = 5.0 / 56.0;
61 at(m - 1, m + 1) = -16.0 / 105.0;
62 at(m - 1, m) = 31.0 / 24.0;
63 at(m - 1, m - 1) = -29.0 / 24.0;
64 at(m - 1, m - 2) = 3.0 / 40.0;
65 at(m - 1, m - 3) = -1.0 / 168.0;
66 // Middle
67 for (u32 i = 2; i < m - 1; i++) {
68 at(i, i - 1) = 1.0 / 24.0;
69 at(i, i) = -9.0 / 8.0;
70 at(i, i + 1) = 9.0 / 8.0;
71 at(i, i + 2) = -1.0 / 24.0;
72 }
73 // Weights
74 P = { 1606.0 / 4535.0 , 941.0 / 766.0 , 1384.0 / 1541.0 , 1371.0 / 1346.0
75 , 701.0 / 700.0 , 1371.0 / 1346.0 , 1384.0 / 1541.0 , 941.0 / 766.0
76 , 1606.0 / 4535.0 };
77 break;
78 case 6:
79 // A
80 at(0, 0) = -13016.0 / 3465.0;
81 at(0, 1) = 693.0 / 128.0;
82 at(0, 2) = -385.0 / 128.0;
83 at(0, 3) = 693.0 / 320.0;
84 at(0, 4) = -495.0 / 448.0;
85 at(0, 5) = 385.0 / 1152.0;
86 at(0, 6) = -63.0 / 1408.0;
87 at(1, 0) = 496.0 / 3465.0;
88 at(1, 1) = -811.0 / 640.0;
89 at(1, 2) = 449.0 / 384.0;
90 at(1, 3) = -29.0 / 960.0;
91 at(1, 4) = -11.0 / 448.0;
92 at(1, 5) = 13.0 / 1152.0;
93 at(1, 6) = -37.0 / 21120.0;
94 at(2, 0) = -8.0 / 385.0;
95 at(2, 1) = 179.0 / 1920.0;
96 at(2, 2) = -153.0 / 128.0;
97 at(2, 3) = 381.0 / 320.0;
98 at(2, 4) = -101.0 / 1344.0;
99 at(2, 5) = 1.0 / 128.0;
100 at(2, 6) = -3.0 / 7040.0;
101 // A'
102 at(m, m + 1) = 13016.0 / 3465.0;
103 at(m, m) = -693.0 / 128.0;
104 at(m, m - 1) = 385.0 / 128.0;
105 at(m, m - 2) = -693.0 / 320.0;
106 at(m, m - 3) = 495.0 / 448.0;
107 at(m, m - 4) = -385.0 / 1152.0;
108 at(m, m - 5) = 63.0 / 1408.0;
109 at(m - 1, m + 1) = -496.0 / 3465.0;
110 at(m - 1, m) = 811.0 / 640.0;
111 at(m - 1, m - 1) = -449.0 / 384.0;
112 at(m - 1, m - 2) = 29.0 / 960.0;
113 at(m - 1, m - 3) = 11.0 / 448.0;
114 at(m - 1, m - 4) = -13.0 / 1152.0;
115 at(m - 1, m - 5) = 37.0 / 21120.0;
116 at(m - 2, m + 1) = 8.0 / 385.0;
117 at(m - 2, m) = -179.0 / 1920.0;
118 at(m - 2, m - 1) = 153.0 / 128.0;
119 at(m - 2, m - 2) = -381.0 / 320.0;
120 at(m - 2, m - 3) = 101.0 / 1344.0;
121 at(m - 2, m - 4) = -1.0 / 128.0;
122 at(m - 2, m - 5) = 3.0 / 7040.0;
123 // Middle
124 for (u32 i = 3; i < m - 2; i++) {
125 at(i, i - 2) = -3.0 / 640.0;
126 at(i, i - 1) = 25.0 / 384.0;
127 at(i, i) = -75.0 / 64.0;
128 at(i, i + 1) = 75.0 / 64.0;
129 at(i, i + 2) = -25.0 / 384.0;
130 at(i, i + 3) = 3.0 / 640.0;
131 }
132 // Weights
133 P = { 420249.0 / 1331069.0 , 2590978.0 / 1863105.0 , 882762.0 / 1402249.0
134 , 1677712.0 / 1359311.0 , 239985.0 / 261097.0 , 664189.0 / 657734.0
135 , 756049.0 / 754729.0 , 664189.0 / 657734.0 , 239985.0 / 261097.0
136 , 1677712.0 / 1359311.0 , 882762.0 / 1402249.0 , 2590978.0 / 1863105.0
137 , 420249.0 / 1331069.0 };
138 break;
139 case 8:
140 // A
141 at(0, 0) = -4856215.0 / 1200963.0;
142 at(0, 1) = 45858154.0 / 7297397.0;
143 at(0, 2) = -23409299.0 / 4789435.0;
144 at(0, 3) = 3799178.0 / 719717.0;
145 at(0, 4) = -4892189.0 / 1089890.0;
146 at(0, 5) = 1789111.0 / 658879.0;
147 at(0, 6) = -1406819.0 / 1289899.0;
148 at(0, 7) = 1154863.0 / 4436807.0;
149 at(0, 8) = -2936602.0 / 105142673.0;
150 at(1, 0) = 86048.0 / 675675.0;
151 at(1, 1) = -131093.0 / 107520.0;
152 at(1, 2) = 5503131.0 / 5166017.0;
153 at(1, 3) = 305249.0 / 2136437.0;
154 at(1, 4) = -1763845.0 / 8250973.0;
155 at(1, 5) = 1562032.0 / 10745723.0;
156 at(1, 6) = -270419.0 / 4422611.0;
157 at(1, 7) = 2983.0 / 199680.0;
158 at(1, 8) = -2621.0 / 1612800.0;
159 at(2, 0) = -3776.0 / 225225.0;
160 at(2, 1) = 8707.0 / 107520.0;
161 at(2, 2) = -17947.0 / 15360.0;
162 at(2, 3) = 29319.0 / 25600.0;
163 at(2, 4) = -533.0 / 21504.0;
164 at(2, 5) = -263.0 / 9216.0;
165 at(2, 6) = 903.0 / 56320.0;
166 at(2, 7) = -283.0 / 66560.0;
167 at(2, 8) = 257.0 / 537600.0;
168 at(3, 0) = 32.0 / 9009.0;
169 at(3, 1) = -543.0 / 35840.0;
170 at(3, 2) = 265.0 / 3072.0;
171 at(3, 3) = -1233.0 / 1024.0;
172 at(3, 4) = 8625.0 / 7168.0;
173 at(3, 5) = -775.0 / 9216.0;
174 at(3, 6) = 639.0 / 56320.0;
175 at(3, 7) = -15.0 / 13312.0;
176 at(3, 8) = 1.0 / 21504.0;
177 // A'
178 at(m, m + 1) = 4856215.0 / 1200963.0;
179 at(m, m) = -45858154.0 / 7297397.0;
180 at(m, m - 1) = 23409299.0 / 4789435.0;
181 at(m, m - 2) = -3799178.0 / 719717.0;
182 at(m, m - 3) = 4892189.0 / 1089890.0;
183 at(m, m - 4) = -1789111.0 / 658879.0;
184 at(m, m - 5) = 1406819.0 / 1289899.0;
185 at(m, m - 6) = -1154863.0 / 4436807.0;
186 at(m, m - 7) = 2936602.0 / 105142673.0;
187 at(m - 1, m + 1) = -86048.0 / 675675.0;
188 at(m - 1, m) = 131093.0 / 107520.0;
189 at(m - 1, m - 1) = -5503131.0 / 5166017.0;
190 at(m - 1, m - 2) = -305249.0 / 2136437.0;
191 at(m - 1, m - 3) = 1763845.0 / 8250973.0;
192 at(m - 1, m - 4) = -1562032.0 / 10745723.0;
193 at(m - 1, m - 5) = 270419.0 / 4422611.0;
194 at(m - 1, m - 6) = -2983.0 / 199680.0;
195 at(m - 1, m - 7) = 2621.0 / 1612800.0;
196 at(m - 2, m + 1) = 3776.0 / 225225.0;
197 at(m - 2, m) = -8707.0 / 107520.0;
198 at(m - 2, m - 1) = 17947.0 / 15360.0;
199 at(m - 2, m - 2) = -29319.0 / 25600.0;
200 at(m - 2, m - 3) = 533.0 / 21504.0;
201 at(m - 2, m - 4) = 263.0 / 9216.0;
202 at(m - 2, m - 5) = -903.0 / 56320.0;
203 at(m - 2, m - 6) = 283.0 / 66560.0;
204 at(m - 2, m - 7) = -257.0 / 537600.0;
205 at(m - 3, m + 1) = -32.0 / 9009.0;
206 at(m - 3, m) = 543.0 / 35840.0;
207 at(m - 3, m - 1) = -265.0 / 3072.0;
208 at(m - 3, m - 2) = 1233.0 / 1024.0;
209 at(m - 3, m - 3) = -8625.0 / 7168.0;
210 at(m - 3, m - 4) = 775.0 / 9216.0;
211 at(m - 3, m - 5) = -639.0 / 56320.0;
212 at(m - 3, m - 6) = 15.0 / 13312.0;
213 at(m - 3, m - 7) = -1.0 / 21504.0;
214 // Middle
215 for (u32 i = 4; i < m - 3; i++) {
216 at(i, i - 3) = 5.0 / 7168.0;
217 at(i, i - 2) = -49.0 / 5120.0;
218 at(i, i - 1) = 245.0 / 3072.0;
219 at(i, i) = -1225.0 / 1024.0;
220 at(i, i + 1) = 1225.0 / 1024.0;
221 at(i, i + 2) = -245.0 / 3072.0;
222 at(i, i + 3) = 49.0 / 5120.0;
223 at(i, i + 4) = -5.0 / 7168.0;
224 }
225 // Weights
226 P = { 267425.0 / 904736.0 , 2307435.0 / 1517812.0 , 847667.0 / 3066027.0
227 , 4050911.0 / 2301238.0 , 498943.0 / 1084999.0 , 211042.0 / 170117.0
228 , 2065895.0 / 2191686.0 , 1262499.0 / 1258052.0 , 1314891.0 / 1312727.0
229 , 1262499.0 / 1258052.0 , 2065895.0 / 2191686.0 , 211042.0 / 170117.0
230 , 498943.0 / 1084999.0 , 4050911.0 / 2301238.0 , 847667.0 / 3066027.0
231 , 2307435.0 / 1517812.0 , 267425.0 / 904736.0 };
232 break;
233 }
234
235 // Scaling
236 *this /= dx;
237}

◆ Gradient() [2/3]

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

2-D Mimetic Gradient Constructor

Parameters
kOrder of accuracy
mNumber of cells in x-direction
nNumber of cells in y-direction
dxSpacing between cells in x-direction
dySpacing between cells in y-direction

Definition at line 240 of file gradient.cpp.

240 {
241 Gradient Gx(k, m, dx);
242 Gradient Gy(k, n, dy);
243
244 sp_mat Im = speye(m + 2, m + 2);
245 sp_mat In = speye(n + 2, n + 2);
246
247 Im.shed_row(0);
248 Im.shed_row(m);
249 In.shed_row(0);
250 In.shed_row(n);
251
252 sp_mat G1 = Utils::spkron(In, Gx);
253 sp_mat G2 = Utils::spkron(Gy, Im);
254
255 // Dimensions = 2*m*n+m+n, (m+2)*(n+2)
256 if (m != n)
257 *this = Utils::spjoin_cols(G1, G2);
258 else {
259 sp_mat A1(2, 1);
260 sp_mat A2(2, 1);
261 A1(0, 0) = A2(1, 0) = 1.0;
262 *this = Utils::spkron(A1, G1) + Utils::spkron(A2, G2);
263 }
264}
Gradient::Gradient
Gradient(u16 k, u32 m, Real dx)
1-D Mimetic Gradient Constructor
Definition gradient.cpp:20
Utils::spjoin_cols
static sp_mat spjoin_cols(const sp_mat &A, const sp_mat &B)
An in place operation for joining two matrices by columns.
Definition utils.cpp:138
Utils::spkron
static sp_mat spkron(const sp_mat &A, const sp_mat &B)
A wrappper for implementing a sparse Kroenecker product.
Definition utils.cpp:70

◆ Gradient() [3/3]

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

3-D Mimetic Gradient Constructor

Parameters
kOrder of accuracy
mNumber of cells in x-direction
nNumber of cells in y-direction
oNumber of cells in z-direction
dxSpacing between cells in x-direction
dySpacing between cells in y-direction
dzSpacing between cells in z-direction

Definition at line 267 of file gradient.cpp.

267 {
268 Gradient Gx(k, m, dx);
269 Gradient Gy(k, n, dy);
270 Gradient Gz(k, o, dz);
271
272 sp_mat Im = speye(m + 2, m + 2);
273 sp_mat In = speye(n + 2, n + 2);
274 sp_mat Io = speye(o + 2, o + 2);
275
276 Im.shed_row(0);
277 Im.shed_row(m);
278 In.shed_row(0);
279 In.shed_row(n);
280 Io.shed_row(0);
281 Io.shed_row(o);
282
283 sp_mat G1 = Utils::spkron(Utils::spkron(Io, In), Gx);
284 sp_mat G2 = Utils::spkron(Utils::spkron(Io, Gy), Im);
285 sp_mat G3 = Utils::spkron(Utils::spkron(Gz, In), Im);
286
287 // Dimensions = 3*m*n*o+m*n+m*o+n*o, (m+2)*(n+2)*(o+2)
288 if ((m != n) || (n != o))
289 *this = Utils::spjoin_cols(Utils::spjoin_cols(G1, G2), G3);
290 else {
291 sp_mat A1(3, 1);
292 sp_mat A2(3, 1);
293 sp_mat A3(3, 1);
294 A1(0, 0) = A2(1, 0) = A3(2, 0) = 1.0;
295 *this =
296 Utils::spkron(A1, G1) + Utils::spkron(A2, G2) + Utils::spkron(A3, G3);
297 }
298}

Member Function Documentation

◆ getP()

vec Gradient::getP ( )

Returns the weights used in the Mimeitc Gradient Operators.

Note
for informational purposes only, can be used in constructing new operators.

Definition at line 301 of file gradient.cpp.

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