diff --git a/doc/lua_api.md b/doc/lua_api.md
index decb61fea02e35e505145f590e8b6d89f546ca4f..2d01e98f347e174750129f32c53fdcb31caff73e 100644
--- a/doc/lua_api.md
+++ b/doc/lua_api.md
@@ -32,6 +32,7 @@ This document describes the API available on the version v0.4 of the solver.
### Parallel execution information (available only if MPI support is enabled)
The parallel solver runs as separate MPI processes. As such, each process loads and executes its own instance of the configuration script. This means that you must beware of global variables or other shared state, because modifications are **not** reflected across processes. Even if the configuration script should in general not depend on the rank of the current process, the following variables are exposed:
+
- `parallel.comm_rank` (integer): The MPI rank of the current solver process
- `parallel.comm_size` (integer): The MPI communicator size
diff --git a/include/libgmshdg/gmsh_io.h b/include/libgmshdg/gmsh_io.h
index c33ac53e6bf7304ff5add92c4a8e899ec892a651..e9baad28c4227f147a3d7a61e54852b18411bbff 100644
--- a/include/libgmshdg/gmsh_io.h
+++ b/include/libgmshdg/gmsh_io.h
@@ -73,6 +73,11 @@ struct boundary_descriptor
#endif /* USE_MPI */
{}
+ int gmsh_tag(void) const
+ {
+ return gmsh_entity;
+ }
+
int material_tag(void) const
{
#ifdef USE_MPI
@@ -95,6 +100,15 @@ struct physical_group_descriptor {
}
};
+inline std::ostream&
+operator<<(std::ostream& os, const physical_group_descriptor& pgd)
+{
+ os << "["<< pgd.dim << ", " << pgd.tag << "] ";
+ for (auto& et : pgd.entity_tags)
+ os << et << " ";
+ return os;
+}
+
#ifdef USE_MPI
struct comm_descriptor
{
@@ -107,8 +121,8 @@ struct comm_descriptor
std::vector<element_key> fks;
};
-inline
-std::ostream& operator<<(std::ostream& os, const comm_descriptor& cd)
+inline std::ostream&
+operator<<(std::ostream& os, const comm_descriptor& cd)
{
os << "tag: " << cd.boundary_tag << ", PM: " << cd.partition_mine << ", PO: " << cd.partition_other;
os << ", RM: " << cd.rank_mine << ", RO: " << cd.rank_other << ", map size: " << cd.dof_mapping.size();
@@ -178,6 +192,11 @@ struct surface_descriptor
std::sort(face_numbers.begin(), face_numbers.end(), comp);
}
+ int gmsh_tag(void) const
+ {
+ return surface_entity;
+ }
+
int material_tag(void) const
{
#ifdef USE_MPI
diff --git a/include/maxwell/maxwell_common.h b/include/maxwell/maxwell_common.h
index 1e902e1efef6f48c4ab87530bae16a7a8946607e..199743514a0214b215bacf3658db73bdf4c16871 100644
--- a/include/maxwell/maxwell_common.h
+++ b/include/maxwell/maxwell_common.h
@@ -25,7 +25,7 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
size_t face_num_base = 0;
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iF = 0; iF < e.num_faces(); iF++)
{
auto gfnum = face_num_base + iF;
@@ -45,7 +45,7 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
vec3d n = state.eds[e.number()].normals.row(iF);
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::UNSPECIFIED:
break;
@@ -67,10 +67,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_boundary_source(bd.material_tag(), pt, state.curr_time);
+ return mpl.eval_boundary_source(bd.gmsh_tag(), pt, state.curr_time);
};
auto fH = [&](const point_3d& pt) -> vec3d {
- vec3d H = mpl.eval_boundary_source(bd.material_tag(), pt, state.curr_time);
+ vec3d H = mpl.eval_boundary_source(bd.gmsh_tag(), pt, state.curr_time);
return n.cross(H)/Z;
};
matxd prjE = e.project_on_face(iF, fE);
@@ -149,9 +149,9 @@ eval_boundary_sources_new(const model& mod, const parameter_loader& mpl,
for (auto& [si, sd] : bds)
{
const auto& e = mod.entity_at(sd.adjacent_volume_entity);
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
- auto mtag = sd.material_tag();
+ auto mtag = sd.gmsh_tag();
switch ( mpl.boundary_type( mtag ) )
{
case boundary_condition::UNSPECIFIED:
@@ -219,7 +219,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
vec3d n = state.eds[e.number()].normals.row(iF);
- switch ( mpl.interface_type(bd.material_tag()) )
+ switch ( mpl.interface_type(bd.gmsh_tag()) )
{
case interface_condition::UNSPECIFIED:
break;
@@ -227,7 +227,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
case interface_condition::E_FIELD: {
/* This is BS actually... */
auto fE = [&](const point_3d& pt) -> vec3d {
- return mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
+ return mpl.eval_interface_source(bd.gmsh_tag(), pt, state.curr_time);
};
matxd prjE = e.project_on_face(iF, fE);
@@ -242,7 +242,7 @@ eval_interface_sources(const model& mod, const parameter_loader& mpl,
case interface_condition::SURFACE_CURRENT: {
auto fJ = [&](const point_3d& pt) -> vec3d {
- vec3d J = mpl.eval_interface_source(bd.material_tag(), pt, state.curr_time);
+ vec3d J = mpl.eval_interface_source(bd.gmsh_tag(), pt, state.curr_time);
return n.cross(J);
};
diff --git a/share/validation/params_twomat.lua b/share/validation/params_twomat.lua
index 67b01602f1f96a415bf3cf2788a11e4a21ea18c1..8be05f653ce4dbc40b5f128e4d369f16d487145b 100644
--- a/share/validation/params_twomat.lua
+++ b/share/validation/params_twomat.lua
@@ -13,35 +13,47 @@ postpro.cycle_print_rate = 100 -- console print rate
postpro["E"].mode = "nodal"
-materials[1] = {}
-materials[1].epsilon = 1
-materials[1].mu = 1
-materials[1].sigma = 0
-
-materials[2] = {}
-materials[2].epsilon = 4
-materials[2].mu = 1
-materials[2].sigma = 0
-
--- ** Boundary conditions **
-local PMC_bnds = { 3, 4, 8, 9 }
-for i,v in ipairs(PMC_bnds) do
- bndconds[v] = {}
- bndconds[v].kind = "pmc"
+function setup_materials()
+ local bulk = volume_physical_group_entities(1);
+ for i,v in ipairs(bulk) do
+ materials[v] = {}
+ materials[v].epsilon = 1
+ materials[v].mu = 1
+ materials[v].sigma = 0
+ end
end
-bndconds[7] = {}
-bndconds[7].kind = "impedance"
-
local freq = 3e8
-function source(tag, x, y, z, t)
+function plane_wave(tag, x, y, z, t)
local Ex = 0
local Ey = 0
local Ez = math.sin(2*math.pi*t*freq)
return Ex, Ey, Ez
end
-bndconds[1] = {}
-bndconds[1].kind = "plane_wave_E"
-bndconds[1].source = source
+function setup_boundary_conditions()
+ local source = surface_physical_group_entities(10);
+ for i,v in ipairs(source) do
+ bndconds[v] = {};
+ bndconds[v].kind = "plane_wave_E";
+ bndconds[v].source = plane_wave;
+ end
+
+ local pmc = surface_physical_group_entities(11);
+ for i,v in ipairs(pmc) do
+ bndconds[v] = {};
+ bndconds[v].kind = "pmc";
+ end
+
+ local abc = surface_physical_group_entities(13);
+ for i,v in ipairs(abc) do
+ bndconds[v] = {};
+ bndconds[v].kind = "impedance";
+ end
+end
+
+function on_mesh_loaded()
+ setup_materials();
+ setup_boundary_conditions();
+end
diff --git a/share/validation/twomat.geo b/share/validation/twomat.geo
index 6545506a9a6c791431f505a4a12411f04b440e31..5ca38adf475090dac8a21ddef1b1007ae6a9b574 100644
--- a/share/validation/twomat.geo
+++ b/share/validation/twomat.geo
@@ -3,4 +3,10 @@ Box(1) = {0, 0, 0, 1, 1, 0.1};
Box(2) = {1, 0, 0, 1, 1, 0.1};
Coherence;
+Physical Volume("bulk", 1) = {1,2};
+Physical Surface("source", 10) = {1};
+Physical Surface("pmc", 11) = {3,4,8,9};
+Physical Surface("pec", 12) = {5,6,10,11};
+Physical Surface("abc", 13) = {7};
+
MeshSize{:} = 0.075;
diff --git a/src/libgmshdg/gmsh_io.cpp b/src/libgmshdg/gmsh_io.cpp
index b67a6b4897ebdcfdcf27afda73c38182ef5a7ab8..7fcedcaf870f78ac94fa576dd4f6e066f055f256 100644
--- a/src/libgmshdg/gmsh_io.cpp
+++ b/src/libgmshdg/gmsh_io.cpp
@@ -423,7 +423,7 @@ model::import_physical_groups(void)
for (auto& pg : physical_groups)
{
MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ MPI_Bcast(&pg.tag, 1, MPI_INT, 0, MPI_COMM_WORLD);
size_t esize = pg.entity_tags.size();
MPI_Bcast(&esize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
MPI_Bcast(pg.entity_tags.data(), esize, MPI_INT, 0, MPI_COMM_WORLD);
@@ -435,7 +435,9 @@ model::import_physical_groups(void)
for (auto& pg : physical_groups)
{
MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
- MPI_Bcast(&pg.dim, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ assert(pg.dim != 0);
+ MPI_Bcast(&pg.tag, 1, MPI_INT, 0, MPI_COMM_WORLD);
+ assert(pg.tag != 0);
size_t esize;
MPI_Bcast(&esize, 1, MPI_UNSIGNED_LONG_LONG, 0, MPI_COMM_WORLD);
pg.entity_tags.resize(esize);
@@ -1114,10 +1116,10 @@ model::lookup_physical_group_entities(int dim, int tag)
physical_group_descriptor pg;
pg.dim = dim;
pg.tag = tag;
- auto elem = std::lower_bound(physical_groups.begin() , physical_groups.end(), pg);
- if (!(elem == physical_groups.end()) && !(pg < *physical_groups.begin()))
- return elem->entity_tags;
-
+ auto pg_begin = physical_groups.begin();
+ auto pg_end = physical_groups.end();
+ auto itor = std::lower_bound(pg_begin , pg_end, pg);
+ if (!(itor == pg_end) && !(pg < *pg_begin))
+ return itor->entity_tags;
return std::vector<int>{};
}
-
diff --git a/src/libgmshdg/param_loader.cpp b/src/libgmshdg/param_loader.cpp
index 569a95916488c256294a2fc077a5488ceb67bd48..0adc570fbbb9da101715ef297e3b3fbfd7f0db25 100644
--- a/src/libgmshdg/param_loader.cpp
+++ b/src/libgmshdg/param_loader.cpp
@@ -16,25 +16,6 @@
#include "libgmshdg/param_loader.h"
-static auto
-lua_getEntitiesForPhysicalGroup(int dim, int tag)
-{
- std::vector<int> ret;
-#ifdef USE_MPI
- int rank;
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- if (rank != 0)
- {
- std::cout << "[WARNING] getEntitiesForPhysicalGroup() ";
- std::cout << "can be called only on MPI rank 0. On the other ";
- std::cout << "ranks you will get an empty list." << std::endl;
- return sol::as_table(ret);
- }
-#endif
- gmsh::model::getEntitiesForPhysicalGroup(dim, tag, ret);
- return sol::as_table(ret);
-}
-
parameter_loader_base::parameter_loader_base()
: m_bnd_sources_active(true), m_ifc_sources_active(true),
m_vol_sources_active(true)
@@ -86,8 +67,13 @@ parameter_loader_base::init(void)
¶meter_loader_base::enable_volume_sources,
this);
- lua.set_function("getEntitiesForPhysicalGroup",
- &lua_getEntitiesForPhysicalGroup);
+ auto exit_abc = [](){
+#ifdef USE_MPI
+ MPI_Finalize();
+#endif /* USE_MPI */
+ exit(0);
+ };
+ lua["exit"] = exit_abc;
#ifdef USE_MPI
int comm_rank, comm_size;
diff --git a/src/maxwell/maxwell_common.cpp b/src/maxwell/maxwell_common.cpp
index bda7f6d2a73831c828c72468a134b6ba83b3b636..eec2b36b32ddb94e39c44b83faeeda61c229765d 100644
--- a/src/maxwell/maxwell_common.cpp
+++ b/src/maxwell/maxwell_common.cpp
@@ -48,7 +48,7 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
else
{
#endif /* USE_MPI */
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::UNSPECIFIED:
case boundary_condition::PEC:
@@ -103,7 +103,7 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
{
auto& pe = e.cell(iT);
auto bar = pe.barycenter();
- auto tag = e.material_tag();
+ auto tag = e.gmsh_tag();
auto my_eps = mpl.epsilon(tag, bar);
auto my_mu = mpl.mu(tag, bar);
@@ -125,7 +125,7 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
auto& ne_e = mod.entity_at(neigh.iE);
auto& ne_pe = ne_e.cell(neigh.iT);
auto ne_bar = ne_pe.barycenter();
- auto ne_tag = ne_e.material_tag();
+ auto ne_tag = ne_e.gmsh_tag();
auto ne_eps = mpl.epsilon(ne_tag, ne_bar);
auto ne_mu = mpl.mu(ne_tag, ne_bar);
@@ -228,7 +228,7 @@ export_vector_field_nodal(const model& mod, silo& sdb, const vecxd& Fx,
auto& pe = e.cell(iT);
auto cgofs = e.cell_world_dof_offset(iT);
auto nodetags = pe.node_tags();
- auto mp = mf(e.material_tag(), pe.barycenter());
+ auto mp = mf(e.gmsh_tag(), pe.barycenter());
for (size_t iD = 0; iD < 4; iD++)
{
size_t tag = nodetags[iD];
@@ -282,7 +282,7 @@ export_vector_field_zonal(const model& mod, silo& sdb, const vecxd& Fx,
auto& pe = e.cell(iT);
auto& re = e.cell_refelem(pe);
vecxd phi_bar = re.basis_functions({1./3., 1./3., 1./3.});
- auto mp = mf(e.material_tag(), pe.barycenter());
+ auto mp = mf(e.gmsh_tag(), pe.barycenter());
auto num_bf = re.num_basis_functions();
auto ofs = e.cell_world_dof_offset(iT);
auto gi = e.cell_world_index_by_gmsh(iT);
diff --git a/src/maxwell/maxwell_cpu.cpp b/src/maxwell/maxwell_cpu.cpp
index b54294b229ebe4e05ad9a27697d9c3d5a9b71708..e6d7782339af0ba3ed72b48dc74b06ecbdb3d9f9 100644
--- a/src/maxwell/maxwell_cpu.cpp
+++ b/src/maxwell/maxwell_cpu.cpp
@@ -33,7 +33,7 @@ init_matparams(const model& mod, solver_state& state,
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
@@ -787,7 +787,7 @@ eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state&
{
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
diff --git a/src/maxwell/maxwell_gpu.cpp b/src/maxwell/maxwell_gpu.cpp
index 5f95498583d7babe575fb4dbe47bd80d3024da66..e13aff21bc11f7e49aee6360227941af541aba79 100644
--- a/src/maxwell/maxwell_gpu.cpp
+++ b/src/maxwell/maxwell_gpu.cpp
@@ -32,7 +32,7 @@ init_matparams(const model& mod, solver_state_gpu& state,
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);
@@ -87,7 +87,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
auto bd = bds[gfnum];
if (bd.type == face_type::INTERFACE)
{
- switch ( mpl.interface_type(bd.material_tag()) )
+ switch ( mpl.interface_type(bd.gmsh_tag()) )
{
case interface_condition::E_FIELD:
case interface_condition::SURFACE_CURRENT:
@@ -102,7 +102,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
if (bd.type == face_type::BOUNDARY)
{
- switch ( mpl.boundary_type(bd.material_tag()) )
+ switch ( mpl.boundary_type(bd.gmsh_tag()) )
{
case boundary_condition::PLANE_WAVE_E:
for (size_t i = 0; i < num_fluxes; i++)
@@ -594,7 +594,7 @@ eval_volume_sources(const model& mod, const parameter_loader& mpl, solver_state_
{
for (auto& e : mod)
{
- auto etag = e.material_tag();
+ auto etag = e.gmsh_tag();
if ( not mpl.volume_has_source(etag) )
continue;
diff --git a/src/maxwell/maxwell_interface.cpp b/src/maxwell/maxwell_interface.cpp
index ab6c81ee7cf4a5ec11e35384fe09487f4911190d..32899c4ea9f4cb4a96eff973df1865a2bfda578c 100644
--- a/src/maxwell/maxwell_interface.cpp
+++ b/src/maxwell/maxwell_interface.cpp
@@ -524,6 +524,20 @@ parameter_loader::parameter_loader()
bool
parameter_loader::validate_materials(const std::string& mpn, int tag) const
{
+#ifdef USE_MPI
+ int rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ if (rank == 0) {
+#endif /* USE_MPI */
+ int pdim, ptag;
+ gm::getParent(3,tag,pdim,ptag);
+ if (ptag == -1)
+ return true; /* This is not a computational entity, only the parent of a partition. */
+#ifdef USE_MPI
+ } else {
+ return true; /* The script is the same for all ranks, do validation only on rank 0. */
+ }
+#endif /* USE_MPI */
auto mfun = lua[SEC_MATERIALS][mpn];
if ( mfun.valid() )
return true;
@@ -636,12 +650,12 @@ parameter_loader::validate_conditions(const model& mod) const
break;
case face_type::BOUNDARY:
- if (not validate_boundary_condition(mod, bd.material_tag()) )
+ if (not validate_boundary_condition(mod, bd.gmsh_tag()) )
success = false;
break;
case face_type::INTERFACE:
- if (not validate_interface_condition(mod, bd.material_tag()) )
+ if (not validate_interface_condition(mod, bd.gmsh_tag()) )
success = false;
break;
@@ -659,6 +673,7 @@ bool
parameter_loader::validate_model_params(const model& mod) const
{
bool success = true;
+
for (auto& e : mod)
{
auto tag = e.material_tag();
@@ -671,7 +686,7 @@ parameter_loader::validate_model_params(const model& mod) const
if ( not validate_conditions(mod) )
success = false;
-
+
return success;
}
diff --git a/src/maxwell/maxwell_postpro.cpp b/src/maxwell/maxwell_postpro.cpp
index 8490a59a1b7f02de7fe22136cdac545b9244905b..c84baf15eeb43068c35407cc6c199cfc599e124c 100644
--- a/src/maxwell/maxwell_postpro.cpp
+++ b/src/maxwell/maxwell_postpro.cpp
@@ -92,8 +92,8 @@ compute_energy(const model& mod, const solver_state& state, const parameter_load
const auto cbs = re.num_basis_functions();
const auto bar = pe.barycenter();
- auto eps = mpl.epsilon( e.material_tag(), bar );
- auto mu = mpl.mu( e.material_tag(), bar );
+ auto eps = mpl.epsilon( e.gmsh_tag(), bar );
+ auto mu = mpl.mu( e.gmsh_tag(), bar );
field_values mat(eps, eps, eps, mu, mu, mu);
diff --git a/src/maxwell/maxwell_solver.cpp b/src/maxwell/maxwell_solver.cpp
index 600713fefb77618eedc6af048c5bf8d1629df831..7c17ab60f7a2429498ffdffab0a43510f191dcd1 100644
--- a/src/maxwell/maxwell_solver.cpp
+++ b/src/maxwell/maxwell_solver.cpp
@@ -157,7 +157,7 @@ reinit_materials(const model& mod, maxwell::solver_state& state, maxwell::parame
{
for (auto& e : mod)
{
- auto tag = e.material_tag();
+ auto tag = e.gmsh_tag();
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
auto& pe = e.cell(iT);