[NL] extrapolator implementation in cpp file

* no property tags anymore
* no templates anymore in the implementation
* all template use is in the adaptor class

NumLib/Extrapolation/ExtrapolatableElement.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENT_H
+#define NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENT_H
+
+namespace NumLib
+{
+/*! Interface for providing shape matrices and integration point values for
+ *  extrapolation,
+ *
+ * Local assemblers that want to have some integration point values extrapolated
+ * using an Extrapolator have to implement this interface.
+ */
+class ExtrapolatableElement
+{
+public:
+    //! Provides the shape matrix at the given integration point.
+    virtual Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
+        const unsigned integration_point) const = 0;
+
+    virtual ~ExtrapolatableElement() = default;
+};
+
+}  // namespace NumLib
+
+#endif  // NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENT_H

NumLib/Extrapolation/ExtrapolatableElementCollection.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENTCOLLECTION_H
+#define NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENTCOLLECTION_H
+
+#include "ExtrapolatableElement.h"
+
+namespace NumLib
+{
+//! Adaptor to get information needed by an Extrapolator from an "arbitrary"
+//! collection of elements (e.g., local assemblers).
+class ExtrapolatableElementCollection
+{
+public:
+    //! Returns the shape matrix of the element with the given \c id at the
+    //! given \c integration_point.
+    virtual Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
+        std::size_t const id, unsigned const integration_point) const = 0;
+
+    //! Returns integration point values of some property of the element with
+    //! the given \c id.
+    virtual std::vector<double> const& getIntegrationPointValues(
+        std::size_t const id, std::vector<double>& cache) const = 0;
+
+    //! Returns the number of elements whose properties shall be extrapolated.
+    virtual std::size_t size() const = 0;
+
+    virtual ~ExtrapolatableElementCollection() = default;
+};
+
+template <typename LocalAssemblerCollection>
+class ExtrapolatableLocalAssemblerCollection
+    : public ExtrapolatableElementCollection
+{
+public:
+    //! LocalAssemblerCollection contains many LocalAssembler's.
+    using LocalAssembler = typename std::decay<decltype(
+        *std::declval<LocalAssemblerCollection>()[static_cast<std::size_t>(
+            0)])>::type;
+
+    static_assert(std::is_base_of<ExtrapolatableElement, LocalAssembler>::value,
+                  "Local assemblers used for extrapolation must be derived "
+                  "from Extrapolatable.");
+
+    /*! A method providing integration point values of some property.
+     *
+     * \remark
+     * \parblock
+     * The vector \c cache can be used to store some derived quantities that
+     * should be used as integration point values. They can be written to the
+     * \c cache and a reference to the \c cache can then be returned by the
+     * implementation of this method.
+     * Ordinary integration point values can be returned directly (if they are
+     * stored in the local assembler), and the \c cache can be left untouched.
+     *
+     * For usage examples see the processes already implemented or the unit test
+     * in TestExtrapolation.cpp
+     * \endparblock
+     *
+     * The method must return reference to a vector containing the integration
+     * point values.
+     */
+    using IntegrationPointValuesMethod = std::vector<double> const& (
+        LocalAssembler::*)(std::vector<double>& cache) const;
+
+    /*! Constructs a new instance.
+     *
+     * \param local_assemblers a collection of local assemblers whose
+     * integration point values shall be extrapolated.
+     * \param integration_point_values_method a LocalAssembler method returning
+     * integration point values of some property for that local assembler.
+     */
+    ExtrapolatableLocalAssemblerCollection(
+        LocalAssemblerCollection const& local_assemblers,
+        IntegrationPointValuesMethod integration_point_values_method)
+        : _local_assemblers(local_assemblers)
+        , _integration_point_values_method(integration_point_values_method)
+    {
+    }
+
+    Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
+        std::size_t const id, unsigned const integration_point) const override
+    {
+        ExtrapolatableElement const& loc_asm = *_local_assemblers[id];
+        return loc_asm.getShapeMatrix(integration_point);
+    }
+
+    std::vector<double> const& getIntegrationPointValues(
+        std::size_t const id, std::vector<double>& cache) const override
+    {
+        auto const& loc_asm = *_local_assemblers[id];
+        return (loc_asm.*_integration_point_values_method)(cache);
+    }
+
+    std::size_t size() const override { return _local_assemblers.size(); }
+private:
+    LocalAssemblerCollection const& _local_assemblers;
+    IntegrationPointValuesMethod const _integration_point_values_method;
+};
+
+//! Creates an ExtrapolatableLocalAssemblerCollection, which can be used to
+//! provide information to an Extrapolator.
+template <typename LocalAssemblerCollection,
+          typename IntegrationPointValuesMethod>
+ExtrapolatableLocalAssemblerCollection<LocalAssemblerCollection>
+makeExtrapolatable(LocalAssemblerCollection const& local_assemblers,
+                   IntegrationPointValuesMethod integration_point_values_method)
+{
+    return ExtrapolatableLocalAssemblerCollection<LocalAssemblerCollection>{
+        local_assemblers, integration_point_values_method};
+}
+}  // namespace NumLib
+
+#endif  // NUMLIB_EXTRAPOLATION_EXTRAPOLATABLEELEMENTCOLLECTION_H

NumLib/Extrapolation/Extrapolator.h

@@ -18,60 +18,14 @@
 namespace NumLib
 {
 
-/*! Interface for providing shape matrices and integration point values for
- *  extrapolation,
- *
- * Local assemblers that want to have some integration point values extrapolated
- * using Extrapolator have to implement this interface.
- *
- * \tparam PropertyTag  type of the property used to query a specific kind of
- *                      integration point value, usually an enum.
- */
-template <typename PropertyTag>
-class Extrapolatable
-{
-public:
-    //! Provides the shape matrix at the given integration point.
-    virtual Eigen::Map<const Eigen::RowVectorXd>
-    getShapeMatrix(const unsigned integration_point) const = 0;
-
-    /*! Provides integration point values for the given property.
-     *
-     * \remark
-     * The vector \c cache can be used to store some derived quantities that
-     * should be used as integration point values. They can be written to the
-     * \c cache and a reference to the \c cache can then be returned by the
-     * implementation of this method.
-     * Ordinary integration point values can be returned directly (if they are
-     * stored in the local assembler), and the \c cache cen be left untouched.
-     *
-     * \returns a reference to a vector containing the integration point values.
-     */
-    virtual std::vector<double> const&
-    getIntegrationPointValues(PropertyTag const property,
-                              std::vector<double>& cache) const = 0;
-
-    virtual ~Extrapolatable() = default;
-};
-
-/*! Interface for classes that extrapolate integration point values to nodal
- *  values.
- *
- * \tparam PropertyTag    type of the property used to query a specific kind of
- *                        integration point value, usually an enum.
- * \tparam LocalAssembler type of the local assembler being queried for
- *                        integration point values.
- */
-template<typename PropertyTag, typename LocalAssembler>
+//! Interface for classes that extrapolate integration point values to nodal
+//! values.
 class Extrapolator
 {
 public:
-    //! Vector of local assemblers, the same as is used in the FEM process.
-    using LocalAssemblers = std::vector<std::unique_ptr<LocalAssembler>>;
-
     //! Extrapolates the given \c property from the given local assemblers.
     virtual void extrapolate(
-            LocalAssemblers const& loc_asms, PropertyTag const property) = 0;
+            ExtrapolatableElementCollection const& extrapolatables) = 0;
 
     /*! Computes residuals from the extrapolation of the given \c property.
      *
@@ -80,7 +34,7 @@ public:
      * \pre extrapolate() must have been called before with the same arguments.
      */
     virtual void calculateResiduals(
-            LocalAssemblers const& loc_asms, PropertyTag const property) = 0;
+            ExtrapolatableElementCollection const& extrapolatables) = 0;
 
     //! Returns the extrapolated nodal values.
     //! \todo Maybe write directly to a MeshProperty.

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h → NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp

@@ -7,24 +7,23 @@
  *
  */
 
+#include "LocalLinearLeastSquaresExtrapolator.h"
+
 #include <functional>
 
-#include <logog/include/logog.hpp>
 #include <Eigen/Core>
+#include <logog/include/logog.hpp>
 
 #include "MathLib/LinAlg/LinAlg.h"
-#include "NumLib/Assembler/SerialExecutor.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "NumLib/Assembler/SerialExecutor.h"
 #include "NumLib/Function/Interpolation.h"
-#include "LocalLinearLeastSquaresExtrapolator.h"
+#include "ExtrapolatableElementCollection.h"
 
 namespace NumLib
 {
-
-template<typename PropertyTag, typename LocalAssembler>
-void
-LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
-extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
+inline void LocalLinearLeastSquaresExtrapolator::extrapolate(
+    ExtrapolatableElementCollection const& extrapolatables)
 {
     _nodal_values.setZero();
 
@@ -32,44 +31,39 @@ extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
     // compute the average afterwards
     auto counts =
         MathLib::MatrixVectorTraits<GlobalVector>::newInstance(_nodal_values);
-    counts->setZero(); // TODO BLAS?
-
-    using Self = LocalLinearLeastSquaresExtrapolator<
-        PropertyTag, LocalAssembler>;
+    counts->setZero();  // TODO BLAS?
 
-    NumLib::SerialExecutor::executeMemberDereferenced(
-        *this, &Self::extrapolateElement, local_assemblers, property, *counts);
+    auto const size = extrapolatables.size();
+    for (std::size_t i=0; i<size; ++i) {
+        extrapolateElement(i, extrapolatables, *counts);
+    }
 
     MathLib::LinAlg::componentwiseDivide(_nodal_values, _nodal_values, *counts);
 }
 
-template<typename PropertyTag, typename LocalAssembler>
-void
-LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
-calculateResiduals(LocalAssemblers const& local_assemblers,
-                   PropertyTag const property)
+inline void LocalLinearLeastSquaresExtrapolator::calculateResiduals(
+    ExtrapolatableElementCollection const& extrapolatables)
 {
-    assert(static_cast<std::size_t>(_residuals.size()) == local_assemblers.size());
-
-    using Self = LocalLinearLeastSquaresExtrapolator<
-        PropertyTag, LocalAssembler>;
+    assert(static_cast<std::size_t>(_residuals.size()) ==
+           extrapolatables.size());
 
-    NumLib::SerialExecutor::executeMemberDereferenced(
-        *this, &Self::calculateResiudalElement, local_assemblers, property);
+    auto const size = extrapolatables.size();
+    for (std::size_t i=0; i<size; ++i) {
+        calculateResiudalElement(i, extrapolatables);
+    }
 }
 
-template<typename PropertyTag, typename LocalAssembler>
-void
-LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
-extrapolateElement(std::size_t const element_index,
-                   LocalAssembler const& loc_asm, PropertyTag const property,
-                   GlobalVector& counts)
+inline void LocalLinearLeastSquaresExtrapolator::extrapolateElement(
+    std::size_t const element_index,
+    ExtrapolatableElementCollection const& extrapolatables,
+    GlobalVector& counts)
 {
-    auto const& integration_point_values = loc_asm.getIntegrationPointValues(
-            property, _integration_point_values_cache);
+    auto const& integration_point_values =
+        extrapolatables.getIntegrationPointValues(
+            element_index, _integration_point_values_cache);
 
     // number of nodes in the element
-    const auto nn = loc_asm.getShapeMatrix(0).cols();
+    const auto nn = extrapolatables.getShapeMatrix(element_index, 0).cols();
     // number of integration points in the element
     const auto ni = integration_point_values.size();
 
@@ -80,9 +74,9 @@ extrapolateElement(std::size_t const element_index,
     auto& N = _local_matrix_cache; // TODO make that local?
     N.resize(ni, nn); // TODO: might reallocate very often
 
-    for (auto int_pt=decltype(ni){0}; int_pt<ni; ++int_pt)
-    {
-        auto const& shp_mat = loc_asm.getShapeMatrix(int_pt);
+    for (auto int_pt = decltype(ni){0}; int_pt < ni; ++int_pt) {
+        auto const& shp_mat =
+            extrapolatables.getShapeMatrix(element_index, int_pt);
         assert(shp_mat.cols() == nn);
 
         // copy shape matrix to extrapolation matrix columnwise
@@ -91,7 +85,7 @@ extrapolateElement(std::size_t const element_index,
 
     // TODO make gp_vals an Eigen::VectorXd const& ?
     Eigen::Map<const Eigen::VectorXd> const integration_point_values_vec(
-                integration_point_values.data(), integration_point_values.size());
+        integration_point_values.data(), integration_point_values.size());
 
     // TODO
     // optimization: Store decomposition of N*N^T or N^T for reuse?
@@ -111,19 +105,18 @@ extrapolateElement(std::size_t const element_index,
     // TODO: for now always zeroth component is used
     auto const& global_indices = _local_to_global(element_index, 0).rows;
 
-    _nodal_values.add(global_indices, tmp); // TODO does that give rise to PETSc problems?
+    _nodal_values.add(global_indices,
+                      tmp);  // TODO does that give rise to PETSc problems?
     counts.add(global_indices, std::vector<double>(global_indices.size(), 1.0));
 }
 
-template<typename PropertyTag, typename LocalAssembler>
-void
-LocalLinearLeastSquaresExtrapolator<PropertyTag, LocalAssembler>::
-calculateResiudalElement(std::size_t const element_index,
-                         LocalAssembler const& loc_asm, PropertyTag const property)
+inline void LocalLinearLeastSquaresExtrapolator::calculateResiudalElement(
+    std::size_t const element_index,
+    ExtrapolatableElementCollection const& extrapolatables)
 {
-    auto const& gp_vals = loc_asm.getIntegrationPointValues(
-            property, _integration_point_values_cache);
-    const unsigned ni = gp_vals.size();        // number of gauss points
+    auto const& gp_vals = extrapolatables.getIntegrationPointValues(
+        element_index, _integration_point_values_cache);
+    const unsigned ni = gp_vals.size();  // number of gauss points
 
     // TODO: for now always zeroth component is used
     const auto& global_indices = _local_to_global(element_index, 0).rows;
@@ -131,7 +124,7 @@ calculateResiudalElement(std::size_t const element_index,
     // filter nodal values of the current element
     std::vector<double> nodal_vals_element;
     nodal_vals_element.resize(global_indices.size());
-    for (unsigned i=0; i<global_indices.size(); ++i) {
+    for (unsigned i = 0; i < global_indices.size(); ++i) {
         // TODO PETSc negative indices?
         nodal_vals_element[i] = _nodal_values[global_indices[i]];
     }
@@ -139,10 +132,11 @@ calculateResiudalElement(std::size_t const element_index,
     double residual = 0.0;
     double gp_val_extrapol = 0.0;
 
-    for (unsigned gp=0; gp<ni; ++gp)
-    {
+    for (unsigned gp = 0; gp < ni; ++gp) {
         NumLib::shapeFunctionInterpolate(
-            nodal_vals_element, loc_asm.getShapeMatrix(gp), gp_val_extrapol);
+            nodal_vals_element,
+            extrapolatables.getShapeMatrix(element_index, gp),
+            gp_val_extrapol);
         auto const& ax_m_b = gp_val_extrapol - gp_vals[gp];
         residual += ax_m_b * ax_m_b;
     }
@@ -150,4 +144,4 @@ calculateResiudalElement(std::size_t const element_index,
     _residuals.set(element_index, std::sqrt(residual / ni));
 }
 
-} // namespace ProcessLib
+}  // namespace NumLib

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h

@@ -35,14 +35,9 @@ namespace NumLib
  * system.
  * \endparblock
  */
-template <typename PropertyTag, typename LocalAssembler>
-class LocalLinearLeastSquaresExtrapolator
-    : public Extrapolator<PropertyTag, LocalAssembler>
+class LocalLinearLeastSquaresExtrapolator : public Extrapolator
 {
 public:
-    using LocalAssemblers =
-        typename Extrapolator<PropertyTag, LocalAssembler>::LocalAssemblers;
-
     /*! Constructs a new instance
      *
      * \note
@@ -77,8 +72,8 @@ public:
                "only one component!");
     }
 
-    void extrapolate(LocalAssemblers const& local_assemblers,
-                     PropertyTag const property) override;
+    void extrapolate(
+            ExtrapolatableElementCollection const& extrapolatables) override;
 
     /*! \copydoc Extrapolator::calculateResiduals()
      *
@@ -87,8 +82,8 @@ public:
      * extrapolation results when interpolated back to the integration points
      * again.
      */
-    void calculateResiduals(LocalAssemblers const& local_assemblers,
-                            PropertyTag const property) override;
+    void calculateResiduals(
+            ExtrapolatableElementCollection const& extrapolatables) override;
 
     GlobalVector const& getNodalValues() const override
     {
@@ -108,14 +103,15 @@ public:
 
 private:
     //! Extrapolate one element.
-    void extrapolateElement(std::size_t const element_index,
-                            LocalAssembler const& local_assembler,
-                            PropertyTag const property, GlobalVector& counts);
+    void extrapolateElement(
+        std::size_t const element_index,
+        ExtrapolatableElementCollection const& extrapolatables,
+        GlobalVector& counts);
 
     //! Compute the residuals for one element
-    void calculateResiudalElement(std::size_t const element_index,
-                                  LocalAssembler const& local_assembler,
-                                  PropertyTag const property);
+    void calculateResiudalElement(
+        std::size_t const element_index,
+        ExtrapolatableElementCollection const& extrapolatables);
 
     GlobalVector& _nodal_values;  //!< extrapolated nodal values
     GlobalVector _residuals;      //!< extrapolation residuals
@@ -129,8 +125,7 @@ private:
     //! Avoids frequent reallocations.
     std::vector<double> _integration_point_values_cache;
 };
-}
 
-#include "LocalLinearLeastSquaresExtrapolator-impl.h"
+}  // namespace NumLib
 
 #endif  // NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H