maxwell_solver.cpp

#include "gmsh.h"
#include "silo_output.hpp"
#include "maxwell_interface.h"

static void
make_geometry(int order, double mesh_h)
{
    gm::add("difftest");

    std::vector<std::pair<int,int>> objects;

    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.0, 1.0, 1.0, 0.05) )
        );
    objects.push_back(
        std::make_pair(3, gmo::addBox(0.0, 0.0, 0.05, 1.0, 1.0, 0.05) )
        );

    std::vector<std::pair<int, int>> tools;
    gmsh::vectorpair odt;
    std::vector<gmsh::vectorpair> odtm;
    gmo::fragment(objects, tools, odt, odtm);

    gmo::synchronize();

    gvp_t vp;
    gm::getEntities(vp);
    gmm::setSize(vp, mesh_h);
}

template<typename State>
void test_it(const model& mod, State& state)
{
    auto E = [](const point_3d& pt) -> vec3d {
        vec3d ret;
        ret(0) = 0;
        ret(1) = 0;
        ret(2) = std::sin(M_PI*pt.x()) * std::sin(M_PI*pt.y());
        return ret;
    };

    maxwell::init_from_model(mod, state);
    maxwell::init_E_field(mod, state, E);

    state.delta_t = 0.001;

    for(size_t i = 0; i < 10000; i++)
    {
        timestep(state);
        std::stringstream ss;
        ss << "maxwell_" << i << ".silo";

        if (i%10 == 0)
            maxwell::export_to_silo(mod, state, ss.str());
        
        std::swap(state.emf_curr, state.emf_next);
    }

    
}

int main(int argc, const char *argv[])
{
    gmsh::initialize();
    make_geometry(1, 0.1);

    model mod(1,3);
    
    maxwell::solver_state state_c;
    test_it(mod, state_c);

    //maxwell::solver_state_gpu state_g;
    //test_it(mod, state_g);

    //gmsh::finalize();

    return 0;
}