Verify material specification

* Add 'debug' key in material specification.
* Add general check in parameter lookup, avoiding segfaults.
* Add test for material lookup and sanitizer

CMakeLists.txt

@@ -31,6 +31,7 @@ add_subdirectory(cmake/modules)
 add_subdirectory(dune)
 add_subdirectory(python)
 add_subdirectory(experiments)
+add_subdirectory(test)
 
 # finalize the dune project, e.g. generating config.h etc.
 finalize_dune_project(GENERATE_CONFIG_H_CMAKE)

dune/structures/material.hh

@@ -3,6 +3,7 @@
 
 #include <dune/blocklab/blocks/blockbase.hh>
 #include <dune/blocklab/utilities/stringsplit.hh>
+#include <dune/common/exceptions.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/parametertree.hh>
@@ -216,8 +217,15 @@ private:
     while (lookup_entity.hasFather())
       lookup_entity = lookup_entity.father();
 
-    return materials.find((*physical_entity_mapping)[is->index(lookup_entity)])
-      ->second;
+    const auto it =
+      materials.find((*physical_entity_mapping)[is->index(lookup_entity)]);
+    if (it == materials.end())
+    {
+      DUNE_THROW(Dune::InvalidStateException,
+                 "Grid mapping index not found in defined materials!");
+    }
+
+    return it->second;
   }
 
   const typename GV::IndexSet* is;
@@ -244,6 +252,47 @@ parse_material(const GV& gv,
   return coll;
 }
 
+/// Check if material definitions and mappings match
+/**
+ *  Iterates over all grid elements and checks if parameters for the given
+ *  physical mapping index can be retrieved.
+ *
+ *  \param material Shared pointer to ElasticMaterialBase
+ *  \throw IOError if a physical group on the grid is not given in the material
+ *                 specification
+ */
+template<typename Material>
+void
+sanitize_material(const std::shared_ptr<Material> material)
+{
+  // NOTE: GridView of material is actually a PartitionView
+  const auto gv = material->gridView();
+  const auto& index_set = gv.indexSet();
+  const auto groups = material->getPhysicalGroups();
+
+  // Iterate over grid elements
+  for (auto it = gv.template begin<0>(); it != gv.template end<0>(); ++it)
+  {
+    const auto& element = *it;
+    const auto index = index_set.index(element);
+    const auto center_local =
+      element.geometry().local(element.geometry().center());
+
+    // Try evaluating the parameterization
+    try
+    {
+      material->parameter(element, center_local, 0);
+    }
+    catch (Dune::InvalidStateException& e)
+    {
+      const auto group = groups->at(index);
+      DUNE_THROW(Dune::IOError,
+                 "Material group " + std::to_string(group)
+                   + " is missing in the material configuration!");
+    }
+  }
+}
+
 template<typename P, typename V>
 class MaterialInitializationBlock : public Dune::BlockLab::BlockBase<P, V>
 {
@@ -257,7 +306,7 @@ public:
   MaterialInitializationBlock(Context& ctx, const YAML::Node& config)
     : Dune::BlockLab::BlockBase<P, V>(ctx, config)
     , root_config(ctx.getRootConfig())
-    , material_config(config["materials"])
+    , material_config(config)
   {
   }
 
@@ -279,8 +328,14 @@ public:
       this->solver->template getVector<0>()->gridFunctionSpace().entitySet(),
       this->solver->template param<std::shared_ptr<std::vector<int>>>(
         "physical"),
-      material_config,
+      material_config["materials"],
       root_config);
+
+    if (material_config["debug"].as<bool>(false))
+    {
+      sanitize_material(material);
+    }
+
     this->solver->introduce_parameter("material",
                                       typename Traits::Parameter(material));
   }
@@ -362,7 +417,10 @@ public:
         "          meta:                                             \n"
         "            title: Prestress                                \n"
         "    meta:                                                   \n"
-        "      title: Materials                                      \n");
+        "      title: Materials                                      \n"
+        "  debug:                                                    \n"
+        "    type: boolean                                           \n"
+        "    default: false                                          \n");
     return data;
   }
 

test/CMakeLists.txt

+dune_add_test(SOURCES test-parameter-assign.cc)

test/test-parameter-assign.cc

+#include "config.h"
+
+#include <exception>
+#include <memory>
+#include <regex>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include <yaml-cpp/yaml.h>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/blocklab/construction/context.hh>
+#include <dune/blocklab/grids.hh>
+#include <dune/blocklab/init.hh>
+#include <dune/blocklab/vectors.hh>
+
+#include <dune/structures/material.hh>
+
+int
+main(int argc, char** argv)
+{
+  auto init = Dune::BlockLab::initBlockLab(argc, argv);
+  Dune::TestSuite test;
+
+  // Valid parameterization
+  std::string config_str = "solver:\n"
+                           "  grid:\n"
+                           "    refinement: 1\n"
+                           "    _blockname: grid_0\n"
+                           "  blocks:\n"
+                           "    material:\n"
+                           "      _blockname: material_0\n"
+                           "      debug: True\n"
+                           "      materials:\n"
+                           "        -\n"
+                           "          group: 0\n"
+                           "          model: linear\n"
+                           "          poisson_ratio: 0.3333\n"
+                           "          youngs_modulus: 300\n";
+  YAML::Node config;
+  config = YAML::Load(config_str);
+
+  // Grid and Vector
+  using Grid = Dune::UGGrid<2>;
+  using GridProvider = Dune::BlockLab::StructuredSimplexGridProvider<Grid>;
+  using VectorProvider = Dune::BlockLab::PkFemVectorProvider<GridProvider, 1>;
+  auto grid = std::make_shared<GridProvider>(config);
+  auto vector = std::make_shared<VectorProvider>(grid);
+
+  // Construction Context
+  using Material = std::shared_ptr<
+    ElasticMaterialBase<Dune::PDELab::OverlappingEntitySet<
+                          typename GridProvider::Grid::Traits::LeafGridView>,
+                        double>>;
+  using ParameterTuple =
+    std::tuple<std::shared_ptr<std::vector<int>>, Material>;
+
+  using Context =
+    Dune::BlockLab::ConstructionContext<ParameterTuple, VectorProvider>;
+
+  // Valid parameterization
+  try
+  {
+    Context ctx(init.helper, config, vector);
+    ctx.template registerBlock<MaterialInitializationBlock>("material");
+    auto solver = ctx.constructSolver(config["solver"]);
+    solver->apply();
+    test.check(true, "Valid Parameterization");
+  }
+  catch (std::exception& e)
+  {
+    test.check(false, "Valid Parameterization")
+      << "Valid parameterization could not be loaded: " << e.what();
+  }
+
+  // Replace group by "1", invalidating the parameterization
+  std::regex group_re("group: (\\d+)");
+  std::smatch group_match;
+  std::regex_search(config_str, group_match, group_re);
+  config_str.replace(group_match[1].first, group_match[1].second, "1", 1);
+  config = YAML::Load(config_str);
+
+  // Invalid parameterization
+  try
+  {
+    Context ctx(init.helper, config, vector);
+    ctx.template registerBlock<MaterialInitializationBlock>("material");
+    auto solver = ctx.constructSolver(config["solver"]);
+    solver->apply();
+    test.check(false, "Invalid Parameterization Group")
+      << "Invalid parameterization was loaded successfully!";
+  }
+  catch (Dune::IOError& e)
+  {
+    const std::string error_msg(e.what());
+    const auto pos = error_msg.find(
+      "Material group 0 is missing in the material configuration!");
+    test.check(pos != std::string::npos, "Invalid Parameterization Group")
+      << "Invalid parameterization led to wrong error message: " << e.what();
+  }
+  catch (...)
+  {
+    test.check(false, "Invalid Parameterization Group")
+      << "Invalid parameterization led to wrong error type!";
+  }
+
+  // Replace debug with "False"
+  std::regex debug_re("debug: (\\w+)");
+  std::smatch debug_match;
+  std::regex_search(config_str, debug_match, debug_re);
+  config_str.replace(debug_match[1].first, debug_match[1].second, "False");
+  config = YAML::Load(config_str);
+
+  // Invalid parameterization without debugging
+  try
+  {
+    Context ctx(init.helper, config, vector);
+    ctx.template registerBlock<MaterialInitializationBlock>("material");
+    auto solver = ctx.constructSolver(config["solver"]);
+    solver->apply();
+    test.check(true, "Invalid Parameterization Group, Debug off")
+      << "Invalid parameterization was not loaded successfully, although debug "
+         "was disabled!";
+  }
+  catch (...)
+  {
+    test.check(false, "Invalid Parameterization Group, Debug off");
+  }
+
+  return test.exit();
+}