diff --git a/include/maxwell_interface.h b/include/maxwell_interface.h
index 0c6a57a66078443657ab31645185f66a267cf1e1..27d62b495f59c4b5ff92e2052e34ef5707f8100f 100644
--- a/include/maxwell_interface.h
+++ b/include/maxwell_interface.h
@@ -34,6 +34,14 @@ enum class boundary_condition
SURFACE_CURRENT
};
+enum class interface_condition
+{
+ UNSPECIFIED,
+ E_FIELD,
+ H_FIELD,
+ SURFACE_CURRENT
+};
+
class parameter_loader : public parameter_loader_base
{
void init(void);
@@ -54,9 +62,11 @@ public:
bool initial_Hfield_defined(void) const;
vec3d initial_Hfield(const point_3d&) const;
- vec3d eval_plane_wave_E(int, const point_3d&, double) const;
+ vec3d eval_boundary_source(int, const point_3d&, double) const;
+ vec3d eval_interface_source(int, const point_3d&, double) const;
- boundary_condition boundary_type(int tag) const;
+ boundary_condition boundary_type(int tag) const;
+ interface_condition interface_type(int tag) const;
};
#endif /* HIDE_THIS_FROM_NVCC */
@@ -243,7 +253,6 @@ struct solver_state_gpu
#endif /* ENABLE_GPU_SOLVER */
#ifndef HIDE_THIS_FROM_NVCC
-
void init_from_model(const model&, solver_state&);
void init_matparams(const model&, solver_state&, const parameter_loader&);
void apply_operator(solver_state&, const field&, field&);
@@ -261,6 +270,7 @@ void timestep(solver_state_gpu&);
void init_boundary_conditions(const model&, const parameter_loader&, vecxd&);
void init_lax_milgram(const model&, const parameter_loader&, vecxd&, vecxd&, vecxd&, vecxd&);
void eval_boundary_sources(const model&, const parameter_loader&, solver_state&, field&);
+void eval_interface_sources(const model&, const parameter_loader&, solver_state&, field&);
#endif /* HIDE_THIS_FROM_NVCC */
template<typename Function>
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index 7a4ff6924c83ff6a8bd8c865c83c45fa4770556d..9f5dff79513e6f7546f8a86479cfe0c04240347e 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -297,8 +297,39 @@ parameter_loader::boundary_type(int tag) const
return boundary_condition::UNSPECIFIED;
}
+interface_condition
+parameter_loader::interface_type(int tag) const
+{
+ auto bnd_data = lua["ifaceconds"][tag];
+ if (not bnd_data.valid())
+ return interface_condition::UNSPECIFIED;
+
+ auto kind = bnd_data["kind"];
+ if (not kind.valid())
+ {
+ std::cout << "[CONFIG] warning: 'kind' not specified on interface ";
+ std::cout << tag << std::endl;
+ return interface_condition::UNSPECIFIED;
+ }
+
+ std::string kind_str = kind;
+
+ if (kind_str == "E_field")
+ return interface_condition::E_FIELD;
+
+ if (kind_str == "H_field")
+ return interface_condition::H_FIELD;
+
+ if (kind_str == "surface_current")
+ return interface_condition::SURFACE_CURRENT;
+
+ std::cout << "[CONFIG] warning: 'kind' invalid on boundary ";
+ std::cout << tag << std::endl;
+ return interface_condition::UNSPECIFIED;
+}
+
vec3d
-parameter_loader::eval_plane_wave_E(int tag, const point_3d& pt, double t) const
+parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) const
{
vec3d ret;
sol::tie(ret(0), ret(1), ret(2)) =
@@ -306,4 +337,13 @@ parameter_loader::eval_plane_wave_E(int tag, const point_3d& pt, double t) const
return ret;
}
+vec3d
+parameter_loader::eval_interface_source(int tag, const point_3d& pt, double t) const
+{
+ vec3d ret;
+ sol::tie(ret(0), ret(1), ret(2)) =
+ lua["ifaceconds"][tag]["source"](tag, pt.x(), pt.y(), pt.z(), t);
+ return ret;
+}
+
} // namespace maxwell
\ No newline at end of file
diff --git a/src/maxwell_cpu.cpp b/src/maxwell_cpu.cpp
index 3ce0b419e730f915811fc697888f3cba91ccf195..c8ffddb7d6eebf63fab84eb599cbce73d3669e1e 100644
--- a/src/maxwell_cpu.cpp
+++ b/src/maxwell_cpu.cpp
@@ -80,10 +80,10 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
case boundary_condition::PLANE_WAVE_E: {
auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_plane_wave_E(bd.gmsh_entity, pt, state.curr_time);
+ return mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
};
auto fH = [&](const point_3d& pt) -> vec3d {
- vec3d E = mpl.eval_plane_wave_E(bd.gmsh_entity, pt, state.curr_time);
+ vec3d E = mpl.eval_boundary_source(bd.gmsh_entity, pt, state.curr_time);
return n.cross(E)/Z;
};
matxd prjE = e.project_on_face(iF, fE);
@@ -116,6 +116,70 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
}
}
+void
+eval_interface_sources(const model& mod, const parameter_loader& mpl,
+ solver_state& state, field& srcfield)
+{
+ auto& bds = mod.boundary_descriptors();
+ assert( srcfield.num_dofs == mod.num_fluxes() );
+
+ size_t face_num_base = 0;
+ for (auto& e : mod)
+ {
+ for (size_t iF = 0; iF < e.num_faces(); iF++)
+ {
+ auto& pf = e.face(iF);
+ auto bar = pf.barycenter();
+ auto eps = mpl.epsilon(e.gmsh_tag(), bar);
+ auto mu = mpl.mu(e.gmsh_tag(), bar);
+ auto Z = std::sqrt(mu/eps);
+
+ vec3d n = state.eds[e.number()].normals.row(iF);
+
+ auto gfnum = face_num_base + iF;
+ auto bd = bds[gfnum];
+ if (bd.type == face_type::NONE or bd.type == face_type::BOUNDARY)
+ continue;
+
+ switch ( mpl.interface_type(bd.gmsh_entity) )
+ {
+ case interface_condition::E_FIELD: {
+ /* This is BS actually... */
+ auto fE = [&](const point_3d& pt) -> vec3d {
+ return mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ };
+
+ matxd prjE = e.project_on_face(iF, fE);
+ auto num_bf = prjE.rows();
+ srcfield.Ex.segment(gfnum*num_bf, num_bf) = prjE.col(0);
+ srcfield.Ey.segment(gfnum*num_bf, num_bf) = prjE.col(1);
+ srcfield.Ez.segment(gfnum*num_bf, num_bf) = prjE.col(2);
+ } break;
+
+ case interface_condition::SURFACE_CURRENT: {
+ auto fJ = [&](const point_3d& pt) -> vec3d {
+ vec3d J = mpl.eval_interface_source(bd.gmsh_entity, pt, state.curr_time);
+ return n.cross(J);
+ };
+
+ matxd prjJ = e.project_on_face(iF, fJ);
+ auto num_bf = prjJ.rows();
+ srcfield.Hx.segment(gfnum*num_bf, num_bf) = prjJ.col(0);
+ srcfield.Hy.segment(gfnum*num_bf, num_bf) = prjJ.col(1);
+ srcfield.Hz.segment(gfnum*num_bf, num_bf) = prjJ.col(2);
+ } break;
+
+ case interface_condition::H_FIELD:
+ break;
+
+ default:
+ break;
+ }
+ }
+
+ face_num_base += e.num_faces();
+ }
+}
void init_from_model(const model& mod, solver_state& state)
{
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 68e43718658239dbdba7bfd003cd84810f765b27..5e3fdfe0e5c6565bde87f20e80c568a1eb7dd521 100644
--- a/src/maxwell_solver.cpp
+++ b/src/maxwell_solver.cpp
@@ -65,6 +65,7 @@ void test_it(const model& mod, State& state, const maxwell::parameter_loader& mp
for(size_t i = 0; i < num_timesteps; i++)
{
maxwell::eval_boundary_sources(mod, mpl, state, state.bndsrcs);
+ maxwell::eval_interface_sources(mod, mpl, state, state.bndsrcs);
timestep(state);
std::stringstream ss;
ss << "maxwell_" << i << ".silo";