diff --git a/include/maxwell/maxwell_common.h b/include/maxwell/maxwell_common.h index 199743514a0214b215bacf3658db73bdf4c16871..d4cfc6be3889b81c37d22f5f671dcbb67ae94c0c 100644 --- a/include/maxwell/maxwell_common.h +++ b/include/maxwell/maxwell_common.h @@ -265,6 +265,47 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl, } } +/* Field must be either 'field' or 'pinned_field'. Too early to use a concept. */ +template<typename State, typename Field> +void +eval_interface_sources_new(const model& mod, const parameter_loader& mpl, + State& state, Field& srcfield) +{ + auto& is = mod.get_interfaces(); + assert( srcfield.num_dofs == mod.num_fluxes() ); + + for (auto& [si, sd] : is) + { + const auto& e = mod.entity_at(sd.adjacent_volume_entity); + //auto etag = e.gmsh_tag(); + auto mtag = sd.gmsh_tag(); + switch ( mpl.interface_type( mtag ) ) + { + case interface_condition::UNSPECIFIED: + case interface_condition::E_FIELD: + case interface_condition::H_FIELD: + break; + + case interface_condition::SURFACE_CURRENT: { + for (auto& [local, global] : sd.face_numbers) + { + vec3d n = state.eds[e.number()].normals.row(local); + auto fJ = [&](const point_3d& pt) -> vec3d { + vec3d J = mpl.eval_interface_source(mtag, pt, state.curr_time); + return n.cross(J); + }; + + matxd prjJ = e.project_on_face(local, fJ); + auto num_bf = prjJ.rows(); + srcfield.Hx.segment(global*num_bf, num_bf) = prjJ.col(0); + srcfield.Hy.segment(global*num_bf, num_bf) = prjJ.col(1); + srcfield.Hz.segment(global*num_bf, num_bf) = prjJ.col(2); + } + } break; /* case interface_condition::SURFACE_CURRENT */ + } + } +} + #ifdef USE_MPI std::vector<double> compute_cell_ranks(const model&); #endif /* USE_MPI */ diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp index 7fcedcaf870f78ac94fa576dd4f6e066f055f256..80119b1e74212e087750a3ce42cc4cf2c73ca3a1 100644 --- a/src/libgmshdg/gmsh_io.cpp +++ b/src/libgmshdg/gmsh_io.cpp @@ -576,10 +576,10 @@ model::map_boundaries_new(void) sd.sort(); } - //std::cout << "Interfaces:" << std::endl; + std::cout << "Interfaces:" << std::endl; for (auto& [si, sd] : interfaces) { - //std::cout << si << " " << sd.face_numbers.size() << std::endl; + std::cout << si << " " << sd.face_numbers.size() << std::endl; sd.sort(); } }