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