From caec3f240dec50efb692907c2f3787d23bb590ff Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Mon, 4 Apr 2022 23:15:57 +0200
Subject: [PATCH] Added EH boundary condition.
---
include/maxwell/maxwell_common.h | 21 +++++++++++++++++++++
include/maxwell/maxwell_interface.h | 2 ++
src/maxwell/maxwell_gpu.cpp | 1 +
src/maxwell/maxwell_interface.cpp | 15 ++++++++++++++-
4 files changed, 38 insertions(+), 1 deletion(-)
diff --git a/include/maxwell/maxwell_common.h b/include/maxwell/maxwell_common.h
index d4cfc6b..4634817 100644
--- a/include/maxwell/maxwell_common.h
+++ b/include/maxwell/maxwell_common.h
@@ -62,6 +62,9 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
case boundary_condition::H_FIELD:
break;
+ case boundary_condition::EH_FIELD:
+ break;
+
case boundary_condition::IMPEDANCE:
break;
@@ -164,6 +167,24 @@ eval_boundary_sources_new(const model& mod, const parameter_loader& mpl,
case boundary_condition::SURFACE_CURRENT:
break;
+ case boundary_condition::EH_FIELD: {
+ for (auto& [local, global] : sd.face_numbers)
+ {
+ auto fEH = [&](const point_3d& pt) -> std::pair<vec3d, vec3d> {
+ return mpl.eval_boundary_source_eh(mtag, pt, state.curr_time);
+ };
+
+ matxd prj = project_EH_to_face(e, local, fEH);
+ auto num_bf = prj.rows();
+ srcfield.Ex.segment(global*num_bf, num_bf) = prj.col(0);
+ srcfield.Ey.segment(global*num_bf, num_bf) = prj.col(1);
+ srcfield.Ez.segment(global*num_bf, num_bf) = prj.col(2);
+ srcfield.Hx.segment(global*num_bf, num_bf) = prj.col(3);
+ srcfield.Hy.segment(global*num_bf, num_bf) = prj.col(4);
+ srcfield.Hz.segment(global*num_bf, num_bf) = prj.col(5);
+ }
+ } break; /* case boundary_condition::EH_FIELD */
+
case boundary_condition::PLANE_WAVE_E: {
for (auto& [local, global] : sd.face_numbers)
{
diff --git a/include/maxwell/maxwell_interface.h b/include/maxwell/maxwell_interface.h
index af0379e..e9bcfde 100644
--- a/include/maxwell/maxwell_interface.h
+++ b/include/maxwell/maxwell_interface.h
@@ -43,6 +43,7 @@ enum class boundary_condition
PLANE_WAVE_H,
E_FIELD,
H_FIELD,
+ EH_FIELD,
SURFACE_CURRENT
};
@@ -80,6 +81,7 @@ public:
bool volume_has_source(int) const;
vec3d eval_volume_source(int, const point_3d&, double) const;
vec3d eval_boundary_source(int, const point_3d&, double) const;
+ std::pair<vec3d,vec3d> eval_boundary_source_eh(int, const point_3d&, double) const;
vec3d eval_interface_source(int, const point_3d&, double) const;
bool has_analytical_solution(void) const;
diff --git a/src/maxwell/maxwell_gpu.cpp b/src/maxwell/maxwell_gpu.cpp
index e13aff2..25b7a28 100644
--- a/src/maxwell/maxwell_gpu.cpp
+++ b/src/maxwell/maxwell_gpu.cpp
@@ -105,6 +105,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::PLANE_WAVE_E:
+ case boundary_condition::EH_FIELD:
for (size_t i = 0; i < num_fluxes; i++)
state.bndsrcs_decomp_table_cpu.push_back(gfnum*num_fluxes+i);
break;
diff --git a/src/maxwell/maxwell_interface.cpp b/src/maxwell/maxwell_interface.cpp
index 32899c4..5f3921d 100644
--- a/src/maxwell/maxwell_interface.cpp
+++ b/src/maxwell/maxwell_interface.cpp
@@ -489,6 +489,7 @@ material_params_gpu::copyin(const material_params& mp)
#define BCOND_TYPE_HPLANEW "plane_wave_H"
#define BCOND_TYPE_EFIELD "E_field"
#define BCOND_TYPE_HFIELD "H_field"
+#define BCOND_TYPE_EHFIELD "EH_field"
#define BCOND_TYPE_SURFCURR "surface_current"
#define SEC_IFCONDS "ifaceconds"
@@ -581,7 +582,7 @@ parameter_loader::validate_boundary_condition(const model&, int tag) const
if (bc_kind == BCOND_TYPE_IMPEDANCE)
return true;
- if (bc_kind == BCOND_TYPE_EPLANEW)
+ if (bc_kind == BCOND_TYPE_EPLANEW || bc_kind == BCOND_TYPE_EHFIELD)
{
if (bc[BC_SOURCE].valid())
return true;
@@ -791,6 +792,9 @@ parameter_loader::boundary_type(int tag) const
if (kind_str == BCOND_TYPE_HFIELD)
return boundary_condition::H_FIELD;
+
+ if (kind_str == BCOND_TYPE_EHFIELD)
+ return boundary_condition::EH_FIELD;
if (kind_str == BCOND_TYPE_SURFCURR)
return boundary_condition::SURFACE_CURRENT;
@@ -846,6 +850,15 @@ parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) co
return ret;
}
+std::pair<vec3d, vec3d>
+parameter_loader::eval_boundary_source_eh(int tag, const point_3d& pt, double t) const
+{
+ vec3d retE, retH;
+ sol::tie(retE(0), retE(1), retE(2), retH(0), retH(1), retH(2)) =
+ lua[SEC_BCONDS][tag][BC_SOURCE](tag, pt.x(), pt.y(), pt.z(), t);
+ return std::make_pair(retE, retH);
+}
+
vec3d
parameter_loader::eval_interface_source(int tag, const point_3d& pt, double t) const
{
--
GitLab