Mesh_To_FEM_Test.cpp 42.8 KB
Newer Older
1
#include "Library_Integration_Test.hpp"
2
3
4

std::string SAVE_DIR = "./test/output/Integration/";

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
TEST(Full_Circle, Magnetostatic_Uniform_Current_Density) {
    // Create Sketch
    Sketch sk;

    auto v0 = sk.new_element<Vertex>(0.0, 0.0);
    auto v1 = sk.new_element<Vertex>(1.0, 0.0);
    auto v2 = sk.new_element<Vertex>(1.0 * std::cos(M_PI * 2.0 / 3.0), 1.0 * std::sin(M_PI * 2.0 / 3.0));
    auto v3 = sk.new_element<Vertex>(1.0 * std::cos(-M_PI * 2.0 / 3.0), 1.0 * std::sin(-M_PI * 2.0 / 3.0));

    auto c0 = sk.new_element<CircularArc>(v1, v2, v0, 1.0);
    auto c1 = sk.new_element<CircularArc>(v2, v3, v0, 1.0);
    auto c2 = sk.new_element<CircularArc>(v3, v1, v0, 1.0);

    auto f0 = sk.new_element<Fixation>(v0);
    auto f1 = sk.new_element<Fixation>(v1);
    auto f2 = sk.new_element<Fixation>(v2);
    auto f3 = sk.new_element<Fixation>(v3);

    double_t tol = sk.solve();
    EXPECT_LE(tol, FLT_EPSILON);

    bool result = sk.build();
    EXPECT_TRUE(result);

    EXPECT_EQ(sk.size_contours(), 1);

    sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, std::string("full_circle_sketch"));

    // Create Refineable Mesh
    Mesh me{sk};
    me.create();

    me.MinimumElementSize = 0.01;
    me.MaximumElementSize = 0.1;
    me.MinimumElementQuality = 0.5;

    me.refine();

43
    me.save_as(SAVE_DIR, std::string("full_circle_mesh"));
44
45

    // Convert to FiniteElementMesh
46
    auto fem = std::make_shared<FiniteElementMesh<1,1>>(me);
47

48
49
50
    for (auto iter : fem->boundaries()) {
        auto b = iter.second;

51
52
        double_t a0{-2 * M_PI};
        for (size_t i : b->nodes()) {
53
            XY const &p = fem->node(i);
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69

            // Test radius of boundary
            EXPECT_NEAR(1.0, std::hypot(p.x(), p.y()), FLT_EPSILON);

            // Test ordering of boundary nodes
            double_t a1{std::atan2(p.y(), p.x())};
            if (a1 < 0 || (a1 == 0 && a0 > M_PI)) {
                a1 += 2 * M_PI;
            }

            EXPECT_TRUE(a1 > a0);
            a0 = a1;
        }
    }

    // Create magnetostatic physics
70
    Magnetostatic msph{fem};
71

72
73
    FunctionArguments fargs;

74
75
    auto cdr = me.select_contour(Point{0.0,0.0});
    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {fem->region(cdr)});
76

77
78
79
    msph.add_magnetic_insulation(fem->boundary(c0));
    msph.add_magnetic_insulation(fem->boundary(c1));
    msph.add_magnetic_insulation(fem->boundary(c2));
80
81
82
83
84
85

    msph.build_matrices();

    msph.apply_conditions();

    // Initialize matrices
86
87
    auto solution = msph.new_solution();

88
    // Set time
89
    msph.calculate_forcing(solution->f, fargs);
90
91

    // Linearize
92
93
    solution->v.setZero();
    msph.linearize(solution);
94
95
96

    // Factor
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
97
    ldlt.compute(solution->J);
98
99
100
    ASSERT_EQ(ldlt.info(), Eigen::Success);

    // Solve
101
    solution->v -= ldlt.solve(solution->r);
102
103

    // Save
104
    Eigen::VectorXd vv = msph.recover_boundary(solution->v);
105
    fem->write_scalar(vv, SAVE_DIR, std::string("full_circle_magnetostatic_uniform_current_density"));
106
107

    // Test
108
    msph.linearize(solution);
109

110
111
    for (size_t i = 0; i != solution->r.size(); ++i) {
        EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
112
113
114
    }

    // Test flux-density values
115
116
    Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * solution->v;
    Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * solution->v;
117
118
119
120
121
122
123
124
125

    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());

    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
    }

126
    std::vector<std::vector<XY>> qp = fem->quadrature_points();
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141

    for (size_t i = 0; i != qp.size(); ++i) {
        double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());

        EXPECT_NEAR(Bmag(i), r, 1.0 * 0.05);

        double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;

        if (a > 180.0) { a -= 360.0; }
        if (a - Bang(i) > 180.0) { a -= 360.0; }
        if (a - Bang(i) < -180.0) { a += 360.0; }
        EXPECT_NEAR(Bang(i), a, 360 * 0.05);
    }
}

142
143
144
145
146
147
148
149
150
151
152
class Sixth_Circle_Mirror_Copy : public ::testing::Test {
public:
    virtual void SetUp() {
        SetUp_Sketch();
        SetUp_Mesh();
    }

    void SetUp_Sketch() {
        v0 = sk.new_element<Vertex>(0.0, 0.0);
        v1 = sk.new_element<Vertex>(1.0, 0.0);
        v2 = sk.new_element<Vertex>(2.0, 0.0);
JasonPries's avatar
JasonPries committed
153
        v3 = sk.new_element<Vertex>(sqrt(3.0) / 2.0, 0.5);
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
        v4 = sk.new_element<Vertex>(sqrt(3.0), 1.0);

        f0 = sk.new_element<Fixation>(v0);

        l01 = sk.new_element<LineSegment>(v0, v1);
        l12 = sk.new_element<LineSegment>(v1, v2);

        l04 = sk.new_element<LineSegment>(v0, v4, true);

        h01 = sk.new_element<Horizontal>(l01);
        h12 = sk.new_element<Horizontal>(l12);

        a0103 = sk.new_element<Angle>(l01, l04, 30.0);

        c013 = sk.new_element<CircularArc>(v1, v3, v0, 1.0);
        c024 = sk.new_element<CircularArc>(v2, v4, v0, 2.0);

        r013 = sk.new_element<Radius>(c013, 1.0);
        r024 = sk.new_element<Radius>(c024, 2.0);

        std::vector<std::shared_ptr<Curve const>> vec{l01, l12, c013, c024};
        m04 = sk.new_element<MirrorCopy>(vec, l04);

        double residual = sk.solve();
        EXPECT_LE(residual, FLT_EPSILON);

        bool result = sk.build();
        ASSERT_TRUE(result);

        sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, "sixth_circle_mirror_copy_sketch");

        periodic_boundary = sk.select_periodic_boundary_pairs(v0, 60.0);
        {
            for (auto &bp : periodic_boundary) {
                if (bp.curve0().get() == l01.get()) {
                    EXPECT_EQ(bp.curve1()->is_identical(l01, v0, 60.0), MatchOrientation::Reverse);
                } else if (bp.curve0().get() == l12.get()) {
                    EXPECT_EQ(bp.curve1()->is_identical(l12, v0, 60.0), MatchOrientation::Reverse);
                } else {
                    GTEST_NONFATAL_FAILURE_("No matching boundary found");
                }
            }
        }

        outer_boundary = sk.select_radial_boundary(v0, 2.0);
        {
            for (auto &c : outer_boundary) {
                auto cc = std::dynamic_pointer_cast<CircularArc const>(c);
                if (cc) {
                    EXPECT_NEAR(cc->radius(), 2.0, FLT_EPSILON);
                    EXPECT_EQ(cc->center().get(), v0.get());
                } else {
                    GTEST_NONFATAL_FAILURE_("dynamic_cast of Curve to CircularArc failed");
                }
            }
        }

        v1p = sk.select_periodic_vertex(v1, v0, 60.0);
        {
            ASSERT_TRUE(v1p); // Make sure that returned vertex pointer is not null
            double_t x0 = v0->x();
            double_t y0 = v0->y();
            double_t x1 = v1->x();
            double_t y1 = v1->y();
            double_t x2 = v1p->x();
            double_t y2 = v1p->y();

            double_t dx = x1 - x0;
            double_t dy = y1 - y0;
            double_t dr = std::hypot(dx, dy);
            double_t da = std::atan2(dy, dx);

            da += M_PI * 1.0 / 3.0;
            x1 = x0 + dr * std::cos(da);
            y1 = y0 + dr * std::sin(da);

            EXPECT_NEAR(x1, x2, dr * FLT_EPSILON);
            EXPECT_NEAR(y1, y2, dr * FLT_EPSILON);
        }
    }

    void SetUp_Mesh() {
        // Create Refineable Mesh
        me = Mesh(sk);
238
239
240

        me.add_mapped_boundary_pair(periodic_boundary);

241
242
243
244
245
246
247
248
        me.create();

        me.MaximumElementSize = 0.1;
        me.MinimumElementSize = 0.01;
        me.MinimumElementQuality = 0.5;

        for (auto &pb : periodic_boundary) {
            auto bc0 = me.boundary_constraint(pb.curve0());
249
            EXPECT_TRUE(bc0);
250
251

            auto bc1 = me.boundary_constraint(pb.curve1());
252
            EXPECT_TRUE(bc1);
253
254
255
256
257
258
259
260
        }

        bool result = me.refine();
        ASSERT_TRUE(result);

        me.save_as(SAVE_DIR, std::string("sixth_circle_mirror_copy_mesh"));

        // Create FiniteElementMesh
261
        fem = std::make_shared<FiniteElementMesh<1, 1>>(me);
262

JasonPries's avatar
JasonPries committed
263
        for (auto &pb : periodic_boundary) {
264
265
            EXPECT_EQ(fem->boundary(pb.curve0()).size(), 1);
            EXPECT_EQ(fem->boundary(pb.curve1()).size(), 1);
266

267
268
            auto nodes0 = fem->boundary(pb.curve0())[0]->nodes();
            auto nodes1 = fem->boundary(pb.curve1())[0]->nodes();
269

270
271
            EXPECT_EQ(nodes0.size(), nodes1.size());

272
            for (size_t i = 0; i != nodes0.size(); ++i) {
273
274
                auto const &n0 = fem->node(nodes0[i]);
                auto const &n1 = fem->node(nodes1[nodes0.size() - 1 - i]);
275
276
277
278
279
280
281
282
283
284
285
286
287

                EXPECT_NEAR(atan2(n0.y(), n0.x()), 0.0, FLT_EPSILON);
                if (&n0 != &n1) {
                    EXPECT_NEAR(atan2(n1.y(), n1.x()), M_PI / 3.0, FLT_EPSILON);
                }
            }
        }

        current_density_contour = me.select_contour(Point{0.5 * cos(M_PI / 6.0), 0.5 * sin(M_PI / 6.0)});
    }

    Sketch sk;
    Mesh me;
288
289

    std::shared_ptr<FiniteElementMeshInterface> fem;
290
291
292
293
294
295
296
297
298
299

    // Sketch Variables
    std::shared_ptr<Vertex> v0;
    std::shared_ptr<Vertex> v1;
    std::shared_ptr<Vertex> v2;
    std::shared_ptr<Vertex> v3;
    std::shared_ptr<Vertex> v4;

    std::shared_ptr<Fixation> f0;

JasonPries's avatar
JasonPries committed
300
    std::shared_ptr<LineSegment> l01;
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
    std::shared_ptr<LineSegment> l12;

    std::shared_ptr<LineSegment> l04;

    std::shared_ptr<Horizontal> h01;
    std::shared_ptr<Horizontal> h12;

    std::shared_ptr<Angle> a0103;

    std::shared_ptr<CircularArc> c013;
    std::shared_ptr<CircularArc> c024;

    std::shared_ptr<Radius> r013;
    std::shared_ptr<Radius> r024;

    std::shared_ptr<MirrorCopy> m04;

    std::shared_ptr<Contour const> current_density_contour;

    decltype(sk.select_periodic_boundary_pairs(v0, 90.0)) periodic_boundary;
    decltype(sk.select_radial_boundary(v0, 2.0)) outer_boundary;
    decltype(sk.select_periodic_vertex(v1, v0, 90.0)) v1p;
};

TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Periodic) {
    // Create Physics
327
    Magnetostatic msph{fem};
328

329
330
331
    FunctionArguments fargs;

    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
332
333
334
335
336
337
338
339
340

    msph.add_magnetic_insulation(outer_boundary);

    msph.add_periodic_boundary(periodic_boundary, false);

    msph.build_matrices();

    msph.apply_conditions();

341
    auto solution = msph.new_solution();
342

343
    solution->v.setZero();
344

345
    msph.calculate_forcing(solution->f, fargs);
346

347
    msph.linearize(solution);
348
349

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
350
    ldlt.compute(solution->J);
351
    ASSERT_EQ(ldlt.info(), Eigen::Success);
352
    solution->v -= ldlt.solve(solution->r);
353
354

    // Verify equation is solved
355
356
357
    msph.linearize(solution);
    for (size_t i = 0; i != solution->r.size(); ++i) {
        EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
358
359
360
361
362
    }

    // Test solution:
    // A = a_f * r^2 + c_f in forced inner region
    // A = a_h * log(r) + c_h in homogeneous outer region
363
    Eigen::VectorXd vv = msph.recover_boundary(solution->v);
364
    fem->write_scalar(vv, SAVE_DIR, std::string("sixth_circle_uniform_current_density_periodic"));
365
366
367
368
369

    double_t a_f = -0.5;
    double_t c_f = -2.0 * a_f * std::log(2.0) - a_f;
    double_t a_h = 2 * a_f;
    double_t c_h = -2.0 * a_f * log(2.0);
370
371
    for (size_t i = 0; i != fem->nodes_size(); ++i) {
        XY const &n = fem->node(i);
372
373
374
375
376
377
378
379
380
381
382
383
384
385
        double_t r = std::hypot(n.x(), n.y());
        if (r <= 1.0) {
            double_t val = a_f * r * r + c_f;
            EXPECT_NEAR(val, vv(i), c_f * 0.01);
        } else {
            double_t val = a_h * log(r) + c_h;
            EXPECT_NEAR(val, vv(i), c_h * 0.01);
        }
    }

    // Test flux-density values
    // B_r = 2 * a_f * r in forced inner region
    // B_r = a_h / r in homogenous outer region
    // TODO: The correct way to test may be to calculate average flux in the element from analytical solution
386
387
    Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * solution->v;
    Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * solution->v;
388
389
390
391
392
393
394
395
396

    Eigen::ArrayXd Bmag(By.size());
    Eigen::ArrayXd Bang(By.size());

    for (size_t i = 0; i != By.size(); ++i) {
        Bmag(i) = hypot(Bx(i), By(i));
        Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
    }

397
    std::vector<std::vector<XY>> qp = fem->quadrature_points();
398
399
400
401
402
403
404
405
406

    for (size_t i = 0; i != qp.size(); ++i) {
        double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());

        if (r <= 1) {
            double_t val = std::abs(2.0 * a_f * r);
            EXPECT_NEAR(Bmag(i), val, std::abs(2 * a_f * 0.05));
        } else {
            double_t val = std::abs(a_h / r);
JasonPries's avatar
JasonPries committed
407
            EXPECT_NEAR(Bmag(i), val, std::abs(a_h * 0.05));
408
409
410
411
412
413
414
415
416
417
418
        }

        double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;

        if (Bang(i) < 0.0) { Bang(i) += 360.0; };
        EXPECT_NEAR(Bang(i), a, 90 * 0.05);
    }
}

TEST_F(Sixth_Circle_Mirror_Copy, Magnetostatic_Uniform_Current_Density_Antieriodic) {
    // Create Physics
419
    Magnetostatic msph{fem};
420

421
422
423
    FunctionArguments fargs;

    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
424
425
426
427
428
429
430
431
432

    msph.add_magnetic_insulation(outer_boundary);

    msph.add_periodic_boundary(periodic_boundary, true);

    msph.build_matrices();

    msph.apply_conditions();

433
    auto solution = msph.new_solution();
434

435
    solution->v.setZero();
436

437
    msph.calculate_forcing(solution->f, fargs);
438

439
    msph.linearize(solution);
440
441

    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
442
    ldlt.compute(solution->J);
443
    ASSERT_EQ(ldlt.info(), Eigen::Success);
444
    solution->v -= ldlt.solve(solution->r);
445
446

    // Verify equation is solved
447
448
449
    msph.linearize(solution);
    for (size_t i = 0; i != solution->r.size(); ++i) {
        EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
450
451
452
453
454
455
456
    }

    // Test solution:
    // The actual analytic solution to this problem is somewhat complicated
    // The test is based on a single harmonic approximation
    // Bounds are obtained from the Fourier series approximation of the square wave current density

457
    Eigen::VectorXd vv = msph.recover_boundary(solution->v);
458
    fem->write_scalar(vv, SAVE_DIR, std::string("sixth_circle_uniform_current_density_antiperiodic"));
459
460
461
462
463
464

    double_t a_f = -2.0 / (5.0);
    double_t b_f = -384.0 / 321.0 * a_f;
    double_t c_h = a_f / 321.0;
    double_t d_h = -64.0 / 321.0 * a_f;

JasonPries's avatar
JasonPries committed
465
    double_t r_max = -(2.0 * b_f) / (3.0 * a_f);
466
    double_t tol_a = (a_f * r_max + b_f) * r_max * r_max * 0.02;
467
468
    for (size_t i = 0; i != fem->nodes_size(); ++i) {
        XY const &n = fem->node(i);
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
        double_t r = std::hypot(n.x(), n.y());
        double_t a = std::atan2(n.y(), n.x());
        if (r <= 1.0) {
            double_t val = ((a_f * r + b_f) * r * r) * sin(3.0 * a);
            EXPECT_LT(val - tol_a, vv(i));
            val = val * 4.0 / M_PI;
            EXPECT_GT(val + tol_a, vv(i));
        } else {
            double_t val = (c_h * (r * r * r) + d_h / (r * r * r)) * sin(3.0 * a);
            EXPECT_LT(val - tol_a, vv(i));
            val = val * 4.0 / M_PI;
            EXPECT_GT(val + tol_a, vv(i));
        }
    }
}

class Sixth_Circle : public ::testing::Test {
public:
    virtual void SetUp() {
        SetUp_Sketch();
        SetUp_Mesh();
    }

    void SetUp_Sketch() {
        v0 = sk.new_element<Vertex>(0.0, 0.0);
        v1 = sk.new_element<Vertex>(1.0, 0.0);
        v2 = sk.new_element<Vertex>(2.0, 0.0);
        v3 = sk.new_element<Vertex>(1.0 * cos(M_PI / 3.0), 1.0 * sin(M_PI / 3.0));
        v4 = sk.new_element<Vertex>(2.0 * cos(M_PI / 3.0), 2.0 * sin(M_PI / 3.0));

        f0 = sk.new_element<Fixation>(v0);

        l01 = sk.new_element<LineSegment>(v0, v1);
        l12 = sk.new_element<LineSegment>(v1, v2);

        l03 = sk.new_element<LineSegment>(v0, v3);
        l34 = sk.new_element<LineSegment>(v3, v4);

        h01 = sk.new_element<Horizontal>(l01);
        h12 = sk.new_element<Horizontal>(l12);

        a0103 = sk.new_element<Angle>(l01, l03, 60.0);
        a0134 = sk.new_element<Angle>(l12, l34, 60.0);

        c013 = sk.new_element<CircularArc>(v1, v3, v0, 1.0);
        c024 = sk.new_element<CircularArc>(v2, v4, v0, 2.0);

        r013 = sk.new_element<Radius>(c013, 1.0);
        r024 = sk.new_element<Radius>(c024, 2.0);

        double residual = sk.solve();
        EXPECT_LE(residual, FLT_EPSILON);

        bool result = sk.build();
        ASSERT_TRUE(result);

        sk.save_as<SaveMethod::Rasterize>(SAVE_DIR, "sixth_circle_sketch");

        periodic_boundary = sk.select_periodic_boundary_pairs(v0, 60.0);
        {
            for (auto &bp : periodic_boundary) {
                if (bp.curve0().get() == l01.get()) {
                    EXPECT_EQ(bp.curve1()->is_identical(l01, v0, 60.0), MatchOrientation::Forward);
                } else if (bp.curve0().get() == l12.get()) {
                    EXPECT_EQ(bp.curve1()->is_identical(l12, v0, 60.0), MatchOrientation::Forward);
                } else {
                    GTEST_NONFATAL_FAILURE_("No matching boundary found");
                }
            }
        }

        outer_boundary = sk.select_radial_boundary(v0, 2.0);
        {
            for (auto &c : outer_boundary) {
                auto cc = std::dynamic_pointer_cast<CircularArc const>(c);
                if (cc) {
                    EXPECT_NEAR(cc->radius(), 2.0, FLT_EPSILON);
                    EXPECT_EQ(cc->center().get(), v0.get());
                } else {
                    GTEST_NONFATAL_FAILURE_("dynamic_cast of Curve to CircularArc failed");
                }
            }
        }

        v1p = sk.select_periodic_vertex(v1, v0, 60.0);
        {
            ASSERT_TRUE(v1p); // Make sure that returned vertex pointer is not null
            double_t x0 = v0->x();
            double_t y0 = v0->y();
            double_t x1 = v1->x();
            double_t y1 = v1->y();
            double_t x2 = v1p->x();
            double_t y2 = v1p->y();

            double_t dx = x1 - x0;
            double_t dy = y1 - y0;
            double_t dr = std::hypot(dx, dy);
            double_t da = std::atan2(dy, dx);

            da += M_PI * 1.0 / 3.0;
            x1 = x0 + dr * std::cos(da);
            y1 = y0 + dr * std::sin(da);

            EXPECT_NEAR(x1, x2, dr * FLT_EPSILON);
            EXPECT_NEAR(y1, y2, dr * FLT_EPSILON);
        }
    }

    void SetUp_Mesh() {
        // Create Refineable Mesh
        me = Mesh(sk);
580
581
582

        me.add_mapped_boundary_pair(periodic_boundary);

583
584
585
586
587
588
589
590
        me.create();

        me.MaximumElementSize = 0.1;
        me.MinimumElementSize = 0.01;
        me.MinimumElementQuality = 0.5;

        for (auto &pb : periodic_boundary) {
            auto bc0 = me.boundary_constraint(pb.curve0());
591
            EXPECT_TRUE(bc0);
592
593

            auto bc1 = me.boundary_constraint(pb.curve1());
594
            EXPECT_TRUE(bc1);
595
596
597
598
599
600
601
602
        }

        bool result = me.refine();
        ASSERT_TRUE(result);

        me.save_as(SAVE_DIR, std::string("sixth_circle_mesh"));

        // Create FiniteElementMesh
603
        fem = std::make_shared<FiniteElementMesh<1, 1>>(me);
604

JasonPries's avatar
JasonPries committed
605
        for (auto &pb : periodic_boundary) {
606
607
            EXPECT_EQ(fem->boundary(pb.curve0()).size(), 1);
            EXPECT_EQ(fem->boundary(pb.curve1()).size(), 1);
608

609
610
            auto nodes0 = fem->boundary(pb.curve0())[0]->nodes();
            auto nodes1 = fem->boundary(pb.curve1())[0]->nodes();
611
612

            for (size_t i = 0; i != nodes0.size(); ++i) {
613
614
                auto const &n0 = fem->node(nodes0[i]);
                auto const &n1 = fem->node(nodes1[i]);
615
616
617
618
619
620
621
622
623
624
625
626
627

                EXPECT_NEAR(atan2(n0.y(), n0.x()), 0.0, FLT_EPSILON);
                if (&n0 != &n1) {
                    EXPECT_NEAR(atan2(n1.y(), n1.x()), M_PI / 3.0, FLT_EPSILON);
                }
            }
        }

        current_density_contour = me.select_contour(Point{0.5 * cos(M_PI / 6.0), 0.5 * sin(M_PI / 6.0)});
    }

    Sketch sk;
    Mesh me;
628
    std::shared_ptr<FiniteElementMesh<1,1>> fem;
629
630
631
632
633
634
635
636
637
638

    // Sketch Variables
    std::shared_ptr<Vertex> v0;
    std::shared_ptr<Vertex> v1;
    std::shared_ptr<Vertex> v2;
    std::shared_ptr<Vertex> v3;
    std::shared_ptr<Vertex> v4;

    std::shared_ptr<Fixation> f0;

JasonPries's avatar
JasonPries committed
639
    std::shared_ptr<LineSegment> l01;
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
    std::shared_ptr<LineSegment> l12;

    std::shared_ptr<LineSegment> l03;
    std::shared_ptr<LineSegment> l34;

    std::shared_ptr<Horizontal> h01;
    std::shared_ptr<Horizontal> h12;

    std::shared_ptr<Angle> a0103;
    std::shared_ptr<Angle> a0134;

    std::shared_ptr<CircularArc> c013;
    std::shared_ptr<CircularArc> c024;

    std::shared_ptr<Radius> r013;
    std::shared_ptr<Radius> r024;

    std::shared_ptr<Contour const> current_density_contour;

    decltype(sk.select_periodic_boundary_pairs(v0, 90.0)) periodic_boundary;
    decltype(sk.select_radial_boundary(v0, 2.0)) outer_boundary;
    decltype(sk.select_periodic_vertex(v1, v0, 90.0)) v1p;
};

TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Periodic_Rotation) {
    // Create Physics
666
    Magnetostatic msph{fem};
667
668
669
670

    auto position = msph.add_sliding_interface(c013, false);

    { // Test node renumbering of triangles
671
672
        auto b_vec = fem->boundary(c013);
        auto nodes = fem->nodes();
673

674
        for (auto t : fem->triangles()) {
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
            double_t xc = (nodes[t.node(0)].x() + nodes[t.node(1)].x() + nodes[t.node(2)].x()) / 3.0;
            double_t yc = (nodes[t.node(0)].y() + nodes[t.node(1)].y() + nodes[t.node(2)].y()) / 3.0;
            double_t rc = std::hypot(xc, yc);
            if (rc > 1.0) {
                for (size_t i = 0; i != 3; ++i) {
                    for (size_t j : b_vec[0]->nodes()) {
                        EXPECT_NE(t.node(i), j);
                        if (t.node(i) == j) {
                            std::cout << nodes[t.node(i)].x() << ',' << nodes[t.node(i)].y() << std::endl;
                        }
                    }
                }
            }
        }
    }

691
692
693
    FunctionArguments fargs;

    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
694
695
696
697
698
699
700
701
702
703

    msph.add_magnetic_insulation(outer_boundary);

    msph.add_periodic_boundary(periodic_boundary, false);

    for (size_t iter = 0; iter != 17; ++iter) {
        msph.build_matrices();

        msph.apply_conditions();

704
        auto solution = msph.new_solution();
705

706
        solution->v.setZero();
707

708
        msph.calculate_forcing(solution->f, fargs);
709

710
        msph.linearize(solution);
711
712

        Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
713
        ldlt.compute(solution->J);
714
        ASSERT_EQ(ldlt.info(), Eigen::Success);
715
        solution->v -= ldlt.solve(solution->r);
716
717

        // Verify equation is solved
718
719
720
        msph.linearize(solution);
        for (size_t i = 0; i != solution->r.size(); ++i) {
            EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
721
722
723
724
725
        }

        // Test solution:
        // A = a_f * r^2 + c_f in forced inner region
        // A = a_h * log(r) + c_h in homogeneous outer region
726
        Eigen::VectorXd vv = msph.recover_boundary(solution->v);
727
        fem->write_scalar(vv, SAVE_DIR, std::string("one_sixth_circle_uniform_current_density_periodic_rotation_") + std::to_string(iter));
728
729
730
731

        double_t a_f = -0.5;
        double_t c_f = -2.0 * a_f * std::log(2.0) - a_f;
        double_t a_h = 2 * a_f;
732
        double_t c_h = -2.0 * a_f * std::log(2.0);
733
734
        for (size_t i = 0; i != fem->nodes_size(); ++i) {
            XY const &n = fem->node(i);
735
736
737
738
739
740
741
742
743
744
745
746
747
748
            double_t r = std::hypot(n.x(), n.y());
            if (r <= 1.0) {
                double_t val = a_f * r * r + c_f;
                EXPECT_NEAR(val, vv(i), c_f * 0.01);
            } else {
                double_t val = a_h * log(r) + c_h;
                EXPECT_NEAR(val, vv(i), c_h * 0.01);
            }
        }

        // Test flux-density values
        // B_r = 2 * a_f * r in forced inner region
        // B_r = a_h / r in homogenous outer region
        // TODO: The correct way to test may be to calculate average flux in the element from analytical solution
749
750
        Eigen::ArrayXd Bx = msph.derivatives().dy(0).transpose() * solution->v;
        Eigen::ArrayXd By = -msph.derivatives().dx(0).transpose() * solution->v;
751
752
753
754
755
756
757
758
759

        Eigen::ArrayXd Bmag(By.size());
        Eigen::ArrayXd Bang(By.size());

        for (size_t i = 0; i != By.size(); ++i) {
            Bmag(i) = hypot(Bx(i), By(i));
            Bang(i) = atan2(By(i), Bx(i)) * 180.0 / M_PI;
        }

760
        std::vector<std::vector<XY>> qp = fem->quadrature_points();
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784

        for (size_t i = 0; i != qp.size(); ++i) {
            double_t r = std::hypot(qp[i][0].x(), qp[i][0].y());

            if (r <= 1) {
                double_t val = std::abs(2.0 * a_f * r);
                EXPECT_NEAR(Bmag(i), val, std::abs(2 * a_f * 0.05));
            } else {
                double_t val = std::abs(a_h / r);
                EXPECT_NEAR(Bmag(i), val, std::abs(a_h * 0.05));
            }

            double_t a = std::atan2(qp[i][0].y(), qp[i][0].x()) * 180 / M_PI + 90.0;

            if (Bang(i) < 0.0) { Bang(i) += 360.0; };
            EXPECT_NEAR(Bang(i), a, 90 * 0.05);
        }

        ++(*position);
    }
}

TEST_F(Sixth_Circle, Magnetostatic_Uniform_Current_Density_Antiperiodic_Rotation) {
    // Create Physics
785
    Magnetostatic msph{fem};
786
787
788
789
790

    auto position = msph.add_sliding_interface(c013, true);
    double_t angle = 0.0;
    double_t delta_angle = (2 * M_PI) / (6.0 * position->size());

791
792
793
    FunctionArguments fargs;

    msph.add_current_density([](FunctionArguments args) { return 1.0 / (2.0 * M_PI * 1e-7); }, {current_density_contour});
794
795
796
797
798
799
800
801
802
803

    msph.add_magnetic_insulation(outer_boundary);

    msph.add_periodic_boundary(periodic_boundary, true);

    for (size_t iter = 0; iter != 17; ++iter) {
        msph.build_matrices();

        msph.apply_conditions();

804
        auto solution = msph.new_solution();
805

806
        solution->v.setZero();
807

808
        msph.calculate_forcing(solution->f, fargs);
809

810
        msph.linearize(solution);
811
812

        Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
813
        ldlt.compute(solution->J);
814
        ASSERT_EQ(ldlt.info(), Eigen::Success);
815
        solution->v -= ldlt.solve(solution->r);
816
817

        // Verify equation is solved
818
819
820
        msph.linearize(solution);
        for (size_t i = 0; i != solution->r.size(); ++i) {
            EXPECT_NEAR(solution->r(i), 0.0, FLT_EPSILON);
821
822
823
824
825
826
827
        }

        // Test solution:
        // The actual analytic solution to this problem is somewhat complicated
        // The test is based on a single harmonic approximation
        // Bounds are obtained from the Fourier series approximation of the square wave current density

828
        Eigen::VectorXd vv = msph.recover_boundary(solution->v);
829
        fem->write_scalar(vv, SAVE_DIR, std::string("sixth_circle_uniform_current_density_antiperiodic_") + std::to_string(iter));
830
831
832
833
834
835
836
837

        double_t a_f = -2.0 / (5.0);
        double_t b_f = -384.0 / 321.0 * a_f;
        double_t c_h = a_f / 321.0;
        double_t d_h = -64.0 / 321.0 * a_f;

        double_t r_max = -(2.0 * b_f) / (3.0 * a_f);
        double_t tol_a = (a_f * r_max + b_f) * r_max * r_max * 0.02;
838
        for (size_t i = 0; i != fem->nodes_size(); ++i) {
839
            XY const &n = fem->node(i);
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
            double_t r = std::hypot(n.x(), n.y());
            double_t a = std::atan2(n.y(), n.x());
            if (r < 1.0 - FLT_EPSILON) {
                double_t val = ((a_f * r + b_f) * r * r) * sin(3.0 * a);
                EXPECT_LT(val - tol_a, vv(i));
                val = val * 4.0 / M_PI;
                EXPECT_GT(val + tol_a, vv(i));
            } else if (r > 1.0 + FLT_EPSILON) {
                double_t val = (c_h * (r * r * r) + d_h / (r * r * r)) * sin(3.0 * (a + angle));
                if (val > 0) {
                    EXPECT_LT(val - tol_a, vv(i));
                    val = val * 4.0 / M_PI;
                    EXPECT_GT(val + tol_a, vv(i));
                } else {
                    EXPECT_GT(val + tol_a, vv(i));
                    val = val * 4.0 / M_PI;
                    EXPECT_LT(val - tol_a, vv(i));
                }
            } else {

            }
        }

        // Rotate
        ++(*position);
        angle -= delta_angle;
    }
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
}

class Salient_Pole_Synchrel : public ::testing::Test {
public:
    virtual void SetUp() {
        SetUp_Sketch();
        SetUp_Mesh();
    }

    void SetUp_Sketch() {
        sk = Sketch();

        // General Dimensions
        poles = 4.0;
        airgap = 30e-3 * 25.4e-3;

        origin = sk.new_element<Vertex>(0.0, 0.0);
        auto fo = sk.new_element<Fixation>(origin);

        // Rotor
        rro = 80e-3;

        auto vr0 = sk.new_element<Vertex>(rro * cos(2.0 * M_PI / poles / 4.0), 0.0);
        auto vr1 = sk.new_element<Vertex>(rro, 0.0);
        auto vr2 = sk.new_element<Vertex>(rro * cos(2.0 * M_PI / poles / 4.0), rro * sin(2.0 * M_PI / poles / 4.0));
        auto vr3 = sk.new_element<Vertex>(rro * cos(2.0 * M_PI / poles / 2.0), rro * sin(2.0 * M_PI / poles / 2.0));

        auto lr0 = sk.new_element<LineSegment>(origin, vr0);
        auto lr1 = sk.new_element<LineSegment>(vr0, vr1);
        auto lr2 = sk.new_element<LineSegment>(vr0, vr2);
        auto lr3 = sk.new_element<LineSegment>(origin, vr3, true);

        auto dr0 = sk.new_element<Distance<LineSegment>>(lr3, lr2, 0.5 * rro * sin(M_PI / poles));

        auto hr0 = sk.new_element<Horizontal>(lr0);
        auto hr1 = sk.new_element<Horizontal>(lr1);

        auto ar0 = sk.new_element<Angle>(lr0, lr3, 90.0 / 2.0);

        auto cr0 = sk.new_element<CircularArc>(vr1, vr2, origin, rro);
        auto cr1 = sk.new_element<CircularArc>(vr2, vr3, origin, rro);

        auto rr0 = sk.new_element<Radius>(cr0, rro);
        auto rr1 = sk.new_element<Radius>(cr1, rro);

        std::vector<std::shared_ptr<Curve const>> mirror_vec_r{lr0, lr1, lr2, cr0, cr1};
        auto mr0 = sk.new_element<MirrorCopy>(mirror_vec_r, lr3);

        // Rotor part of airgap
916
        /*
917
918
919
920
921
922
923
924
925
926
927
        auto vra0 = sk.new_element<Vertex>(rro + airgap / 2.0, 0.0);
        auto vra1 = sk.new_element<Vertex>(0.0, rro + airgap / 2.0);

        auto fra0 = sk.new_element<Fixation>(vra0);
        auto fra1 = sk.new_element<Fixation>(vra1);

        auto vr1p = sk.select_periodic_vertex(vr1, origin, 90.0);
        auto lra0 = sk.new_element<LineSegment>(vr1, vra0);
        auto lra1 = sk.new_element<LineSegment>(vr1p, vra1);

        continuous_airgap = sk.new_element<CircularArc>(vra0, vra1, origin, rro + airgap / 2.0);
928
        */
929
930

        // Stator
931
932
        rsi = rro + airgap;
        rso = 132e-3;
JasonPries's avatar
JasonPries committed
933
934
        rsb = (1.0 * rso + 1.0 * rsi) / 2.0;
        tooth_width = rsi * (2.0 * M_PI / (4.0 * 3.0 * 2.0)) * (2.0 / 3.0);
935
936

        auto vs0 = sk.new_element<Vertex>(rsi, 0.0);
937
938
939
940
941
942
943
944
945
946
947
948
949
        auto vs1 = sk.new_element<Vertex>(rso, 0.0);
        auto vs2 = sk.new_element<Vertex>(rso * cos(2.0 * M_PI / 24.0), rso * sin(2.0 * M_PI / 24.0));
        auto vs3 = sk.new_element<Vertex>(rsb * cos(2.0 * M_PI / 24.0), rsb * sin(2.0 * M_PI / 24.0));
        auto vs4 = sk.new_element<Vertex>(rsb * cos(2.0 * M_PI / 24.0), rsi * sin(2.0 * M_PI / 48.0));
        auto vs5 = sk.new_element<Vertex>(rsi * cos(2.0 * M_PI / 48.0), rsi * sin(2.0 * M_PI / 48.0));
        auto vs6 = sk.new_element<Vertex>(rsi * cos(2.0 * M_PI / 24.0), rsi * sin(2.0 * M_PI / 24.0));

        auto ls01 = sk.new_element<LineSegment>(vs0, vs1);
        auto ls63 = sk.new_element<LineSegment>(vs6, vs3, true);
        auto ls32 = sk.new_element<LineSegment>(vs3, vs2, true);
        auto ls54 = sk.new_element<LineSegment>(vs5, vs4);

        auto h01 = sk.new_element<Horizontal>(ls01);
JasonPries's avatar
JasonPries committed
950
951
952
953

        auto d01_54 = sk.new_element<Distance<LineSegment>>(ls01, ls54, tooth_width);
        //auto h54 = sk.new_element<Horizontal>(ls54);
        //std::cout << "//TODO: Fix distance between l54 and l01. l54 will not need the horizontal constraint" << std::endl;
954
955
956
957
958
959
960
961
962
963
964
965
966

        auto a63 = sk.new_element<Angle>(ls01, ls63, 360.0 / 24.0);
        auto a32 = sk.new_element<Angle>(ls01, ls32, 360.0 / 24.0);

        auto co12 = sk.new_element<CircularArc>(vs1, vs2, origin, rso);
        auto co43 = sk.new_element<CircularArc>(vs4, vs3, origin, rsb);
        auto co05 = sk.new_element<CircularArc>(vs0, vs5, origin, rsi);
        auto co56 = sk.new_element<CircularArc>(vs5, vs6, origin, rsi);

        auto r12 = sk.new_element<Radius>(co12, rso);
        auto r43 = sk.new_element<Radius>(co43, rsb);
        auto r05 = sk.new_element<Radius>(co05, rsi);
        auto r56 = sk.new_element<Radius>(co56, rsi);
967

968
969
        double_t residual = sk.solve();
        EXPECT_LE(residual, FLT_EPSILON * rro);
970

971
972
        std::vector<std::shared_ptr<Curve const>> copy_items{ls01, ls54, co12, co43, co05, co56};
        auto mirror0 = sk.new_element<MirrorCopy>(copy_items, ls63, true);
973

974
975
976
977
978
979
980
981
982
983
        residual = sk.solve();
        EXPECT_LE(residual, FLT_EPSILON * rro);

        copy_items.insert(copy_items.end(), mirror0->curves().begin(), mirror0->curves().end());

        auto rotation0 = sk.new_element<RotateCopy>(copy_items, origin, 360.0 / 12.0, 2.0, true);

        residual = sk.solve();
        EXPECT_LE(residual, FLT_EPSILON * rro);

984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
        // Rotor part of airgap
        std::cout << vr1->x() << "," << vr1->y() << std::endl;
        std::cout << vs0->x() << "," << vs0->y() << std::endl;

        CylindricalAirgap<1> cyag{sk,origin,vr1,vs0,90.0};

        //auto vra0 = sk.new_element<Vertex>(rro + airgap / 2.0, 0.0);
        //auto vra1 = sk.new_element<Vertex>(0.0, rro + airgap / 2.0);

        auto fra0 = sk.new_element<Fixation>(cyag.vertex0(0));
        auto fra1 = sk.new_element<Fixation>(cyag.vertex1(0));

        //auto vr1p = sk.select_periodic_vertex(vr1, origin, 90.0);
        //auto lra0 = sk.new_element<LineSegment>(vr1, vra0);
        //auto lra1 = sk.new_element<LineSegment>(vr1p, vra1);

        //continuous_airgap = sk.new_element<CircularArc>(vra0, vra1, origin, rro + airgap / 2.0);
For faster browsing, not all history is shown. View entire blame