LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
hodge_and_dirac_test.cc
1
7#include <hodge_and_dirac_experiment.h>
8
9using complex = std::complex<double>;
10
23int main(int argc, char *argv[]) {
24 if (argc != 4 && argc != 5) {
25 std::cerr << "Usage: " << argv[0]
26 << " max_refinement_level min_k max_k step=0.1 " << std::endl;
27 exit(1);
28 }
29
30 double min_k;
31 double max_k;
32 double step = 0.1;
33 const unsigned refinement_level = atoi(argv[1]);
34 sscanf(argv[2], "%lf", &min_k);
35 sscanf(argv[3], "%lf", &max_k);
36 if (argc == 5) sscanf(argv[4], "%lf", &step);
37 std::cout << "max_refinement_level : " << refinement_level << std::endl;
38 std::cout << "min_k : " << min_k << ", max_k = " << max_k
39 << " step = " << step << std::endl;
40
41 // needs to be a reference such that it can be reassigned
42 double &k = min_k;
43
44 std::vector<unsigned> refinement_levels(refinement_level + 1);
45 for (int i = 0; i < refinement_level + 1; i++) {
46 refinement_levels[i] = i;
47 }
48
49 std::vector<double> ks;
50 for (double i = min_k; i <= max_k; i += step) {
51 ks.push_back(i);
52 }
53
54 // mathematica function output requries the following helpers
55 auto Power = [](double a, double b) -> double { return std::pow(a, b); };
56 auto Sin = [](double a) -> double { return std::sin(a); };
57 auto Cos = [](double a) -> double { return std::cos(a); };
58 auto Sqrt = [](double a) -> double { return std::sqrt(a); };
59 auto Complex = [](double a, double b) -> complex { return complex(a, b); };
60
61 double r = 1.;
62
63 // righthandside for the zero and two form
64 auto f_zero = [&](const Eigen::Vector3d &x_vec) -> complex {
65 // first scale to the circle
66 Eigen::Vector3d x_ = x_vec / x_vec.norm();
67 double x = x_(0);
68 double y = x_(1);
69 double z = x_(2);
70
71 // autogenerated by mathematica
72 complex ret = Complex(0, 1) * k * (Cos(y) + Sin(x) + Sin(z)) +
73 (x * z * Cos(x) + x * y * Cos(y) + y * z * Cos(z) +
74 2 * z * Sin(x) + 2 * x * Sin(y) + 2 * y * Sin(z)) /
75 (Power(x, 2) + Power(y, 2) + Power(z, 2));
76 return ret;
77 };
78
79 // righthandside for the zero and two form
80 auto f_zero_dirac = [&](const Eigen::Vector3d &x_vec) -> complex {
81 // first scale to the circle
82 Eigen::Vector3d x_ = x_vec / x_vec.norm();
83 double x = x_(0);
84 double y = x_(1);
85 double z = x_(2);
86
87 // autogenerated by mathematica
88 complex ret = (2 * x * Cos(x) +
89 ((1 + Power(k, 2)) * Power(x, 2) + Power(z, 2) +
90 Power(k, 2) * (Power(y, 2) + Power(z, 2))) *
91 Cos(y) +
92 2 * z * Cos(z) + Power(k, 2) * Power(x, 2) * Sin(x) +
93 Power(y, 2) * Sin(x) + Power(k, 2) * Power(y, 2) * Sin(x) +
94 Power(z, 2) * Sin(x) + Power(k, 2) * Power(z, 2) * Sin(x) -
95 2 * y * Sin(y) + Power(x, 2) * Sin(z) +
96 Power(k, 2) * Power(x, 2) * Sin(z) + Power(y, 2) * Sin(z) +
97 Power(k, 2) * Power(y, 2) * Sin(z) +
98 Power(k, 2) * Power(z, 2) * Sin(z)) /
99 (Power(x, 2) + Power(y, 2) + Power(z, 2));
100 return ret;
101 };
102
103 // righthandside for the one form
104 auto f_one = [&](const Eigen::Vector3d &x_vec) -> Eigen::VectorXcd {
105 // first scale to the circle
106 Eigen::Vector3d x_ = x_vec / x_vec.norm();
107 double x = x_(0);
108 double y = x_(1);
109 double z = x_(2);
110
111 // autogenerated by mathematica
112 Eigen::VectorXcd ret(3);
113 ret << ((Power(y, 2) + Power(z, 2)) * Cos(x) - x * z * Cos(z) +
114 x * y * Sin(y) +
115 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
116 (y * Cos(z) + z * Sin(y)) -
117 Complex(0, 1) * k *
118 (x * z * Sin(x) - (Power(y, 2) + Power(z, 2)) * Sin(y) +
119 x * y * Sin(z))) /
120 (Power(x, 2) + Power(y, 2) + Power(z, 2)),
121 -((x * y * Cos(x) + y * z * Cos(z) +
122 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
123 (-(z * Cos(x)) + x * Cos(z)) +
124 (Power(x, 2) + Power(z, 2)) * Sin(y) +
125 Complex(0, 1) * k *
126 (y * z * Sin(x) + x * y * Sin(y) -
127 (Power(x, 2) + Power(z, 2)) * Sin(z))) /
128 (Power(x, 2) + Power(y, 2) + Power(z, 2))),
129 (-(x * z * Cos(x)) + (Power(x, 2) + Power(y, 2)) * Cos(z) +
130 y * z * Sin(y) -
131 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
132 (y * Cos(x) + x * Sin(y)) +
133 Complex(0, 1) * k *
134 ((Power(x, 2) + Power(y, 2)) * Sin(x) -
135 z * (x * Sin(y) + y * Sin(z)))) /
136 (Power(x, 2) + Power(y, 2) + Power(z, 2));
137 return ret;
138 };
139
140 // righthandside for the one form
141 auto f_one_dirac = [&](const Eigen::Vector3d &x_vec) -> Eigen::VectorXcd {
142 // first scale to the circle
143 Eigen::Vector3d x_ = x_vec / x_vec.norm();
144 double x = x_(0);
145 double y = x_(1);
146 double z = x_(2);
147
148 // autogenerated by mathematica
149 Eigen::VectorXcd ret(3);
150 ret << ((-(Power(x, 2) * z) + 3 * z * (Power(y, 2) + Power(z, 2)) +
151 Complex(0, 1) * k * (Power(y, 2) + Power(z, 2)) *
152 (Power(x, 2) + Power(y, 2) + Power(z, 2))) *
153 Cos(x) +
154 y * (-3 * Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(y) -
155 Complex(0, 1) * k * Power(x, 3) * z * Cos(z) -
156 4 * x * y * z * Cos(z) -
157 Complex(0, 1) * k * x * Power(y, 2) * z * Cos(z) -
158 Complex(0, 1) * k * x * Power(z, 3) * Cos(z) - 2 * x * z * Sin(x) -
159 x * Power(y, 2) * z * Sin(x) - x * Power(z, 3) * Sin(x) +
160 Complex(0, 1) * k * Power(x, 3) * y * Sin(y) +
161 2 * Power(y, 2) * Sin(y) + Power(x, 2) * Power(y, 2) * Sin(y) +
162 Complex(0, 1) * k * x * Power(y, 3) * Sin(y) +
163 2 * Power(z, 2) * Sin(y) +
164 Complex(0, 1) * k * x * y * Power(z, 2) * Sin(y) -
165 2 * x * y * Sin(z) + x * y * Power(z, 2) * Sin(z)) /
166 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2) +
167 (-(z * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(x)) +
168 y * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(y) +
169 Complex(0, 1) * k * Power(x, 2) * y * Cos(z) +
170 Complex(0, 1) * k * Power(y, 3) * Cos(z) +
171 Complex(0, 1) * k * y * Power(z, 2) * Cos(z) +
172 Complex(0, 1) * k * Power(x, 2) * z * Sin(y) +
173 Complex(0, 1) * k * Power(y, 2) * z * Sin(y) +
174 Complex(0, 1) * k * Power(z, 3) * Sin(y) +
175 Power(z, 2) * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
176 Sin(y) -
177 x * y * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
178 Sin(z)) /
179 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 1.5) -
180 Complex(0, 1) * k *
181 (((Power(y, 2) + Power(z, 2)) * Cos(x) - x * z * Cos(z) +
182 x * y * Sin(y)) /
183 (Power(x, 2) + Power(y, 2) + Power(z, 2)) +
184 (y * Cos(z) + z * Sin(y)) /
185 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) +
186 (Complex(0, 1) * k *
187 (-(x * z * Sin(x)) + (Power(y, 2) + Power(z, 2)) * Sin(y) -
188 x * y * Sin(z))) /
189 (Power(x, 2) + Power(y, 2) + Power(z, 2))),
190 (Complex(0, -1) *
191 (x * y *
192 (Complex(0, -4) * z +
193 k * (Power(x, 2) + Power(y, 2) + Power(z, 2))) *
194 Cos(x) +
195 Complex(0, 1) * x *
196 (3 * Power(x, 2) - Power(y, 2) + 3 * Power(z, 2)) * Cos(y) +
197 Complex(0, 1) * Power(x, 2) * z * Cos(z) +
198 k * Power(x, 2) * y * z * Cos(z) -
199 Complex(0, 3) * Power(y, 2) * z * Cos(z) +
200 k * Power(y, 3) * z * Cos(z) + Complex(0, 1) * Power(z, 3) * Cos(z) +
201 k * y * Power(z, 3) * Cos(z) - Complex(0, 2) * y * z * Sin(x) +
202 Complex(0, 1) * Power(x, 2) * y * z * Sin(x) +
203 k * Power(x, 4) * Sin(y) - Complex(0, 2) * x * y * Sin(y) -
204 Complex(0, 1) * Power(x, 3) * y * Sin(y) +
205 k * Power(x, 2) * Power(y, 2) * Sin(y) +
206 2 * k * Power(x, 2) * Power(z, 2) * Sin(y) -
207 Complex(0, 1) * x * y * Power(z, 2) * Sin(y) +
208 k * Power(y, 2) * Power(z, 2) * Sin(y) + k * Power(z, 4) * Sin(y) +
209 Complex(0, 2) * Power(x, 2) * Sin(z) +
210 Complex(0, 2) * Power(z, 2) * Sin(z) +
211 Complex(0, 1) * Power(y, 2) * Power(z, 2) * Sin(z))) /
212 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2) +
213 (Complex(0, 1) * k * z * (Power(x, 2) + Power(y, 2) + Power(z, 2)) *
214 Cos(x) -
215 x * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(y) -
216 Complex(0, 1) * k * Power(x, 3) * Cos(z) -
217 Complex(0, 1) * k * x * Power(y, 2) * Cos(z) -
218 Complex(0, 1) * k * x * Power(z, 2) * Cos(z) +
219 z * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(z) -
220 y * z * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Sin(x) +
221 Power(x, 2) * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) *
222 Sin(z)) /
223 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 1.5) -
224 Complex(0, 1) * k *
225 ((z * Cos(x) - x * Cos(z)) /
226 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) -
227 (x * y * Cos(x) + y * z * Cos(z) +
228 (Power(x, 2) + Power(z, 2)) * Sin(y)) /
229 (Power(x, 2) + Power(y, 2) + Power(z, 2)) +
230 (Complex(0, 1) * k *
231 (-(y * z * Sin(x)) - x * y * Sin(y) +
232 (Power(x, 2) + Power(z, 2)) * Sin(z))) /
233 (Power(x, 2) + Power(y, 2) + Power(z, 2))),
234 ((x * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) -
235 Complex(0, 1) * k * y * (Power(x, 2) + Power(y, 2) + Power(z, 2))) *
236 Cos(x) -
237 y * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Cos(z) +
238 Power(y, 2) * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Sin(x) -
239 Complex(0, 1) * k * Power(x, 3) * Sin(y) -
240 Complex(0, 1) * k * x * Power(y, 2) * Sin(y) -
241 Complex(0, 1) * k * x * Power(z, 2) * Sin(y) -
242 x * z * Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) * Sin(y)) /
243 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 1.5) +
244 (x *
245 (Power(z, 2) * (-3. - Complex(0, 1) * k * z) +
246 Power(x, 2) * (1. - Complex(0, 1) * k * z) +
247 Power(y, 2) * (1. - Complex(0, 1) * k * z)) *
248 Cos(x) -
249 4 * x * y * z * Cos(y) + Complex(0, 1) * k * Power(x, 4) * Cos(z) +
250 3 * Power(x, 2) * y * Cos(z) +
251 Complex(0, 2) * k * Power(x, 2) * Power(y, 2) * Cos(z) +
252 3 * Power(y, 3) * Cos(z) +
253 Complex(0, 1) * k * Power(y, 4) * Cos(z) +
254 Complex(0, 1) * k * Power(x, 2) * Power(z, 2) * Cos(z) -
255 y * Power(z, 2) * Cos(z) +
256 Complex(0, 1) * k * Power(y, 2) * Power(z, 2) * Cos(z) +
257 2 * Power(x, 2) * Sin(x) + 2 * Power(y, 2) * Sin(x) +
258 Power(x, 2) * Power(z, 2) * Sin(x) - 2 * x * z * Sin(y) +
259 Complex(0, 1) * k * Power(x, 2) * y * z * Sin(y) +
260 x * Power(y, 2) * z * Sin(y) +
261 Complex(0, 1) * k * Power(y, 3) * z * Sin(y) +
262 Complex(0, 1) * k * y * Power(z, 3) * Sin(y) - 2 * y * z * Sin(z) -
263 Power(x, 2) * y * z * Sin(z) - Power(y, 3) * z * Sin(z)) /
264 Power(Power(x, 2) + Power(y, 2) + Power(z, 2), 2) -
265 Complex(0, 1) * k *
266 (-((y * Cos(x) + x * Sin(y)) /
267 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2))) +
268 (-(x * z * Cos(x)) + (Power(x, 2) + Power(y, 2)) * Cos(z) +
269 y * z * Sin(y)) /
270 (Power(x, 2) + Power(y, 2) + Power(z, 2)) +
271 (Complex(0, 1) * k *
272 ((Power(x, 2) + Power(y, 2)) * Sin(x) -
273 z * (x * Sin(y) + y * Sin(z)))) /
274 (Power(x, 2) + Power(y, 2) + Power(z, 2)));
275 return ret;
276 };
277
278 // righthandside for the two form
279 auto f_two = [&](const Eigen::Vector3d &x_vec) -> complex {
280 // first scale to the circle
281 Eigen::Vector3d x_ = x_vec / x_vec.norm();
282 double x = x_(0);
283 double y = x_(1);
284 double z = x_(2);
285
286 // autogenerated by mathematica
287 complex ret = (y * Cos(x) + z * Cos(y) + x * Cos(z)) /
288 Sqrt(Power(x, 2) + Power(y, 2) + Power(z, 2)) +
289 Complex(0, 1) * k * (Cos(y) + Sin(x) + Sin(z));
290 return ret;
291 };
292
293 // righthandside for the two form
294 auto f_two_dirac = [&](const Eigen::Vector3d &x_vec) -> complex {
295 // first scale to the circle
296 Eigen::Vector3d x_ = x_vec / x_vec.norm();
297 double x = x_(0);
298 double y = x_(1);
299 double z = x_(2);
300
301 // autogenerated by mathematica
302 complex ret = (2 * x * Cos(x) +
303 ((1 + Power(k, 2)) * Power(x, 2) + Power(z, 2) +
304 Power(k, 2) * (Power(y, 2) + Power(z, 2))) *
305 Cos(y) +
306 2 * z * Cos(z) + Power(k, 2) * Power(x, 2) * Sin(x) +
307 Power(y, 2) * Sin(x) + Power(k, 2) * Power(y, 2) * Sin(x) +
308 Power(z, 2) * Sin(x) + Power(k, 2) * Power(z, 2) * Sin(x) -
309 2 * y * Sin(y) + Power(x, 2) * Sin(z) +
310 Power(k, 2) * Power(x, 2) * Sin(z) + Power(y, 2) * Sin(z) +
311 Power(k, 2) * Power(y, 2) * Sin(z) +
312 Power(k, 2) * Power(z, 2) * Sin(z)) /
313 (Power(x, 2) + Power(y, 2) + Power(z, 2));
314 return ret;
315 };
316
317 projects::hldo_sphere::experiments::HodgeAndDiracExperiment experiment(
318 f_zero_dirac, f_one_dirac, f_two_dirac, f_zero, f_one, f_two, k, "test");
319
320 experiment.Compute(refinement_levels, ks);
321
322 return 0;
323}
projects::hldo_sphere::experiments::HodgeAndDiracExperiment
Creates and solves the discretised Hodge Laplacian source problems and the Dirac operator source Prob...
Definition: hodge_and_dirac_experiment.h:51
Generated by doxygen 1.9.3