A lot of cleanup, corrected DirichletBC, removed boost

dune/amdis/AMDiS.cpp

 #include "AMDiS.hpp"
 
-// std c++ headers
-#include <string>
-
-#include <boost/program_options.hpp>
-
 // AMDiS includes
 #include <dune/amdis/Initfile.hpp>
 #include <dune/amdis/Output.hpp>
 
 namespace AMDiS
 {
-//   using namespace std;
-
   Dune::MPIHelper& init(int& argc, char**& argv, std::string initFileName)
   {
     // Maybe initialize MPI
@@ -20,73 +13,19 @@ namespace AMDiS
 
     Parameters::clearData();
 
-    // read commandline arguments
-    namespace po = boost::program_options;
-
-    // Declare the supported options.
-    po::options_description desc("Usage: " + std::string(argv[0]) + " init-file [options]\nAllowed options");
-    desc.add_options()
-    ("help", "produce help message")
-    ("init-file", po::value<std::string>(), "set init file")
-    ("parameters", po::value<std::string>(), "set parameter in init file\nsyntax: \"key1: value1; key2: value2...\"");
-
-    po::options_description hidden("Hidden options");
-    hidden.add_options()
-    ("unknown", po::value<std::vector<std::string>>(), "unknown options");
-
-    po::options_description cmdline_options;
-    cmdline_options.add(desc).add(hidden);
-
-    // first argument is init-filename
-    po::positional_options_description p;
-    p.add("init-file", 1);
-    p.add("unknown", -1);
-
-    // parse comandline
-    po::variables_map vm;
-    po::store(po::command_line_parser(argc, argv).options(cmdline_options).positional(p).allow_unregistered().run(), vm);
-    po::notify(vm);
-
-    // print help message
-    if (vm.count("help"))
-    {
-      std::cout << desc << "\n";
-      exit(1);
-    }
-
-    // set parameters before reading the initfile
-    // if (vm.count("parameters"))
-    //   Parameters::readArgv(vm["parameters"].as<std::string>());
-
-    if (initFileName == "")
-    {
-      if (vm.count("init-file"))
-        Parameters::init(vm["init-file"].as<std::string>());
+    if (initFileName == "") {
+      if (argc > 1)
+        Parameters::init(argv[1]);
       else
         error_exit("No init file specified!\n");
-    }
-    else
-    {
+    } else {
       Parameters::init(initFileName);
     }
 
-    // reset parameters from command line
-    bool ignoreCommandline = false;
-    Parameters::get("ignore commandline options", ignoreCommandline);
-    // if (vm.count("parameters") && !ignoreCommandline)
-    //   Parameters::readArgv(vm["parameters"].as<std::string>(),0);
-    // TODO: add command-line arguments again.
-
     return mpiHelper;
   }
 
 
-  void init(std::string initFileName)
-  {
-    Parameters::init(initFileName);
-  }
-
-
   void finalize()
   {}
 

dune/amdis/AMDiS.hpp

@@ -3,6 +3,7 @@
 #ifdef HAVE_CONFIG_H
 # include "config.h"
 #endif
+
 // std c++ headers
 #include <string>
 
@@ -13,8 +14,6 @@ namespace AMDiS
 {
   Dune::MPIHelper& init(int& argc, char**& argv, std::string initFileName = "");
 
-  void init(std::string initFileName);
-
   void finalize();
 
 } // end namespace AMDiS

dune/amdis/AdaptBase.hpp

@@ -9,7 +9,7 @@ namespace AMDiS
   class AdaptInfo;
   class ProblemIterationInterface;
   class ProblemTimeInterface;
-  
+
   /// Interface for adaption loops.
   class AdaptBase
   {
@@ -18,8 +18,8 @@ namespace AMDiS
     AdaptBase(std::string sname,
               ProblemIterationInterface* problemIteration,
               AdaptInfo& adapt,
-              ProblemTimeInterface* problemTime = NULL,
-              AdaptInfo* initialAdaptInfo = NULL)
+              ProblemTimeInterface* problemTime = nullptr,
+              AdaptInfo* initialAdaptInfo = nullptr)
       : name(sname)
       , problemIteration(problemIteration)
       , adaptInfo(adapt)

dune/amdis/AdaptInfo.cpp

@@ -27,13 +27,13 @@ namespace AMDiS
   }
 
 
-  AdaptInfo::AdaptInfo(std::string name_, int size)
-    : name(name_),
-      scalContents(size)
+  AdaptInfo::AdaptInfo(std::string name, int size)
+    : name(name)
   {
     // init();
     Parameters::get(name + "->start time", startTime);
     time = startTime;
+
     Parameters::get(name + "->timestep", timestep);
     Parameters::get(name + "->end time", endTime);
     Parameters::get(name + "->max iteration", maxSpaceIteration);
@@ -45,23 +45,17 @@ namespace AMDiS
     Parameters::get(name + "->time tolerance", globalTimeTolerance);
 
     for (int i = 0; i < size; i++)
-    {
-      scalContents[i] = new ScalContent(name + "[" + std::to_string(i) + "]");
-    }
+      scalContents.emplace_back(new ScalContent(name + "[" + std::to_string(i) + "]"));
   }
 
 
   void AdaptInfo::setScalContents(int newSize)
   {
-    int oldSize = static_cast<int>(scalContents.size());
-
-    if (newSize > oldSize)
-    {
-      scalContents.resize(newSize);
+    int oldSize = int(scalContents.size());
 
-      for (int i = oldSize; i < newSize; i++)
-        scalContents[i] =
-          new ScalContent(name + "[" + std::to_string(i) + "]");
+    if (newSize > oldSize) {
+      for (int i = oldSize; i < newSize; ++i)
+        scalContents.emplace_back(new ScalContent(name + "[" + std::to_string(i) + "]"));
     }
   }
 

dune/amdis/AdaptInfo.hpp

@@ -82,14 +82,10 @@ namespace AMDiS
 
   public:
     /// Constructor.
-    AdaptInfo(std::string name_, int size = 1);
+    AdaptInfo(std::string name, int size = 1);
 
     /// Destructor.
-    virtual ~AdaptInfo()
-    {
-      for (std::size_t i = 0;  i < scalContents.size(); i++)
-        delete scalContents[i];
-    }
+    virtual ~AdaptInfo() {}
 
     /// Resets all variables to zero (or something equivalent)
     void reset();
@@ -304,8 +300,7 @@ namespace AMDiS
     double getEstSum(int index) const
     {
       AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
-      test_exit_dbg(static_cast<std::size_t>(index) < scalContents.size(),
-                    "Wrong index for adaptInfo!\n");
+      test_exit_dbg(size_t(index) < scalContents.size(), "Wrong index for adaptInfo!\n");
 
       return scalContents[index]->est_sum;
     }
@@ -320,8 +315,7 @@ namespace AMDiS
     double getEstMax(int index) const
     {
       AMDIS_FUNCNAME_DBG("AdaptInfo::getEstSum()");
-      test_exit_dbg(static_cast<std::size_t>(index) < scalContents.size(),
-                    "Wrong index for adaptInfo!\n");
+      test_exit_dbg(size_t(index) < scalContents.size(), "Wrong index for adaptInfo!\n");
 
       return scalContents[index]->est_max;
     }
@@ -535,7 +529,7 @@ namespace AMDiS
 
     int getSize() const
     {
-      return static_cast<int>(scalContents.size());
+      return int(scalContents.size());
     }
 
     // TODO: remove from AdaptInfo
@@ -690,7 +684,7 @@ namespace AMDiS
     double globalTimeTolerance = 1.0;
 
     /// Scalar adapt infos.
-    std::vector<ScalContent*> scalContents;
+    std::vector<std::unique_ptr<ScalContent>> scalContents;
 
     /// Is true, if the adaptive procedure was deserialized from a file. TODO: remove deserialization
     bool deserialized = false;

dune/amdis/AdaptInstationary.cpp

@@ -15,13 +15,8 @@ namespace AMDiS
                                        AdaptInfo& adaptInfo,
                                        ProblemTimeInterface& problemInstat,
                                        AdaptInfo& initialInfo)
-    : AdaptBase(name, &problemStat, adaptInfo, &problemInstat, &initialInfo),
-      breakWhenStable(0)
+    : AdaptBase(name, &problemStat, adaptInfo, &problemInstat, &initialInfo)
   {
-    strategy = 0;
-    timeDelta1 = 0.7071;
-    timeDelta2 = 1.4142;
-
     Parameters::get(name + "->strategy", strategy);
     Parameters::get(name + "->time delta 1", timeDelta1);
     Parameters::get(name + "->time delta 2", timeDelta2);

dune/amdis/AdaptInstationary.hpp

@@ -19,7 +19,8 @@ namespace AMDiS
    * problems (see ProblemInstat). It contains a pointer to a ProblemInstat
    * object.
    */
-  class AdaptInstationary : public AdaptBase
+  class AdaptInstationary
+      : public AdaptBase
   {
   public:
     /// Creates a AdaptInstationary object with the given name for the time
@@ -75,19 +76,19 @@ namespace AMDiS
 
   protected:
     /// Strategy for choosing one timestep
-    int strategy;
+    int strategy = 0;
 
     /// Parameter \f$ \delta_1 \f$ used in time step reduction
-    double timeDelta1;
+    double timeDelta1 = 0.7071;
 
     /// Parameter \f$ \delta_2 \f$ used in time step enlargement
-    double timeDelta2;
+    double timeDelta2 = 1.4142;
 
     /// If this parameter is 1 and the instationary problem is stable, hence the number
     /// of solver iterations to solve the problem is zero, the adaption loop will stop.
-    int breakWhenStable;
+    bool breakWhenStable = false;
 
-    ///
+    /// min-timestep == max-timestep
     bool fixedTimestep;
   };
 

dune/amdis/AdaptStationary.cpp

@@ -6,46 +6,45 @@
 #include "Initfile.hpp"
 #include "ProblemIterationInterface.hpp"
 
-namespace AMDiS
+namespace AMDiS {
+
+AdaptStationary::AdaptStationary(std::string name,
+                                 ProblemIterationInterface& prob,
+                                 AdaptInfo& adaptInfo)
+  : AdaptBase(name, &prob, adaptInfo)
 {
+  Parameters::get(name + "->info", info);
+}
+
 
-  AdaptStationary::AdaptStationary(std::string name,
-                                   ProblemIterationInterface& prob,
-                                   AdaptInfo& adaptInfo)
-    : AdaptBase(name, &prob, adaptInfo)
+int AdaptStationary::adapt()
+{
+  // initial iteration
+  if (adaptInfo.getSpaceIteration() == -1)
   {
-    Parameters::get(name + "->info", info);
+    problemIteration->beginIteration(adaptInfo);
+    problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
+    problemIteration->endIteration(adaptInfo);
+    adaptInfo.incSpaceIteration();
   }
 
-
-  int AdaptStationary::adapt()
+  // adaption loop
+  while (!adaptInfo.spaceToleranceReached() &&
+          (adaptInfo.getSpaceIteration() < adaptInfo.getMaxSpaceIteration() ||
+          adaptInfo.getMaxSpaceIteration() < 0) )
   {
-    // initial iteration
-    if (adaptInfo.getSpaceIteration() == -1)
-    {
-      problemIteration->beginIteration(adaptInfo);
-      problemIteration->oneIteration(adaptInfo, NO_ADAPTION);
-      problemIteration->endIteration(adaptInfo);
-      adaptInfo.incSpaceIteration();
-    }
-
-    // adaption loop
-    while (!adaptInfo.spaceToleranceReached() &&
-           (adaptInfo.getSpaceIteration() < adaptInfo.getMaxSpaceIteration() ||
-            adaptInfo.getMaxSpaceIteration() < 0) )
-    {
-
-      problemIteration->beginIteration(adaptInfo);
-      Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
-      problemIteration->endIteration(adaptInfo);
-
-      if (adapted == Flag{0})
-        break;
-
-      adaptInfo.incSpaceIteration();
-    }
-
-    return 0;
+
+    problemIteration->beginIteration(adaptInfo);
+    Flag adapted = problemIteration->oneIteration(adaptInfo, FULL_ITERATION);
+    problemIteration->endIteration(adaptInfo);
+
+    if (adapted == Flag{0})
+      break;
+
+    adaptInfo.incSpaceIteration();
   }
 
+  return 0;
+}
+
 } // end namespace AMDiS

dune/amdis/AdaptStationary.hpp

-/** \defgroup Adaption Adaption module
- * @{ <img src="adaption.png"> @}
- *
- * \brief
- * Contains all classes needed for adaption.
- */
-
 #pragma once
 
 // std c++ headers
@@ -19,12 +12,21 @@ namespace AMDiS
   class AdaptInfo;
   class ProblemIterationInterface;
 
+  /** \defgroup Adaption Adaption module
+  * @{ <img src="adaption.png"> @}
+  *
+  * \brief
+  * Contains all classes needed for space and time adaption.
+  */
+
+
   /** \ingroup Adaption
    * \brief
    * AdaptStationary contains information about the adaptive procedure and the
    * adapt procedure itself
    */
-  class AdaptStationary : public AdaptBase
+  class AdaptStationary
+      : public AdaptBase
   {
   public:
     /// Creates a AdaptStationary object with given name.

dune/amdis/Assembler.hpp

 #pragma once
 
+// std c++ headers
 #include <memory>
+
+// AMDiS includes
 #include <dune/amdis/QuadratureRules.hpp>
 #include <dune/amdis/common/TypeDefs.hpp>
 
@@ -13,6 +16,7 @@ namespace AMDiS
     static constexpr int dim = GridView::dimension;
     using ctype = typename GridView::ctype;
 
+    using LocalQuadratureRules = Dune::QuadratureRules<ctype, LocalContext::mydimension>;
     using QuadratureRule = QuadratureRuleFactory_t<LocalContext, ctype, dim>;
     using Geometry = typename Impl::Get<LocalContext>::Geometry;
 
@@ -24,7 +28,7 @@ namespace AMDiS
       : geometry(get_geometry(element))
     {
       int order = op.getQuadratureDegree(geometry.type(), geometry, degree, type);
-      auto const& quad = Dune::QuadratureRules<ctype, LocalContext::mydimension>::rule(geometry.type(), order);
+      auto const& quad = LocalQuadratureRules::rule(geometry.type(), order);
 
       quadrature.reset(new QuadratureRule(element, quad));
     }
@@ -36,7 +40,7 @@ namespace AMDiS
      * elements this is a wrapper around the classical element quadrature rule.
      * Access the underlying dune-quadrature rule, with `quadrature->getRule()`.
      **/
-    auto const& getQuadrature() const
+    QuadratureRule const& getQuadrature() const
     {
       return *quadrature;
     }

dune/amdis/Concepts.hpp

-#pragma once
-
-#include <dune/common/concept.hh>
-
-namespace AMDiS
-{
-  namespace Concepts
-  {
-    namespace Definition
-    {
-      // Concept for a Dune::Entity
-      struct Entity
-      {
-        template<class E>
-        auto require(const E& e) -> decltype(
-          e.level(),
-          e.partitionType(),
-          e.geometry(),
-          e.type(),
-          e.subEntities((unsigned int)(0)),
-          e.seed()
-        );
-      };
-
-      // Concept for a Dune::Intersection
-      struct Intersection
-      {
-        template<class I>
-        auto require(const I& i) -> decltype(
-          i.boundary(),
-          i.neighbor(),
-          i.inside(),
-          i.outside(),
-          i.geometryInInside(),
-          i.geometryInOutside(),
-          i.geometry(),
-          i.type(),
-          i.indexInInside(),
-          i.indexInOutside()
-        );
-      };
-
-    } // end namespace Definition
-
-    template <class E>
-    constexpr bool Entity = Dune::models<Definition::Entity, E>();
-
-    template <class I>
-    constexpr bool Intersection = Dune::models<Definition::Intersection, I>();
-
-  } // end namespace Concepts
-
-} // end namespace AMDiS

dune/amdis/CreatorInterface.hpp

@@ -8,7 +8,6 @@
 
 namespace AMDiS
 {
-
   /** \ingroup Common
    * \brief
    * Interface for the implementation of the factory method pattern.
@@ -41,11 +40,12 @@ namespace AMDiS
    * Interface for creators with name.
    */
   template <class BaseClass>
-  class CreatorInterfaceName : public CreatorInterface<BaseClass>
+  class CreatorInterfaceName
+      : public CreatorInterface<BaseClass>
   {
   public:
 
-    virtual shared_ptr<BaseClass> create() override
+    virtual std::shared_ptr<BaseClass> create() final
     {
       error_exit("Should not be called. Call create(string) instead!");
       return {};
@@ -56,7 +56,7 @@ namespace AMDiS
      * Creates a new instance of the sub class of BaseClass by passing a
      * string to the constructor.
      */
-    virtual shared_ptr<BaseClass> create(std::string) = 0;
+    virtual std::shared_ptr<BaseClass> create(std::string) = 0;
   };
 
 

dune/amdis/CreatorMap.hpp

@@ -35,8 +35,7 @@ namespace AMDiS
     static void addCreator(std::string key, CreatorInterface<BaseClass>* creator)
     {
       init();
-      test_exit(!creatorMap[key],
-                "There is already a creator for key ", key);
+      test_exit(!creatorMap[key], "There is already a creator for key ", key);
       creatorMap[key] = creator;
     }
 
@@ -50,8 +49,7 @@ namespace AMDiS
         key = "default";
 
       auto creator = creatorMap[key];
-      test_exit(creator,
-                "No creator for key \"", key, "\" defined in init file for parameter \"", initFileStr, "\"");
+      test_exit(creator, "No creator for key `", key, "` defined in init file for parameter `", initFileStr, "`");
 
       return creator;
     }

dune/amdis/DirichletBC.hpp

@@ -14,7 +14,7 @@ namespace AMDiS
   /**
    * By calling the methods \ref init() and \ref finish before and after
    * assembling the system-matrix, respectively, dirichlet boundary conditions
-   * can be applied to the matrix and system vector. Therefore a predicate
+   * can be applied to the matrix and system vector. Therefore, a predicate
    * functions indicates the DOFs where values should be enforced and a second
    * functor provided in the constructor is responsible for determining the
    * values to be set at the DOFs.
@@ -38,6 +38,9 @@ namespace AMDiS
     {}
 
 
+    /// Prepare the matrix, rhs and solution block for assembling dirichlet
+    /// boundary conditions, e.g. collect the corresponding indices of boundary
+    /// DOFS for later elimination in \ref finish.
     template <class Matrix, class VectorX, class VectorB>
     void init(bool apply,
               Matrix& matrix,
@@ -45,8 +48,16 @@ namespace AMDiS
               VectorB& solution);
 
 
+    /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row
+    /// to the matrix and optionally delete the corresponding matrix-column.
+    /** If DBC is set for matrix block {R,C}, then \p apply_row means that
+      * the block-row R is assembled, and \p apply_col means that the block-col C
+      * is assembled. The \p matrix, \p rhs, and \p solution correspond to the
+      * currently visited blocks of system-matrix and system-vector.
+      **/
     template <class Matrix, class VectorX, class VectorB>
-    void finish(bool apply,
+    void finish(bool apply_row,
+                bool apply_col,
                 Matrix& matrix,
                 VectorX& rhs,
                 VectorB& solution);

dune/amdis/DirichletBC.inc.hpp

@@ -7,10 +7,10 @@ namespace AMDiS
 {
   template <class WorldVector>
     template <class Matrix, class VectorX, class VectorB>
-  void DirichletBC<WorldVector>::init(bool apply,
+  void DirichletBC<WorldVector>::init(bool /*apply*/,
 				                              Matrix& matrix,
-                                      VectorX& solution,
-				                              VectorB& rhs)
+                                      VectorX& /*solution*/,
+				                              VectorB& /*rhs*/)
   {
     using Dune::Functions::interpolate;
 
@@ -23,7 +23,8 @@ namespace AMDiS
 
   template <class WorldVector>
     template <class Matrix, class VectorX, class VectorB>
-  void DirichletBC<WorldVector>::finish(bool apply,
+  void DirichletBC<WorldVector>::finish(bool apply_row,
+                                        bool apply_col,
 					                              Matrix& matrix,
                                         VectorX& solution,
 					                              VectorB& rhs)
@@ -31,20 +32,16 @@ namespace AMDiS
     using Dune::Functions::interpolate;
 
     test_exit_dbg(initialized, "Boundary condition not initialized!");
-    auto columns = matrix.applyDirichletBC(dirichletNodes, apply);
+    auto columns = matrix.applyDirichletBC(dirichletNodes, apply_row, apply_col);
 
-
-    if (apply) {
+    if (apply_col) {
       interpolate(matrix.getRowFeSpace(), wrapper(rhs.getVector()), values, dirichletNodes);
-      interpolate(matrix.getColFeSpace(), wrapper(solution.getVector()), values, dirichletNodes);
+      if (apply_row)
+        interpolate(matrix.getColFeSpace(), wrapper(solution.getVector()), values, dirichletNodes);
 
       // subtract columns of dirichlet nodes from rhs
-      for (size_t i = 0; i < columns.size(); ++i) {
-        for (auto const& aij : columns[i]) {
-          rhs[i] -= aij.second * solution[aij.first];
-        }
-      }
-