diff --git a/include/libgmshdg/gmsh_io.h b/include/libgmshdg/gmsh_io.h
index e9baad28c4227f147a3d7a61e54852b18411bbff..441b7190c4452c13203939c11f652b5f65550f7b 100644
--- a/include/libgmshdg/gmsh_io.h
+++ b/include/libgmshdg/gmsh_io.h
@@ -132,25 +132,6 @@ operator<<(std::ostream& os, const comm_descriptor& cd)
#define IMPLIES(a, b) ((not(a)) or (b))
#define IFF(a,b) (IMPLIES(a,b) and IMPLIES(b,a))
-struct surface_identifier
-{
- int surface_entity;
- int adjacent_volume_entity;
-
- bool operator<(const surface_identifier& other) const
- {
- return (surface_entity < other.surface_entity) or
- ((surface_entity == other.surface_entity) and
- (adjacent_volume_entity < other.adjacent_volume_entity));
- }
-};
-
-inline std::ostream&
-operator<<(std::ostream& os, const surface_identifier& si)
-{
- os << "(" << si.surface_entity << ", " << si.adjacent_volume_entity << ")";
- return os;
-}
struct surface_descriptor
{
@@ -244,8 +225,8 @@ class model
std::vector<boundary_descriptor> bnd_descriptors;
element_connectivity<element_key> conn;
- std::map<surface_identifier, surface_descriptor> boundaries;
- std::map<surface_identifier, surface_descriptor> interfaces;
+ std::map<int, surface_descriptor> boundaries;
+ std::map<int, surface_descriptor> interfaces;
#ifdef USE_MPI
interprocess_connectivity<element_key> ipconn;
@@ -284,7 +265,6 @@ class model
void populate_from_gmsh_rank0(void);
#ifdef USE_MPI
- bool is_interprocess_boundary(int) const;
void populate_from_gmsh_rankN(void);
void import_gmsh_entities_rankN(void);
int my_partition(void) const;
@@ -315,13 +295,15 @@ public:
void map_boundaries_new(void);
- const std::map<surface_identifier, surface_descriptor>&
+ const std::map<int, surface_descriptor>&
get_boundaries() const;
- const std::map<surface_identifier, surface_descriptor>&
+ const std::map<int, surface_descriptor>&
get_interfaces() const;
#ifdef USE_MPI
+ bool is_interprocess_boundary(int) const;
+
const std::map<int, comm_descriptor>& comm_descriptors(void) const {
return ipc_boundary_comm_table;
}
diff --git a/share/validation/params_twomat.lua b/share/validation/params_twomat.lua
index 8be05f653ce4dbc40b5f128e4d369f16d487145b..7b640ee1a5b34e19cfc53d6a35a35dba9ad8a9d2 100644
--- a/share/validation/params_twomat.lua
+++ b/share/validation/params_twomat.lua
@@ -13,6 +13,9 @@ postpro.cycle_print_rate = 100 -- console print rate
postpro["E"].mode = "nodal"
+debug = {};
+debug.dump_cell_ranks = true;
+
function setup_materials()
local bulk = volume_physical_group_entities(1);
for i,v in ipairs(bulk) do
diff --git a/share/validation/twomat.geo b/share/validation/twomat.geo
index 5ca38adf475090dac8a21ddef1b1007ae6a9b574..c37a400240b677633590270f1a20f5dd8972a8c5 100644
--- a/share/validation/twomat.geo
+++ b/share/validation/twomat.geo
@@ -9,4 +9,4 @@ Physical Surface("pmc", 11) = {3,4,8,9};
Physical Surface("pec", 12) = {5,6,10,11};
Physical Surface("abc", 13) = {7};
-MeshSize{:} = 0.075;
+MeshSize{:} = 0.05;
diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp
index 80119b1e74212e087750a3ce42cc4cf2c73ca3a1..3212b70eba93b57e965f7cb8622eccf4cf88736d 100644
--- a/src/libgmshdg/gmsh_io.cpp
+++ b/src/libgmshdg/gmsh_io.cpp
@@ -528,37 +528,27 @@ model::map_boundaries_new(void)
if ( (itor == bfk.end()) or not ((*itor).first == fk) )
continue;
- surface_identifier si;
- si.surface_entity = (*itor).second;
- si.adjacent_volume_entity = iE;
+ int si = (*itor).second;
- using sm_t = std::map<surface_identifier, surface_descriptor>;
+ using sm_t = std::map<int, surface_descriptor>;
auto add_face = [&](sm_t& m, size_t se, size_t ave) -> void {
if ( m.find(si) == m.end() )
{
m[si].surface_entity = se;
m[si].adjacent_volume_entity = ave;
-#ifdef USE_MPI
- if (num_partitions > 1)
- {
- auto sitor = surface_to_parent_map.find( se );
- if ( sitor != surface_to_parent_map.end() )
- m[si].parent_entity = surface_to_parent_map.at( se );
- }
-#endif /* USE_MPI */
}
m[si].add_face(iF, fbase+iF);
};
- switch ( face_type(si.surface_entity, fk) )
+ switch ( face_type(si, fk) )
{
case face_type::BOUNDARY:
- add_face(boundaries, si.surface_entity, si.adjacent_volume_entity);
+ add_face(boundaries, si, iE);
break;
case face_type::INTERFACE:
- add_face(interfaces, si.surface_entity, si.adjacent_volume_entity);
+ add_face(interfaces, si, iE);
break;
default:
@@ -584,13 +574,13 @@ model::map_boundaries_new(void)
}
}
-const std::map<surface_identifier, surface_descriptor>&
+const std::map<int, surface_descriptor>&
model::get_boundaries() const
{
return boundaries;
}
-const std::map<surface_identifier, surface_descriptor>&
+const std::map<int, surface_descriptor>&
model::get_interfaces() const
{
return interfaces;
@@ -661,14 +651,7 @@ model::map_boundaries(void)
}
/* and store the gmsh tag in the descriptor. */
bd.gmsh_entity = (*itor).second;
-#ifdef USE_MPI
- if (num_partitions > 1)
- {
- auto sitor = surface_to_parent_map.find( (*itor).second );
- if ( sitor != surface_to_parent_map.end() )
- bd.parent_entity = surface_to_parent_map.at( (*itor).second );
- }
-#endif /* USE_MPI */
+
/* Finally, put the descriptor in the global array of faces. */
bnd_descriptors.at(fbase + iF) = bd;
}
diff --git a/src/maxwell/maxwell_common.cpp b/src/maxwell/maxwell_common.cpp
index eec2b36b32ddb94e39c44b83faeeda61c229765d..ea7cba2aefbd2d1ac5ba51710e5c8408bc6a352b 100644
--- a/src/maxwell/maxwell_common.cpp
+++ b/src/maxwell/maxwell_common.cpp
@@ -23,7 +23,7 @@ namespace maxwell {
* coefficients 0, 1 and 2 to select the appropriate boundary condition.
* See the jump kernels for details about the 0, 1, and 2. */
void
-init_boundary_conditions(const model& mod, const parameter_loader& mpl,
+init_boundary_conditions_old(const model& mod, const parameter_loader& mpl,
vecxd& bcc)
{
auto bds = mod.boundary_descriptors();
@@ -82,6 +82,56 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
}
}
+/* Take a model and a parameter loader and return a vector with the
+ * coefficients 0, 1 and 2 to select the appropriate boundary condition.
+ * See the jump kernels for details about the 0, 1, and 2. */
+void
+init_boundary_conditions(const model& mod, const parameter_loader& mpl,
+ vecxd& bcc)
+{
+ bcc = vecxd::Ones( mod.num_fluxes() );
+
+ auto& bds = mod.get_boundaries();
+
+ for (auto& [si, sd] : bds)
+ {
+ const auto& e = mod.entity_at(sd.adjacent_volume_entity);
+ auto mtag = sd.gmsh_tag();
+
+ for (auto& [local, global] : sd.face_numbers)
+ {
+ double bc_coeff = DIRICHLET;
+
+ switch ( mpl.boundary_type(mtag) )
+ {
+ case boundary_condition::UNSPECIFIED:
+ case boundary_condition::PEC:
+ case boundary_condition::E_FIELD:
+ bc_coeff = DIRICHLET;
+ break;
+
+ case boundary_condition::IMPEDANCE:
+ case boundary_condition::PLANE_WAVE_E:
+ case boundary_condition::PLANE_WAVE_H:
+ bc_coeff = ROBIN;
+ break;
+
+ case boundary_condition::PMC:
+ case boundary_condition::SURFACE_CURRENT:
+ case boundary_condition::H_FIELD:
+ bc_coeff = NEUMANN;
+ break;
+ }
+
+ auto& pf = e.face(local);
+ auto& rf = e.face_refelem(pf);
+ auto num_bf = rf.num_basis_functions();
+ for (size_t i = 0; i < num_bf; i++)
+ bcc(global*num_bf + i) = bc_coeff;
+ }
+ }
+}
+
/* Take a model and a parameter loader and return the four vectors with
* the Lax-Milgram flux coefficients */
void