diff --git a/NumLib/Extrapolation/ExtrapolatableElement.h b/NumLib/Extrapolation/ExtrapolatableElement.h
new file mode 100644
index 0000000000000000000000000000000000000000..3a24dc3fdbcf6a653e19325bec9f2adaa8f5ed53
--- /dev/null
+++ b/NumLib/Extrapolation/ExtrapolatableElement.h
@@ -0,0 +1,33 @@
+/**
+ * \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
diff --git a/NumLib/Extrapolation/ExtrapolatableElementCollection.h b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
new file mode 100644
index 0000000000000000000000000000000000000000..7be326d3faa881d6018cf2d57cb4bb35b62ab9b5
--- /dev/null
+++ b/NumLib/Extrapolation/ExtrapolatableElementCollection.h
@@ -0,0 +1,121 @@
+/**
+ * \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
diff --git a/NumLib/Extrapolation/Extrapolator.h b/NumLib/Extrapolation/Extrapolator.h
index 1ac340fa0a75ce0e1498522463c9b05d0e9418b7..ee21daa77840a18042fd9e446290d235ed275fa4 100644
--- a/NumLib/Extrapolation/Extrapolator.h
+++ b/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.
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
similarity index 56%
rename from NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
rename to NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.cpp
index 5e03eece5e3ba05e9c49ec5d82bf2ad31b270eba..2ed625bec24e12eac69166671cc375ba5452d272 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
+++ b/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
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index df3f0ffb093eb3ec788929846ff4e95cf72f7c85..db1154f717cb9134c7bccb7eb6cee06201b51af3 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/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