Parameter validation.

CMakeLists.txt

@@ -280,17 +280,19 @@ target_link_libraries(test_mass gmshdg)
 add_executable(test_differentiation tests/test_differentiation.cpp)
 target_link_libraries(test_differentiation gmshdg)
 
-add_executable(test_differentiation_gpu tests/test_differentiation_gpu.cpp)
-target_link_libraries(test_differentiation_gpu gmshdg)
-
-add_executable(profile_differentiation_gpu tests/profile_differentiation_gpu.cpp)
-target_link_libraries(profile_differentiation_gpu gmshdg)
-
 add_executable(test_lifting tests/test_lifting.cpp)
 target_link_libraries(test_lifting gmshdg)
 
-add_executable(test_lifting_gpu tests/test_lifting_gpu.cpp)
-target_link_libraries(test_lifting_gpu gmshdg)
+if (OPT_ENABLE_GPU_SOLVER)
+    add_executable(test_differentiation_gpu tests/test_differentiation_gpu.cpp)
+    target_link_libraries(test_differentiation_gpu gmshdg)
+
+    add_executable(profile_differentiation_gpu tests/profile_differentiation_gpu.cpp)
+    target_link_libraries(profile_differentiation_gpu gmshdg)
+
+    add_executable(test_lifting_gpu tests/test_lifting_gpu.cpp)
+    target_link_libraries(test_lifting_gpu gmshdg)
+endif()
 
 #add_executable(test test.cpp test_basics.cpp test_mass.cpp test_differentiation.cpp test_differentiation_gpu.cpp test_lifting.cpp ${COMMON_SOURCES})
 #target_link_libraries(test ${LINK_LIBS})

include/gmsh_io.h

@@ -38,6 +38,7 @@ struct boundary_descriptor
     face_type   type;
     int         gmsh_entity;
 
+    auto operator<=>(const boundary_descriptor&) const = default;
     boundary_descriptor()
         : type(face_type::NONE), gmsh_entity(0)
     {}

include/maxwell_interface.h

@@ -49,12 +49,15 @@ class parameter_loader : public parameter_loader_base
 {
     void init(void);
 
-    bool    validate_materials(const std::string&, int);
+    bool    validate_materials(const std::string&, int) const;
+    bool    validate_boundary_condition(const model&, int) const;
+    bool    validate_interface_condition(const model&, int) const;
+    bool    validate_conditions(const model&) const;
 
 public:
     parameter_loader();
 
-    bool    validate_model_params(const model&);
+    bool    validate_model_params(const model&) const;
     double  epsilon(int, const point_3d&) const;
     double  mu(int, const point_3d&) const;
     double  sigma(int, const point_3d&) const;

include/physical_element.h

@@ -93,7 +93,12 @@ class physical_elements_factory
     std::vector<double>     cellDets;
     std::vector<int>        cellOrientations;
     std::vector<double>     barycenters;
+#ifdef USE_INITIAL_4_8_4_API
     gmsh::vectorpair        keypairs;
+#else
+    std::vector<int>        tagKeys;
+    std::vector<size_t>     entityKeys;
+#endif
     int                     keys_per_elem;
 
 public:

src/maxwell_common.cpp

@@ -5,6 +5,33 @@
 #define ROBIN       1.0
 #define NEUMANN     0.0
 
+#define SEC_BCONDS              "bndconds"
+#define BC_KIND                 "kind"
+#define BC_SOURCE               "source"
+#define BCOND_TYPE_PEC          "pec"
+#define BCOND_TYPE_PMC          "pmc"
+#define BCOND_TYPE_IMPEDANCE    "impedance"
+#define BCOND_TYPE_EPLANEW      "plane_wave_E"
+#define BCOND_TYPE_HPLANEW      "plane_wave_H"
+#define BCOND_TYPE_EFIELD       "E_field"
+#define BCOND_TYPE_HFIELD       "H_field"
+#define BCOND_TYPE_SURFCURR     "surface_current"
+
+#define SEC_IFCONDS             "ifaceconds"
+#define IFC_KIND                "kind"
+#define IFC_SOURCE              "source"
+#define IFCOND_TYPE_EFIELD      "E_field"
+#define IFCOND_TYPE_HFIELD      "H_field"
+#define IFCOND_TYPE_SURFCURR    "surface_current"
+
+
+#define SEC_MATERIALS           "materials"
+#define MAT_EPS                 "epsilon"
+#define MAT_EPS0                "eps0"
+#define MAT_MU                  "mu"
+#define MAT_MU0                 "mu0"
+#define MAT_SIGMA               "sigma"
+
 namespace maxwell {
 
 /* Take a model and a parameter loader and return a vector with the
@@ -138,44 +165,156 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
 
 parameter_loader::parameter_loader()
 {
-    lua["const"]["eps0"] = 8.8541878128e-12;
-    lua["const"]["mu0"] = 4e-7*M_PI;
+    lua["const"][MAT_EPS0] = 8.8541878128e-12;
+    lua["const"][MAT_MU0]  = 4e-7*M_PI;
 }
 
 bool
-parameter_loader::validate_materials(const std::string& mpn, int tag)
+parameter_loader::validate_materials(const std::string& mpn, int tag) const
 {
-    auto mfun = lua["materials"][mpn];
+    auto mfun = lua[SEC_MATERIALS][mpn];
     if ( mfun.valid() )
         return true;
 
-    auto mparams = lua["materials"][tag];
+    auto mparams = lua[SEC_MATERIALS][tag];
     if ( mparams.valid() )
     {
-        auto matparams_mpn = lua["materials"][tag][mpn];
+        auto matparams_mpn = lua[SEC_MATERIALS][tag][mpn];
         if ( matparams_mpn.valid() )
             return true;
     }
 
-    std::cout << "[CONFIG] 'materials[" << tag << "]." << mpn << "' not defined and ";
-    std::cout << "'materials." << mpn << "(tag,x,y,z)' not present." << std::endl;
+    std::cout << "[CONFIG] '"SEC_MATERIALS".[" << tag << "]." << mpn;
+    std::cout << "' not defined and '"SEC_MATERIALS"." << mpn;
+    std::cout << "(tag,x,y,z)' not present." << std::endl;
+    return false;
+}
+
+bool
+parameter_loader::validate_boundary_condition(const model& mod, int tag) const
+{
+    auto bc = lua[SEC_BCONDS][tag];
+    if (not bc.valid())
+        return true;
+
+    if ( not bc[BC_KIND].valid() )
+    {
+        std::cout << "[CONFIG] '"BC_KIND"' not specified on interface ";
+        std::cout << tag << std::endl;
+        return false;
+    }
+
+    std::string bc_kind = bc[BC_KIND];
+
+    if (bc_kind == BCOND_TYPE_PEC)
+        return true;
+
+    if (bc_kind == BCOND_TYPE_PMC)
+        return true;
+
+    if (bc_kind == BCOND_TYPE_IMPEDANCE)
+        return true;
+
+    if (bc_kind == BCOND_TYPE_EPLANEW)
+    {
+        if (bc[BC_SOURCE].valid())
+            return true;
+        
+        std::cout << "[CONFIG] '"BC_SOURCE"' not specified for plane ";
+        std::cout << "wave condition on surface " << tag << std::endl;
+        return false;
+    }
+
+    std::cout << "[CONFIG] boundary condition not implemented on ";
+    std::cout << "surface " << tag << std::endl;
+    return false;
+}
+
+bool
+parameter_loader::validate_interface_condition(const model& mod, int tag) const
+{
+    auto ic = lua[SEC_IFCONDS][tag];
+    if (not ic.valid())
+        return true;
+
+    if ( not ic[IFC_KIND].valid() )
+    {
+        std::cout << "[CONFIG '"IFC_KIND"' not specified on interface ";
+        std::cout << tag << std::endl;
+        return false;
+    }
+
+    std::string ic_kind = ic[IFC_KIND];
+
+    if (ic_kind == IFCOND_TYPE_SURFCURR)
+    {
+        if (ic[IFC_SOURCE].valid())
+            return true;
+        
+        std::cout << "[CONFIG] '"IFC_SOURCE"' not specified for surface ";
+        std::cout << "current condition on surface " << tag << std::endl;
+        return false;
+    }
+
+    std::cout << "[CONFIG] interface condition not implemented on ";
+    std::cout << "surface " << tag << std::endl;
     return false;
 }
 
 bool
-parameter_loader::validate_model_params(const model& mod)
+parameter_loader::validate_conditions(const model& mod) const
+{
+    bool success = true;
+    auto bds = mod.boundary_descriptors();
+    std::sort(bds.begin(), bds.end());
+    bds.erase( std::unique(bds.begin(), bds.end()), bds.end() );
+    
+    auto ft_none = [](const boundary_descriptor& bd) -> bool {
+        return bd.type == face_type::NONE;
+    };
+
+    bds.erase( std::remove_if(bds.begin(), bds.end(), ft_none), bds.end() );
+
+    for (auto& bd : bds)
+    {
+        switch (bd.type)
+        {
+            case face_type::NONE:
+                assert(false);
+                break;
+
+            case face_type::BOUNDARY:
+                if (not validate_boundary_condition(mod, bd.gmsh_entity) )
+                    success = false;
+                break;
+
+            case face_type::INTERFACE:
+                if (not validate_interface_condition(mod, bd.gmsh_entity) )
+                    success = false;
+                break;
+        }
+    }
+
+    return success;
+}
+
+bool
+parameter_loader::validate_model_params(const model& mod) const
 {
     bool success = true;
     for (auto& e : mod)
     {
         auto tag = e.gmsh_tag();
-        bool eps_valid = validate_materials("epsilon", tag);
-        bool mu_valid = validate_materials("mu", tag);
-        bool sigma_valid = validate_materials("sigma", tag);
+        bool eps_valid = validate_materials(MAT_EPS, tag);
+        bool mu_valid = validate_materials(MAT_MU, tag);
+        bool sigma_valid = validate_materials(MAT_SIGMA, tag);
         if ( (not eps_valid) or (not mu_valid) or (not sigma_valid) )
             success = false;
     }
 
+    if ( not validate_conditions(mod) )
+        success = false;
+
     return success;
 }
 
@@ -215,14 +354,14 @@ double
 parameter_loader::epsilon(int tag, const point_3d& pt) const
 {
     double eps0 = lua["const"]["eps0"];
-    auto eps_dom = lua["materials"][tag]["epsilon"];
+    auto eps_dom = lua[SEC_MATERIALS][tag][MAT_EPS];
     if (eps_dom.valid())
     {
         double ret = eps_dom;
         return eps0*ret;
     }
 
-    double ret = lua["materials"]["epsilon"](tag, pt.x(), pt.y(), pt.z());
+    double ret = lua[SEC_MATERIALS][MAT_EPS](tag, pt.x(), pt.y(), pt.z());
     return eps0*ret;
 }
 
@@ -230,101 +369,83 @@ double
 parameter_loader::mu(int tag, const point_3d& pt) const
 {
     double mu0 = lua["const"]["mu0"];
-    auto mu_dom = lua["materials"][tag]["mu"];
+    auto mu_dom = lua[SEC_MATERIALS][tag][MAT_MU];
     if (mu_dom.valid())
     {
         double ret = mu_dom;
         return mu0*ret;
     }
 
-    double ret = lua["materials"]["mu"](tag, pt.x(), pt.y(), pt.z());
+    double ret = lua[SEC_MATERIALS][MAT_MU](tag, pt.x(), pt.y(), pt.z());
     return mu0*ret;
 }
 
 double
 parameter_loader::sigma(int tag, const point_3d& pt) const
 {
-    auto sigma_dom = lua["materials"][tag]["sigma"];
+    auto sigma_dom = lua[SEC_MATERIALS][tag][MAT_SIGMA];
     if (sigma_dom.valid())
         return sigma_dom;
 
-    return lua["materials"]["sigma"](tag, pt.x(), pt.y(), pt.z());
+    return lua[SEC_MATERIALS][MAT_SIGMA](tag, pt.x(), pt.y(), pt.z());
 }
 
 boundary_condition
 parameter_loader::boundary_type(int tag) const
 {
-    auto bnd_data = lua["bndconds"][tag];
+    auto bnd_data = lua[SEC_BCONDS][tag];
     if (not bnd_data.valid())
         return boundary_condition::UNSPECIFIED;
-    
-    auto kind = bnd_data["kind"];
-    if (not kind.valid())
-    {
-        std::cout << "[CONFIG] warning: 'kind' not specified on boundary ";
-        std::cout << tag << std::endl;
-        return boundary_condition::UNSPECIFIED;
-    }
 
-    std::string kind_str = kind;
+    std::string kind_str = bnd_data[BC_KIND];
 
-    if (kind_str == "pec")
+    if (kind_str == BCOND_TYPE_PEC)
         return boundary_condition::PEC;
 
-    if (kind_str == "pmc")
+    if (kind_str == BCOND_TYPE_PMC)
         return boundary_condition::PMC;
 
-    if (kind_str == "impedance")
+    if (kind_str == BCOND_TYPE_IMPEDANCE)
         return boundary_condition::IMPEDANCE;
 
-    if (kind_str == "plane_wave_E")
+    if (kind_str == BCOND_TYPE_EPLANEW)
         return boundary_condition::PLANE_WAVE_E;
 
-    if (kind_str == "plane_wave_H")
+    if (kind_str == BCOND_TYPE_HPLANEW)
         return boundary_condition::PLANE_WAVE_H;
 
-    if (kind_str == "E_field")
+    if (kind_str == BCOND_TYPE_EFIELD)
         return boundary_condition::E_FIELD;
     
-    if (kind_str == "H_field")
+    if (kind_str == BCOND_TYPE_HFIELD)
         return boundary_condition::H_FIELD;
 
-    if (kind_str == "surface_current")
+    if (kind_str == BCOND_TYPE_SURFCURR)
         return boundary_condition::SURFACE_CURRENT;
 
-    std::cout << "[CONFIG] warning: 'kind' invalid on boundary ";
-    std::cout << tag << std::endl;
     return boundary_condition::UNSPECIFIED;
 }
 
 interface_condition
 parameter_loader::interface_type(int tag) const
 {
-    auto bnd_data = lua["ifaceconds"][tag];
-    if (not bnd_data.valid())
+    auto if_data = lua[SEC_IFCONDS][tag];
+    if (not if_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;
-    }
+    auto kind = if_data[IFC_KIND];
 
     std::string kind_str = kind;
 
-    if (kind_str == "E_field")
+    if (kind_str == IFCOND_TYPE_EFIELD)
         return interface_condition::E_FIELD;
     
-    if (kind_str == "H_field")
+    if (kind_str == IFCOND_TYPE_HFIELD)
         return interface_condition::H_FIELD;
 
-    if (kind_str == "surface_current")
+    if (kind_str == IFCOND_TYPE_SURFCURR)
         return interface_condition::SURFACE_CURRENT;
 
-    std::cout << "[CONFIG] warning: 'kind' invalid on boundary ";
-    std::cout << tag << std::endl;
     return interface_condition::UNSPECIFIED;
 }
 
@@ -333,7 +454,7 @@ parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) co
 {
     vec3d ret;
     sol::tie(ret(0), ret(1), ret(2)) =
-        lua["bndconds"][tag]["source"](tag, pt.x(), pt.y(), pt.z(), t);
+        lua[SEC_BCONDS][tag][BC_SOURCE](tag, pt.x(), pt.y(), pt.z(), t);
     return ret;
 }
 
@@ -342,8 +463,8 @@ parameter_loader::eval_interface_source(int tag, const point_3d& pt, double t) c
 {
     vec3d ret;
     sol::tie(ret(0), ret(1), ret(2)) =
-        lua["ifaceconds"][tag]["source"](tag, pt.x(), pt.y(), pt.z(), t);
+        lua[SEC_IFCONDS][tag][IFC_SOURCE](tag, pt.x(), pt.y(), pt.z(), t);
     return ret;
 }
 
-} // namespace maxwell
\ No newline at end of file
+} // namespace maxwell

src/physical_element.cpp

@@ -114,10 +114,6 @@ physical_element::barycenter(void) const
     return m_barycenter;
 }
 
-
-
-
-
 physical_elements_factory::physical_elements_factory(const entity_params& ep)
     : dim(ep.dim), tag(ep.tag), elemType(ep.etype), geom_order(ep.gorder),
       approx_order(ep.aorder)
@@ -131,8 +127,15 @@ physical_elements_factory::physical_elements_factory(const entity_params& ep)
     gmm::getBarycenters(elemType, tag, false, false, barycenters);
 
     std::vector<double>     coord;
+
+#ifdef USE_INITIAL_4_8_4_API
     gmm::getKeysForElements(elemType, basis_func_name(approx_order),
         keypairs, coord, tag, false);
+#else
+    gmm::getKeysForElements(elemType, basis_func_name(approx_order),
+        tagKeys, entityKeys, coord, tag, false);
+    assert(tk.size() == ek.size());
+#endif
 
     keys_per_elem = gmm::getNumberOfKeysForElements(elemType, basis_func_name(approx_order));
     assert(keys_per_elem*cellTags.size() == keypairs.size());
@@ -170,7 +173,12 @@ physical_elements_factory::get_elements()
         new_pe.m_bf_keys.resize(keys_per_elem);
         for (size_t i = 0; i < keys_per_elem; i++)
         {
+#ifdef USE_INITIAL_4_8_4_API
             auto [vi, vu] = keypairs[keys_per_elem*elem + i]; 
+#else
+            auto vi = tagKeys[keys_per_elem*elem + i];
+            auto vu = entityKeys[keys_per_elem*elem + i];
+#endif
             new_pe.m_bf_keys[i] = bf_key(vi,vu);
         }