samurai/demos/tutorial/proj_on_mesh.cpp at main · hpc-maths/samurai · GitHub

1
2
3
4
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
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

#include <iostream>

#include <samurai/field.hpp>
#include <samurai/io/hdf5.hpp>
#include <samurai/mr/adapt.hpp>
#include <samurai/mr/mesh.hpp>
#include <samurai/reconstruction.hpp>
#include <samurai/samurai.hpp>

template <class Mesh>
auto init_field(Mesh& mesh, double dec)
{
    static constexpr std::size_t dim = Mesh::dim;
    auto f                           = samurai::make_scalar_field<double>("u", mesh);
    f.fill(0.);
    samurai::for_each_cell(
        mesh,
        [&](const auto& cell)
        {
            auto x = cell.center(0);
            if constexpr (dim == 1)
            {
                if (x + dec > 0.4 && x + dec < 0.6)
                {
                    f[cell] = 1;
                }
            }
            else if constexpr (dim == 2)
            {
                auto y              = cell.center(1);
                const double radius = .2;
                if ((x - dec - .5) * (x - dec - .5) + (y - dec - .5) * (y - dec - .5) < radius * radius)
                {
                    f[cell] = 1;
                }
            }
            else if constexpr (dim == 3)
            {
                auto y              = cell.center(1);
                auto z              = cell.center(2);
                const double radius = .2;
                if ((x - dec - .5) * (x - dec - .5) + (y - dec - .5) * (y - dec - .5) + (z - dec - .5) * (z - dec - .5) < radius * radius)
                {
                    f[cell] = 1;
                }
            }
        });
    return f;
}

int main()
{
    samurai::initialize();

    constexpr std::size_t dim = 3;
    auto box                  = samurai::Box<double, dim>({0., 0., 0.}, {1., 1., 1.});

    auto mesh_cfg = samurai::mesh_config<dim>().min_level(2).max_level(8);
    auto mesh1    = samurai::mra::make_mesh(box, mesh_cfg);
    auto f1       = init_field(mesh1, 0);

    auto mesh2 = samurai::mra::make_mesh(box, mesh_cfg);
    auto f2    = init_field(mesh2, 0.1);

    auto adapt_1 = samurai::make_MRAdapt(f1);
    auto adapt_2 = samurai::make_MRAdapt(f2);

    auto mra_config = samurai::mra_config().epsilon(1e-3).regularity(2.);
    adapt_1(mra_config);
    adapt_2(mra_config);

    samurai::update_ghost_mr(f1);
    samurai::transfer(f1, f2);

    samurai::save("solution_src", mesh1, f1);
    samurai::save("solution_dst", mesh2, f2);

    samurai::finalize();
    return 0;
}